!*************************************************************************** ! Radiant's arrays start at the toa and go down to the surface !*************************************************************************** program read_first_guess use NETCDF use retrieval_subs use interpolation use radiant2s_driver real,dimension(1) :: vap,lat,lon double precision, dimension(1) :: tsurf,modis_va integer :: status,ncid,varid,i integer :: cldy_gates,numheights,profilenum,profiles,numrad integer :: year,month,day,doy,hour real :: az,el,soldia,soldst double precision :: sza,zenith_angle,albedo,toa_rad_550ch,refl_550ch,toa_rad_21ch,refl_21ch character(len=100) :: filename integer, dimension(:), allocatable :: vert_index_vector real, dimension(:), allocatable :: cl_nx_l,cl_nx_s,cl_lx_l,cl_lx_s,cl_alpha_l,cl_alpha_s,& cl_am_s,cl_bm_s,cl_am_l,cl_bm_l,cl_ba_s,cl_ba_l,shum real, dimension(:), allocatable :: nx_l,nx_s,lx_l,lx_s,alpha_l,alpha_s,& am_s,bm_s,am_l,bm_l,ba_s,ba_l double precision, dimension(:), allocatable :: height,height_km,temp,pres ! Variables for interpolating the radiant text files integer :: mod_numheights double precision, dimension(12) :: mod_height,mod_h2o,mod_co2,mod_o3,ray_scat,& kabs_co2,kabs_h2o double precision, dimension(:), allocatable :: mod_h2o_int,mod_co2_int,mod_o3_int,ray_scat_int,& kabs_h2o_int,kabs_co2_int ! Variables for calculating re and lwc real, dimension(:), allocatable :: re_large_mode,re_small_mode,lwc_large_mode,lwc_small_mode ! Variables from single_scatter_lookup subroutine integer :: n_lwc,n_re,n_alpha real, dimension(:), allocatable :: mie_lwc_0p5,mie_re_0p5,mie_alpha_0p5_1mode real, dimension(:,:,:), allocatable :: mie_sca_0p5_1mode,mie_ext_0p5_1mode,& mie_g_0p5_1mode,mie_ssa_0p5_1mode real, dimension(:), allocatable :: mie_lwc_2p1,mie_re_2p1,mie_alpha_2p1_1mode real, dimension(:,:,:), allocatable :: mie_sca_2p1_1mode,mie_ext_2p1_1mode,& mie_g_2p1_1mode,mie_ssa_2p1_1mode ! Variables for interpolating the lookup table real, dimension(:), allocatable :: ext_lookup_s,ssa_lookup_s,g_lookup_s,& ext_lookup_l,ssa_lookup_l,g_lookup_l ! Variables to be sent into radiant double precision, dimension(:), allocatable :: liq_tau_solar,liq_omega_solar,liq_g_solar,& top_rad,cos_ang double precision,external :: gammam ! Constants - these need to come from netcdf file year=2007 month=01 day=06 !hour=19 !*************************************************************************** ! Read in the profile data from the netcdf file !*************************************************************************** !filename='inversion_inputs_atrain_retrieval_south_pacific_20070106_03691_pro15090.cdf' filename='inversion_inputs_atrain_retrieval_south_pacific_20070106_03691_pro15150.cdf' ! Open the netcdf file status=nf90_open(path=filename,mode=nf90_nowrite,ncid=ncid) ! Read in the dimensions status=nf90_inq_dimid(ncid,"cldy_gates",varid) status=nf90_inquire_dimension(ncid,varid,len=cldy_gates) print*,'cldy_gates',cldy_gates status=nf90_inq_dimid(ncid,"cs_heights",varid) status=nf90_inquire_dimension(ncid,varid,len=numheights) print*,'numheights',numheights status=nf90_inq_dimid(ncid,"profiles",varid) status=nf90_inquire_dimension(ncid,varid,len=profiles) print*,'profiles',profiles ! Allocate the arrays to read in data from netcdf file allocate ( height(numheights)) allocate ( temp(numheights)) allocate ( pres(numheights)) allocate ( shum(numheights)) allocate ( height_km(numheights)) !allocate ( vap(profiles)) allocate ( vert_index_vector(cldy_gates)) allocate ( cl_nx_l(cldy_gates)) allocate ( cl_nx_s(cldy_gates)) allocate ( cl_lx_l(cldy_gates)) allocate ( cl_lx_s(cldy_gates)) allocate ( cl_alpha_l(cldy_gates)) allocate ( cl_alpha_s(cldy_gates)) allocate ( cl_am_s(cldy_gates)) allocate ( cl_bm_s(cldy_gates)) allocate ( cl_am_l(cldy_gates)) allocate ( cl_bm_l(cldy_gates)) allocate ( cl_ba_s(cldy_gates)) allocate ( cl_ba_l(cldy_gates)) ! Read in variables status=nf90_inq_varid(ncid,"orbit_pro_idx",varid) status=nf90_get_var(ncid,varid,profilenum) profilenum=profilenum+1 !indexes from idl start at 0 !print*,'profilenum',profilenum !status=nf90_inq_varid(ncid,"amsre_vap",varid) !status=nf90_get_var(ncid,varid,vap) !print*,'vap',vap(profilenum) status=nf90_inq_varid(ncid,"hour",varid) status=nf90_get_var(ncid,varid,hour) !print*,'hour',hour status=nf90_inq_varid(ncid,"zenith_angle",varid) status=nf90_get_var(ncid,varid,zenith_angle) !print*,'zenith_angle',zenith_angle status=nf90_inq_varid(ncid,"amsre_vap",varid) status=nf90_get_var(ncid,varid,vap,start=(/profilenum/),count=(/1/)) !print*,'vap',vap status=nf90_inq_varid(ncid,"geo_lat",varid) status=nf90_get_var(ncid,varid,lat,start=(/profilenum/),count=(/1/)) !print*,'lat',lat status=nf90_inq_varid(ncid,"geo_lon",varid) status=nf90_get_var(ncid,varid,lon,start=(/profilenum/),count=(/1/)) !print*,'lon',lon status=nf90_inq_varid(ncid,"sktemp",varid) !surface temp (K) status=nf90_get_var(ncid,varid,tsurf,start=(/profilenum/),count=(/1/)) !print*,'tsurf',tsurf status=nf90_inq_varid(ncid,"modis_va",varid) !view angle (degrees from zenith) status=nf90_get_var(ncid,varid,modis_va,start=(/profilenum/),count=(/1/)) !print*,'modis_va',modis_Va status=nf90_inq_varid(ncid,"temp_profile",varid) !Kelvin status=nf90_get_var(ncid,varid,temp) !print*,'temp',(temp(i),i=1,numheights) status=nf90_inq_varid(ncid,"press_profile",varid) !mb status=nf90_get_var(ncid,varid,pres) !print*,'pres',(pres(i),i=1,numheights) status=nf90_inq_varid(ncid,"spec_hum_pro_trk",varid) status=nf90_get_var(ncid,varid,shum,start=(/profilenum,1/),count=(/1,numheights/)) !print*,'shum',(shum(i),i=1,numheights) status=nf90_inq_varid(ncid,"cs_heights",varid) status=nf90_get_var(ncid,varid,height) !print*,'height',(height(i),i=1,numheights) status=nf90_inq_varid(ncid,"vert_index_vector",varid) status=nf90_get_var(ncid,varid,vert_index_vector) !print*,'vert_index_vector',(vert_index_vector(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"Nx_out_precip",varid) status=nf90_get_var(ncid,varid,cl_nx_l) !print*,'cl_nx_l',(cl_nx_l(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"Nx_out_re",varid) status=nf90_get_var(ncid,varid,cl_nx_s) !print*,'cl_nx_s',(cl_nx_s(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"Lx_out_precip",varid) status=nf90_get_var(ncid,varid,cl_lx_l) !print*,'cl_lx_l',(cl_lx_l(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"Lx_out_re",varid) status=nf90_get_var(ncid,varid,cl_lx_s) !print*,'cl_lx_s',(cl_lx_s(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"alpha_out_precip",varid) status=nf90_get_var(ncid,varid,cl_alpha_l) !print*,'cl_alpha_l',(cl_alpha_l(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"alpha_out_re",varid) status=nf90_get_var(ncid,varid,cl_alpha_s) !print*,'cl_alpha_s',(cl_alpha_s(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"am_s",varid) status=nf90_get_var(ncid,varid,cl_am_s) !print*,'cl_am_s',(cl_am_s(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"bm_s",varid) status=nf90_get_var(ncid,varid,cl_bm_s) !print*,'cl_bm_s',(cl_bm_s(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"am_l",varid) status=nf90_get_var(ncid,varid,cl_am_l) !print*,'cl_am_l',(cl_am_l(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"bm_l",varid) status=nf90_get_var(ncid,varid,cl_bm_l) !print*,'cl_bm_l',(cl_bm_l(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"ba_s",varid) status=nf90_get_var(ncid,varid,cl_ba_s) !print*,'cl_ba_s',(cl_ba_s(i),i=1,cldy_gates) status=nf90_inq_varid(ncid,"ba_l",varid) status=nf90_get_var(ncid,varid,cl_ba_l) !print*,'cl_ba_l',(cl_ba_l(i),i=1,cldy_gates) ! Close the netcdf file status = nf90_close(ncid) !*************************************************************************** ! Place the cloudy values into the correct location of the complete profile !*************************************************************************** ! Allocate the arrays for a complete profile of values for the ! variable dimensioned by cldy_gates allocate ( nx_l(numheights)) ; nx_l=0 allocate ( nx_s(numheights)) ; nx_s=0 allocate ( lx_l(numheights)) ; lx_l=0 allocate ( lx_s(numheights)) ; lx_s=0 allocate ( alpha_l(numheights)) ; alpha_l=0 allocate ( alpha_s(numheights)) ; alpha_s=0 allocate ( am_s(numheights)) ; am_s=0 allocate ( bm_s(numheights)) ;bm_s=0 allocate ( am_l(numheights)) ;am_l=0 allocate ( bm_l(numheights)) ;bm_l=0 allocate ( ba_s(numheights)) ;ba_s=0 allocate ( ba_l(numheights)) ;ba_l=0 do i=1,cldy_gates j=vert_index_vector(i)+1 !indexes from idl start at 0 !print*,'i',i,'j',j nx_l(j)=cl_nx_l(i) nx_s(j)=cl_nx_s(i) lx_l(j)=cl_lx_l(i) lx_s(j)=cl_lx_s(i) alpha_l(j)=cl_alpha_l(i) alpha_s(j)=cl_alpha_s(i) am_s(j)=cl_am_s(i) bm_s(j)=cl_bm_s(i) am_l(j)=cl_am_l(i) bm_l(j)=cl_bm_l(i) ba_s(j)=cl_ba_s(i) ba_l(j)=cl_ba_l(i) enddo !print*,(nx_l(i),i=1,numheights) !************************************************ ! Calculate the julian date (day of year) !************************************************ call julian_date(year,month,day,doy) print*,'day of year',doy !*********************************************** ! Calculate the solar zenith angle for radiant ! The sunae value was not the same as the one in the cdf file ! so I am using the one in the cdf file !*********************************************** rhour=real(hour) print*,year,doy,rhour,lat,lon call sunae(year,doy,rhour,lat,lon,az,el,soldia,soldst) print*,az,el,soldia,soldst sza=90.0-el print*,'sza',sza,zenith_angle sza=zenith_angle !*********************************************** ! Calculate the range gate from the height array (cm) !*********************************************** dz=(height(5)-height(4))*100.0 !*********************************************** ! Calculate height in km !*********************************************** height_km=height/1000.0 !heights in km !************************************************ !converting pres from mb to Pa !************************************************ pres = pres*100.0 !converting to Pa !*********************************************** ! Calculate effective radius (cm) and liquid water ! content (g/cm3) for both modes !*********************************************** allocate ( re_large_mode(numheights)) ;re_large_mode=0 allocate ( re_small_mode(numheights)) ;re_small_mode=0 allocate ( lwc_large_mode(numheights)) ;lwc_large_mode=0 allocate ( lwc_small_mode(numheights)) ;lwc_small_mode=0 do i=1,cldy_gates j=vert_index_vector(i)+1 !indexes from idl start at 0 ! cm re_large_mode(j)=(lx_l(j)/2.)*((gammam(dble(alpha_l(j)+ba_l(j)+2.)))/(gammam(dble(alpha_l(j)+ba_l(j)+1.)))) re_small_mode(j)=(lx_s(j)/2.)*((gammam(dble(alpha_s(j)+ba_s(j)+2.)))/(gammam(dble(alpha_s(j)+ba_s(j)+1.)))) ! g/cm3 lwc_large_mode(j)=nx_l(j)*am_l(j)*(lx_l(j)**(bm_l(j)+1.))*(gammam(dble(alpha_l(j)+bm_l(j)+1.))) lwc_small_mode(j)=nx_s(j)*am_s(j)*(lx_s(j)**(bm_s(j)+1.))*(gammam(dble(alpha_s(j)+bm_s(j)+1.))) enddo !*********************************************** ! Read in the single scattering lookup tables for 0.5 and 2.1 microns ! dimensions of the arrays: n_lwc, n_re, n_alpha !*********************************************** call single_scatter_lookup(mie_lwc_0p5,mie_re_0p5,mie_alpha_0p5_1mode,& mie_ext_0p5_1mode,mie_ssa_0p5_1mode,mie_g_0p5_1mode,& mie_lwc_2p1,mie_re_2p1,mie_alpha_2p1_1mode,& mie_ext_2p1_1mode,mie_ssa_2p1_1mode,mie_g_2p1_1mode,n_lwc,n_re,n_alpha) !*********************************************** ! Interpolate the 0.5 look up table scattering values based on lwc,re,alpha ! Create the tau, omega, g profiles !*********************************************** allocate ( ext_lookup_s(numheights)) ;ext_lookup_s=0 allocate ( ssa_lookup_s(numheights)) ;ssa_lookup_s=0 allocate ( g_lookup_s(numheights)) ;g_lookup_s=0 allocate ( ext_lookup_l(numheights)) ;ext_lookup_l=0 allocate ( ssa_lookup_l(numheights)) ;ssa_lookup_l=0 allocate ( g_lookup_l(numheights)) ;g_lookup_l=0 allocate ( liq_tau_solar(numheights)) ;liq_tau_solar=0 allocate ( liq_omega_solar(numheights)) ;liq_omega_solar=0 allocate ( liq_g_solar(numheights)) ;liq_g_solar=0 do i=1,cldy_gates j=vert_index_vector(i)+1 !indexes from idl start at 0 ! Scattering values for the small mode call get_liq_single_scatter(lwc_small_mode(j),re_small_mode(j),alpha_s(j),& ext_lookup_s(j),ssa_lookup_s(j),g_lookup_s(j),& mie_lwc_0p5,mie_re_0p5,mie_alpha_0p5_1mode,& mie_ext_0p5_1mode,mie_ssa_0p5_1mode,mie_g_0p5_1mode,n_lwc,n_re,n_alpha) ! Scattering values for the large mode call get_liq_single_scatter(lwc_large_mode(j),re_large_mode(j),alpha_l(j),& ext_lookup_l(j),ssa_lookup_l(j),g_lookup_l(j),& mie_lwc_0p5,mie_re_0p5,mie_alpha_0p5_1mode,& mie_ext_0p5_1mode,mie_ssa_0p5_1mode,mie_g_0p5_1mode,n_lwc,n_re,n_alpha) !print*,'exts',ext_lookup_s(j),'ssas',ssa_lookup_s(j),'gs',g_lookup_s(j) !print*,'extl',ext_lookup_l(j),'ssal',ssa_lookup_l(j),'gl',g_lookup_l(j) liq_tau_solar(j)=(ext_lookup_s(j)+ext_lookup_l(j))*dz liq_omega_solar(j)=((ssa_lookup_s(j)*ext_lookup_s(j))+& (ssa_lookup_l(j)*ext_lookup_l(j)))/& (ext_lookup_s(j)+ext_lookup_l(j)) liq_g_solar(j)=((g_lookup_s(j)*ext_lookup_s(j))+& (g_lookup_l(j)*ext_lookup_l(j)))/& (ext_lookup_s(j)+ext_lookup_l(j)) enddo !print*,'liq_tau_solar',(liq_tau_solar(i),i=vert_index_vector(1)+1,vert_index_vector(cldy_gates)+1) !print*,'liq_omega_solar',(liq_omega_solar(i),i=vert_index_vector(1)+1,vert_index_vector(cldy_gates)+1) !print*,'liq_g_solar',(liq_g_solar(i),i=vert_index_vector(1)+1,vert_index_vector(cldy_gates)+1) !*********************************************** ! Check for some bad values in the tau,omega,g profiles !*********************************************** ! where liq_tau_solar goes inf or nan, set to 0 ! dont have a way to do this yet ! where liq_omega_solar gt 1, set to 1 where (liq_omega_solar > 1.0) liq_omega_solar = 1.0 end where ! where liq_g_solar gt 0.95, set to 0.95 where (liq_g_solar > 0.95) liq_g_solar = 0.95 end where ! Overiding with exact values from idl to test !liq_tau_solar(7:11:1)=(/0.979723,13.559127,11.495794,1.293457,1.829391/) !liq_omega_solar(7:11:1)=(/1.0,1.0,1.0,1.0,1.0/) !liq_g_solar(7:11:1)=(/0.863588,0.844672,0.846318,0.868940,0.867581/) !**************************************************************** ! This is the file 'modelatm.dat' provided with radiant !**************************************************************** !mass mixing ratios of absorbing gases: (DC standard, 1=h2o, 2=co2, 3=o3) ! alt temp pressure gas mixing ratios ! (km) (k) (bars) 1 2 3 4 !1 50.0 255.1 6.7315E-04 3.56E-06 5.04E-04 3.78E-06 5.24E-09 !2 40.0 251.4 2.5148E-03 3.12E-06 5.04E-04 1.01E-05 2.10E-08 !3 30.0 224.5 1.0352E-02 2.76E-06 5.04E-04 1.41E-05 1.23E-07 !4 20.0 205.5 5.0441E-02 2.19E-06 5.04E-04 4.17E-06 3.32E-07 !5 12.0 211.2 1.9028E-01 5.10E-05 5.04E-04 3.18E-07 4.66E-07 !6 10.0 220.5 2.5907E-01 1.55E-04 5.04E-04 1.91E-07 4.68E-07 !7 8.0 233.8 3.5312E-01 4.61E-04 5.04E-04 1.15E-07 4.69E-07 !8 6.0 247.8 4.6857E-01 1.12E-03 5.04E-04 7.19E-08 4.71E-07 !9 4.0 261.3 6.1124E-01 2.36E-03 5.04E-04 4.61E-08 4.72E-07 !10 2.0 272.8 7.8252E-01 4.30E-03 5.04E-04 3.02E-08 4.73E-07 !11 1.0 276.0 8.9126E-01 5.27E-03 5.04E-04 2.50E-08 4.74E-07 !12 0.0 279.3 1.0000E+00 6.24E-03 5.04E-04 1.97E-08 4.75E-07 mod_numheights=12 ! Instead of reading in a text file, I define the arrays mod_height=(/0.0, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 20.0, 30.0, 40.0, 50.0/) mod_h2o=(/6.24E-03,5.27E-03,4.30E-03,2.36E-03,1.12E-03,4.61E-04,1.55E-04,5.10E-05,2.19E-06,2.76E-06,3.12E-06,3.56E-06/) mod_co2=(/5.04E-04,5.04E-04,5.04E-04,5.04E-04,5.04E-04,5.04E-04,5.04E-04,5.04E-04,5.04E-04,5.04E-04,5.04E-04,5.04E-04/) mod_o3=(/1.97E-08,2.50E-08,3.02E-08,4.61E-08,7.19E-08,1.15E-07,1.91E-07,3.18E-07,4.17E-06,1.41E-05,1.01E-05,3.78E-06/) !**************************************************************** ! This is the file 'sigs_ray_sample.dat' provided with radiant !**************************************************************** !1 !0.17419E-06 0.17420E-06 0.17420E-06 !0.12467E-05 0.12467E-05 0.12468E-05 !0.70061E-05 0.70063E-05 0.70065E-05 !0.26201E-04 0.26201E-04 0.26202E-04 !0.49202E-04 0.49203E-04 0.49205E-04 !0.63246E-04 0.63248E-04 0.63250E-04 !0.80460E-04 0.80462E-04 0.80464E-04 !0.10129E-03 0.10130E-03 0.10130E-03 !0.12619E-03 0.12619E-03 0.12619E-03 !0.14761E-03 0.14762E-03 0.14762E-03 !0.16345E-03 0.16346E-03 0.16346E-03 ray_scat=(/0.16345E-03,0.14761E-03,0.12619E-03,0.10129E-03,0.80460E-04,0.63246E-04,0.49202E-04,& 0.26201E-04,0.70061E-05,0.12467E-05,0.17419E-06,0.17419E-06/) !**************************************************************** ! This is the file 'kabs_co2_sample.dat' !**************************************************************** !1 !0.55235E-03 0.20187E-02 0.12509E+02 !0.20899E-02 0.76306E-02 0.12672E+02 !0.92240E-02 0.33635E-01 0.12344E+02 !0.44683E-01 0.16216E+00 0.91816E+01 !0.15008E+00 0.50638E+00 0.42721E+01 !0.18648E+00 0.59639E+00 0.33736E+01 !0.22591E+00 0.67403E+00 0.26625E+01 !0.26673E+00 0.72868E+00 0.21026E+01 !0.30607E+00 0.75156E+00 0.16720E+01 !0.34039E+00 0.74038E+00 0.13369E+01 !0.35513E+00 0.72449E+00 0.11990E+01 kabs_co2=(/0.35513E+00,0.34039E+00,0.30607E+00,0.26673E+00,0.22591E+00,0.18648E+00,0.15008E+00,0.44683E-01,& 0.92240E-02,0.20899E-02,0.55235E-03,0.55235E-03/) !**************************************************************** ! This is the file 'kabs_h2o_sample.dat' !**************************************************************** !1 !0.52274E-05 0.48387E-05 0.46884E-05 !0.19327E-04 0.17908E-04 0.17359E-04 !0.82683E-04 0.76786E-04 0.74505E-04 !0.39306E-03 0.36591E-03 0.35541E-03 !0.13570E-02 0.12658E-02 0.12300E-02 !0.17422E-02 0.16259E-02 0.15794E-02 !0.22013E-02 0.20590E-02 0.20006E-02 !0.27474E-02 0.25809E-02 0.25097E-02 !0.33839E-02 0.31981E-02 0.31146E-02 !0.41220E-02 0.39250E-02 0.38311E-02 !0.45383E-02 0.43400E-02 0.42420E-02 kabs_h2o=(/0.45383E-02,0.41220E-02,0.33839E-02,0.27474E-02,0.22013E-02,0.17422E-02,0.13570E-02,0.39306E-03,& 0.82683E-04,0.19327E-04,0.52274E-05,0.52274E-05/) !*********************************************** ! Interpolate the provided model atmospere to the profile heights !printf,1, mod_o3_pro[i], mod_co2_pro[i], mod_h2o_pro[i], ray_scat_pro[i], kabs_h2o_pro[i], kabs_co2_pro[i] !*********************************************** allocate ( mod_h2o_int(numheights)) ;mod_h2o_int=0 allocate ( mod_co2_int(numheights)) ;mod_co2_int=0 allocate ( mod_o3_int(numheights)) ;mod_o3_int=0 allocate ( ray_scat_int(numheights)) ;ray_scat_int=0 allocate ( kabs_h2o_int(numheights)) ;kabs_h2o_int=0 allocate ( kabs_co2_int(numheights)) ;kabs_co2_int=0 do i=1,numheights j=1 do while(mod_height(j) .le. height_km(i) .and. j .lt. mod_numheights) j=j+1 enddo mod_h2o_int(i)=mod_h2o(j-1)+(height_km(i)-mod_height(j-1))*((mod_h2o(j)-mod_h2o(j-1))/(mod_height(j)-mod_height(j-1))) mod_co2_int(i)=mod_co2(j-1)+(height_km(i)-mod_height(j-1))*((mod_co2(j)-mod_co2(j-1))/(mod_height(j)-mod_height(j-1))) mod_o3_int(i)=mod_o3(j-1)+(height_km(i)-mod_height(j-1))*((mod_o3(j)-mod_o3(j-1))/(mod_height(j)-mod_height(j-1))) ray_scat_int(i)=ray_scat(j-1)+(height_km(i)-mod_height(j-1))*((ray_scat(j)-ray_scat(j-1))/(mod_height(j)-mod_height(j-1))) kabs_h2o_int(i)=kabs_h2o(j-1)+(height_km(i)-mod_height(j-1))*((kabs_h2o(j)-kabs_h2o(j-1))/(mod_height(j)-mod_height(j-1))) kabs_co2_int(i)=kabs_co2(j-1)+(height_km(i)-mod_height(j-1))*((kabs_co2(j)-kabs_co2(j-1))/(mod_height(j)-mod_height(j-1))) enddo !*********************************************** ! Choose albedo based on wavelength ! from Sun-Mack et al., Proc SPIE !*********************************************** albedo=0.08 !0.550 !albedo=0.02 !1.6 broken on aqua !albedo=0.015 !2.1 !*********************************************** ! Run radiant for the 550 channel case ! Radiant uses double precision instead of real values !************************************************ ! All these arrays need to be reversed because radiant expects toa to surface ! array = array(n:1:-1) height_km=height_km(numheights:1:-1) pres=pres(numheights:1:-1) temp=temp(numheights:1:-1) liq_tau_solar=liq_tau_solar(numheights:1:-1) liq_omega_solar=liq_omega_solar(numheights:1:-1) liq_g_solar=liq_g_solar(numheights:1:-1) mod_o3_int=mod_o3_int(numheights:1:-1) mod_co2_int=mod_co2_int(numheights:1:-1) mod_h2o_int=mod_h2o_int(numheights:1:-1) ray_scat_int=ray_scat_int(numheights:1:-1) kabs_h2o_int=kabs_h2o_int(numheights:1:-1) kabs_co2_int=kabs_co2_int(numheights:1:-1) ! Prints profiles to the screen and writes a file !10 FORMAT(I2,1X,F9.6,1X,F7.3,1X,F7.4,1X,F8.5,1X,F10.3,1X,F7.3,1X,F9.6,1X,F9.6,1X,F9.6) !open(unit=10,file='cloud_pro_sb.txt') !do i=1,numheights !write(*,10)i,sza,tsurf(1),modis_va(1),height_km(i),pres(i),temp(i),liq_tau_solar(i),& ! liq_omega_solar(i),liq_g_solar(i) ! write(10,10)i,sza,tsurf(1),modis_va(1),height_km(i),pres(i),temp(i),liq_tau_solar(i),& ! liq_omega_solar(i),liq_g_solar(i) !enddo !close(10) !20 FORMAT(E11.5,1X,F11.9,1X,E11.5,1X,E11.5,1X,E11.5,1X,F11.9) !open(unit=11,file='mod_atm_int.txt') !do i=1,numheights !write(*,20)mod_o3_int(i),mod_co2_int(i),mod_h2o_int(i),ray_scat_int(i),kabs_h2o_int(i),& ! kabs_co2_int(i) ! write(11,20)mod_o3_int(i),mod_co2_int(i),mod_h2o_int(i),ray_scat_int(i),kabs_h2o_int(i),& ! kabs_co2_int(i) !enddo !close(11) ! This runs radiant. The output is top_rad call radiant_driver(sza,tsurf(1),modis_va(1),height_km,pres,temp,liq_tau_solar,& liq_omega_solar,liq_g_solar,numheights,mod_o3_int,mod_co2_int,mod_h2o_int,& ray_scat_int,kabs_h2o_int,kabs_co2_int,albedo,top_rad) ! Get the size of the top_rad array since it was allocated in radiant numrad=size(top_rad) ! Radiance for the channel, using the last stream toa_rad_550ch=top_rad(numrad) ! Calculate reflectance refl_550ch=0.0 if (toa_rad_550ch .gt. 0.0) then refl_550ch=(toa_rad_550ch*3.14)/COS(sza*(3.14/180.0)) end if ! Results for 550 channel print*,toa_rad_550ch,refl_550ch,'rad,refl for 550' ! Prints profiles !10 FORMAT(I2,1X,F9.6,1X,F7.3,1X,F7.4,1X,F8.5,1X,F7.3,1X,F7.3,1X,F9.6,1X,F9.6,1X,F9.6) !open(unit=12,file='cloud_pro_sb_after.txt') !do i=1,numheights !write(*,10)i,sza,tsurf(1),modis_va(1),height_km(i),pres(i),temp(i),liq_tau_solar(i),& !liq_omega_solar(i),liq_g_solar(i) ! write(12,10)i,sza,tsurf(1),modis_va(1),height_km(i),pres(i),temp(i),liq_tau_solar(i),& ! liq_omega_solar(i),liq_g_solar(i) !enddo !close(12) !20 FORMAT(E11.5,1X,F11.9,1X,E11.5,1X,E11.5,1X,E11.5,1X,F11.9) !open(unit=13,file='mod_atm_int_after.txt') !do i=1,numheights !write(*,20)mod_o3_int(i),mod_co2_int(i),mod_h2o_int(i),ray_scat_int(i),kabs_h2o_int(i),& !kabs_co2_int(i) ! write(13,20)mod_o3_int(i),mod_co2_int(i),mod_h2o_int(i),ray_scat_int(i),kabs_h2o_int(i),& ! kabs_co2_int(i) !enddo !close(13) !************************************************************** !************************************************************** ! Now run for 2.1 case !************************************************************** !************************************************************** ext_lookup_s=0 ssa_lookup_s=0 g_lookup_s=0 ext_lookup_l=0 ssa_lookup_l=0 g_lookup_l=0 liq_tau_solar=0 liq_omega_solar=0 liq_g_solar=0 do i=1,cldy_gates j=vert_index_vector(i)+1 !indexes from idl start at 0 ! Scattering values for the small mode call get_liq_single_scatter(lwc_small_mode(j),re_small_mode(j),alpha_s(j),& ext_lookup_s(j),ssa_lookup_s(j),g_lookup_s(j),& mie_lwc_2p1,mie_re_2p1,mie_alpha_2p1_1mode,& mie_ext_2p1_1mode,mie_ssa_2p1_1mode,mie_g_2p1_1mode,n_lwc,n_re,n_alpha) ! Scattering values for the large mode call get_liq_single_scatter(lwc_large_mode(j),re_large_mode(j),alpha_l(j),& ext_lookup_l(j),ssa_lookup_l(j),g_lookup_l(j),& mie_lwc_2p1,mie_re_2p1,mie_alpha_2p1_1mode,& mie_ext_2p1_1mode,mie_ssa_2p1_1mode,mie_g_2p1_1mode,n_lwc,n_re,n_alpha) liq_tau_solar(j)=(ext_lookup_s(j)+ext_lookup_l(j))*dz liq_omega_solar(j)=((ssa_lookup_s(j)*ext_lookup_s(j))+& (ssa_lookup_l(j)*ext_lookup_l(j)))/& (ext_lookup_s(j)+ext_lookup_l(j)) liq_g_solar(j)=((g_lookup_s(j)*ext_lookup_s(j))+& (g_lookup_l(j)*ext_lookup_l(j)))/& (ext_lookup_s(j)+ext_lookup_l(j)) enddo ! where liq_tau_solar goes inf or nan, set to 0 ! dont have a way to do this yet ! where liq_omega_solar gt 1, set to 1 where (liq_omega_solar > 1.0) liq_omega_solar = 1.0 end where ! where liq_g_solar gt 0.95, set to 0.95 where (liq_g_solar > 0.95) liq_g_solar = 0.95 end where ! Overiding with exact values from idl to test !liq_tau_solar(7:11:1)=(/0.979723,13.559127,11.495794,1.293457,1.829391/) !liq_g_solar(7:11:1)=(/0.863588,0.844672,0.846318,0.868940,0.867581/) albedo=0.015 !2.1 ! Reverse these arrays liq_tau_solar=liq_tau_solar(numheights:1:-1) liq_omega_solar=liq_omega_solar(numheights:1:-1) liq_g_solar=liq_g_solar(numheights:1:-1) ! Prints profiles to the screen and writes a file !10 FORMAT(I2,1X,F9.6,1X,F7.3,1X,F7.4,1X,F8.5,1X,F10.3,1X,F7.3,1X,F9.6,1X,F9.6,1X,F9.6) !open(unit=10,file='cloud_pro_sb2p1.txt') !do i=1,numheights !write(*,10)i,sza,tsurf(1),modis_va(1),height_km(i),pres(i),temp(i),liq_tau_solar(i),& !liq_omega_solar(i),liq_g_solar(i) !write(10,10)i,sza,tsurf(1),modis_va(1),height_km(i),pres(i),temp(i),liq_tau_solar(i),& !liq_omega_solar(i),liq_g_solar(i) !enddo !close(10) !20 FORMAT(E11.5,1X,F11.9,1X,E11.5,1X,E11.5,1X,E11.5,1X,F11.9) !open(unit=11,file='mod_atm_int2p1.txt') !do i=1,numheights !write(*,20)mod_o3_int(i),mod_co2_int(i),mod_h2o_int(i),ray_scat_int(i),kabs_h2o_int(i),& !kabs_co2_int(i) !write(11,20)mod_o3_int(i),mod_co2_int(i),mod_h2o_int(i),ray_scat_int(i),kabs_h2o_int(i),& !kabs_co2_int(i) !enddo !close(11) ! This runs radiant. The output is top_rad call radiant_driver(sza,tsurf(1),modis_va(1),height_km,pres,temp,liq_tau_solar,& liq_omega_solar,liq_g_solar,numheights,mod_o3_int,mod_co2_int,mod_h2o_int,& ray_scat_int,kabs_h2o_int,kabs_co2_int,albedo,top_rad) ! Get the size of the top_rad array since it was allocated in radiant numrad=size(top_rad) ! Radiance for the channel, using the last stream toa_rad_21ch=top_rad(numrad) ! Calculate reflectance refl_21ch=0.0 if (toa_rad_21ch .gt. 0.0) then refl_21ch=(toa_rad_21ch*3.14)/COS(sza*(3.14/180.0)) end if ! Results for 550 channel print*,toa_rad_21ch,refl_21ch,'rad,refl for 2.1' end program read_first_guess