source: TOOLS/CMIP6_FORCING/AER_OPTICS/mie_6bands_SW_LW_dust.f @ 3385

Last change on this file since 3385 was 3385, checked in by oboucher, 4 years ago

Uploading package to prepare tropospheric aerosol optics LUT in LMDz fo rCMIP6

File size: 23.2 KB
Line 
1        PROGRAM mie
2        IMPLICIT NONE
3C
4C-------Mie computations for a size distribution 
5C       of homogeneous spheres.
6c
7C==========================================================
8C--Ref : Toon and Ackerman, Applied Optics, 1981
9C        Stephens, CSIRO, 1979
10C Attention : surdimensionement des tableaux
11C to be compiled with double precision option (-r8 on Sun)
12C AUTHOR: Olivier Boucher, Christoph Kleinschmitt
13C-------SIZE distribution properties----------------
14C--sigma_g : geometric standard deviation 
15C--r_0     : geometric number mean radius (um)/modal radius
16C--Ntot    : total concentration in m-3 
17C--rho     : dry density in kg/m3
18c--mmd=2.5 um from Yves
19c
20        INTEGER Ndis, dis
21        PARAMETER (Ndis=1)
22        REAL sigma_g(Ndis), r_0(Ndis), Ntot(Ndis), rho
23        DATA r_0    /0.277E-6/ 
24        DATA sigma_g/2.0/ 
25        DATA Ntot   /1.0/
26        PARAMETER (rho=2.650E3)  !--dry density kg/m3
27
28        CHARACTER*1  :: run
29        CHARACTER(*), PARAMETER :: mainfolder="./"
30        CHARACTER(*), PARAMETER :: version="./"
31        CHARACTER(*), PARAMETER :: output="./" 
32        CHARACTER(*), PARAMETER :: source="./"
33c
34        REAL masse,volume,surface
35c
36        REAL rmin, rmax    !----integral bounds in  m
37        PARAMETER (rmin=0.002E-6,rmax=100.E-6)
38c
39c-------------------------------------
40c
41        COMPLEX m          !----refractive index m=n_r-i*n_i
42        INTEGER Nmax,Nstart !--number of iterations
43        COMPLEX k2, k3, z1, z2
44        COMPLEX u1,u5,u6,u8
45        COMPLEX a(1:21000), b(1:21000)
46        COMPLEX I 
47        INTEGER n  !--loop index
48        REAL pi, nnn
49        COMPLEX nn
50        REAL Q_ext, Q_abs, Q_sca, g, omega   !--parameters for radius r
51        REAL x  !--size parameter
52        REAL r  !--radius
53        REAL sigma_sca, sigma_ext, sigma_abs
54        REAL omegatot,  gtot !--averaged parameters
55        COMPLEX ksiz2(-1:21000), psiz2(1:21000)
56        COMPLEX nu1z1(1:21010), nu1z2(1:21010)
57        COMPLEX nu3z2(0:21000)
58        REAL number, deltar
59        INTEGER bin, Nbin, k
60        PARAMETER (Nbin=941)
61c
62        INTEGER nb_lambda_dust_lw
63        PARAMETER (nb_lambda_dust_lw=601)
64c
65C---wavelengths STREAMER
66        INTEGER Nwv, NwvmaxSW
67        PARAMETER (NwvmaxSW=24)
68        REAL lambda(1:NwvmaxSW+1)
69        DATA lambda/0.28E-6, 0.30E-6, 0.33E-6, 0.36E-6, 0.40E-6,
70     .              0.44E-6, 0.48E-6, 0.52E-6, 0.57E-6, 0.64E-6,
71     .              0.69E-6, 0.75E-6, 0.78E-6, 0.87E-6, 1.00E-6,
72     .              1.10E-6, 1.19E-6, 1.28E-6, 1.53E-6, 1.64E-6,
73     .              2.13E-6, 2.38E-6, 2.91E-6, 3.42E-6, 4.00E-6/
74c
75c---wavelengths de references dans le SW
76        INTEGER nb, nb_lambda_ref
77        PARAMETER (nb_lambda_ref=5)
78        REAL lambda_ref(nb_lambda_ref)
79        DATA lambda_ref /0.443E-6,0.550E-6,0.670E-6,0.765E-6,0.865E-6/ 
80c
81c--LW 
82c--Tb=representative blackbody temperature of aerosol (in K)
83        INTEGER NwvmaxLW
84        PARAMETER (NwvmaxLW=500)
85        REAL Tb, hh, cc, kb
86        PARAMETER (Tb=270.0, hh=6.62607e-34)
87        PARAMETER (cc=2.99792e8, kb=1.38065e-23)
88c
89C---TOA fluxes - Streamer Cs
90        REAL weight(1:NwvmaxSW), weightLW
91c        DATA weight/0.839920E1, 0.231208E2, 0.322393E2, 0.465058E2,
92c     .              0.678199E2, 0.798964E2, 0.771359E2, 0.888472E2,
93c     .              0.115281E3, 0.727565E2, 0.816992E2, 0.336172E2,
94c     .              0.914603E2, 0.112706E3, 0.658840E2, 0.524470E2,
95c     .              0.391067E2, 0.883864E2, 0.276672E2, 0.681812E2,
96c     .              0.190966E2, 0.250766E2, 0.128704E2, 0.698720E1/
97C---TOA fluxes - Tad
98        DATA weight/ 4.20, 11.56, 16.12, 23.25, 33.91, 39.95, 38.57,
99     .              44.42, 57.64, 29.36, 47.87, 16.81, 45.74, 56.35,
100     .              32.94, 26.22, 19.55, 44.19, 13.83, 34.09,  9.55,
101     .              12.54,  6.44,  3.49/
102C---BOA fluxes - Tad
103c        DATA weight/ 0.01,  4.05, 9.51,  15.99, 26.07, 33.10, 33.07,
104c     .              39.91, 52.67, 27.89, 43.60, 13.67, 42.22, 40.12,
105c     .              32.70, 14.44, 19.48, 14.23, 13.43, 16.42,  8.33,
106c     .               0.95,  0.65, 2.76/
107c
108        REAL lambda_int(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW), ll
109        REAL dlambda_int(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW), dl
110c
111        REAL n_r(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW)
112        REAL n_i(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW)
113c
114        REAL final_a(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW)
115        REAL final_g(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW)
116        REAL final_w(1:NwvmaxSW+nb_lambda_ref+NwvmaxLW)
117c
118        INTEGER band, bandSW, bandLW, NbandSW, NbandLW
119        PARAMETER (NbandSW=6, NbandLW=16)
120c
121c---wavelengths SW RRTM
122        REAL wv_rrtm_SW(NbandSW+1)
123        DATA wv_rrtm_SW/  0.185E-6, 0.25E-6, 0.44E-6, 0.69E-6,
124     .                     1.19E-6, 2.38E-6, 4.00E-6/
125c
126c---wavenumbers (wn) and wavelengths (wv) LW RRTM
127        REAL wn_rrtm(NbandLW+1), wv_rrtm(NbandLW+1)
128        DATA wn_rrtm/  10.,  250.,  500.,  630.,  700.,  820., 
129     .                980., 1080., 1180., 1390., 1480., 1800., 
130     .               2080., 2250., 2380., 2600., 3000./
131c
132c--GCM results
133        REAL gcm_a(NbandSW+NbandLW)
134        REAL gcm_g(NbandSW+NbandLW)
135        REAL gcm_w(NbandSW+NbandLW)
136        REAL gcm_weight_a(NbandSW+NbandLW)
137        REAL gcm_weight_g(NbandSW+NbandLW)
138        REAL gcm_weight_w(NbandSW+NbandLW)
139c
140c--all results
141        REAL ss_a(NbandSW+NbandLW+nb_lambda_ref)
142        REAL ss_w(NbandSW+NbandLW+nb_lambda_ref)
143        REAL ss_g(NbandSW+NbandLW+nb_lambda_ref)
144c
145c--index for SW 
146        INTEGER wv, nb_wv_r, nb_wv_i
147        PARAMETER (nb_wv_r=78, nb_wv_i=78)
148        REAL wv_r(1:nb_wv_r), index_r(1:nb_wv_r)
149        REAL wv_i(1:nb_wv_i), index_i(1:nb_wv_i)
150        REAL count_n_r, count_n_i
151c--index for LW 
152        REAL wavenumber
153        REAL, DIMENSION(nb_lambda_dust_LW)   :: ref_wv
154        REAL, DIMENSION(nb_lambda_dust_LW,4) :: ref_ind_r, ref_ind_i
155c
156c---------------preferred choice of Claudia's model
157c--change setting of imod if you wish
158        INTEGER, PARAMETER :: imin=1, imax=2, imean=3, imedian=4
159        INTEGER, PARAMETER :: imod=imean
160        INTEGER :: ilambda1, ilambda2
161c
162        CHARACTER(*), PARAMETER :: fileplace=mainfolder//version//output
163        CHARACTER(*), PARAMETER :: sourcefile=mainfolder//source
164        CHARACTER*100 :: dummy
165c
166c---------------------------------------------------------
167        WRITE(run,'(i1)') imod
168c
169c--set up SW refractive index
170c
171        DO Nwv=1, NwvmaxSW
172          lambda_int(Nwv)=( lambda(Nwv)+lambda(Nwv+1) ) /2.
173        ENDDO
174c
175        DO nb=1, nb_lambda_ref
176          lambda_int(NwvmaxSW+nb)=lambda_ref(nb)
177        ENDDO
178c
179c--set up LW refractive index
180c--conversion wavenumber in cm-1 to wavelength in m
181c--be careful wavelengths are in decreasing order
182        DO Nwv=1, NbandLW+1
183          wv_rrtm(Nwv)=10000./wn_rrtm(Nwv)*1.e-6
184        ENDDO
185c
186c--spread lambda_int logarithmically in the LW
187c--still in decreasing order
188        DO Nwv=1, NwvmaxLW
189          lambda_int(NwvmaxSW+nb_lambda_ref+Nwv)=
190     .    exp( log(wv_rrtm(1))+float(Nwv)/float(NwvmaxLW+1)*
191     .    (log(wv_rrtm(NbandLW+1))-log(wv_rrtm(1))) )
192        ENDDO
193c--computing the dlamdas
194        Nwv=1
195        dlambda_int(NwvmaxSW+nb_lambda_ref+Nwv)=
196     .        lambda_int(NwvmaxSW+nb_lambda_ref+Nwv)-
197     .        lambda_int(NwvmaxSW+nb_lambda_ref+Nwv+1)
198        DO Nwv=2, NwvmaxLW-1
199        dlambda_int(NwvmaxSW+nb_lambda_ref+Nwv)=
200     .        (lambda_int(NwvmaxSW+nb_lambda_ref+Nwv-1)-
201     .         lambda_int(NwvmaxSW+nb_lambda_ref+Nwv+1))/2.
202        ENDDO
203        Nwv=NwvmaxLW
204        dlambda_int(NwvmaxSW+nb_lambda_ref+Nwv)=
205     .        lambda_int(NwvmaxSW+nb_lambda_ref+Nwv-1)-
206     .        lambda_int(NwvmaxSW+nb_lambda_ref+Nwv)
207c        DO Nwv=1, NwvmaxLW
208c          print *,'ll dl=', lambda_int(NwvmaxSW+nb_lambda_ref+Nwv),
209c     .                      dlambda_int(NwvmaxSW+nb_lambda_ref+Nwv)
210c        ENDDO
211c
212c--read Yves hematite's data
213        OPEN(unit=10,file='r_1v5_hematite.dat')
214        DO wv=1, nb_wv_r
215          READ (10,*) wv_r(wv), index_r(wv)
216        ENDDO
217        CLOSE(10)
218c
219        OPEN(unit=10,file='i_1v5_hematite.dat')
220        DO wv=1, nb_wv_i
221          READ (10,*) wv_i(wv), index_i(wv)
222        ENDDO
223        CLOSE(10)
224c
225c--interpolating
226        DO Nwv=1, NwvmaxSW+nb_lambda_ref
227           n_r(Nwv)=0.0
228           n_i(Nwv)=0.0
229        ENDDO
230c
231c--first the 24 Streamer wavelengths
232        DO Nwv=1, NwvmaxSW
233c
234          count_n_r=0.0
235          DO wv=1, nb_wv_r
236            IF (wv_r(wv)/1.e9.GT.lambda(Nwv).AND.
237     .          wv_r(wv)/1.e9.LT.lambda(Nwv+1)) THEN
238              n_r(Nwv)=n_r(Nwv)+index_r(wv)
239              count_n_r=count_n_r+1.0
240            ENDIF
241          ENDDO
242c
243          IF (count_n_r.GT.0.5) THEN
244c--on moyenne
245            n_r(Nwv)=n_r(Nwv)/count_n_r
246          ELSE
247c--sinon plus proche voisin
248          DO wv=1, nb_wv_r
249            IF (wv_r(wv)/1.e9.LT.lambda_int(Nwv)) THEN
250              n_r(Nwv)=index_r(wv)
251            ENDIF
252          ENDDO
253          ENDIF
254c
255          count_n_i=0.0
256          DO wv=1, nb_wv_i
257            IF (wv_i(wv)/1.e9.GT.lambda(Nwv).AND.
258     .          wv_i(wv)/1.e9.LT.lambda(Nwv+1)) THEN
259              n_i(Nwv)=n_i(Nwv)+index_i(wv)
260              count_n_i=count_n_i+1.0
261            ENDIF
262          ENDDO
263c
264          IF (count_n_i.GT.0.5) THEN
265c--on moyenne
266            n_i(Nwv)=n_i(Nwv)/count_n_i
267          ELSE
268c--sinon plus proche voisin
269          DO wv=1, nb_wv_i
270            IF (wv_i(wv)/1.e9.LT.lambda_int(Nwv)) THEN
271              n_i(Nwv)=index_i(wv)
272            ENDIF
273          ENDDO
274          ENDIF
275c
276        ENDDO
277c
278c--then the reference wavelengths
279        DO nb=1, nb_lambda_ref
280c
281          DO wv=1, nb_wv_r
282            IF (wv_r(wv)/1.e9.LT.lambda_ref(nb)) THEN
283              n_r(NwvmaxSW+nb)=index_r(wv)
284            ENDIF
285          ENDDO
286c
287          DO wv=1, nb_wv_i
288            IF (wv_i(wv)/1.e9.LT.lambda_ref(nb)) THEN
289              n_i(NwvmaxSW+nb)=index_i(wv)
290            ENDIF
291          ENDDO
292c
293        ENDDO
294c
295c--read Claudia's LW refractive index data
296c--format 4x imaginary parts + 4 real parts
297        OPEN (unit=11,file=sourcefile//'Refractive_Index_Dust_LW.txt')
298        READ(11,*) dummy
299        READ(11,*) dummy
300        DO nb=1,nb_lambda_dust_LW
301          READ(11,*) ref_wv(nb),ref_ind_i(nb,:),
302     .               ref_wv(nb),ref_ind_r(nb,:)
303          ref_wv(nb)=ref_wv(nb)*1.e-6  !--convert in m
304        ENDDO
305        CLOSE(11)
306c
307c--extrapolate LW refractive index data
308       DO Nwv=NwvmaxSW+nb_lambda_ref+1, NwvmaxSW+nb_lambda_ref+NwvmaxLW
309          ! for lambda larger than largest value, take largest value
310          IF (lambda_int(Nwv).GE.ref_wv(nb_lambda_dust_LW)) THEN
311           n_r(Nwv)=ref_ind_r(nb_lambda_dust_LW,imod)
312           n_i(Nwv)=ref_ind_i(nb_lambda_dust_LW,imod)
313          ELSEIF (lambda_int(Nwv).LE.ref_wv(1)) THEN
314           n_r(Nwv)=ref_ind_r(1,imod)
315           n_i(Nwv)=ref_ind_i(1,imod)
316          ELSE
317            DO nb=1,nb_lambda_dust_LW
318             IF (ref_wv(nb).LT.lambda_int(Nwv)) ilambda1=nb
319            ENDDO
320            DO nb=nb_lambda_dust_LW,1,-1
321             IF (ref_wv(nb).GT.lambda_int(Nwv)) ilambda2=nb
322            ENDDO
323            n_r(Nwv)=ref_ind_r(ilambda1,imod)+
324     .              (lambda_int(Nwv)-ref_wv(ilambda1))/
325     .              (ref_wv(ilambda2)-ref_wv(ilambda1))*
326     .              (ref_ind_r(ilambda2,imod)-ref_ind_r(ilambda1,imod))
327            n_i(Nwv)=ref_ind_i(ilambda1,imod)+
328     .              (lambda_int(Nwv)-ref_wv(ilambda1))/
329     .              (ref_wv(ilambda2)-ref_wv(ilambda1))*
330     .              (ref_ind_i(ilambda2,imod)-ref_ind_i(ilambda1,imod))
331          ENDIF
332       ENDDO
333c
334c--save extrapolated refr index
335       OPEN (unit=11,file=sourcefile//'Refractive_Index_Dust_SW_ex.txt')
336       DO Nwv=1, NwvmaxSW+nb_lambda_ref
337         WRITE(11,*) lambda_int(Nwv), n_r(Nwv), n_i(Nwv)
338       ENDDO
339       CLOSE(11)
340c
341       OPEN (unit=11,file=sourcefile//'Refractive_Index_Dust_LW_ex.txt')
342       DO Nwv=NwvmaxSW+nb_lambda_ref+1, NwvmaxSW+nb_lambda_ref+NwvmaxLW
343         WRITE(11,*) lambda_int(Nwv), n_r(Nwv), n_i(Nwv)
344       ENDDO
345       CLOSE(11)
346c
347c---Loop on wavelengths
348        DO Nwv=1, NwvmaxSW+nb_lambda_ref+NwvmaxLW
349c
350        m=CMPLX(n_r(Nwv),-n_i(Nwv))
351c
352        pi=4.*ATAN(1.)
353c
354        I=CMPLX(0.,1.)
355c
356        sigma_sca=0.0
357        sigma_ext=0.0
358        sigma_abs=0.0
359        gtot=0.0
360        omegatot=0.0
361        masse = 0.0
362        volume=0.0
363        surface=0.0
364c
365        DO bin=0, Nbin !---loop on size bins
366 
367        r=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
368!        PRINT *, 'r =', r
369        x=2.*pi*r/lambda_int(Nwv)
370        deltar=1./FLOAT(Nbin)*(log(rmax)-log(rmin))
371c
372        number=0
373        DO dis=1, Ndis
374!        PRINT *, 'Ntot(dis) =', Ntot(dis)
375!        PRINT *, 'r_0(dis) =', r_0(dis)
376!        PRINT *, 'sigma_g(dis) =', sigma_g(dis)
377          number=number+Ntot(dis)/SQRT(2.*pi)/log(sigma_g(dis))*
378     .           exp(-0.5*(log(r/r_0(dis))/log(sigma_g(dis)))**2)
379        masse=masse  +4./3.*pi*(r**3)*number*
380     .                          deltar*rho*1.E3          !--g/m3
381        volume=volume+4./3.*pi*(r**3)*number*deltar
382        surface=surface+4.*pi*r**2*number*deltar
383        ENDDO
384
385!        PRINT *, 'number =', number
386c
387        k2=m
388        k3=CMPLX(1.0,0.0)
389
390        z2=CMPLX(x,0.0)
391        z1=m*z2
392 
393        IF (0.0.LE.x.AND.x.LE.8.) THEN
394        Nmax=INT(x+4*x**(1./3.)+1.)+2
395        ELSEIF (8..LT.x.AND.x.LT.4200.) THEN
396        Nmax=INT(x+4.05*x**(1./3.)+2.)+1
397        ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN
398        Nmax=INT(x+4*x**(1./3.)+2.)+1
399        ELSE
400        PRINT *, 'x out of bound, x=', x
401        STOP
402        ENDIF
403
404        Nstart=Nmax+10
405
406C-----------loop for nu1z1, nu1z2
407
408       nu1z1(Nstart)=CMPLX(0.0,0.0)
409       nu1z2(Nstart)=CMPLX(0.0,0.0)
410       DO n=Nstart-1, 1 , -1
411       nn=CMPLX(FLOAT(n),0.0)
412       nu1z1(n)=(nn+1.)/z1 - 1./( (nn+1.)/z1 + nu1z1(n+1) )
413       nu1z2(n)=(nn+1.)/z2 - 1./( (nn+1.)/z2 + nu1z2(n+1) )
414       ENDDO
415
416C------------loop for nu3z2
417
418       nu3z2(0)=-I
419       DO n=1, Nmax
420       nn=CMPLX(FLOAT(n),0.0)
421       nu3z2(n)=-nn/z2 + 1./ (nn/z2 - nu3z2(n-1) )
422       ENDDO
423
424C-----------loop for psiz2 and ksiz2 (z2)
425       ksiz2(-1)=COS(REAL(z2))-I*SIN(REAL(z2))
426       ksiz2(0)=SIN(REAL(z2))+I*COS(REAL(z2))
427       DO n=1,Nmax
428       nn=CMPLX(FLOAT(n),0.0)
429       ksiz2(n)=(2.*nn-1.)/z2 * ksiz2(n-1) - ksiz2(n-2)
430       psiz2(n)=CMPLX(REAL(ksiz2(n)),0.0)
431       ENDDO
432
433C-----------loop for a(n) and b(n)
434   
435       DO n=1, Nmax
436        u1=k3*nu1z1(n) - k2*nu1z2(n)
437        u5=k3*nu1z1(n) - k2*nu3z2(n)
438        u6=k2*nu1z1(n) - k3*nu1z2(n)
439        u8=k2*nu1z1(n) - k3*nu3z2(n)
440        a(n)=psiz2(n)/ksiz2(n) * u1/u5
441        b(n)=psiz2(n)/ksiz2(n) * u6/u8
442       ENDDO
443
444C-----------------final loop--------------
445        Q_ext=0.0
446        Q_sca=0.0
447        g=0.0
448        DO n=Nmax-1,1,-1
449          nnn=FLOAT(n)
450          Q_ext=Q_ext+ (2.*nnn+1.) * REAL( a(n)+b(n) )
451          Q_sca=Q_sca+ (2.*nnn+1.) *
452     .               REAL( a(n)*CONJG(a(n)) + b(n)*CONJG(b(n)) )
453          g=g + nnn*(nnn+2.)/(nnn+1.) *
454     .       REAL( a(n)*CONJG(a(n+1))+b(n)*CONJG(b(n+1)) )  +
455     .          (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n)))
456        ENDDO
457        Q_ext=2./x**2 * Q_ext
458        Q_sca=2./x**2 * Q_sca
459        Q_abs=Q_ext-Q_sca
460        IF (AIMAG(m).EQ.0.0) Q_abs=0.0
461        omega=Q_sca/Q_ext
462        g=g*4./x**2/Q_sca
463c
464        sigma_sca=sigma_sca+r**2*Q_sca*number*deltar
465        sigma_abs=sigma_abs+r**2*Q_abs*number*deltar
466        sigma_ext=sigma_ext+r**2*Q_ext*number*deltar
467        omegatot=omegatot+r**2*Q_ext*omega*number*deltar
468        gtot    =gtot+r**2*Q_sca*g*number*deltar
469c
470        ENDDO   !---bin
471C------------------------------------------------------------------
472
473        sigma_sca=pi*sigma_sca
474        sigma_abs=pi*sigma_abs
475        sigma_ext=pi*sigma_ext
476        gtot=pi*gtot/sigma_sca
477        omegatot=pi*omegatot/sigma_ext
478c
479        final_g(Nwv)=gtot
480        final_w(Nwv)=omegatot
481        final_a(Nwv)=sigma_ext/masse
482c
483        ENDDO  !--loop on wavelength
484c
485c---averaging over LMDZ wavebands
486c
487        DO band=1, NbandSW+NbandLW
488          gcm_a(band)=0.0
489          gcm_g(band)=0.0
490          gcm_w(band)=0.0
491          gcm_weight_a(band)=0.0
492          gcm_weight_g(band)=0.0
493          gcm_weight_w(band)=0.0
494        ENDDO
495c
496c---band 1 is now in the UV, so we take the first wavelength as being representative
497       DO Nwv=1,1
498          bandSW=1
499          gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv)
500          gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv)
501          gcm_w(bandSW)=gcm_w(bandSW)+
502     .                final_w(Nwv)*final_a(Nwv)*weight(Nwv)
503          gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+
504     .                       final_a(Nwv)*weight(Nwv)
505          gcm_g(bandSW)=gcm_g(bandSW)+
506     .                final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
507          gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+
508     .                       final_a(Nwv)*final_w(Nwv)*weight(Nwv)
509        ENDDO
510c
511       DO Nwv=1,NwvmaxSW
512c
513        IF (lambda_int(Nwv).LE.wv_rrtm_SW(3)) THEN      !--RRTM spectral interval 2
514          bandSW=2
515        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(4)) THEN  !--RRTM spectral interval 3
516          bandSW=3
517        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(5)) THEN  !--RRTM spectral interval 4
518          bandSW=4
519        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(6)) THEN  !--RRTM spectral interval 5
520          bandSW=5
521        ELSE                                            !--RRTM spectral interval 6
522          bandSW=6
523        ENDIF
524c
525        gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv)
526        gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv)
527        gcm_w(bandSW)=gcm_w(bandSW)+
528     .              final_w(Nwv)*final_a(Nwv)*weight(Nwv)
529        gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+
530     .                     final_a(Nwv)*weight(Nwv)
531        gcm_g(bandSW)=gcm_g(bandSW)+
532     .              final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
533        gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+
534     .                     final_a(Nwv)*final_w(Nwv)*weight(Nwv)
535
536        ENDDO
537c
538        DO Nwv=1,NwvmaxLW
539          ll=lambda_int(NwvmaxSW+nb_lambda_ref+Nwv)
540          dl=dlambda_int(NwvmaxSW+nb_lambda_ref+Nwv)
541          weightLW=1./ll**5./(exp(hh*cc/kb/Tb/ll)-1.)*dl
542          DO band=1, NbandLW
543            IF (wv_rrtm(band+1).LE.ll.AND.ll.LE.wv_rrtm(band)) THEN
544             bandLW=band
545            ENDIF
546          ENDDO
547c          print *,'Nwv =',Nwv,lambda_int(NwvmaxSW+nb_lambda_ref+Nwv),
548c     .                    wv_rrtm(bandLW),wv_rrtm(bandLW+1),bandLW
549c--extinction
550          gcm_a(NbandSW+bandLW)=gcm_a(NbandSW+bandLW)+
551     .                  final_a(NwvmaxSW+nb_lambda_ref+Nwv)*weightLW
552          gcm_weight_a(NbandSW+bandLW)=
553     .                  gcm_weight_a(NbandSW+bandLW)+weightLW
554c--single scattering albedo
555          gcm_w(NbandSW+bandLW)=gcm_w(NbandSW+bandLW)+
556     .         final_w(NwvmaxSW+nb_lambda_ref+Nwv)*
557     .         final_a(NwvmaxSW+nb_lambda_ref+Nwv)*weightLW
558          gcm_weight_w(NbandSW+bandLW)=gcm_weight_w(NbandSW+bandLW)+
559     .         final_a(NwvmaxSW+nb_lambda_ref+Nwv)*weightLW
560c--asymetry parameter
561          gcm_g(NbandSW+bandLW)=gcm_g(NbandSW+bandLW)+
562     .         final_g(NwvmaxSW+nb_lambda_ref+Nwv)*
563     .         final_a(NwvmaxSW+nb_lambda_ref+Nwv)*
564     .         final_w(NwvmaxSW+nb_lambda_ref+Nwv)*weightLW
565          gcm_weight_g(NbandSW+bandLW)=gcm_weight_g(NbandSW+bandLW)+
566     .         final_a(NwvmaxSW+nb_lambda_ref+Nwv)*
567     .         final_w(NwvmaxSW+nb_lambda_ref+Nwv)*weightLW
568        ENDDO
569c
570        DO band=1, NbandSW
571          gcm_a(band)=gcm_a(band)/gcm_weight_a(band)
572          gcm_w(band)=gcm_w(band)/gcm_weight_w(band)
573          gcm_g(band)=gcm_g(band)/gcm_weight_g(band)
574          ss_a(band)=gcm_a(band)
575          ss_w(band)=gcm_w(band)
576          ss_g(band)=gcm_g(band)
577        ENDDO
578c
579        DO band=NbandSW+1, NbandSW+NbandLW
580          gcm_a(band)=gcm_a(band)/gcm_weight_a(band)
581          gcm_w(band)=gcm_w(band)/gcm_weight_w(band)
582          gcm_g(band)=gcm_g(band)/gcm_weight_g(band)
583          ss_a(band)=gcm_a(band)
584          ss_w(band)=gcm_w(band)
585          ss_g(band)=gcm_g(band)
586        ENDDO
587c
588        DO nb=1, nb_lambda_ref
589          ss_a(NbandSW+NbandLW+nb)=final_a(NwvmaxSW+nb)
590          ss_w(NbandSW+NbandLW+nb)=final_w(NwvmaxSW+nb)
591          ss_g(NbandSW+NbandLW+nb)=final_g(NwvmaxSW+nb)
592        ENDDO
593c--------------------------------------------------------------
594c--saving output data
595c--------------------------------------------------------------
596c        OPEN (unit=14, file=fileplace//
597c     .  'aer_optical_properties_streamer_imod'//run//'.txt')
598c
599c        DO nwv=1,NwvmaxSW
600c        WRITE(14,*)lambda_int(Nwv),
601c     .  final_a(nwv),final_w(nwv),final_g(nwv)
602c        ENDDO
603c
604c        CLOSE(14)
605c--------------------------------------------------------------
606        OPEN (unit=14, file=fileplace//'SEXT_dust_6bands.txt')
607        WRITE(14,*) '  ! Dust insoluble'
608        WRITE(14,952) (ss_a(k),k=1,NbandSW) 
609        CLOSE(14)
610        OPEN (unit=14, file=fileplace//'SSA_dust_6bands.txt')
611        WRITE(14,*) '  ! Dust insoluble'
612        WRITE(14,952) (ss_w(k),k=1,NbandSW) 
613        CLOSE(14)
614        OPEN (unit=14, file=fileplace//'G_dust_6bands.txt')
615        WRITE(14,*) '  ! Dust insoluble'
616        WRITE(14,952) (ss_g(k),k=1,NbandSW) 
617        CLOSE(14)
618952     FORMAT(1X,6(F6.3,','),' &')
619c--------------------------------------------------------------
620        OPEN (unit=14, file=fileplace//'SEXT_dust_5wave.txt')
621        WRITE(14,*) '  ! extinction Dust insoluble'
622        WRITE(14,953) (ss_a(k),
623     .      k=NbandSW+NbandLW+1,NbandSW+NbandLW+nb_lambda_ref)
624        CLOSE(14)
625        OPEN (unit=14, file=fileplace//'SABS_dust_5wave.txt')
626        WRITE(14,*) '  ! absorption Dust insoluble'
627        WRITE(14,953) ((1.0-ss_w(k))*ss_a(k),
628     .      k=NbandSW+NbandLW+1,NbandSW+NbandLW+nb_lambda_ref)
629        CLOSE(14)
630953     FORMAT(1X,5(F6.3,','),' &')
631c--------------------------------------------------------------
632c        OPEN (unit=14, file=fileplace//
633c     .  'aer_optical_properties_LW_16bands_imod'//run//'.txt')
634c
635c        DO k=NbandSW+1,NbandSW+NbandLW
636c        WRITE(14,*)wv_rrtm(k-NbandSW),wv_rrtm(k-NbandSW+1),
637c     .             ss_a(k),ss_w(k),ss_g(k)
638c        ENDDO
639c
640c        CLOSE(14)
641c--------------------------------------------------------------
642c-LW alpha for RRtm: only absorption, no scattering 
643c decreasing wavelength or increasing wavenumber
644        OPEN (unit=14, file=fileplace//
645     .    'SABS_dust_16bands_imod'//run//'.txt') 
646        WRITE(14,*) '  ! Dust insoluble'
647        WRITE(14,954) (ss_a(NbandSW+k)*(1.-ss_w(NbandSW+k)),
648     .                  k=1,NbandLW/2)
649        WRITE(14,955) (ss_a(NbandSW+k)*(1.-ss_w(NbandSW+k)),
650     .                  k=NbandLW/2+1,NbandLW)
651        CLOSE(14)
652954     FORMAT(1X,8(F6.3,','),' &')
653955     FORMAT(1X,7(F6.3,','),1(F6.3),' /')
654c
655c---------------------------------------------------------------
656        END
Note: See TracBrowser for help on using the repository browser.