source: TOOLS/CMIP6_FORCING/AER_OPTICS/mie_6bands_sul_bc_for_internal_mixture.f @ 4226

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

Uploading package to prepare tropospheric aerosol optics LUT in LMDz fo rCMIP6
File size: 29.0 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
13C-------SIZE distribution properties----------------
14C--sigma_g : geometric standard deviation 
15C--r_0     : geometric number mean dry radius (um)=modal dry radius
16C--rho     : density in kg/m3
17C--Ntot    : total concentration in m-3 (does not matter)
18c
19        REAL rmin, rmax    !----integral bounds in  m
20        PARAMETER (rmin=0.002E-6, rmax=30.E-6)
21c
22        INTEGER Ndis, dis
23        PARAMETER (Ndis=1)
24        REAL sigma_g(Ndis), r_0(Ndis), Ntot(Ndis)
25        DATA r_0    /0.1E-6/ 
26        DATA sigma_g/1.8/
27        DATA Ntot   /1.0/
28c
29        INTEGER irh,Nrh
30        PARAMETER(Nrh=12)
31        REAL RH_tab(Nrh),RH
32        DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./
33        REAL rh_growth(Nrh), rho_SUL(Nrh), rhgrowth
34c
35        REAL nombre, masse, massebc, massesul, volume,rho, rho_BC
36c
37        PARAMETER (rho_BC=1.8E3)  !--dry density kg/m3
38c
39c--growth factor for radius for sulfate (=non-BC) aerosols 
40c--output from exact calculations as a function of RH
41        DATA rh_growth/1.000, 1.000, 1.000, 1.000, 1.169, 1.220, 
42     .                 1.282, 1.363, 1.485, 1.580, 1.735, 2.085 /
43c
44c--density for sulfate (=non-BC) aerosols
45c--output from exact calculations as a function of RH
46        DATA rho_SUL/ 1769., 1769., 1769., 1769., 1430., 1390., 
47     .                1349., 1302., 1245., 1210., 1165., 1101./ !--kg/m3
48c
49        INTEGER class,Nclass
50        PARAMETER (Nclass=7)
51c
52c--dry mass content of BC
53        REAL bc_content(Nclass) 
54        DATA bc_content/0.0, 0.001, 0.01, 0.02, 0.05, 0.1, 0.2/
55c
56        REAL bccontentbymass, sulcontentbymass, drycontentbymass
57        REAL bccontentbyvol, sulcontentbyvol
58c
59c-------------------------------------
60c
61        COMPLEX m          !----refractive index m=n_r-i*n_i
62c--the following is added by Rong Wang
63        COMPLEX zmax1, zmax2, zmax3
64c--end
65        INTEGER Nmax,Nstart !--number of iterations
66        COMPLEX k2, k3, z1, z2
67        COMPLEX u1,u5,u6,u8
68        COMPLEX a(1:21000), b(1:21000)
69        COMPLEX I 
70        INTEGER n  !--loop index
71        REAL pi, nnn
72        COMPLEX nn
73        REAL Q_ext, Q_abs, Q_sca, g, omega   !--parameters for radius r
74        REAL x  !--size parameter
75        REAL r  !--radius
76        REAL sigma_sca, sigma_ext, sigma_abs
77        REAL omegatot,  gtot !--averaged parameters
78        COMPLEX ksiz2(-1:21000), psiz2(1:21000)
79        COMPLEX nu1z1(1:21010), nu1z2(1:21010)
80        COMPLEX nu3z2(0:21000)
81        REAL dnumber, deltar
82        INTEGER bin, Nbin, k
83        PARAMETER (Nbin=700)
84c
85C---wavelengths STREAMER
86        INTEGER Nwv, NwvmaxSW
87        PARAMETER (NwvmaxSW=25)
88        REAL lambda(1:NwvmaxSW)
89        DATA lambda/0.28E-6, 0.30E-6, 0.33E-6, 0.36E-6, 0.40E-6,
90     .              0.44E-6, 0.48E-6, 0.52E-6, 0.57E-6, 0.64E-6,
91     .              0.69E-6, 0.75E-6, 0.78E-6, 0.87E-6, 1.00E-6,
92     .              1.10E-6, 1.19E-6, 1.28E-6, 1.53E-6, 1.64E-6,
93     .              2.13E-6, 2.38E-6, 2.91E-6, 3.42E-6, 4.00E-6/
94c
95        INTEGER nb, nb_lambda
96        PARAMETER (nb_lambda=5)
97        REAL lambda_ref(nb_lambda)
98        DATA lambda_ref /0.443E-6,0.550E-6,0.670E-6,0.765E-6,0.865E-6/ 
99c
100C---BOA fluxes from typical STREAMER run - Tad
101        REAL weight(1:NwvmaxSW-1)
102        DATA weight/ 0.01,  4.05, 9.51,  15.99, 26.07, 33.10, 33.07,
103     .              39.91, 52.67, 27.89, 43.60, 13.67, 42.22, 40.12,
104     .              32.70, 14.44, 19.48, 14.23, 13.43, 16.42,  8.33,
105     .               0.95,  0.65, 2.76/
106c
107        REAL lambda_int(1:NwvmaxSW-1+nb_lambda)
108c
109        REAL n_r_BC(1:NwvmaxSW-1+nb_lambda)
110        REAL n_i_BC(1:NwvmaxSW-1+nb_lambda)
111c
112        REAL n_r_SUL(1:Nrh,1:NwvmaxSW-1+nb_lambda)
113        REAL n_i_SUL(1:Nrh,1:NwvmaxSW-1+nb_lambda)
114c
115        REAL final_a(1:NwvmaxSW-1+nb_lambda)
116        REAL final_g(1:NwvmaxSW-1+nb_lambda)
117        REAL final_w(1:NwvmaxSW-1+nb_lambda)
118c
119        INTEGER band, NbandSW
120        PARAMETER (NbandSW=6)
121c
122        REAL gcm_a(NbandSW)
123        REAL gcm_g(NbandSW)
124        REAL gcm_w(NbandSW)
125        REAL gcm_weight_a(NbandSW)
126        REAL gcm_weight_g(NbandSW)
127        REAL gcm_weight_w(NbandSW)
128c
129        REAL ss_a(Nclass,Nrh,NbandSW+nb_lambda)
130        REAL ss_w(Nclass,Nrh,NbandSW+nb_lambda)
131        REAL ss_g(Nclass,Nrh,NbandSW+nb_lambda)
132c
133        REAL ss_a_bc(Nclass,Nrh,NbandSW+nb_lambda)
134        REAL ss_w_bc(Nclass,Nrh,NbandSW+nb_lambda)
135        REAL ss_g_bc(Nclass,Nrh,NbandSW+nb_lambda)
136c
137        INTEGER wv, nb_wv_BCr, nb_wv_BCi, nb_wv_SUL
138        REAL count_n_r, count_n_i
139c
140        REAL n_r, n_i
141c
142c--refractive index for soluble stuff other than BC
143        PARAMETER (nb_wv_SUL=100)
144        REAL wv_SUL(1:nb_wv_SUL)
145        REAL index_r_SUL(1:Nrh,1:nb_wv_SUL)
146        REAL index_i_SUL(1:Nrh,1:nb_wv_SUL)
147
148c--refractive index for BC 
149c--Greg Schuster data
150c        PARAMETER (nb_wv_BCr=21, nb_wv_BCi=21)
151c--Sheffield
152        PARAMETER (nb_wv_BCr=61, nb_wv_BCi=61)
153        REAL wv_BCr(1:nb_wv_BCr), wv_BCi(1:nb_wv_BCi)
154        REAL index_r_BC(1:nb_wv_BCr), index_i_BC(1:nb_wv_BCi)
155        REAL dummy
156c--this is used by Rong Wang to test the square root
157c         zmax1=CMPLX(0,2)
158c         zmax1=zmax1**(1./2.)
159c         n_r=REAL(zmax1)
160c         n_i=AIMAG(zmax1)
161c         print *, n_r,n_i
162c       
163c         zmax1=CMPLX(-3,4) 
164c         zmax1=zmax1**(1./2.)
165c         n_r=REAL(zmax1)
166c         n_i=AIMAG(zmax1)
167c         print *, n_r,n_i
168c--end
169c--reading real part of refractive index
170c        OPEN(unit=10,file='r_bc_greg.dat')
171        OPEN(unit=10,file='r_bcsoot_she.dat')
172        DO wv=1, nb_wv_BCr
173          READ (10,*) wv_BCr(wv), index_r_BC(wv)
174        ENDDO
175        CLOSE(10)
176c
177c--reading imaginary part of refractive index
178c        OPEN(unit=10,file='i_bc_greg.dat')
179        OPEN(unit=10,file='i_bcsoot_she.dat')
180        DO wv=1, nb_wv_BCi
181          READ (10,*) wv_BCi(wv), index_i_BC(wv)
182        ENDDO
183        CLOSE(10)
184c
185c--reading imaginary part of refractive index
186c--for sulfate (=non-BC) aerosols
187        OPEN(unit=10,file='ri_sul_v2')
188        DO irh=1, Nrh
189        DO wv=1, nb_wv_SUL
190          READ (10,*) dummy, wv_SUL(wv), 
191     .                index_r_SUL(irh,wv), index_i_SUL(irh,wv)
192        ENDDO
193        ENDDO
194        CLOSE(10)
195c
196c------------------------------------------------------------
197c--initialising BC refractice indices before interpolation
198        DO Nwv=1, NwvmaxSW-1+nb_lambda
199           n_r_BC(Nwv)=0.0
200           n_i_BC(Nwv)=0.0
201        ENDDO
202c
203c--interpolating on our wavelengths
204        DO Nwv=1, NwvmaxSW-1
205c
206          lambda_int(Nwv)=( lambda(Nwv)+lambda(Nwv+1) ) /2.
207c
208c--first real part
209          count_n_r=0.0
210          DO wv=1, nb_wv_BCr
211            IF (wv_BCr(wv)/1.e9.GT.lambda(Nwv).AND.
212     .          wv_BCr(wv)/1.e9.LT.lambda(Nwv+1)) THEN
213              n_r_BC(Nwv)=n_r_BC(Nwv)+index_r_BC(wv)
214              count_n_r=count_n_r+1.0
215            ENDIF 
216          ENDDO
217c
218          IF (count_n_r.GT.0.5) THEN
219c--averaging
220            n_r_BC(Nwv)=n_r_BC(Nwv)/count_n_r
221          ELSE
222c--otherwise closest neighbour
223          DO wv=1, nb_wv_BCr
224            IF (wv_BCr(wv)/1.e9.LT.lambda_int(Nwv)) THEN
225              n_r_BC(Nwv)=index_r_BC(wv)
226            ENDIF 
227          ENDDO
228          ENDIF
229c
230c--then imaginary part
231          count_n_i=0.0
232          DO wv=1, nb_wv_BCi
233            IF (wv_BCi(wv)/1.e9.GT.lambda(Nwv).AND.
234     .          wv_BCi(wv)/1.e9.LT.lambda(Nwv+1)) THEN
235              n_i_BC(Nwv)=n_i_BC(Nwv)+index_i_BC(wv)
236              count_n_i=count_n_i+1.0
237            ENDIF 
238          ENDDO
239c
240          IF (count_n_i.GT.0.5) THEN
241c--averaging
242            n_i_BC(Nwv)=n_i_BC(Nwv)/count_n_i
243          ELSE
244c--otherwise closest neighbour
245          DO wv=1, nb_wv_BCi
246            IF (wv_BCi(wv)/1.e9.LT.lambda_int(Nwv)) THEN
247              n_i_BC(Nwv)=index_i_BC(wv)
248            ENDIF 
249          ENDDO
250          ENDIF
251c
252        ENDDO
253c
254c------------------------------------------------------------
255c--initialising SUL refractice indices before interpolation
256        DO irh=1, Nrh
257c
258        DO Nwv=1, NwvmaxSW-1+nb_lambda
259           n_r_SUL(irh,Nwv)=0.0
260           n_i_SUL(irh,Nwv)=0.0
261        ENDDO
262c
263c--interpolating on our wavelengths
264c
265        DO Nwv=1, NwvmaxSW-1
266c
267          lambda_int(Nwv)=( lambda(Nwv)+lambda(Nwv+1) ) /2.
268c
269c--first real part
270          count_n_r=0.0
271          DO wv=1, nb_wv_SUL
272            IF (wv_SUL(wv).GT.lambda(Nwv).AND.
273     .          wv_SUL(wv).LT.lambda(Nwv+1)) THEN
274              n_r_SUL(irh,Nwv)=n_r_SUL(irh,Nwv)+index_r_SUL(irh,wv)
275              count_n_r=count_n_r+1.0
276            ENDIF 
277          ENDDO
278c
279          IF (count_n_r.GT.0.5) THEN
280c--averaging
281            n_r_SUL(irh,Nwv)=n_r_SUL(irh,Nwv)/count_n_r
282          ELSE
283c--otherwise closest neighbour
284          DO wv=1, nb_wv_SUL
285            IF (wv_SUL(wv).LT.lambda_int(Nwv)) THEN
286              n_r_SUL(irh,Nwv)=index_r_SUL(irh,wv)
287            ENDIF 
288          ENDDO
289          ENDIF
290c
291c--then imaginary part
292          count_n_i=0.0
293          DO wv=1, nb_wv_SUL
294            IF (wv_SUL(wv).GT.lambda(Nwv).AND.
295     .          wv_SUL(wv).LT.lambda(Nwv+1)) THEN
296              n_i_SUL(irh,Nwv)=n_i_SUL(irh,Nwv)+index_i_SUL(irh,wv)
297              count_n_i=count_n_i+1.0
298            ENDIF 
299          ENDDO
300c
301          IF (count_n_i.GT.0.5) THEN
302c--averaging
303            n_i_SUL(irh,Nwv)=n_i_SUL(irh,Nwv)/count_n_i
304          ELSE
305c--otherwise closest neighbour
306          DO wv=1, nb_wv_SUL
307            IF (wv_SUL(wv).LT.lambda_int(Nwv)) THEN
308              n_i_SUL(irh,Nwv)=index_i_SUL(irh,wv)
309            ENDIF 
310          ENDDO
311          ENDIF
312c
313        ENDDO
314c
315        ENDDO
316c
317c-----------------------------------------------------------
318c
319c--now defining nr and ni for my set of reference wavelengths
320        DO nb=1, nb_lambda
321c
322          lambda_int(NwvmaxSW-1+nb)=lambda_ref(nb)
323c
324c--for BC
325          DO wv=1, nb_wv_BCr
326            IF (wv_BCr(wv)/1.e9.LT.lambda_ref(nb)) THEN
327              n_r_BC(NwvmaxSW-1+nb)=index_r_BC(wv)
328            ENDIF 
329          ENDDO
330
331          DO wv=1, nb_wv_BCi
332            IF (wv_BCi(wv)/1.e9.LT.lambda_ref(nb)) THEN
333              n_i_BC(NwvmaxSW-1+nb)=index_i_BC(wv)
334            ENDIF
335          ENDDO
336c
337c--for SUL
338          DO irh=1,Nrh
339          DO wv=1, nb_wv_SUL
340            IF (wv_SUL(wv).LT.lambda_ref(nb)) THEN
341              n_r_SUL(irh,NwvmaxSW-1+nb)=index_r_SUL(irh,wv)
342              n_i_SUL(irh,NwvmaxSW-1+nb)=index_i_SUL(irh,wv)
343            ENDIF 
344          ENDDO
345          ENDDO
346c
347        ENDDO
348c
349c--n_r and n_i are defined manually for first band
350        n_r_BC(1)=index_r_BC(1)
351        n_i_BC(1)=index_i_BC(1)
352c
353c--here check that n_r and n_i are fine after interpolation scheme
354c        PRINT *,'n_r_BC=',   n_r_BC
355c        PRINT *,'n_i_BC=',   n_i_BC
356c        PRINT *,'n_r_SUL=', n_r_SUL
357c        PRINT *,'n_i_SUL=', n_i_SUL
358c
359c
360c---------------------------------------------------------------------
361c------MASTER LOOPS 
362c
363        DO class=1, Nclass
364c
365        DO irh=1, Nrh
366c
367c        print *,'class=', class, ' irh=', irh
368c
369c--Mie calculation starts here
370        DO Nwv=1, NwvmaxSW-1+nb_lambda
371c
372c--BC and non-BC mass fraction
373c--for the dry aerosol
374        bccontentbymass=bc_content(class)      !--dry mass fraction
375        sulcontentbymass=1.0-bccontentbymass  !--dry mass fraction
376c--rhgrowth for the mixed particle
377        rhgrowth=(( bccontentbymass/rho_BC + 
378     .              sulcontentbymass/rho_SUL(1)*rh_growth(irh)**3.) / 
379     .            ( bccontentbymass/rho_BC + 
380     .              sulcontentbymass/rho_SUL(1) ))**(1./3.)
381c--added by Rong Wang to compute dry mass content
382        drycontentbymass= ( bccontentbymass + sulcontentbymass ) /
383     .            (bccontentbymass +
384     .            sulcontentbymass*rh_growth(irh)**3.0*
385     .                  rho_SUL(irh)/rho_SUL(1) ) !--wet mass fraction
386c--BC and non-BC mass fraction
387c--correcting for the wet aerosol
388c--sulcontent is for sulfate + water
389        bccontentbymass=bccontentbymass / ( bccontentbymass + 
390     .                  sulcontentbymass*rh_growth(irh)**3.0*
391     .                  rho_SUL(irh)/rho_SUL(1) ) !--wet mass fraction
392        sulcontentbymass=1.0-bccontentbymass     !--wet mass fraction
393c--BC and non-BC volume fraction
394c--for the wet aerosol
395        bccontentbyvol=(bccontentbymass/rho_BC) / 
396     .                 (bccontentbymass/rho_BC + 
397     .                  sulcontentbymass/rho_SUL(irh)) !--wet vol fraction
398        sulcontentbyvol=1.0-bccontentbyvol             !--wet vol fraction
399c
400c
401c        print *, rhgrowth, rh_growth(irh)
402c        print *,'BC=',bc_content(class), bccontentbymass, bccontentbyvol
403c
404c--density of mixture rho = (Mbc+Msul)/(Vbc+Vsul)
405c                         = (Mbc+Msul)/Msul*Msul/Vsul*Vsul/(Vbc+Vsul)
406c                         = rhoSUL*sulcontentbyvol/sulcontentbymass
407        rho=rho_SUL(irh)*sulcontentbyvol/sulcontentbymass
408c        print *,'rho=',rho_BC, rho_SUL(irh), rho, 
409c
410c--volume weighed refractive index
411c--the following is modified by Rong Wang
412c        n_r=n_r_BC(Nwv)*bccontentbyvol+n_r_SUL(irh,Nwv)*sulcontentbyvol
413c        n_i=n_i_BC(Nwv)*bccontentbyvol+n_i_SUL(irh,Nwv)*sulcontentbyvol
414         zmax1=CMPLX(n_r_BC(Nwv),-n_i_BC(Nwv))
415         zmax2=CMPLX(n_r_SUL(irh,Nwv),-n_i_SUL(irh,Nwv))
416         zmax3=(zmax2**2)  *
417     .         (zmax1**2+2.*zmax2**2+2.*bccontentbyvol*(zmax1**2
418     .         -zmax2**2)) /
419     .         (zmax1**2+2.*zmax2**2-bccontentbyvol*(zmax1**2
420     .         -zmax2**2)) 
421         zmax3=zmax3**(1./2.)
422         n_r=REAL(zmax3)
423         n_i=-AIMAG(zmax3)
424c
425c----test print 550 nm refractive index
426c        if (Nwv.eq.NwvmaxSW-1+2) then
427c        print *, 'n_r=',n_r, n_r_BC(Nwv), n_r_SUL(irh,Nwv)
428c        print *, 'n_i=',n_i, n_i_BC(Nwv), n_i_SUL(irh,Nwv)
429c        endif
430c
431        m=CMPLX(n_r,-n_i)
432c
433        pi=4.*ATAN(1.)
434c
435        I=CMPLX(0.,1.)
436c
437        sigma_sca=0.0
438        sigma_ext=0.0
439        sigma_abs=0.0
440        gtot=0.0
441        omegatot=0.0
442        masse = 0.0
443        massebc = 0.0
444        massesul = 0.0
445        volume=0.0
446        nombre=0.0
447c
448        DO bin=0, Nbin !---loop on size bins
449 
450c--this radius is the wet aerosol radius
451        r=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
452        x=2.*pi*r/lambda_int(Nwv)
453        deltar=1./FLOAT(Nbin)*(log(rmax)-log(rmin))
454c
455c--computing number concentration using log-normal for the dry aerosols
456        dnumber=0
457        DO dis=1, Ndis
458         dnumber=dnumber+Ntot(dis)/SQRT(2.*pi)/log(sigma_g(dis))*
459     .            exp(-0.5*(log(r/(r_0(dis)*rhgrowth))
460     .                /log(sigma_g(dis)))**2)
461        ENDDO
462c--be careful here we compute the mass of BC in the mixed particle
463c--using the rho of the mixed particle gives total mass 
464c--then multiply by bc content by mass for the wet aerosol
465        massebc = massebc  +4./3.*pi*(r**3)*dnumber*
466     .                          deltar*rho*1.E3*bccontentbymass      !--g/m3
467c--added by Rong Wang
468        masse = masse  + 4./3.*pi*(r**3)*dnumber*
469     .                          deltar*rho*1.E3*drycontentbymass     !--g/m3
470        massesul = massesul  + 4./3.*pi*(r**3)*dnumber*
471     .      deltar*rho*1.E3*(drycontentbymass - bccontentbymass)     !--g/m3
472c--be careful here we compute the volume of BC in the mixed particle
473        volume=volume+4./3.*pi*(r**3)*dnumber*deltar*bccontentbyvol
474c--total number 
475        nombre=nombre+dnumber*deltar
476c
477        k2=m
478        k3=CMPLX(1.0,0.0)
479
480        z2=CMPLX(x,0.0)
481        z1=m*z2
482 
483        IF (0.0.LE.x.AND.x.LE.8.) THEN
484        Nmax=INT(x+4*x**(1./3.)+1.)+2
485        ELSEIF (8..LT.x.AND.x.LT.4200.) THEN
486        Nmax=INT(x+4.05*x**(1./3.)+2.)+1
487        ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN
488        Nmax=INT(x+4*x**(1./3.)+2.)+1
489        ELSE
490        WRITE(10,*) 'x out of bound, x=', x
491        STOP
492        ENDIF
493
494        Nstart=Nmax+10
495
496C-----------loop for nu1z1, nu1z2
497
498       nu1z1(Nstart)=CMPLX(0.0,0.0)
499       nu1z2(Nstart)=CMPLX(0.0,0.0)
500       DO n=Nstart-1, 1 , -1
501       nn=CMPLX(FLOAT(n),0.0)
502       nu1z1(n)=(nn+1.)/z1 - 1./( (nn+1.)/z1 + nu1z1(n+1) ) 
503       nu1z2(n)=(nn+1.)/z2 - 1./( (nn+1.)/z2 + nu1z2(n+1) )
504       ENDDO
505
506C------------loop for nu3z2
507
508       nu3z2(0)=-I
509       DO n=1, Nmax
510       nn=CMPLX(FLOAT(n),0.0)
511       nu3z2(n)=-nn/z2 + 1./ (nn/z2 - nu3z2(n-1) ) 
512       ENDDO
513
514C-----------loop for psiz2 and ksiz2 (z2)
515       ksiz2(-1)=COS(REAL(z2))-I*SIN(REAL(z2))
516       ksiz2(0)=SIN(REAL(z2))+I*COS(REAL(z2))
517       DO n=1,Nmax
518       nn=CMPLX(FLOAT(n),0.0)
519       ksiz2(n)=(2.*nn-1.)/z2 * ksiz2(n-1) - ksiz2(n-2)
520       psiz2(n)=CMPLX(REAL(ksiz2(n)),0.0)
521       ENDDO
522
523C-----------loop for a(n) and b(n)
524   
525       DO n=1, Nmax
526        u1=k3*nu1z1(n) - k2*nu1z2(n)
527        u5=k3*nu1z1(n) - k2*nu3z2(n)
528        u6=k2*nu1z1(n) - k3*nu1z2(n)
529        u8=k2*nu1z1(n) - k3*nu3z2(n)
530        a(n)=psiz2(n)/ksiz2(n) * u1/u5
531        b(n)=psiz2(n)/ksiz2(n) * u6/u8
532       ENDDO
533
534C-----------------final loop--------------
535        Q_ext=0.0
536        Q_sca=0.0
537        g=0.0
538        DO n=Nmax-1,1,-1
539          nnn=FLOAT(n)
540          Q_ext=Q_ext+ (2.*nnn+1.) * REAL( a(n)+b(n) ) 
541          Q_sca=Q_sca+ (2.*nnn+1.) * 
542     .               REAL( a(n)*CONJG(a(n)) + b(n)*CONJG(b(n)) )
543          g=g + nnn*(nnn+2.)/(nnn+1.) * 
544     .       REAL( a(n)*CONJG(a(n+1))+b(n)*CONJG(b(n+1)) )  + 
545     .          (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n)))
546        ENDDO
547        Q_ext=2./x**2 * Q_ext
548        Q_sca=2./x**2 * Q_sca
549        Q_abs=Q_ext-Q_sca
550        IF (AIMAG(m).EQ.0.0) Q_abs=0.0
551        omega=Q_sca/Q_ext
552        g=g*4./x**2/Q_sca
553c
554        sigma_sca=sigma_sca+r**2*Q_sca*dnumber*deltar
555        sigma_abs=sigma_abs+r**2*Q_abs*dnumber*deltar
556        sigma_ext=sigma_ext+r**2*Q_ext*dnumber*deltar
557        omegatot=omegatot+r**2*Q_ext*omega*dnumber*deltar
558        gtot    =gtot+r**2*Q_sca*g*dnumber*deltar 
559c
560        ENDDO   !---bin
561C------------------------------------------------------------------
562
563        sigma_sca=pi*sigma_sca
564        sigma_abs=pi*sigma_abs
565        sigma_ext=pi*sigma_ext
566        gtot=pi*gtot/sigma_sca
567        omegatot=pi*omegatot/sigma_ext
568c
569        final_g(Nwv)=gtot
570        final_w(Nwv)=min(1.,omegatot)
571c-- Rong Wang modify this to compute the alpha based on the dry mass of
572c-- whole mixture (black carbon + sulfate)
573         final_a(Nwv)=sigma_ext/masse
574c
575        ENDDO  !--loop on wavelength
576c
577c---averaging over LMDZ wavebands
578c
579c        print *, final_a(6), final_a(NwvmaxSW), final_a(NwvmaxSW+1)
580        DO band=1, NbandSW
581          gcm_a(band)=0.0
582          gcm_g(band)=0.0
583          gcm_w(band)=0.0
584          gcm_weight_a(band)=0.0
585          gcm_weight_g(band)=0.0
586          gcm_weight_w(band)=0.0
587        ENDDO
588c
589c---band 1 is now in the UV, so we take the first wavelength as being representative
590c---it doesn't matter anyway because all radiation is absorbed in the stratosphere
591       DO Nwv=1,1
592          band=1
593          gcm_a(band)=gcm_a(band)+final_a(Nwv)*weight(Nwv)
594          gcm_weight_a(band)=gcm_weight_a(band)+weight(Nwv)
595          gcm_w(band)=gcm_w(band)+
596     .                final_w(Nwv)*final_a(Nwv)*weight(Nwv)
597          gcm_weight_w(band)=gcm_weight_w(band)+
598     .                       final_a(Nwv)*weight(Nwv)
599          gcm_g(band)=gcm_g(band)+
600     .                final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
601          gcm_weight_g(band)=gcm_weight_g(band)+
602     .                       final_a(Nwv)*final_w(Nwv)*weight(Nwv)
603        ENDDO
604c
605       DO Nwv=1,NwvmaxSW-1
606c
607        IF (Nwv.LE.5) THEN          !--RRTM spectral interval 2
608          band=2
609        ELSEIF (Nwv.LE.10) THEN     !--RRTM spectral interval 3
610          band=3
611        ELSEIF (Nwv.LE.16) THEN     !--RRTM spectral interval 4
612          band=4
613        ELSEIF (Nwv.LE.21) THEN     !--RRTM spectral interval 5
614          band=5
615        ELSE                        !--RRTM spectral interval 6
616          band=6
617        ENDIF
618c
619        gcm_a(band)=gcm_a(band)+final_a(Nwv)*weight(Nwv)
620        gcm_weight_a(band)=gcm_weight_a(band)+weight(Nwv)
621        gcm_w(band)=gcm_w(band)+
622     .              final_w(Nwv)*final_a(Nwv)*weight(Nwv)
623        gcm_weight_w(band)=gcm_weight_w(band)+
624     .                     final_a(Nwv)*weight(Nwv)
625        gcm_g(band)=gcm_g(band)+
626     .              final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
627        gcm_weight_g(band)=gcm_weight_g(band)+
628     .                     final_a(Nwv)*final_w(Nwv)*weight(Nwv)
629        ENDDO
630c
631        DO band=1, NbandSW
632          gcm_a(band)=gcm_a(band)/gcm_weight_a(band)
633          gcm_w(band)=gcm_w(band)/gcm_weight_w(band)
634          gcm_g(band)=gcm_g(band)/gcm_weight_g(band)
635          ss_a(class,irh,band)=gcm_a(band)
636          ss_w(class,irh,band)=gcm_w(band)
637          ss_g(class,irh,band)=gcm_g(band)
638        ENDDO
639c
640        DO nb=NbandSW+1, NbandSW+nb_lambda
641          ss_a(class,irh,nb)=final_a(NwvmaxSW-1+nb-NbandSW)
642          ss_w(class,irh,nb)=final_w(NwvmaxSW-1+nb-NbandSW)
643          ss_g(class,irh,nb)=final_g(NwvmaxSW-1+nb-NbandSW)
644        ENDDO
645c
646c-- Rong Wang (July 7, 2016))
647c-- From this line, it should be very careful
648c-- Here, I compute the ext, abs, sca for dry BC mass
649c-- I have checked these variables for mass:
650c--      masse, dry total mass; massebc, dry BC; massesul, dry sulfate
651c--         masse = massebc + massesul
652c--
653c-- I derive the ext, abs, sca for dry sulfate mass by setting class=1
654c--                                      where BCcontent=0
655c-- So, I compute the ext since class=2 here
656
657       IF (class.EQ.1) THEN ! For class=1, we do nothing (no BC)
658
659          DO nb=1, NbandSW+nb_lambda
660
661            ss_a_bc(class,irh,nb) = 0.0
662            ss_w_bc(class,irh,nb) = 1.0
663            ss_g_bc(class,irh,nb) = 0.0
664           
665          ENDDO
666
667       ELSE
668
669c--here we compute the properties for an equivalent BC that would be in
670c--external mixture. The properties do not make sense as such but have
671c--to be recombined with that of SUL (classe=1) to be those of the mixed
672c--particle
673
674          DO nb=1, NbandSW+nb_lambda
675
676            ! this ss_a is for dry BC only hereafter
677            ss_a_bc(class,irh,nb) = ( ss_a(class,irh,nb)*masse -
678     .        ss_a(1,irh,nb)*massesul ) / massebc
679
680            ! this ss_w is for dry BC only hereafter
681            ss_w_bc(class,irh,nb) = (ss_w(class,irh,nb)*
682     .        ss_a(class,irh,nb)*masse -
683     .        ss_w(1,irh,nb)*ss_a(1,irh,nb)*massesul ) /
684     .        (ss_a_bc(class,irh,nb)*massebc)
685
686            ! this ss_g is for dry BC only hereafter
687            ss_g_bc(class,irh,nb) = ( ss_g(class,irh,nb)*
688     .        ss_w(class,irh,nb)*ss_a(class,irh,nb)*masse -
689     .        ss_g(1,irh,nb)*ss_w(1,irh,nb)*ss_a(1,irh,nb)*massesul) /
690     .        (ss_w_bc(class,irh,nb)*ss_a_bc(class,irh,nb)*massebc)
691
692          ENDDO
693       ENDIF
694c-- End (Rong Wang)
695
696        ENDDO !--irh
697c
698        ENDDO !--Nclass
699c
700c--saving MEC, g and SSA for SUL for the model two bands
701        OPEN (unit=14,file='SEXT_sulfate_internal_mixture_6bands.txt')
702        WRITE(14,*) '  !-- Sulfate Accumulation (BC content=0)'
703        DO k=1,NbandSW
704        WRITE(14,950) (ss_a(1,irh,k),irh=1, Nrh)
705        ENDDO
706        CLOSE(14)
707
708        OPEN (unit=14,file='SSA_sulfate_internal_mixture_6bands.txt')
709        WRITE(14,*) '  !-- Sulfate Accumulation (BC content=0)'
710        DO k=1,NbandSW
711        WRITE(14,950) (ss_w(1,irh,k),irh=1, Nrh)
712        ENDDO
713        CLOSE(14)
714
715        OPEN (unit=14,file='G_sulfate_internal_mixture_6bands.txt')
716        WRITE(14,*) '  !-- Sulfate Accumulation (BC content=0)'
717        DO k=1,NbandSW
718        WRITE(14,950) (ss_g(1,irh,k),irh=1, Nrh)
719        ENDDO
720        CLOSE(14)
721
722        OPEN 
723     .    (unit=14,file='SEXT_sulfate+2%bc_internal_mixture_6bands.txt')
724        WRITE(14,*) '  !-- Sulfate Accumulation (BC content=2%)'
725        DO k=1,NbandSW
726        WRITE(14,950) (ss_a(4,irh,k),irh=1, Nrh)
727        ENDDO
728        CLOSE(14)
729
730        OPEN 
731     .    (unit=14,file='SSA_sulfate+2%bc_internal_mixture_6bands.txt')
732        WRITE(14,*) '  !-- Sulfate Accumulation (BC content=2%)'
733        DO k=1,NbandSW
734        WRITE(14,950) (ss_w(4,irh,k),irh=1, Nrh)
735        ENDDO
736        CLOSE(14)
737
738        OPEN 
739     .    (unit=14,file='G_sulfate+2%bc_internal_mixture_6bands.txt')
740        WRITE(14,*) '  !-- Sulfate Accumulation (BC content=2%)'
741        DO k=1,NbandSW
742        WRITE(14,950) (ss_g(4,irh,k),irh=1, Nrh)
743        ENDDO
744        CLOSE(14)
745
746c--saving MEC, g and SSA for equivalent BC the model two bands
747        OPEN (unit=14,file='DATA_BC_internal_mixture_6bands.txt')
748        WRITE(14,*) '  DATA alpha_MG_6bands/ &'
749        DO class=2, Nclass
750        WRITE(14,900) bc_content(class)
751        DO k=1,NbandSW
752        WRITE(14,951) (ss_a_bc(class,irh,k),irh=1, Nrh)
753        ENDDO
754        ENDDO
755
756        WRITE(14,*) '  '
757        WRITE(14,*) '  DATA piz_MG_6bands/ &'
758        DO class=2, Nclass
759        WRITE(14,900) bc_content(class)
760        DO k=1,NbandSW
761        WRITE(14,951) (ss_w_bc(class,irh,k),irh=1,Nrh)
762        ENDDO
763        ENDDO
764
765        WRITE(14,*) '  '
766        WRITE(14,*) '  DATA cg_MG_6bands/ &'
767        DO class=2, Nclass
768        WRITE(14,900) bc_content(class)
769        DO k=1,NbandSW
770        WRITE(14,951) (ss_g_bc(class,irh,k),irh=1,Nrh)
771        ENDDO
772        ENDDO
773        CLOSE(14)
774c
775c--saving MEC and MAC for SUL for our reference wavelengths
776        OPEN (unit=14, file='SEXT_sulfate_internal_mixture_5wave.txt')
777        WRITE(14,*) 
778     .      '  !-- Extinction Sulfate Accumulation (BC content=0)'
779        DO k=NbandSW+1,NbandSW+nb_lambda
780        WRITE(14,950) (ss_a(1,irh,k), irh=1,Nrh)
781        ENDDO
782        CLOSE(14)
783        OPEN (unit=14, file='SABS_sulfate_internal_mixture_5wave.txt')
784        WRITE(14,*) 
785     .      '  !-- Absorption Sulfate Accumulation (BC content=0)'
786        DO k=NbandSW+1,NbandSW+nb_lambda
787        WRITE(14,950) ((1.0-ss_w(1,irh,k))*ss_a(1,irh,k), irh=1,Nrh)
788        ENDDO
789        CLOSE(14)
790
791        OPEN 
792     .    (unit=14, file='SEXT_sulfate+2%bc_internal_mixture_5wave.txt')
793        WRITE(14,*) 
794     .      '  !-- Extinction Sulfate Accumulation (BC content=2%)'
795        DO k=NbandSW+1,NbandSW+nb_lambda
796        WRITE(14,950) (ss_a(4,irh,k), irh=1,Nrh)
797        ENDDO
798        CLOSE(14)
799        OPEN 
800     .    (unit=14, file='SABS_sulfate+2%bc_internal_mixture_5wave.txt')
801        WRITE(14,*) 
802     .      '  !-- Absorption Sulfate Accumulation (BC content=2%)'
803        DO k=NbandSW+1,NbandSW+nb_lambda
804        WRITE(14,950) ((1.0-ss_w(4,irh,k))*ss_a(4,irh,k), irh=1,Nrh)
805        ENDDO
806        CLOSE(14)
807
808c        OPEN (unit=14, file='SSA_sulfate_internal_mixture_5wave.txt')
809c        WRITE(14,*) '  !-- Sulfate Accumulation (BC content=0)'
810c        DO k=NbandSW+1,NbandSW+nb_lambda
811c        WRITE(14,951) (ss_w(1,irh,k), irh=1,Nrh)
812c        ENDDO
813c        CLOSE(14)
814c
815c        OPEN (unit=14, file='CG_sulfate_internal_mixture_5wave.txt')
816c        WRITE(14,*) '  !-- Sulfate Accumulation (BC content=0)'
817c        DO k=NbandSW+1,NbandSW+nb_lambda
818c        WRITE(14,951) (ss_g(1,irh,k), irh=1,Nrh)
819c        ENDDO
820c        CLOSE(14)
821c
822c--saving MEC and MAC for equivalent BC for our reference wavelengths
823        OPEN (unit=14, file='DATA_BC_internal_mixture_5wave.txt')
824        WRITE(14,*) '  DATA alpha_MG_5wv/ &' 
825        DO class=2, Nclass
826        WRITE(14,900) bc_content(class)
827        DO k=NbandSW+1,NbandSW+nb_lambda
828        WRITE(14,951) (ss_a_bc(class,irh,k),irh=1,Nrh)
829        ENDDO
830        ENDDO
831        WRITE(14,*) '  ' 
832        WRITE(14,*) '  DATA abs_MG_5wv/ &' 
833        DO class=2, Nclass
834        WRITE(14,900) bc_content(class)
835        DO k=NbandSW+1,NbandSW+nb_lambda
836        WRITE(14,951) 
837     .     ((1.0-ss_w_bc(class,irh,k))*ss_a_bc(class,irh,k),irh=1,Nrh)
838        ENDDO
839        ENDDO
840        CLOSE(14)
841
842c        OPEN (unit=14, file='SSA_BC_internal_mixture_5wave.txt')
843c        WRITE(14,*) '  DATA piz_MG_5wv/ &' 
844c        DO class=2, Nclass
845c        WRITE(14,900) bc_content(class)
846c        DO k=NbandSW+1,NbandSW+nb_lambda
847c        WRITE(14,951) (ss_w_bc(class,irh,k),irh=1,Nrh)
848c        ENDDO
849c        ENDDO
850c        CLOSE(14)
851c
852c        OPEN (unit=14, file='G_BC_internal_mixture_5wave.txt')
853c        WRITE(14,*) '  DATA cg_MG_5wv/ &' 
854c        DO class=2, Nclass
855c        WRITE(14,900) bc_content(class)
856c        DO k=NbandSW+1,NbandSW+nb_lambda
857c        WRITE(14,951) (ss_g_bc(class,irh,k),irh=1,Nrh)
858c        ENDDO
859c        ENDDO
860c        CLOSE(14)
861c
862900     FORMAT(1X,'!--BC content=',F5.3)
863950     FORMAT(1X,12(F6.3,','),' &')
864951     FORMAT(1X,12(F7.3,','),' &')
865c
866        END
Note: See TracBrowser for help on using the repository browser.