program FRONT_AERONET C*********************************************************************** C AUTHOR: Nickolay Krotkov and Alexander P. Vasilkov C C FUNCTION: c To fit imaginary refr index C using AERONET PSD, real refr index and MFRSR diffuse and total irradiance measurenets C C WARNING c 8) MFRSR AOT need to be corrected at the calibration step for the difference between AERONET assumed NO2 and actual measured NO2 c 7) Temperature and NO2 profiles are not correct c 6) The program can not handle SZA=??.?5 because of filename trucation c temporal fix: put fill values in the output: SSA=-999 c 5) Surface albedo at all wavelengths = AERONET V2 albedo at 440nm !! (underestimate SSA) C 1) Only one wavelength version (lambda(1) and icount=1 ONLY !!! ) C re-calculated at each wavelength ! C 2) only bi-modal log-normal PSD C 3) Refractive index is the same for both modes (one component model) C 4) surface pressure should be 1 atm=1013.25 in both flat.in and pressure.prf C (otherwise tau_aer is changed)! Pressure in flat.in can not be larger than C surface pressure in pressure.prf file. C C HISTORY: c modified Aug 12 2015 [nk] to use box aerosol profile for Bolivia c modified Mar 25 2014 to add new input with mfrsr_dump(35) = f_vbias [mV] c modified dec 20 2014 to change output NO2 format, dd_fit=.true. cimelaot=.true. c modified dec 8 2013 to process using total irradiance fitting t_fit=.true. c modified Dec 5-7 2013 to read new input with mfrsr_dump(34) = water_aeronet [cm] and Jacobians from gs_fit(...dd_k,dss_dk,das_dk) c modified August9 2013 to process new UV MFRSR527 2005-2013 min_aot440>0.25 using MFRSR AOT cimelaot=.true. c modified July 28 2013 to process new UV MFRSR527 2010-2013 AOT440>0.25 using MFRSR AOT cimelaot=.true. c modified July 12 2013 to process new UV MFRSR582 2010-2013 AOT440>0.25 using MFRSR AOT cimelaot=.true. c modified July 9 2013 to process new VIS MFRSR579 2012 data AOT440>0.25 using MFRSR AOT cimelaot=.false. c modified July 8 2013 to process new VIS MFRSR579 2012 data AOT440>0.25 using cimelaot=.true. c modified June 26 2013 to process new MFRSR582 data c modified June 7 2013 to use MFRSR AOT in inversion : cimelaot=.false. c modified June 6 2013 to use AOT fliter AOT440 > 0.3 c modified may 15 2013 to use AERONET extrapolated AOT in inversion : cimelaot=.true. c modified may 7 2013 to read re-calibrated 527 match file c modified jan 18 2013 to read new 582 match file c modified Dec 09 2012 to read data from new calibration format c mfrsr_dump(31) = tau_no2chan c modified Dec 07 2012 to read data from UV-MFRSR582 c modified Sep 21 2012 to use relative fit_tolerance=ABS( calc/meas -1.) c modified Sep 21 2012 to remove AOT correction due to mismatch between AERONET NO2 climatology and measured NO2 (should be done during calibration) c modified aug 07 2012 to process new match file: match21.440nm.nofit.lev15.txt : using dd fit and tolerance=0.001 c modified jul 08 2012 to process new match file: match21 440nm : using dd fit and tolerance=0.001 c modified jul 03 2012 to constrain maximal imaginary refractive index k_real< max_k_real=1.0 [max BC value] : tolerance = 0.0002 c modified jun 19 2012 change max number of iteraqtion to 5 to process new match file: match21 305nm : tolerance=0.0001 c modified jun 17 2012 to process new match file: match21 440nm : tolerance=0.0001 c modified feb 01 2012 to process new match file: match20.440nm.gsfc.test c modified oct 01 2011 to process new match file: match20.368nm.gsfc.test and new AERONET inversion file with 156 columns c modified sep 27 2011 to read new match file: match20.440nm.gsfc.test and new AERONET inversion file with 156 columns c Modified aug 25 , New AERONET retrieval on 09/27 c Modified aug 09 , surf albedo 0.01 at 317nm 2008 Santa Cruz , Bolivia c Modified aug 04 , 2008 Santa Cruz , Bolivia c Modified jul 15 , 2008 Decrease AERONET surf. albedo by 0.02 for sensitivty study (need removal) c Modified jul 11 , 2008 Remove unnecessary pressure scaling for tau_rau c Modified may 22 , 2008 To process Mexico campaign at 368nm c Modified Feb 13 , 2008 To re-process Santa Cruz campaign at 368nm c Modified Dec 13 , 2007 To re-process Santa Cruz campaign at 325nm c Modified Oct 15 , 2007 include AERONET clim NO2 corrections c Modified Oct 13 , 2007 to re-process aeronet version 2 inversions at 440nm in Santa Cruz without NO2 corrections c Modified Oct 13 , 2007 scale t_ray with surface pressure c Modified Oct 13 , 2007 set NO2=0 for fill values c Modified May 11 , 2007 to re-process all aeronet version 2 inversions at 368nm with NO2 corrections c Modified May 8 , 2007 to re-process all aeronet version 2 inversions at 440nm with NO2 corrections c Modified July 15 , 2005 to re-process all inversions with NO2 corrections c Modified January 9, 2005 to include ozone absorption and replace O2-O2 by NO2 in water vapor absorption c Modified January 5, 2005 get rid of coe.dat file ( read lambda and x-sections from MFRSR input file) c Modified January 4, 2005 reduce AOT according to tau_NO2 at 368nm c Modified January 4, 2005 change ozone (NO2) profile to 1st Umkehr layer only c Modified January 4, 2005 change total ozone to total NO2 (still ozone profile) c Modified January 4, 2005 to add NO2 correction (using Brewer direct sun NO2 measurements) c Modified March 19, 2004 to add interpolated AERONET AOT to the output c Modified March 18, 2004 to compensate for ARIZONA code diffuse flux overestimation by 0.3% c which is equivalent to enhancing measured diffuse by the same amount *1.003 c Modified March 16, 2004 to fix bug with new format input file reading c Modified March 15, 2004 to re-calculate total after diffuse correction c Modified March 14, 2004 to change diffuse due to additional diffuse cosine correction! c Modified March 11, 2004 to use CT fit: t_fit=.true. c Modified March 10, 2004 to use dt fit tolerance=0.0001 and ipseudo=1 c Modified March 10, 2004 to use diffuse/total fraction fit: dt_fit=.true. c Modified March 2, 2004 to use absolute DD fit tolerance= 0.001 c Modified March 2, 2004 to read new input data format (AERONET+20mfrsr) c Modified Sep 28, 2003 to double AERONET k c Modified Sep 26, 2003 to use only AERONET data and no fitting c Modified July 30, 2003 to use MFRSR AOT correction at 368nm only !!! c Modified July 30, 2003 to fix output format for RVC>10micron c Modified July 29, 2003 to use AERONET interpolated AOT at 368nm for SSA at 368nm c Modified July 18, 2003 to read new format AERONET/MFRSR data c Modified July 10, 2003 restrict fitting to SZA<70deg c Modified July 8 , 2003 to use absolute DD fit tolerance= 0.01 c Modified July 6 , 2003 to constrain number of k iterations: max_k_iter c Modified July 6 , 2003 to process missing MFRSR AOT data for zero direct irradiance) c Developed June - July, 2003 C template program: front_end_TOMRAD_GS_V6_440nm.f Developed by AV in June, 2000-June 2003 C*********************************************************************** implicit none integer :: imfrsr_ldir integer aeronet_dump_length parameter (aeronet_dump_length=156) ! added new solar zenith average column (next to second last) c parameter (aeronet_dump_length=155) ! V2 AERONET inversions: c:\abstrats\AGU_JA07\memo_v2 real aeronet_dump(aeronet_dump_length) ! AERONET web data only: ! it's also affects get_aer_mf_v2.f and get_aeronet subs !!! real c0b,c1b,c2b,rhob ! Ozone, and depolarization values are defined here and passed to gs_fit to get rid of coe.dat gas2.coe files ! Rayleigh optical thickness is red from MFRSR input file real xno2_ap,xno2_used ! xno2_ap a priori NO2 abs coeff [atm-cm-1] averaged over MFRSR 2nm bandwidth (10nm for 440nm channel) ! xno2_used NO2 abs coeff [atm-cm-1] used as input to RT calculations logical real_refr_const ! .TRUE. to use constant real refractive index (ignore AERONET retrieval) real real_refr_value parameter (real_refr_const=.false.,real_refr_value=1.4) ! set constant value for real refractive index logical dd_fit ! .T. use diffuse/direct fit of k logical dt_fit ! .T. use diffuse/total fit of k logical t_fit ! .T. use total fit of k logical d_fit ! .T. use diffuse transmittance fit of k logical log_fit ! .T. use logarithms of dd and k logical no_radiance_output c--- choose only one true fitting parameter and corresponding fit_tolerance : d_fit,dd_fit,dt_fit,t_fit real fit_tolerance ! relative fit tolerance: 0.01= 1% c23456789012345678901234567890123456789012345678901234567890123456789012 parameter (d_fit=.false.,dt_fit=.false., *dd_fit=.true.,t_fit=.false.) parameter (fit_tolerance=0.001) ! 0.1% use with dd fit: absolute parameter fit tolerance c *dd_fit=.false.,t_fit=.true.) c parameter (fit_tolerance=0.0002) ! smaller tolerace use with fit: absolute parameter fit tolerance for 305nm c--- choose type of fitting (log vs linear) parameter (log_fit=.false.) parameter (no_radiance_output=.true.) c--- choose min CIMEL AOT for MFRSR inversion processing real min_aot440 parameter (min_aot440 = 0.15) c--- choose first multilpier for the first iteration real firstmult_k parameter (firstmult_k=2.0) ! 332nm and 368nm Initial iAERONET k multiplyer to start 1st iteration c parameter (firstmult_k=20.0) ! 305nm Initial iAERONET k multiplyer to start 1st iteration c parameter (firstmult_k=0.1) ! 440nm and 368nm Initial iAERONET k multiplyer to start 1st iteration c--- choose constraint on Imaginary part of retrieved refractive index max_k_real<1.0 (pure BC value) real max_k_real parameter (max_k_real=1.0) c--- choose what corrections (if any) are applied to MFRSR measurements logical mtau_corr,mdif_corr parameter (mtau_corr=.false.,mdif_corr=.false.) ! do not correct MFRSR AOT and/or measured diffuse transmittance DIF c parameter (mtau_corr=.false.,mdif_corr=.true.) ! correct MFRSR AOT and/or measured diffuse transmittance DIF ( for sensitivity study) real mean_lnk,mdif_corr_cos parameter (mean_lnk=+0.08,mdif_corr_cos=1.30) ! dtau=myu*, increase measured diffuse irradiance to account for uncorrected exces sky radiance blocking by shadowband c--- choose whether to use AERONET extrapolated AOT or MFRSR AOT logical cimelaot c parameter (cimelaot=.false.) ! use MFRSR AOT in the flux forward calculation parameter (cimelaot=.true.) ! use CIMEL AOT in DD fit real alpha c--- corection for gaseous absorption real gas2maxh,gas2width,gas2tot ! gas2tot (other than ozone =nO2 [atm-cm]) is obtained from the input UV-MFRSR-CIMEL merged file parameter (gas2maxh=0.0) ! Gas2 (NO2) constant mixing profile in the boundary layer parameter (gas2width=800.0) ! top of the boundary layer [mb], constant mixing profile logical mie_only parameter (mie_only=.false.) logical gsfc_psd_clim,gsfc_refr_clim,ref_440 parameter (gsfc_psd_clim=.false.,gsfc_refr_clim=.false.) ! .T. to use AERONET climatological values for PSD and refr. index c--- if gsfc_refr_clim=.true. then ref_440 does not matter parameter (ref_440 =.true. ) ! Use constant value of ref_440 (real) to extrapolate into UV wavelengths ! use linear exrapolation (440,670) c---- Fitting parameters integer iskip,max_k_iter,n_iter(3) parameter (iskip=0,max_k_iter=11) ! resume retrievals skipping X raws ; c parameter (iskip=0,max_k_iter=0) ! No mfrsr fitting real err_dd(3),dd_new(3),jacobian_k_dd,jacobian_ssa_dk real jacobian_asym_dk integer i,j,k,nstr,maxnum,iaertype,nlayer integer iaermod,isurfmodel,icount,iozpr,nozpr integer jj,jheader,nvza,nphi c vert azimuth parameter (nvza=6,nphi=7) ! max dimensions of vza and phi arrays parameter (icount=1) ! just one wavelength calculation (determined by MFRSR input file) real sza,szab(10),vza(nvza),phi(nphi),vzatom,phitom,vzar,phir real lambda(5000),wave real aermaxh,aerwidthh,lat ! Define aerosol layer height (peak conc for Gaussian ormiddle for box) real wavea,lamalb(12),surfalb(12),salb440v2 real*8 surfpres,ozone(11),ozoneb(11),temp(11) real refrr1(12),refri1(12),refrr2(12),refri2(12),lam(12) real sigma1,rmean1,Nc1,sigma2,rmean2,Nc2 real tau440,totoz,taumd,taur,tauc,pi,dt real tau_real(1000),om_real(1000),asym_real(1000) real n_real(1000),k_real(1000) real ratio(5000,nvza,nphi),taum(12),omega(12) real rad_rayl(5000,nvza,nphi) real e_dir_ray(50),e_dir_ratio(50),e_dif_ray(50),e_dif_ratio(50) real e_dir(50),e_dif(50),dd_calc(50) c--- Aeronet Water measurements real water_aeronet c--- Nighttime voltage bias real f_vbias c--- aeronet AOT measurements integer iaeronet real aeronet_aotchan real aot340m,aot380m,aot440m,aot500m,aot670m,aot870m,aot1020m c--- aeronet retrievals real oleg_aot340,oleg_aot380,oleg_aot440,oleg_aot670,oleg_aot870 real oleg_aot1020 real oleg_w440,oleg_w670,oleg_w870,oleg_w1020 real oleg_dayofyear integer oleg_dd,oleg_month,oleg_yyyy,oleg_hh,oleg_mm,oleg_ss c--- aeronet size distribution parameters appoximated as bi-modal log-normal distribution real lnsf,rvf,vf,lnsc,rvc,vc ! Volumne size distribution parameters (AERONET) real rnf,nf, rnc,nc ! Number size distribution parameters (ARIZONA RT) real rr440,rr670,rr870,rr1020 ! Oleg's refractive index real ri440,ri670,ri870,ri1020 ! Oleg's refractive index real g440,g670,g870,g1020 ! Oleg's asymetry parameter (total) c--- aeronet diagnostic real sphericity,sky_error,sun_error c--- MFRSR measurements integer imfrsr,n_meas ! number of mfrsr measurements for a given AERONET inversion real mfrsr_dump(35) ! mfrsr_dump with added Brewer NO2 and AERONET NO2 and O3 (col 31) plus water (col 32) & nighttime bias real mfrsr_aot(3), *mfrsr_tot(3),mfrsr_dif(3),mfrsr_dirn,mfrsr_v0raw real mfrsr_dd(3),mfrsr_dt(3) real meas_par c123456789012345678901234567890123456789012345678901234567890123456789012 real mfrsr_tr,mfrsr_sza,mfrsr_ozone,md12_pressure real aeronet_no2 real mfrsr_no2_atmcm,mfrsr_errno2DU,mfrsr_dtno2, * mfrsr_tno2,mfrsr_toz,mfrsr_ldir real mfrsr_hour,mfrsr_min,mfrsr_sec character*256 aeronethead character*580 header character*90 fin,fout character*90 str(24),buf(8) character*15 fdata character*25 buff,faerosol character*37 buftom character*4 ozprof character*7 aers character*1 surfs,ozlat character*14 environm parameter (nstr=18,nozpr=26) ! number of dump rows in arizona output files data vza/0.0,15.0,30.0,45.0,60.0,75.0/ data phi/0.0,30.0,60.0,90.0,120.0,150.0,180.0/ data szab/0.0,30.0,45.0,60.0,70.0,77.0,81.0,84.0,86.0,88.0/ c-- V7 325 ozone & Temp profile c data ozoneb/16.0,14.0,26.0,45.0,74.7,66.9,41.7,24.5,11.1,3.7,1.4/ c data temp/273.0,239.0,219.1,216.6,217.0,220.8,227.6,239.4, c &253.6,263.9,262.6/ c-- V7 275 ozone & Temp profile data ozoneb/16.0,12.0,15.0,29.0,58.0,63.7,40.6,24.5,11.1,3.7,1.4/ ! default TOMS V7 clim ozone profile c data ozoneb/265.0,1.,1.,1., 1.,1.,1., 1.,1.,1.,1. / data temp/273.0,239.0,217.1,212.2,214.9,220.4,227.6,239.4, ! default TOMS V7 clim temp profile &253.6,263.9,262.6/ jheader=1 c -- open output file open(28,file='./conv527.dd440_bolivia2.txt', *status='unknown') open (25,file='DATA/out527.dd440_bolivia2.txt', *status='unknown') c23456789012345678901234567890123456789012345678901234567890123456789012 header='oleg_dd,oleg_month,oleg_yyyy,oleg_hh,oleg_mm,oleg_ss,oleg_! V2 AERONET retrieval *dayofyear,aot1020m,aot870m,aot670m,aot500m,aot440m,aot380m,aot340m *,oleg_aot440,oleg_aot670,oleg_aot870,oleg_aot1020,oleg_w440,oleg_w *670,oleg_w870,oleg_w1020,rr440,rr670,rr870,rr1020,ri440,ri670,ri87 *0,ri1020,vf,vc, rvf,rvc, lnsf,lnsc,hh,mm,sec,sza,ozone,no2,surfpr, *lam,n_real,k_real,aot_mfrsr,ssa_mfrsr,aot_aer,tau_ray,tau_ozn,lam_ *dir,trans_tot,trans_dif,trans_dirn,v0_raw, zenith_ray, zenith_rati *o,e_dir_ray,e_dir_ratio,e_dif_ray,e_dif_ratio,calc_par,meas_par,ni *t,errno2DU,no2_delta_time,asymetry,jacobian_k_dd, jacobian_ssa_dk, *jacobian_asym_dk,water,sphericity,sky_error,sun_error,nt_bias' write(25,*) header C----------set up orbit and environment models iaertype=4 ! only bi-modal log-normal PSD accepted with AERONET data open(17,file='./matchfile527.440_bolivia2.txt', c open(17,file='./matchfile527.440_greece.txt', ! 440: MFRSR527 & form='formatted',status='old',err=900) do iaeronet=1,iskip read(17,"(1a256)",ERR=901) aeronethead enddo c write(6,*) aeronethead c write(6,*) 'Reading AERONET-MFRSR retrievals' ! reading header 100 continue ! MAIN LOOP THOUGH ALL AERONET-MFRSRS DATA write(*,*) write(*,*) write(*,*) write(6,*) 'start reading new AERONET retrieval' read(17,*,err=901,end=210) aeronet_dump oleg_dd=aeronet_dump(1) oleg_month=aeronet_dump(2) oleg_yyyy=aeronet_dump(3) oleg_hh=aeronet_dump(4) oleg_mm=aeronet_dump(5) oleg_ss=aeronet_dump(6) oleg_dayofyear=aeronet_dump(7) C aot1020m =aeronet_dump(9) aot870m =aeronet_dump(10) aot670m =aeronet_dump(11) aot500m =aeronet_dump(17) aot440m =aeronet_dump(20) aot380m =aeronet_dump(22) aot340m =aeronet_dump(23) oleg_aot440 =aeronet_dump(25) oleg_aot670 =aeronet_dump(26) oleg_aot870 =aeronet_dump(27) oleg_aot1020 =aeronet_dump(28) oleg_w440 =aeronet_dump(38) oleg_w670 =aeronet_dump(39) oleg_w870 =aeronet_dump(40) oleg_w1020 =aeronet_dump(41) rr440 =aeronet_dump(47) rr670 =aeronet_dump(48) rr870 =aeronet_dump(49) rr1020 =aeronet_dump(50) ri440 =aeronet_dump(51) ri670 =aeronet_dump(52) ri870 =aeronet_dump(53) ri1020 =aeronet_dump(54) c reading phase function asymmetry parameter now g440 =aeronet_dump(55) g670 =aeronet_dump(56) g870 =aeronet_dump(57) g1020 =aeronet_dump(58) c Parameters of fitted bi-modal log-normal size distribution vf =aeronet_dump(94) vc =aeronet_dump(98) rvf =aeronet_dump(96) rvc =aeronet_dump(100) lnsf =aeronet_dump(97) lnsc =aeronet_dump(101) call get_aer_mf_v2(real_refr_const,real_refr_value, *gsfc_psd_clim,gsfc_refr_clim,ref_440, *rr440,rr670,ri440,ri670, * vf,vc,rvf,rvc,lnsf,lnsc, * aot440m, *lam,refrr1,refri1,refrr2,refri2, *sigma1,rmean1,Nc1,sigma2,rmean2,Nc2) write(*,*) 'Nc1 Nc2',Nc1,Nc2 write(*,*) 'sigma1,rmean1,sigma2,rmean2',sigma1, *rmean1,sigma2,rmean2 write(*,*) 'lam , real, imaginary index' write(*,'(5f10.4)') (lam(i),refrr1(i),refri1(i),refrr2(i), &refri2(i),i=1,12) c--- Surface albedo model salb440v2 = aeronet_dump(151) ! AERONET V2 albedo ( as function of SZA) isurfmodel=6 ! this default surface albedp is modified using AERONET V2 albedo ( as function of SZA) call surfmodel(salb440v2,isurfmodel,lamalb,surfalb) ! too low surf albedo at GSFC write(*,*) 'lamalb, spec albedo' write(*,'(2f8.3)') (lamalb(i),surfalb(i),i=1,12) c--- Read Aeronet diagnostic sky_error = aeronet_dump(129) sun_error = aeronet_dump(130) sphericity = aeronet_dump(133) c--- mfrsr measuremets read(17,*) n_meas ! get number of MFRSR measurements for a given almucantar retrieval write(6,*) n_meas do imfrsr=1,n_meas ! main cycle of MFRSR retrievals read(17,*,err=902) mfrsr_dump c write(6,*) mfrsr_dump c Aeronet water [cm] water_aeronet = mfrsr_dump(34) c nighttime voltage bias [mV] f_vbias = mfrsr_dump(35) c Aerosol optical Thickness mfrsr_aot(1) =mfrsr_dump(21) ! Using MFRSR AOT for selected wavelength from match file if(cimelaot) then ! Using interpolated CIMEL AOT instead of MFRSR AOT write(6,*) 'Using CIMEL AOT !!!' c write(28,*) 'Using CIMEL AOT !!!' mfrsr_aot(1) =mfrsr_dump(18) ! Using interpolated AERONET AOT for selected channel wavelength c mfrsr_aot(2) =aot340m*1.02410**alpha ! 332nm c mfrsr_aot(3) =aot340m*1.04620**alpha ! 325nm endif c if (mfrsr_aot(1) .lt. 0.0) then c write(6,*) 'front_mfrsr: NEGATIVE MFRSR AOT' c write(6,*) 'Skip to the next measurement:' c goto 113 c endif pi=4.0*atan(1.0) dt=pi/180.0 mfrsr_hour =mfrsr_dump(5) ! hour of the 3 min averaging interval mfrsr_min =mfrsr_dump(6) ! minute of 3 min averaging interval mfrsr_sec =mfrsr_dump(7) ! sec of 3 min averaging interval mfrsr_sza =mfrsr_dump(13) ! SZA at the end of 3 min averaging interval, [deg] mfrsr_v0raw =mfrsr_dump(11) ! V0_raw at selected channel c--- Calculate MFRSR total and diffuse transmittances mfrsr_tot(1) =mfrsr_dump(8)/mfrsr_v0raw ! USDA Total Transmittance mfrsr_dif(1) =mfrsr_dump(9)/mfrsr_v0raw ! Diffuse Transmittance mfrsr_dirn =mfrsr_dump(10)/mfrsr_v0raw ! Direc_normal transmittance if(mdif_corr) then ! Using corrected MFRSR diffuse write(28,*) 'Diff transmittance correction: tdif*', * mdif_corr_cos mfrsr_dif(1)=mdif_corr_cos*mfrsr_dif(1) ! increase diffuse Voltage mfrsr_tot(1) =mfrsr_dirn*cos(dt*mfrsr_sza)+mfrsr_dif(1) ! re-calculate total transmittance mfrsr_dif correction !!! endif c Effective monochromatic wavelength and spectral constants for the channel (former coe.dat and gas2.coe) mfrsr_ldir =mfrsr_dump(25) ! mfrsr_ldir is radiatively effective wavelength (direct irradiance) lambda(1) = 10.*mfrsr_ldir ! convert to Angstroms dynamic mfrsr_ldir from mfrsr_dump(25) write(6,*) 'Effective channel wavelength: L_dir=',mfrsr_ldir c Default gaseous absorption coefficients SELECT CASE (INT(mfrsr_ldir)) CASE (495:505) ! 497.0 channel FWHM=10nm c0b=0.027037199 ! 10nm slit avergaed Daumont C0 temp coefficient [atm-cm-1] for ~497 center wavelength c1b=-2.9762499e-06 c2b=-2.3297379e-07 xno2_ap=6.2315e+00 ! NO2 x-section at 497nm [atm_cm-1] rectangular slit FWHM=10nm, T=294K rhob= 2.8459e-02 ! Rayleigh depolarization factor at 497nm CASE (438:442) ! 440.4 channel FWHM=10nm c0b=0.0036 ! 10nm slit avergaed Daumont C0 temp coefficient [atm-cm-1] for ~440 center wavelength c1b=2.4264951e-06 c2b=9.3288679e-09 xno2_ap=1.32200+01 ! NO2 x-section at 440nm [atm_cm-1] rectangular slit FWHM=10nm, T=294K rhob= 2.9096e-02 ! Rayleigh depolarization factor at 440nm CASE (410:420) ! 415 channel FWHM=10nm c0b=0.00081212989 ! 10nm slit avergaed Daumont C0 temp coefficient [atm-cm-1] for ~415 center wavelength c1b=7.1728601e-07 c2b=-2.5112200e-08 xno2_ap=1.5656e+01 ! NO2 x-section at 415nm [atm_cm-1] rectangular slit FWHM=10nm, T=294K rhob= 2.9364e-02 ! Rayleigh depolarization factor at 415nm CASE (378:382) ! 378.8 channel FWHM=2nm c0b=0.0001 ! 2nm slit avergaed Daumont C0 temp coefficient [atm-cm-1] for ~380 center wavelength c1b=2.0891280e-06 c2b=9.6004180e-09 xno2_ap=1.4891e+01 ! NO2 x-section at 380nm [atm_cm-1] rectangular slit FWHM=2nm, T=294K rhob=2.9514e-02 ! Rayleigh depolarization factor at 380nm CASE (364:370) ! 368 channel FWHM=2nm c0b=0.00026483981 ! 2nm Daumont C0 temp coefficient [atm-cm-1] for ~368 center channel wavelength c1b=5.9200480e-06 ! c2b=5.1175269e-08 xno2_ap=1.4765e+01 ! NO2 x-section at 368nm [atm_cm-1] rectangular slit FWHM=2nm, T=294K rhob= 3.0100e-02 ! Rayleigh depolarization factor at 368nm CASE (338:342) ! 339.1 channel FWHM=2nm c0b=0.032 ! 2nm Daumont C0 temp coefficient [atm-cm-1] for ~340 center wavelength c1b=0.00032550268 c2b=2.1126916e-06 xno2_ap=1.0058e+01 ! NO2 x-section at 340nm [atm_cm-1] rectangular slit FWHM=2nm, T=294K rhob=3.0588e-02 ! Rayleigh depolarization factor at 340nm CASE (330:334) ! 332.6 channel FWHM=2nm c0b=0.099 ! 2nm Daumont C0 temp coefficient [atm-cm-1] for ~332nmc enter channel wavelength c1b=0.00072243119 ! c2b=4.0252022e-06 xno2_ap=8.6540e+00 ! NO2 x-section at 332nm [atm_cm-1] rectangular slit FWHM=2nm, T=294K rhob= 3.0600e-02 ! Rayleigh depolarization factor at 332nm CASE (323:327) ! 325.7 channel FWHM=2nm c0b=0.30299431 ! 2nm Daumont C0 temp coefficient [atm-cm-1] for ~325nm center channel wavelength c1b=0.00150827425 ! c2b=9.3629095e-06 xno2_ap=7.6300e+00 ! NO2 x-section at 325nm [atm_cm-1] rectangular slit FWHM=2nm, T=294K rhob= 3.0897e-02 ! Rayleigh depolarization factor at 325nm CASE (315:319) ! 317 channel FWHM=2nm c0b=1.0087953 ! 2nm Daumont C0 temp coefficient [atm-cm-1] for 317nm center channel wavelength c1b=0.0038389232 ! c2b=2.6983636e-05 xno2_ap=6.0148e+00 ! NO2 x-section at 317nm [atm_cm-1] rectangular slit FWHM=2nm, T=294K rhob= 3.1347e-02 ! Rayleigh depolarization factor at 317nm CASE (309:313) ! 311 channel FWHM=2nm c0b=2.13 ! 2nm Daumont C0 temp coefficient [atm-cm-1] for 311nm center channel wavelength c1b=0.0065817842 ! c2b=4.4819110e-05 xno2_ap=5.2945e+00 ! NO2 x-section at 311nm [atm_cm-1] rectangular slit FWHM=2nm, T=294K rhob= 3.1648e-02 ! Rayleigh depolarization factor at 311nm CASE (303:307) ! 305 channel FWHM=2nm c0b=5.04 ! 2nm Daumont C0 temp coefficient [atm-cm-1] for ~305nm center channel wavelength c1b=0.012160761 ! c2b=9.8356421e-05 xno2_ap=4.2412e+00 ! NO2 x-section at 305nm [atm_cm-1] rectangular slit FWHM=2nm, T=294K rhob= 3.1700e-02 ! Rayleigh depolarization factor at 305nm CASE (298:302) ! 300 channel FWHM=2nm c0b=10.354233 ! 2nm B&P C0 temp coefficient [atm-cm-1] for given center channel wavelength c1b=0.016812323 ! c2b=8.2238028e-05 xno2_ap=3.4158e+00 ! NO2 x-section at 300nm [atm_cm-1] rectangular slit FWHM=2nm, T=294K rhob= 3.1718e-02 ! Rayleigh depolarization factor at 300nm CASE DEFAULT write(*,*) 'mfrsr_ldir is outside UV-MFRSR spectral channels: ', * mfrsr_ldir stop END SELECT write(6,*) 'Default Ozone abs for T=228K',c0b +c1b*(-45.0) + * c2b*(-45.0)*(-45.0) c measured O3 optical thickness and effective dynamic x-section mfrsr_ozone =mfrsr_dump(16) ! True measured column ozone,[DU]: Brewer or PANDORA fill in with OMI data or AERONET clim: !! no fill in values allowed !! mfrsr_toz =mfrsr_dump(23) ! tau_ozone subtraction used in calibration based on true measured ozone and T=228K c measured NO2 optical thickness mfrsr_tno2 =mfrsr_dump(31) ! tau_no2 is measured NO2 optical thickness determined during calibration c re-define c0b,c1b=0,c2b=0,rhob here and pass them to gs_fit routine c Define effective ozone x-sec (c0b) that will produce exactly the same tau_ozone (mfrsr_toz) used during calibration of the channel c0b=1000.*mfrsr_toz/mfrsr_ozone ! effective dynamic ozone x-section [atm-cm-1] for given channel and measurement ! using the same effective temperature and true ozone during calibration of this channel c1b=0.0000 ! neglect Temperature dependence of ozone absorption c2b=0.000 ! neglect Temperature dependence of ozone absorption write(6,*) 'Effective ozone abs coeff [atm-cm-1]',c0b c Define NO2 optical thickness c mfrsr_no2_atmcm and xno2_used are passed to gs_fit to add as 2nd (water) absorber mfrsr_no2_atmcm =0.001*mfrsr_dump(27) ! Measured NO2 total amount in [DU] and convert to[atm*cm] aeronet_no2 =0.001*mfrsr_dump(30) ! AERONET assumed NO2 climatology (SCIAMACHY) in [DU] converted to [atm-cm] if (mfrsr_no2_atmcm .le. 0.0) then ! If measured NO2 has fill value: -9.999 revert to AERONET climatology mfrsr_no2_atmcm=aeronet_no2 ! treat fill values as AERONET NO2 climatology c mfrsr_no2_atmcm=0.001*0.14 ! AERONET clim NO2 value [atm*cm] at Santa Cruz, Bolivia, Sep 14 endif mfrsr_errno2DU =0.001*mfrsr_dump(28) ! NO2 expanded uncertainty [DU] converted to [atm*cm] mfrsr_dtno2 =mfrsr_dump(29) ! minimal time difference between NO2 and MFRSR measurements [min] c tau_no2= mfrsr_no2_atmcm*xno2 ! NO2 effective optical thickness (rectangular, FWHM=2nm, T=294K ) c Version 2 AERONET AOT is already corrected for climatological tau NO2 c !!! correction for the difference between measured true and AERONET NO2 is done before extrapolation during calibration step write(6,*) 'Default NO2 abs coeff for T=294K [atm-cm]',xno2_ap xno2_used = mfrsr_tno2 / mfrsr_no2_atmcm write(6,*) 'Used NO2 abs coeff [atm-cm-1]',xno2_used if (mfrsr_aot(1) .lt. 0.0) then write(6,*) 'mfrsr_aot < 0: ',mfrsr_aot(1) write(6,*) 'Skip to the next measurement:' goto 113 endif c---- skip MFRSR retrieval if oleg_aot440 < min_aot440 if (oleg_aot440 .lt. min_aot440) then write (6,*) 'oleg_aot440=',oleg_aot440,' < min_aot440 = ' write(6,*) 'Skip to the next measurement:' goto 113 endif if (mfrsr_sza .gt. 70.0) then write(6,*) mfrsr_hour,mfrsr_min,mfrsr_sec,mfrsr_sza write(6,*) 'front_mfrsr: mfrsr_sza>70 deg' write(6,*) 'Skip to the next MFRSR measurement:' goto 113 ! skip to the next MFRSR record endif ! end check SZA>70 c--- V0 correction: V0_corr=V0*k at 368nm => dtau=myu*ln(k), ln(k) from CIMEL comparisons if(mtau_corr) then ! Using corrected MFRSR AOT write(6,*) 'Using corrected mfrsr 368 AOT: dtau=myu*ln(k), * ln(k)=',mean_lnk write(28,*) 'Using corrected mfrsr 368 AOT: dtau=myu*ln(k), * ln(k)=',mean_lnk mfrsr_aot(1) =mfrsr_aot(1) +mean_lnk*cos(pi*mfrsr_sza/180.) endif c--- Calculation of fitting parameter do i=1,1 if (abs(mfrsr_dirn ) .lt. 0.0001) then mfrsr_dd(i)= - 9999.999 mfrsr_dt(i)= - 9999.999 mfrsr_tot(i) = - 9999.999 mfrsr_dif(i) = - 9999.999 else mfrsr_dd(i) = mfrsr_dif(i)/mfrsr_dirn ! normalize to dir_normal mfrsr_dt(i) = mfrsr_dif(i)/mfrsr_tot(i) ! normalize to total_hor endif enddo c--- Calculate MFRSR total and diffuse transmittances c mfrsr_v0raw =mfrsr_dump(11) ! V0_raw at 368nm c mfrsr_tot(1) =mfrsr_dump(8)/mfrsr_v0raw ! Total transmittance c mfrsr_dif(1) =mfrsr_dump(9)/mfrsr_v0raw ! Diffuse transmittance c mfrsr_dirn =mfrsr_dump(10)/mfrsr_v0raw ! direct normal transmittance if(mie_only) go to 102 ! no irradiance calculation c--- vertical structure model nlayer=100 c aerwidthh=0.5 ! Gaussian LAyer StDev aerwidthh=3.0 ! Box profile for Bolivia [nk- 2015/08/12] c aermaxh=3.0 ! Layer height above surface [km] (Santa Cruz, Bolivia) aermaxh=0.0 ! Bolivia changed to box profile aeronet_aotchan =mfrsr_dump(18) ! interpolated AERONET AOT to the UV-MFRSR channel md12_pressure =mfrsr_dump(15) ! Pressure measurement [mbars] - not used as input to Arizona code mfrsr_tr =mfrsr_dump(22) ! mfrsr_tr is rayleigh tau for the AERONET surface pressure(md12_pressure) and l_dir, used for determining V0 during calibration c mfrsr_tr =mfrsr_dump(22)*md12_pressure/1013.25 !! No additional scaling !!! mfrsr_tr is already scaled for the actual surface pressure and l_dir c!!! Keep surface pressure to 1 atm as Arizona code input, so not to trucate Rayleigh and aerosol optical thicknesses surfpres=1013.25 c--- ozone vert profile model iozpr=0 !iozpr=0; if iozpr=1 then enhanced trop. ozone totoz= 0.001*mfrsr_ozone ! Arizona input total ozone [atmcm] or NO2 (substitution) do i=1,11 ozone(i)=ozoneb(i) ! ozone profile in DU/Umkehr enddo if(iozpr.eq.1) then c totoz=totoz+0.5*ozone(1)/1000.0 ozone(1)=1.5*ozone(1) endif write(ozprof,'(i3,''M'')') int(totoz*1000.0) write(*,*) 'ozone profile ',ozprof c call ozonemodel(ozprof,ozone,temp) write(*,*) 'umkehr layer, ozone, temperature' write(*,'(i2,2f9.2)') (i,ozone(i),temp(i),i=1,11) C----------call a subrotine interpolating the TOMRAD ozone and temperature C----------profiles call profiles(1,surfpres,ozone,temp) C----------call a subroutine adding an aerosol profile and gas2 profile to be used for C----------Gauss-Seidel computations gas2tot=mfrsr_no2_atmcm ! NO2 column amount [atm-cm] call aerprof(aermaxh,aerwidthh,gas2maxh,gas2width,gas2tot,nlayer) c--- Cloud parameters tauc=0.0 ! NO CLOUDS C-----------end the orbit and environmental models 102 continue ! skip here if Mie only calculations c do jj=1,10 ! loop through SZA do jj=1,1 c sza=szab(jj) sza=mfrsr_sza ! MFRSR SZA tau440=aeronet_dump(20) ! AERONET measurement at 441nm write(*,*) 'tau440=',tau440 call gs_fit(firstmult_k,xno2_used,c0b,c1b,c2b,rhob, *max_k_iter,dd_fit,dt_fit,t_fit,d_fit, *fit_tolerance,log_fit,mie_only, *totoz,gas2tot,surfpres,mfrsr_tr, ! gas2tot is currently = mfrsr_no2_atmcm ( + so2 in the future ) &mfrsr_dd,mfrsr_dt,mfrsr_tot,mfrsr_dif,mfrsr_aot, &tau440,tauc,sza,vza,phi,nvza,nphi, &iaermod,iaertype,lam,refrr1,refri1,refrr2,refri2, &sigma1,rmean1,Nc1,sigma2,rmean2,Nc2,isurfmodel,lamalb,surfalb, &icount,lambda,ratio,taum,omega,tau_real,om_real,asym_real, & rad_rayl, & n_real,k_real,e_dir_ray,e_dir_ratio,e_dif_ray,e_dif_ratio, *dd_new,n_iter,err_dd,meas_par,max_k_real,jacobian_k_dd, *jacobian_ssa_dk,jacobian_asym_dk ) call system ('clean') do i=1,icount ! icount < 50 if (e_dir_ray(i) .eq. 0.0) go to 100 e_dir(i)=e_dir_ray(i)*e_dir_ratio(i) e_dif(i)=e_dif_ray(i)*e_dif_ratio(i) if (e_dir(i) .lt. 0.00001 ) then dd_calc(i)= -9999.999 else dd_calc(i)=e_dif(i)/e_dir(i) endif enddo if(mie_only) go to 101 C------------open radiance output file if ( no_radiance_output ) go to 101 if(jheader.eq.1) then jheader=2 if(iaermod.lt.10) then write(aers,'(a,i1,''t'',f4.2)') '0',iaermod,tau440 else write(aers,'(i2,''t'',f4.2)') iaermod,tau440 endif write(surfs,'(i1)') int(aermaxh) environm=ozprof//'_'//aers//'_'//surfs write(fout,'(''DATA/faer_en'',a,''.dat'')') environm WRITE(*,*) 'output filename',fout open(11,file=fout,status='unknown') C-------------write a header write(11,'(a)') 'ENVIRONMENTAL VARIABLES' write(11,'(a,a)') 'TOMS OZONE & TEMPERATURE PROFILES ',ozprof write(11,'(a)') 'umkehr layer, ozone, temperature' write(11,'(i2,2f9.2)') (i,ozone(i),temp(i),i=1,11) write(11,'(a,i2,a,f6.3)') 'AEROSOL MODEL',iaermod,' tau440=', &tau440 write(11,'(a,f6.2,a,f6.2,a,e8.2,a,f6.2,a,f6.2,a,e8.2)') 'log-normal & distr: rmean1, ',rmean1,' sigma1 ',sigma1,' Nc1 ',Nc1, &' rmean2 ',rmean2,' sigma2 ',sigma2,' Nc2 ',Nc2 write(11,'(a)') 'wavelength (mu) real, imager, optical & thickness, omega0' write(11,'(f8.3,4e12.3,2f12.3)') (lam(i),refrr1(i),refri1(i), &refrr2(i),refri2(i),taum(i),omega(i),i=1,12) write(11,'(f8.3,2f14.5)') (lambda(i), &tau_real(i),om_real(i),i=1,12) write(11,'(a,f5.1,a,f5.1)') 'PROFILE: height=',aermaxh,' scale=', &aerwidthh write(11,'(a,f8.1)') 'SURFACE PRESSURE, mb ',md12_pressure write(11,'(a,i2)') 'SURFACE MODEL ',isurfmodel write(11,'(a)') ' wavelength, mu albedo' write(11,'(f8.3,f12.3)') (lamalb(i),surfalb(i),i=1,12) write(11,*) endif do j=1,nvza do k=1,nphi write(11,'(a)') 'GEOMETRY' write(11,'(a,f7.2,a,f7.2,a,f9.2)') ' sza =',sza,' vza =',vza(j), &' phi =',phi(k) write(11,'(a)') 'WAVE (A), sky_rad_rayl aer/ray ratio' C--------------write data do i=1,icount ! icount < 50 write(11,'(3e13.5)') lambda(i),rad_rayl(i,j,k),ratio(i,j,k) enddo enddo enddo enddo ! end sza loop 101 continue ! skip radiance calculations c--- write flux and MFRSR results c write(6,*) 'DEBUG!!!: icount,mfrsr_tot(i),mfrsr_dif(i),mfrsr_dirn c *,mfrsr_v0raw' c write(6,*)icount,mfrsr_tot(1),mfrsr_dif(1),mfrsr_dirn,mfrsr_v0raw write(28,81) &oleg_dd,oleg_month,oleg_yyyy, ! UT day #, month # &oleg_hh,oleg_mm,oleg_ss,oleg_dayofyear, ! UT hh, min,sec, decimal doy # *int(mfrsr_hour),int(mfrsr_min),int(mfrsr_sec) *,mfrsr_sza,mfrsr_ozone,1000.*mfrsr_no2_atmcm,md12_pressure *,(0.1*lambda(i),n_real(i),k_real(i), * tau_real(i),om_real(i),aeronet_aotchan, * mfrsr_tr,mfrsr_toz,mfrsr_ldir, * mfrsr_tot(i),mfrsr_dif(i),mfrsr_dirn,mfrsr_v0raw, ! mfrsr measurements of tot/V0, dif/V0 and dir_norm/V0 components *rad_rayl(i,1,1),ratio(i,1,1),e_dir_ray(i),e_dir_ratio(i), *e_dif_ray(i),e_dif_ratio(i), *dd_new(i),meas_par,n_iter(i),i=1,icount), *mfrsr_errno2DU,mfrsr_dtno2,asym_real(1), *jacobian_k_dd,jacobian_ssa_dk,jacobian_asym_dk write(25,8) &oleg_dd,oleg_month,oleg_yyyy, ! UT day #, month # &oleg_hh,oleg_mm,oleg_ss,oleg_dayofyear, ! UT hh, min,sec, decimal doy # &aot1020m,aot870m,aot670m,aot500m,aot440m,aot380m,aot340m, ! measured aot &oleg_aot440,oleg_aot670,oleg_aot870, &oleg_aot1020, ! OLeg aot &oleg_w440,oleg_w670,oleg_w870,oleg_w1020, ! OLeg W &rr440,rr670,rr870,rr1020, ! OLeg Real refr &ri440,ri670,ri870,ri1020, ! OLeg Im refr & vf,vc, rvf,rvc, lnsf,lnsc ! PSD: fine and coarse mode *,int(mfrsr_hour),int(mfrsr_min),int(mfrsr_sec) *,mfrsr_sza,mfrsr_ozone,1000.*mfrsr_no2_atmcm, * md12_pressure *,(0.1*lambda(i),n_real(i),k_real(i), * tau_real(i),om_real(i),aeronet_aotchan, * mfrsr_tr,mfrsr_toz,mfrsr_ldir, * mfrsr_tot(i),mfrsr_dif(i),mfrsr_dirn,mfrsr_v0raw, ! mfrsr measurements of tot/V0, dif/V0 and dir_norm/V0 components *rad_rayl(i,1,1),ratio(i,1,1),e_dir_ray(i),e_dir_ratio(i), *e_dif_ray(i),e_dif_ratio(i), *dd_new(i),meas_par,n_iter(i),i=1,icount), *mfrsr_errno2DU,mfrsr_dtno2,asym_real(1), *jacobian_k_dd,jacobian_ssa_dk,jacobian_asym_dk, *water_aeronet,sphericity,sky_error,sun_error,f_vbias c23456789012345678901234567890123456789012345678901234567890123456789012 8 format( & 2(1x,I2), !UT day #, UT month # & 1x, I4, ! UT year # & 3(1x, I2), ! UT hour # & 1x, f10.6, ! decimal UT day of year & 7(1x,f12.6),2x, ! aeronet measured aot @ 1020, 870nm,670nm,500nm,440nm,380nm,340nm c & 2(1x,f5.3),2x, ! OLeg aot @ [340, 380,] & 4(1x,f5.3),2x, ! OLeg aot @ 440, 670, 870 and 1020 & 4(1x,f5.3),2x, ! OLeg W @ 1020, 870nm, 670,440 & 4(1x,f8.6),2x, ! OLeg n @ 1020, 870nm, 670,440 & 4(1x,f8.6),2x, ! OLeg k @ 1020, 870nm, 670,440 & 6(1x,f9.6),2x, ! Log-Norm PSD parameters & 3(i2,1x), ! int(mfrsr_hour),int(mfrsr_min),int(mfrsr_sec) & f5.2,1x,f5.1,1x,f5.3,1x,f7.2,2x, ! mfrsr_sza, mfrsr_ozone, 1000.*mfrsr_no2_atmcm, surfpres c lam (n,k) (mt,w,aer_t368) & f5.1,1x,f7.4,1x,f9.6,1x,3(f7.3,1x),2x,2(f8.6,1x),f8.4,2x, c t_tot,t_dif,t_dir V0 & 3(f7.5,1x),1(f9.2,1x),3x, ! mfrsr_tot(i),mfrsr_dif(i),mfrsr_dirn,mfrsr_v0raw & 6(e13.5,1x), 2(f9.4,1x),3x,i2,7(1x,f13.6),2x,f8.4) 81 format( & 2(1x,I2), !UT day #, UT month # & 1x, I4, ! UT year # & 3(1x, I2), ! UT hour # & 1x, f10.6,5x, ! decimal UT day of year & 3(i2,1x), ! int(mfrsr_hour),int(mfrsr_min),int(mfrsr_sec) & f5.2,1x,2(f6.3,1x),f7.2,2x, ! mfrsr_sza, mfrsa_ozone, 1000.*mfrsr_no2_atmcm, surfpres & f6.1,1x,2(f9.6,1x),3(f6.3,1x),2x,2(f8.6,1x),f8.4,2x, & 3(f7.5,1x),1(f9.2,1x),3x, ! mfrsr_tot(i),mfrsr_dif(i),mfrsr_dirn,mfrsr_v0raw & 6(e13.5,1x), 2(f9.4,1x),3x,i2,8(1x,f14.6)) write(6,*) 'mfrsr measurement processed: ' 113 continue ! mfrsr measurement not processed because mfrsr_sza> 78deg enddo ! continue mfrsr data loop for single AERONET inversion write(6,*) 'one aeronet inversion done !' go to 100 ! return to read new AERONET retrieval 210 continue ! end of AERONET file encountered write(6,*) 'end of AERONET file encountered !!' C------------close all files close(17) close(11) close(25) close(28) stop c-- error handling 900 continue write(6,*) 'front_end: could not open merged AERONET/MFRSR file ' stop 901 continue write(6,*) 'front_end: error while reading data in & AERONET file ' stop 902 continue write(6,*) 'front_end: error while reading data in & MFRSR file ' stop end