C program: AERONET database........Sagnik Dey (11th August 2004) c modified om Oct 17 2014 to fit koc332 assuming fixed BC=0.021 (670/870/1020) and OC=0.028 (from 368/440) c modified om Sep 24 2014 C modified on 23rd Aug 2005 C C modified on 16th Feb 2009 by Antti Arola C modified on 27th May 2014 by Antti Arola to add MFRSR k measurements c modified on 13th Oct 2014 by NK to retrieve BC and OC fitting only AERONET wavelengths: 440,670,870,1020 c modified on 13th Oct 2014 by NK to retrieve BC and OC fitting only MFRSR wavelengths: 332,368 and 440 C real constants are replaced by double precision Double precision f3from368440,koc332 Double precision f1,f2,f3,c1,c2,c3,df1,df3,mrmix1,mrmix2,c4,c5 Double precision df4,f1_select,f2_select,f3_select,f4_select Double precision mrmix1_select,mrmix2_select,mrmix3_select Double precision mrmix4_select,mrmix5_select,mrmix6_select Double precision mrmix7_select,mrmix8_select,mrmix9_select Double precision fe1_select,fe2_select Double precision mrrtrv5i,mrrtrv6i,mrrtrv7i,mrrtrv8i,mrrtrv9i Double precision mrrtrv10i Double precision mrrtrv2i,mrrtrv3i,mrrtrv4i,mrrtrv1r,mrrtrv2r Double precision mrrtrv3r,mrrtrv4r,mrmix3,xsq1,xsq2,xsq3 Double precision minxsq1,minxsq2,minxsq3,mrrtrv1i,f4,fe1,dfe1 Double precision fe2,dfe2,localminxsq1 Double precision localfe1,localfe2,localf1,localf2,localf3 Double precision diff1,diff2,diff3,diff4 Double precision iEmg1,iEmg2,iEmg3,iEmg4,iEmg5,mrmix4,mrmix5 Double precision iEmg6,iEmg7,iEmg8,iEmg9 Double precision mrmix7,mrmix8,mrmix9 Double precision adiff1,adiff2,adiff3,adiff4,adiffsum,mrmix6 Double precision dreal,dimag,leftvol Double complex em,e1,e2,e3,e4,e5,ans1,ans2,ans3,ans4,ans5 Double complex e6,e7,e8,e9,e10,e11,e12,e13 Double complex Emg1,Emg2,Emg3,ans6,ans7,ans8,ans9,ans10,ans11 Double complex Emg4,Emg5,ans12,ans13,ans14,ans15,ans16,ans17 Double complex ans18,ans19,ans20,ans21,ans22,ans23,ans24,ans25 Double complex ans26,ans27,ans28,ans29,ans30 Double complex Emg6,Emg7,Emg8,Emg9 integer i,j,k,l,m,nf1,nf2,nf3,nf4,nfe1,nfe2,nf2left c open(9,file='real.dat',form='formatted') c Here we initialize constants such as complex dielectric constants : e(real part, imaginary part) C em=water, e1=BC, e2=nss-SO4, e3,e4,e5,...,e10 = OC at 305,311,317,325,332,368,440,670 ,870,1020 c Water em=(1.7689,0.0) c BC c e1=(2.99,1.8) c BC Bond and Bergstrom 2006 e1=(3.1784,3.081) c Ammonium Sulfate e2=(2.3409,0.0) c f1 - initial values for volume fractions in MG inside the loop f1=0.0 f2=0.0 f3=0.0 f4=0.0 fe1=0.0 fe2=0.0 c ####################### determine loops through volume fractions c BC steps nf1=100 df1=0.001 c OC steps nf3=200 df3=0.004 c ammonium sulfate steps nfe2=150 dfe2=0.001 c ####################### nf4=1000 df4=0.001 read(*,*)mrrtrv1i read(*,*)mrrtrv2i read(*,*)mrrtrv3i read(*,*)mrrtrv4i read(*,*)mrrtrv5i read(*,*)mrrtrv6i read(*,*)mrrtrv7i read(*,*)mrrtrv8i read(*,*)mrrtrv9i read(*,*)mrrtrv10i c f1 is fixed BC fraction from AERONET fitting: 0.021 read(*,*)f1 c f3from368440 is fixed OC fraction from previous MFRSR fitting using 368 and 440: 0.028 read(*,*) f3from368440 C#### HERE we start looping to find f3, f1 is fixed ############## c write(*,*) f1 c f1=f1*1.1 !consider uncertainty of BC c write(*,*) f1 minxsq1=1.e7 localminxsq1=1.e7 f3=0.0 do 20 j=1,nf3 fe2=0.0 c do 40 l=1,nfe2 xsq1=0.0 C Calculate mixture dielectric function Emg c e1 - k_BC = const c e2 - k_host=0 c k_OC at 305 fe2=0.2184 dreal=1.53**2-fe2**2 dimag=2.0*1.53*fe2 e3=cmplx(dreal,dimag) c k_OC at 311 fe2=0.21168 dreal=1.53**2-fe2**2 dimag=2.0*1.53*fe2 e4=cmplx(dreal,dimag) c k_OC at 317 fe2=0.20496 dreal=1.53**2-fe2**2 dimag=2.0*1.53*fe2 e5=cmplx(dreal,dimag) c k_OC at 325 fe2=0.196 dreal=1.53**2-fe2**2 dimag=2.0*1.53*fe2 e6=cmplx(dreal,dimag) c k_OC at 332 c koc332=0.18816 - Kirchetter's extrapolation from 350nm to 332nm koc332=0.48000 dreal=1.53**2-koc332**2 dimag=2.0*1.53*koc332 e7=cmplx(dreal,dimag) c k_OC at 368 fe2=0.1478 dreal=1.53**2-fe2**2 dimag=2.0*1.53*fe2 e8=cmplx(dreal,dimag) c k_OC at 440 fe2=0.073 dreal=1.53**2-fe2**2 dimag=2.0*1.53*fe2 e9=cmplx(dreal,dimag) c k_OC at 670 fe2=0.0034 dreal=1.53**2-fe2**2 dimag=2.0*1.53*fe2 e10=cmplx(dreal,dimag) c k_OC at 870 fe2=0.0 dreal=1.53**2-fe2**2 dimag=2.0*1.53*fe2 e11=cmplx(dreal,dimag) c e12=e11: k_OC at 1020 =0 is the same as OC at 870 ans1 = ((e1-em) / (e1+(2.*em))) ans2 = ((e2-em) / (e2+(2.*em))) ans3 = ((e3-em) / (e3+(2.*em))) ans4 = ((e4-em) / (e4+(2.*em))) ans5 = ((e5-em) / (e5+(2.*em))) ans6 = ((e6-em) / (e6+(2.*em))) ans7 = ((e7-em) / (e7+(2.*em))) ans8 = ((e8-em) / (e8+(2.*em))) ans9 = ((e9-em) / (e9+(2.*em))) ans10 = ((e10-em) / (e10+(2.*em))) ans11 = ((e11-em) / (e11+(2.*em))) ans12 = 3 * (f1*ans1 + f2*ans2 + f3*ans3) ans13 = 3 * (f1*ans1 + f2*ans2 + f3*ans4) ans14 = 3 * (f1*ans1 + f2*ans2 + f3*ans5) ans15 = 3 * (f1*ans1 + f2*ans2 + f3*ans6) ans16 = 3 * (f1*ans1 + f2*ans2 + f3*ans7) ans17 = 3 * (f1*ans1 + f2*ans2 + f3*ans8) ans18 = 3 * (f1*ans1 + f2*ans2 + f3*ans9) ans19 = 3 * (f1*ans1 + f2*ans2 + f3*ans10) ans20 = 3 * (f1*ans1 + f2*ans2 + f3*ans11) ans21 = 1-(f1*ans1) - (f2*ans2) - (f3*ans3) ans22 = 1-(f1*ans1) - (f2*ans2) - (f3*ans4) ans23 = 1-(f1*ans1) - (f2*ans2) - (f3*ans5) ans24 = 1-(f1*ans1) - (f2*ans2) - (f3*ans6) ans25 = 1-(f1*ans1) - (f2*ans2) - (f3*ans7) ans26 = 1-(f1*ans1) - (f2*ans2) - (f3*ans8) ans27 = 1-(f1*ans1) - (f2*ans2) - (f3*ans9) ans28 = 1-(f1*ans1) - (f2*ans2) - (f3*ans10) ans29 = 1-(f1*ans1) - (f2*ans2) - (f3*ans11) Emg1 =em *(1 + ans12/ans21) Emg2 =em *(1 + ans13/ans22) Emg3 =em *(1 + ans14/ans23) Emg4 =em *(1 + ans15/ans24) Emg5 =em *(1 + ans16/ans25) Emg6 =em *(1 + ans17/ans26) Emg7 =em *(1 + ans18/ans27) Emg8 =em *(1 + ans19/ans28) Emg9 =em *(1 + ans20/ans29) c3 =sqrt(REAL(Emg1)**2 + AIMAG(Emg1)**2) - REAL(Emg1) iEmg1 = sqrt(c3/2.0) c4 =sqrt(REAL(Emg2)**2 + AIMAG(Emg2)**2) - REAL(Emg2) iEmg2 = sqrt(c4/2.0) c5 =sqrt(REAL(Emg3)**2 + AIMAG(Emg3)**2) - REAL(Emg3) iEmg3 = sqrt(c5/2.0) c6 =sqrt(REAL(Emg4)**2 + AIMAG(Emg4)**2) - REAL(Emg4) iEmg4 = sqrt(c6/2.0) c7 =sqrt(REAL(Emg5)**2 + AIMAG(Emg5)**2) - REAL(Emg5) iEmg5 = sqrt(c7/2.0) c8 =sqrt(REAL(Emg6)**2 + AIMAG(Emg6)**2) - REAL(Emg6) iEmg6 = sqrt(c8/2.0) c9 =sqrt(REAL(Emg7)**2 + AIMAG(Emg7)**2) - REAL(Emg7) iEmg7 = sqrt(c9/2.0) c10 =sqrt(REAL(Emg8)**2 + AIMAG(Emg8)**2) - REAL(Emg8) iEmg8 = sqrt(c10/2.0) c11 =sqrt(REAL(Emg9)**2 + AIMAG(Emg9)**2) - REAL(Emg9) iEmg9 = sqrt(c11/2.0) mrmix1=iEmg1 mrmix2=iEmg2 mrmix3=iEmg3 mrmix4=iEmg4 mrmix5=iEmg5 mrmix6=iEmg6 mrmix7=iEmg7 mrmix8=iEmg8 mrmix9=iEmg9 c Here we fit all MFRSR(1-6) and AERONET (7-10) channels c xsq1=xsq1+(((mrrtrv1i-mrmix1)/mrrtrv1i)**2)+ c & (((mrrtrv2i-mrmix2)/mrrtrv2i)**2)+ c & (((mrrtrv3i-mrmix3)/mrrtrv3i)**2)+ c & (((mrrtrv4i-mrmix4)/mrrtrv4i)**2)+ c & (((mrrtrv5i-mrmix5)/mrrtrv5i)**2)+ c & (((mrrtrv6i-mrmix6)/mrrtrv6i)**2)+ c & (((mrrtrv7i-mrmix7)/mrrtrv7i)**2)+ c & (((mrrtrv8i-mrmix8)/mrrtrv8i)**2)+ c & (((mrrtrv9i-mrmix9)/mrrtrv9i)**2)+ c & (((mrrtrv10i-mrmix9)/mrrtrv10i)**2) cc Here we fit ONLY AERONET (7-10) channels c xsq1=xsq1+ cc & (((mrrtrv1i-mrmix1)/mrrtrv1i)**2)+ cc & (((mrrtrv2i-mrmix2)/mrrtrv2i)**2)+ cc & (((mrrtrv3i-mrmix3)/mrrtrv3i)**2)+ cc & (((mrrtrv4i-mrmix4)/mrrtrv4i)**2)+ cc & (((mrrtrv5i-mrmix5)/mrrtrv5i)**2)+ cc & (((mrrtrv6i-mrmix6)/mrrtrv6i)**2)+ c & (((mrrtrv7i-mrmix7)/mrrtrv7i)**2)+ c & (((mrrtrv8i-mrmix8)/mrrtrv8i)**2)+ c & (((mrrtrv9i-mrmix9)/mrrtrv9i)**2)+ c & (((mrrtrv10i-mrmix9)/mrrtrv10i)**2) c Here we fit ONLY MFRSR 332 channel xsq1=xsq1 c & +(((mrrtrv1i-mrmix1)/mrrtrv1i)**2) c & +(((mrrtrv2i-mrmix2)/mrrtrv2i)**2) c & +(((mrrtrv3i-mrmix3)/mrrtrv3i)**2) c & +(((mrrtrv4i-mrmix4)/mrrtrv4i)**2) & +(((mrrtrv5i-mrmix5)/mrrtrv5i)**2) c & +(((mrrtrv6i-mrmix6)/mrrtrv6i)**2) c & +(((mrrtrv7i-mrmix7)/mrrtrv7i)**2) c & +(((mrrtrv8i-mrmix8)/mrrtrv8i)**2) c & +(((mrrtrv9i-mrmix9)/mrrtrv9i)**2) c & +(((mrrtrv10i-mrmix9)/mrrtrv10i)**2) if (xsq1.lt.minxsq1) then minxsq1=xsq1 f1_select=f1 f2_select=f2 f3_select=f3 fe1_select=fe1 fe2_select=fe2 mrmix1_select=mrmix1 mrmix2_select=mrmix2 mrmix3_select=mrmix3 mrmix4_select=mrmix4 mrmix5_select=mrmix5 mrmix6_select=mrmix6 mrmix7_select=mrmix7 mrmix8_select=mrmix8 mrmix9_select=mrmix9 diff1=((mrrtrv1i-mrmix4)/mrrtrv1i)**2 diff2=((mrrtrv2i-mrmix5)/mrrtrv2i)**2 diff3=((mrrtrv3i-mrmix6)/mrrtrv3i)**2 diff4=((mrrtrv4i-mrmix6)/mrrtrv4i)**2 c write(*,801)xsq1,f1,f3,fe1,fe2,diff1,diff2,diff3,diff4 c 801 format(1p,9e12.4,2x) diff1=mrmix4 diff2=mrmix5 diff3=mrmix6 diff4=mrmix6 endif if (xsq1.lt.localminxsq1) then localminxsq1=xsq1 localf1=f1 localf2=f2 localf3=f3 localfe1=fe1 localfe2=fe2 endif c fe2=dfe2*l c 40 enddo f3=f3+df3 20 enddo C#### HERE we go to f1xsq to find f2 ########### xsq1=minxsq1 f1=f1_select f2=f2_select f3=f3_select fe1=fe1_select fe2=fe2_select f4=f1 write(*,803) f1,f3, & mrrtrv1i,mrrtrv2i,mrrtrv3i,mrrtrv4i,mrrtrv5i, & mrrtrv6i,mrrtrv7i,mrrtrv8i,mrrtrv9i,mrrtrv10i, & mrmix1_select,mrmix2_select,mrmix3_select, & mrmix4_select,mrmix5_select, & mrmix6_select,mrmix7_select,mrmix8_select, & mrmix9_select,mrmix9_select 803 format(1p,22e12.4,2x) 50 continue STOP END