c************************************************************************ subroutine stalk_extinction(N_layer, tauv,tauh) c************************************************************************ integer nl parameter(nl=250) c consider replacing the above with an include! complex th,tv complex tauh(nl),tauv(nl) real d,del_leaf,z_lc,hole common /layer_struct/ d,del_leaf,z_lc,hole real z_delta,z,z_shift integer N_layer,layer integer s_tensor_flag,s_tau_flag,z0_flag,s_field_flag(2,4) common /stalk_flags/ s_tensor_flag,s_tau_flag & ,z0_flag,s_field_flag real shadow_bound common /shadow/ shadow_bound c s_tau_flag=0 c s_tau_flag serves no purpose ! z_delta=-d/N_layer do layer=1,N_layer z=layer*z_delta z_shift=z+shadow_bound call stalk_tau(z_shift, tv,th) tauh(layer)=th tauv(layer)=tv enddo return end c************************************************************************ subroutine stalk_tau(zz, tauv,tauh) c************************************************************************ integer v,h parameter(v=1) parameter(h=2) complex ii,coef_fs,coef_ks complex tauh,tauv complex tauh_ka,tauv_ka complex tauh_ks,tauv_ks complex stalk_fs_int,stalk_ks_int real pi,mpi real zz integer pol,mech,mech_c common /index/ pol,mech,mech_c real theta_i,phi_i,k_0 common /incident/ theta_i,phi_i,k_0 complex stalk_fs_ave_h,stalk_fs_ave_v complex stalk_ks_ave_h,stalk_ks_ave_v common /store_st/ stalk_fs_ave_h,stalk_fs_ave_v, & stalk_ks_ave_h,stalk_ks_ave_v real N_ppscm common /unfrm_data/ N_ppscm integer s_tensor_flag,s_tau_flag,z0_flag,s_field_flag(2,4) common /stalk_flags/ s_tensor_flag,s_tau_flag & ,z0_flag,s_field_flag real z common /tau_position/ z external stalk_fs_int,stalk_ks_int z=zz ii=cmplx(0.0,1.0) pi=3.14159265 mpi=-pi coef_fs=(ii*2.0*pi*N_ppscm/k_0) coef_ks=(-N_ppscm/2.0) pol=h call integral1(stalk_fs_int,mpi,pi, stalk_fs_ave_h) c write(6,*) 'stalk_fs_ave_h = ',stalk_fs_ave_h call integral1(stalk_ks_int,mpi,pi, stalk_ks_ave_h) pol=v call integral1(stalk_fs_int,mpi,pi, stalk_fs_ave_v) c write(6,*) 'stalk_fs_ave_v = ',stalk_fs_ave_v c write(6,*) 'coef = ',coef call integral1(stalk_ks_int,mpi,pi, stalk_ks_ave_v) tauh_ka=coef_fs*stalk_fs_ave_h tauv_ka=coef_fs*stalk_fs_ave_v tauh_ks=coef_ks*stalk_ks_ave_h tauv_ks=coef_ks*stalk_ks_ave_v tauh=tauh_ka+tauh_ks tauv=tauv_ka+tauv_ks write(6,*) 'tauh,tauh_ka,tauh_ks' write(6,*) tauh,tauh_ka,tauh_ks write(6,*) 'tauv,tauv_ka,tauv_ks' write(6,*) tauv,tauv_ka,tauv_ks c write(6,*)'tauv,tauh ',tauv,tauh return end c************************************************************************ complex function stalk_fs_int(phii) c************************************************************************ complex g_st,result real stalk_alpha_pdf real amin,amax,z0min,z0max real phi_s,a_s,z0 real phii common /stalk_limits/ amin,amax,z0min,z0max common /stalk_state/ phi_s,a_s,z0 external g_st phi_s=phii c print*,'in stalk_fs_int' call integral2(g_st,amin,amax,result) stalk_fs_int=result*stalk_alpha_pdf(phi_s) c print*,'in stalk_fs_int = ',stalk_fs_int return end c************************************************************************ function g_st(aa) c************************************************************************ complex g_st,i_st,result real stalk_a_pdf real amin,amax,z0min,z0max real phi_s,a_s,z0 real aa external i_st common /stalk_limits/ amin,amax,z0min,z0max common /stalk_state/ phi_s,a_s,z0 a_s=aa c print*,'in g_st' call integral3(i_st,z0min,z0max,result) g_st=result*stalk_a_pdf(a_s) return end c************************************************************************ function i_st(zz0) c************************************************************************ complex i_st,stalk_forward_scat,sfs real zz0 real stalk_z0_pdf real z common /tau_position/ z integer pol,mech,mech_c common /index/ pol,mech,mech_c real phi_s,a_s,z0 common /stalk_state/ phi_s,a_s,z0 z0=zz0 c print*,'in i_st' sfs=stalk_forward_scat(phi_s,a_s,z0,pol,z) i_st=sfs*stalk_z0_pdf(z0) c print*,sfs,phi_s,a_s,z0,pol,z return end c************************************************************************ complex function stalk_ks_int(phii) c************************************************************************ complex g_sks,result real stalk_alpha_pdf real amin,amax,z0min,z0max real phi_s,a_s,z0 real phii common /stalk_limits/ amin,amax,z0min,z0max common /stalk_state/ phi_s,a_s,z0 external g_sks phi_s=phii c print*,'in stalk_ks_int' call integral2(g_sks,amin,amax,result) stalk_ks_int=result*stalk_alpha_pdf(phi_s) c print*,'in stalk_ks_int = ',stalk_ks_int return end c************************************************************************ function g_sks(aa) c************************************************************************ complex g_sks,i_sks,result real stalk_a_pdf real amin,amax,z0min,z0max real phi_s,a_s,z0 real aa external i_sks common /stalk_limits/ amin,amax,z0min,z0max common /stalk_state/ phi_s,a_s,z0 a_s=aa c print*,'in g_sks' call integral3(i_sks,z0min,z0max,result) g_sks=result*stalk_a_pdf(a_s) return end c************************************************************************ function i_sks(zz0) c************************************************************************ complex i_sks,stalk_ks,sks real zz0 real stalk_z0_pdf real z common /tau_position/ z integer pol,mech,mech_c common /index/ pol,mech,mech_c real phi_s,a_s,z0 common /stalk_state/ phi_s,a_s,z0 z0=zz0 c print*,'in i_sks' sks=stalk_ks(phi_s,a_s,z0,pol,z) i_sks=sks*stalk_z0_pdf(z0) c print*,sks,phi_s,a_s,z0,pol,z return end