PROGRAM mie IMPLICIT NONE C C-------Mie computations for a size distribution C of homogeneous spheres. c C========================================================== C--Ref : Toon and Ackerman, Applied Optics, 1981 C Stephens, CSIRO, 1979 C Attention : surdimensionement des tableaux C to be compiled with double precision option (-r8 on Sun) C AUTHOR: Olivier Boucher, Christoph Kleinschmitt C-------SIZE distribution properties---------------- C--sigma_g : geometric standard deviation C--r_0 : geometric number mean radius (um)/modal radius C--Ntot : total concentration in m-3 C--rho : dry density in kg/m3 c--mmd=2.5 um from Yves c INTEGER Ndis, dis PARAMETER (Ndis=1) REAL sigma_g(Ndis), r_0(Ndis), Ntot(Ndis), rho DATA r_0 /0.277E-6/ DATA sigma_g/2.0/ DATA Ntot /1.0/ PARAMETER (rho=2.650E3) !--dry density kg/m3 CHARACTER*1 :: run CHARACTER(*), PARAMETER :: mainfolder="./" CHARACTER(*), PARAMETER :: version="./" CHARACTER(*), PARAMETER :: output="./" CHARACTER(*), PARAMETER :: source="./" c REAL masse,volume,surface c REAL rmin, rmax !----integral bounds in m PARAMETER (rmin=0.002E-6,rmax=100.E-6) c c------------------------------------- c COMPLEX m !----refractive index m=n_r-i*n_i INTEGER Nmax,Nstart !--number of iterations COMPLEX k2, k3, z1, z2 COMPLEX u1,u5,u6,u8 COMPLEX a(1:21000), b(1:21000) COMPLEX I INTEGER n !--loop index REAL pi, nnn COMPLEX nn REAL Q_ext, Q_abs, Q_sca, g, omega !--parameters for radius r REAL x !--size parameter REAL r !--radius REAL sigma_sca, sigma_ext, sigma_abs REAL omegatot, gtot !--averaged parameters COMPLEX ksiz2(-1:21000), psiz2(1:21000) COMPLEX nu1z1(1:21010), nu1z2(1:21010) COMPLEX nu3z2(0:21000) REAL number, deltar INTEGER bin, Nbin, k PARAMETER (Nbin=941) c INTEGER nb_lambda_dust_lw PARAMETER (nb_lambda_dust_lw=601) c C---wavelengths STREAMER INTEGER Nwv, NwvmaxSW PARAMETER (NwvmaxSW=24) REAL lambda(1:NwvmaxSW+1) DATA lambda/0.28E-6, 0.30E-6, 0.33E-6, 0.36E-6, 0.40E-6, . 0.44E-6, 0.48E-6, 0.52E-6, 0.57E-6, 0.64E-6, . 0.69E-6, 0.75E-6, 0.78E-6, 0.87E-6, 1.00E-6, . 1.10E-6, 1.19E-6, 1.28E-6, 1.53E-6, 1.64E-6, . 2.13E-6, 2.38E-6, 2.91E-6, 3.42E-6, 4.00E-6/ c c---wavelengths de references dans le SW INTEGER nb, nb_lambda_ref PARAMETER (nb_lambda_ref=5) REAL lambda_ref(nb_lambda_ref) DATA lambda_ref /0.443E-6,0.550E-6,0.670E-6,0.765E-6,0.865E-6/ c c--LW c--Tb=representative blackbody temperature of aerosol (in K) INTEGER NwvmaxLW PARAMETER (NwvmaxLW=500) REAL Tb, hh, cc, kb PARAMETER (Tb=270.0, hh=6.62607e-34) PARAMETER (cc=2.99792e8, kb=1.38065e-23) c C---TOA fluxes - Streamer Cs REAL weight(1:NwvmaxSW), weightLW c DATA weight/0.839920E1, 0.231208E2, 0.322393E2, 0.465058E2, c . 0.678199E2, 0.798964E2, 0.771359E2, 0.888472E2, c . 0.115281E3, 0.727565E2, 0.816992E2, 0.336172E2, c . 0.914603E2, 0.112706E3, 0.658840E2, 0.524470E2, c . 0.391067E2, 0.883864E2, 0.276672E2, 0.681812E2, c . 0.190966E2, 0.250766E2, 0.128704E2, 0.698720E1/ C---TOA fluxes - Tad DATA weight/ 4.20, 11.56, 16.12, 23.25, 33.91, 39.95, 38.57, . 44.42, 57.64, 29.36, 47.87, 16.81, 45.74, 56.35, . 32.94, 26.22, 19.55, 44.19, 13.83, 34.09, 9.55, . 12.54, 6.44, 3.49/ C---BOA fluxes - Tad c DATA weight/ 0.01, 4.05, 9.51, 15.99, 26.07, 33.10, 33.07, c . 39.91, 52.67, 27.89, 43.60, 13.67, 42.22, 40.12, c . 32.70, 14.44, 19.48, 14.23, 13.43, 16.42, 8.33, c . 0.95, 0.65, 2.76/ c REAL lambda_int(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW), ll REAL dlambda_int(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW), dl c REAL n_r(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW) REAL n_i(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW) c REAL final_a(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW) REAL final_g(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW) REAL final_w(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW) c INTEGER band, bandSW, bandLW, NbandSW, NbandLW PARAMETER (NbandSW=6, NbandLW=16) c c---wavelengths SW RRTM REAL wv_rrtm_SW(NbandSW+1) DATA wv_rrtm_SW/ 0.185E-6, 0.25E-6, 0.44E-6, 0.69E-6, . 1.19E-6, 2.38E-6, 4.00E-6/ c c---wavenumbers (wn) and wavelengths (wv) LW RRTM REAL wn_rrtm(NbandLW+1), wv_rrtm(NbandLW+1) DATA wn_rrtm/ 10., 250., 500., 630., 700., 820., . 980., 1080., 1180., 1390., 1480., 1800., . 2080., 2250., 2380., 2600., 3000./ c c--GCM results REAL gcm_a(NbandSW+NbandLW) REAL gcm_g(NbandSW+NbandLW) REAL gcm_w(NbandSW+NbandLW) REAL gcm_weight_a(NbandSW+NbandLW) REAL gcm_weight_g(NbandSW+NbandLW) REAL gcm_weight_w(NbandSW+NbandLW) c c--all results REAL ss_a(NbandSW+NbandLW+nb_lambda_ref) REAL ss_w(NbandSW+NbandLW+nb_lambda_ref) REAL ss_g(NbandSW+NbandLW+nb_lambda_ref) c c--index for SW INTEGER wv, nb_wv_r, nb_wv_i PARAMETER (nb_wv_r=78, nb_wv_i=78) REAL wv_r(1:nb_wv_r), index_r(1:nb_wv_r) REAL wv_i(1:nb_wv_i), index_i(1:nb_wv_i) REAL count_n_r, count_n_i c--index for LW REAL wavenumber REAL, DIMENSION(nb_lambda_dust_LW) :: ref_wv REAL, DIMENSION(nb_lambda_dust_LW,4) :: ref_ind_r, ref_ind_i c c---------------preferred choice of Claudia's model c--change setting of imod if you wish INTEGER, PARAMETER :: imin=1, imax=2, imean=3, imedian=4 INTEGER, PARAMETER :: imod=imean INTEGER :: ilambda1, ilambda2 c CHARACTER(*), PARAMETER :: fileplace=mainfolder//version//output CHARACTER(*), PARAMETER :: sourcefile=mainfolder//source CHARACTER*100 :: dummy c c--------------------------------------------------------- WRITE(run,'(i1)') imod c c--set up SW refractive index c DO Nwv=1, NwvmaxSW lambda_int(Nwv)=( lambda(Nwv)+lambda(Nwv+1) ) /2. ENDDO c DO nb=1, nb_lambda_ref lambda_int(NwvmaxSW+nb)=lambda_ref(nb) ENDDO c c--set up LW refractive index c--conversion wavenumber in cm-1 to wavelength in m c--be careful wavelengths are in decreasing order DO Nwv=1, NbandLW+1 wv_rrtm(Nwv)=10000./wn_rrtm(Nwv)*1.e-6 ENDDO c c--spread lambda_int logarithmically in the LW c--still in decreasing order DO Nwv=1, NwvmaxLW lambda_int(NwvmaxSW+nb_lambda_ref+Nwv)= . exp( log(wv_rrtm(1))+float(Nwv)/float(NwvmaxLW+1)* . (log(wv_rrtm(NbandLW+1))-log(wv_rrtm(1))) ) ENDDO c--computing the dlamdas Nwv=1 dlambda_int(NwvmaxSW+nb_lambda_ref+Nwv)= . lambda_int(NwvmaxSW+nb_lambda_ref+Nwv)- . lambda_int(NwvmaxSW+nb_lambda_ref+Nwv+1) DO Nwv=2, NwvmaxLW-1 dlambda_int(NwvmaxSW+nb_lambda_ref+Nwv)= . (lambda_int(NwvmaxSW+nb_lambda_ref+Nwv-1)- . lambda_int(NwvmaxSW+nb_lambda_ref+Nwv+1))/2. ENDDO Nwv=NwvmaxLW dlambda_int(NwvmaxSW+nb_lambda_ref+Nwv)= . lambda_int(NwvmaxSW+nb_lambda_ref+Nwv-1)- . lambda_int(NwvmaxSW+nb_lambda_ref+Nwv) c DO Nwv=1, NwvmaxLW c print *,'ll dl=', lambda_int(NwvmaxSW+nb_lambda_ref+Nwv), c . dlambda_int(NwvmaxSW+nb_lambda_ref+Nwv) c ENDDO c c--read Yves hematite's data OPEN(unit=10,file='r_1v5_hematite.dat') DO wv=1, nb_wv_r READ (10,*) wv_r(wv), index_r(wv) ENDDO CLOSE(10) c OPEN(unit=10,file='i_1v5_hematite.dat') DO wv=1, nb_wv_i READ (10,*) wv_i(wv), index_i(wv) ENDDO CLOSE(10) c c--interpolating DO Nwv=1, NwvmaxSW+nb_lambda_ref n_r(Nwv)=0.0 n_i(Nwv)=0.0 ENDDO c c--first the 24 Streamer wavelengths DO Nwv=1, NwvmaxSW c count_n_r=0.0 DO wv=1, nb_wv_r IF (wv_r(wv)/1.e9.GT.lambda(Nwv).AND. . wv_r(wv)/1.e9.LT.lambda(Nwv+1)) THEN n_r(Nwv)=n_r(Nwv)+index_r(wv) count_n_r=count_n_r+1.0 ENDIF ENDDO c IF (count_n_r.GT.0.5) THEN c--on moyenne n_r(Nwv)=n_r(Nwv)/count_n_r ELSE c--sinon plus proche voisin DO wv=1, nb_wv_r IF (wv_r(wv)/1.e9.LT.lambda_int(Nwv)) THEN n_r(Nwv)=index_r(wv) ENDIF ENDDO ENDIF c count_n_i=0.0 DO wv=1, nb_wv_i IF (wv_i(wv)/1.e9.GT.lambda(Nwv).AND. . wv_i(wv)/1.e9.LT.lambda(Nwv+1)) THEN n_i(Nwv)=n_i(Nwv)+index_i(wv) count_n_i=count_n_i+1.0 ENDIF ENDDO c IF (count_n_i.GT.0.5) THEN c--on moyenne n_i(Nwv)=n_i(Nwv)/count_n_i ELSE c--sinon plus proche voisin DO wv=1, nb_wv_i IF (wv_i(wv)/1.e9.LT.lambda_int(Nwv)) THEN n_i(Nwv)=index_i(wv) ENDIF ENDDO ENDIF c ENDDO c c--then the reference wavelengths DO nb=1, nb_lambda_ref c DO wv=1, nb_wv_r IF (wv_r(wv)/1.e9.LT.lambda_ref(nb)) THEN n_r(NwvmaxSW+nb)=index_r(wv) ENDIF ENDDO c DO wv=1, nb_wv_i IF (wv_i(wv)/1.e9.LT.lambda_ref(nb)) THEN n_i(NwvmaxSW+nb)=index_i(wv) ENDIF ENDDO c ENDDO c c--read Claudia's LW refractive index data c--format 4x imaginary parts + 4 real parts OPEN (unit=11,file=sourcefile//'Refractive_Index_Dust_LW.txt') READ(11,*) dummy READ(11,*) dummy DO nb=1,nb_lambda_dust_LW READ(11,*) ref_wv(nb),ref_ind_i(nb,:), . ref_wv(nb),ref_ind_r(nb,:) ref_wv(nb)=ref_wv(nb)*1.e-6 !--convert in m ENDDO CLOSE(11) c c--extrapolate LW refractive index data DO Nwv=NwvmaxSW+nb_lambda_ref+1, NwvmaxSW+nb_lambda_ref+NwvmaxLW ! for lambda larger than largest value, take largest value IF (lambda_int(Nwv).GE.ref_wv(nb_lambda_dust_LW)) THEN n_r(Nwv)=ref_ind_r(nb_lambda_dust_LW,imod) n_i(Nwv)=ref_ind_i(nb_lambda_dust_LW,imod) ELSEIF (lambda_int(Nwv).LE.ref_wv(1)) THEN n_r(Nwv)=ref_ind_r(1,imod) n_i(Nwv)=ref_ind_i(1,imod) ELSE DO nb=1,nb_lambda_dust_LW IF (ref_wv(nb).LT.lambda_int(Nwv)) ilambda1=nb ENDDO DO nb=nb_lambda_dust_LW,1,-1 IF (ref_wv(nb).GT.lambda_int(Nwv)) ilambda2=nb ENDDO n_r(Nwv)=ref_ind_r(ilambda1,imod)+ . (lambda_int(Nwv)-ref_wv(ilambda1))/ . (ref_wv(ilambda2)-ref_wv(ilambda1))* . (ref_ind_r(ilambda2,imod)-ref_ind_r(ilambda1,imod)) n_i(Nwv)=ref_ind_i(ilambda1,imod)+ . (lambda_int(Nwv)-ref_wv(ilambda1))/ . (ref_wv(ilambda2)-ref_wv(ilambda1))* . (ref_ind_i(ilambda2,imod)-ref_ind_i(ilambda1,imod)) ENDIF ENDDO c c--save extrapolated refr index OPEN (unit=11,file=sourcefile//'Refractive_Index_Dust_SW_ex.txt') DO Nwv=1, NwvmaxSW+nb_lambda_ref WRITE(11,*) lambda_int(Nwv), n_r(Nwv), n_i(Nwv) ENDDO CLOSE(11) c OPEN (unit=11,file=sourcefile//'Refractive_Index_Dust_LW_ex.txt') DO Nwv=NwvmaxSW+nb_lambda_ref+1, NwvmaxSW+nb_lambda_ref+NwvmaxLW WRITE(11,*) lambda_int(Nwv), n_r(Nwv), n_i(Nwv) ENDDO CLOSE(11) c c---Loop on wavelengths DO Nwv=1, NwvmaxSW+nb_lambda_ref+NwvmaxLW c m=CMPLX(n_r(Nwv),-n_i(Nwv)) c pi=4.*ATAN(1.) c I=CMPLX(0.,1.) c sigma_sca=0.0 sigma_ext=0.0 sigma_abs=0.0 gtot=0.0 omegatot=0.0 masse = 0.0 volume=0.0 surface=0.0 c DO bin=0, Nbin !---loop on size bins r=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin))) ! PRINT *, 'r =', r x=2.*pi*r/lambda_int(Nwv) deltar=1./FLOAT(Nbin)*(log(rmax)-log(rmin)) c number=0 DO dis=1, Ndis ! PRINT *, 'Ntot(dis) =', Ntot(dis) ! PRINT *, 'r_0(dis) =', r_0(dis) ! PRINT *, 'sigma_g(dis) =', sigma_g(dis) number=number+Ntot(dis)/SQRT(2.*pi)/log(sigma_g(dis))* . exp(-0.5*(log(r/r_0(dis))/log(sigma_g(dis)))**2) masse=masse +4./3.*pi*(r**3)*number* . deltar*rho*1.E3 !--g/m3 volume=volume+4./3.*pi*(r**3)*number*deltar surface=surface+4.*pi*r**2*number*deltar ENDDO ! PRINT *, 'number =', number c k2=m k3=CMPLX(1.0,0.0) z2=CMPLX(x,0.0) z1=m*z2 IF (0.0.LE.x.AND.x.LE.8.) THEN Nmax=INT(x+4*x**(1./3.)+1.)+2 ELSEIF (8..LT.x.AND.x.LT.4200.) THEN Nmax=INT(x+4.05*x**(1./3.)+2.)+1 ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN Nmax=INT(x+4*x**(1./3.)+2.)+1 ELSE PRINT *, 'x out of bound, x=', x STOP ENDIF Nstart=Nmax+10 C-----------loop for nu1z1, nu1z2 nu1z1(Nstart)=CMPLX(0.0,0.0) nu1z2(Nstart)=CMPLX(0.0,0.0) DO n=Nstart-1, 1 , -1 nn=CMPLX(FLOAT(n),0.0) nu1z1(n)=(nn+1.)/z1 - 1./( (nn+1.)/z1 + nu1z1(n+1) ) nu1z2(n)=(nn+1.)/z2 - 1./( (nn+1.)/z2 + nu1z2(n+1) ) ENDDO C------------loop for nu3z2 nu3z2(0)=-I DO n=1, Nmax nn=CMPLX(FLOAT(n),0.0) nu3z2(n)=-nn/z2 + 1./ (nn/z2 - nu3z2(n-1) ) ENDDO C-----------loop for psiz2 and ksiz2 (z2) ksiz2(-1)=COS(REAL(z2))-I*SIN(REAL(z2)) ksiz2(0)=SIN(REAL(z2))+I*COS(REAL(z2)) DO n=1,Nmax nn=CMPLX(FLOAT(n),0.0) ksiz2(n)=(2.*nn-1.)/z2 * ksiz2(n-1) - ksiz2(n-2) psiz2(n)=CMPLX(REAL(ksiz2(n)),0.0) ENDDO C-----------loop for a(n) and b(n) DO n=1, Nmax u1=k3*nu1z1(n) - k2*nu1z2(n) u5=k3*nu1z1(n) - k2*nu3z2(n) u6=k2*nu1z1(n) - k3*nu1z2(n) u8=k2*nu1z1(n) - k3*nu3z2(n) a(n)=psiz2(n)/ksiz2(n) * u1/u5 b(n)=psiz2(n)/ksiz2(n) * u6/u8 ENDDO C-----------------final loop-------------- Q_ext=0.0 Q_sca=0.0 g=0.0 DO n=Nmax-1,1,-1 nnn=FLOAT(n) Q_ext=Q_ext+ (2.*nnn+1.) * REAL( a(n)+b(n) ) Q_sca=Q_sca+ (2.*nnn+1.) * . REAL( a(n)*CONJG(a(n)) + b(n)*CONJG(b(n)) ) g=g + nnn*(nnn+2.)/(nnn+1.) * . REAL( a(n)*CONJG(a(n+1))+b(n)*CONJG(b(n+1)) ) + . (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n))) ENDDO Q_ext=2./x**2 * Q_ext Q_sca=2./x**2 * Q_sca Q_abs=Q_ext-Q_sca IF (AIMAG(m).EQ.0.0) Q_abs=0.0 omega=Q_sca/Q_ext g=g*4./x**2/Q_sca c sigma_sca=sigma_sca+r**2*Q_sca*number*deltar sigma_abs=sigma_abs+r**2*Q_abs*number*deltar sigma_ext=sigma_ext+r**2*Q_ext*number*deltar omegatot=omegatot+r**2*Q_ext*omega*number*deltar gtot =gtot+r**2*Q_sca*g*number*deltar c ENDDO !---bin C------------------------------------------------------------------ sigma_sca=pi*sigma_sca sigma_abs=pi*sigma_abs sigma_ext=pi*sigma_ext gtot=pi*gtot/sigma_sca omegatot=pi*omegatot/sigma_ext c final_g(Nwv)=gtot final_w(Nwv)=omegatot final_a(Nwv)=sigma_ext/masse c ENDDO !--loop on wavelength c c---averaging over LMDZ wavebands c DO band=1, NbandSW+NbandLW gcm_a(band)=0.0 gcm_g(band)=0.0 gcm_w(band)=0.0 gcm_weight_a(band)=0.0 gcm_weight_g(band)=0.0 gcm_weight_w(band)=0.0 ENDDO c c---band 1 is now in the UV, so we take the first wavelength as being representative DO Nwv=1,1 bandSW=1 gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv) gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv) gcm_w(bandSW)=gcm_w(bandSW)+ . final_w(Nwv)*final_a(Nwv)*weight(Nwv) gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+ . final_a(Nwv)*weight(Nwv) gcm_g(bandSW)=gcm_g(bandSW)+ . final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv) gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+ . final_a(Nwv)*final_w(Nwv)*weight(Nwv) ENDDO c DO Nwv=1,NwvmaxSW c IF (lambda_int(Nwv).LE.wv_rrtm_SW(3)) THEN !--RRTM spectral interval 2 bandSW=2 ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(4)) THEN !--RRTM spectral interval 3 bandSW=3 ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(5)) THEN !--RRTM spectral interval 4 bandSW=4 ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(6)) THEN !--RRTM spectral interval 5 bandSW=5 ELSE !--RRTM spectral interval 6 bandSW=6 ENDIF c gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv) gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv) gcm_w(bandSW)=gcm_w(bandSW)+ . final_w(Nwv)*final_a(Nwv)*weight(Nwv) gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+ . final_a(Nwv)*weight(Nwv) gcm_g(bandSW)=gcm_g(bandSW)+ . final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv) gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+ . final_a(Nwv)*final_w(Nwv)*weight(Nwv) ENDDO c DO Nwv=1,NwvmaxLW ll=lambda_int(NwvmaxSW+nb_lambda_ref+Nwv) dl=dlambda_int(NwvmaxSW+nb_lambda_ref+Nwv) weightLW=1./ll**5./(exp(hh*cc/kb/Tb/ll)-1.)*dl DO band=1, NbandLW IF (wv_rrtm(band+1).LE.ll.AND.ll.LE.wv_rrtm(band)) THEN bandLW=band ENDIF ENDDO c print *,'Nwv =',Nwv,lambda_int(NwvmaxSW+nb_lambda_ref+Nwv), c . wv_rrtm(bandLW),wv_rrtm(bandLW+1),bandLW c--extinction gcm_a(NbandSW+bandLW)=gcm_a(NbandSW+bandLW)+ . final_a(NwvmaxSW+nb_lambda_ref+Nwv)*weightLW gcm_weight_a(NbandSW+bandLW)= . gcm_weight_a(NbandSW+bandLW)+weightLW c--single scattering albedo gcm_w(NbandSW+bandLW)=gcm_w(NbandSW+bandLW)+ . final_w(NwvmaxSW+nb_lambda_ref+Nwv)* . final_a(NwvmaxSW+nb_lambda_ref+Nwv)*weightLW gcm_weight_w(NbandSW+bandLW)=gcm_weight_w(NbandSW+bandLW)+ . final_a(NwvmaxSW+nb_lambda_ref+Nwv)*weightLW c--asymetry parameter gcm_g(NbandSW+bandLW)=gcm_g(NbandSW+bandLW)+ . final_g(NwvmaxSW+nb_lambda_ref+Nwv)* . final_a(NwvmaxSW+nb_lambda_ref+Nwv)* . final_w(NwvmaxSW+nb_lambda_ref+Nwv)*weightLW gcm_weight_g(NbandSW+bandLW)=gcm_weight_g(NbandSW+bandLW)+ . final_a(NwvmaxSW+nb_lambda_ref+Nwv)* . final_w(NwvmaxSW+nb_lambda_ref+Nwv)*weightLW ENDDO c DO band=1, NbandSW gcm_a(band)=gcm_a(band)/gcm_weight_a(band) gcm_w(band)=gcm_w(band)/gcm_weight_w(band) gcm_g(band)=gcm_g(band)/gcm_weight_g(band) ss_a(band)=gcm_a(band) ss_w(band)=gcm_w(band) ss_g(band)=gcm_g(band) ENDDO c DO band=NbandSW+1, NbandSW+NbandLW gcm_a(band)=gcm_a(band)/gcm_weight_a(band) gcm_w(band)=gcm_w(band)/gcm_weight_w(band) gcm_g(band)=gcm_g(band)/gcm_weight_g(band) ss_a(band)=gcm_a(band) ss_w(band)=gcm_w(band) ss_g(band)=gcm_g(band) ENDDO c DO nb=1, nb_lambda_ref ss_a(NbandSW+NbandLW+nb)=final_a(NwvmaxSW+nb) ss_w(NbandSW+NbandLW+nb)=final_w(NwvmaxSW+nb) ss_g(NbandSW+NbandLW+nb)=final_g(NwvmaxSW+nb) ENDDO c-------------------------------------------------------------- c--saving output data c-------------------------------------------------------------- c OPEN (unit=14, file=fileplace// c . 'aer_optical_properties_streamer_imod'//run//'.txt') c c DO nwv=1,NwvmaxSW c WRITE(14,*)lambda_int(Nwv), c . final_a(nwv),final_w(nwv),final_g(nwv) c ENDDO c c CLOSE(14) c-------------------------------------------------------------- OPEN (unit=14, file=fileplace//'SEXT_dust_6bands.txt') WRITE(14,*) ' ! Dust insoluble' WRITE(14,952) (ss_a(k),k=1,NbandSW) CLOSE(14) OPEN (unit=14, file=fileplace//'SSA_dust_6bands.txt') WRITE(14,*) ' ! Dust insoluble' WRITE(14,952) (ss_w(k),k=1,NbandSW) CLOSE(14) OPEN (unit=14, file=fileplace//'G_dust_6bands.txt') WRITE(14,*) ' ! Dust insoluble' WRITE(14,952) (ss_g(k),k=1,NbandSW) CLOSE(14) 952 FORMAT(1X,6(F6.3,','),' &') c-------------------------------------------------------------- OPEN (unit=14, file=fileplace//'SEXT_dust_5wave.txt') WRITE(14,*) ' ! extinction Dust insoluble' WRITE(14,953) (ss_a(k), . k=NbandSW+NbandLW+1,NbandSW+NbandLW+nb_lambda_ref) CLOSE(14) OPEN (unit=14, file=fileplace//'SABS_dust_5wave.txt') WRITE(14,*) ' ! absorption Dust insoluble' WRITE(14,953) ((1.0-ss_w(k))*ss_a(k), . k=NbandSW+NbandLW+1,NbandSW+NbandLW+nb_lambda_ref) CLOSE(14) 953 FORMAT(1X,5(F6.3,','),' &') c-------------------------------------------------------------- c OPEN (unit=14, file=fileplace// c . 'aer_optical_properties_LW_16bands_imod'//run//'.txt') c c DO k=NbandSW+1,NbandSW+NbandLW c WRITE(14,*)wv_rrtm(k-NbandSW),wv_rrtm(k-NbandSW+1), c . ss_a(k),ss_w(k),ss_g(k) c ENDDO c c CLOSE(14) c-------------------------------------------------------------- c-LW alpha for RRtm: only absorption, no scattering c decreasing wavelength or increasing wavenumber OPEN (unit=14, file=fileplace// . 'SABS_dust_16bands_imod'//run//'.txt') WRITE(14,*) ' ! Dust insoluble' WRITE(14,954) (ss_a(NbandSW+k)*(1.-ss_w(NbandSW+k)), . k=1,NbandLW/2) WRITE(14,955) (ss_a(NbandSW+k)*(1.-ss_w(NbandSW+k)), . k=NbandLW/2+1,NbandLW) CLOSE(14) 954 FORMAT(1X,8(F6.3,','),' &') 955 FORMAT(1X,7(F6.3,','),1(F6.3),' /') c c--------------------------------------------------------------- END