source: TOOLS/CMIP6_FORCING/AER_OPTICS/mie_6bands_bc_v4.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: 14.4 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 radius (um)/modal radius
16C--Ntot    : total concentration in m-3 
17c
18        REAL rmin, rmax    !----integral bounds in  m
19        PARAMETER (rmin=0.002E-6, rmax=30.E-6)
20c
21c--Nmode= 1 mode for external bc
22c--NDis=1 distribution per mode
23        INTEGER Nmode, Ndis, mode, dis
24        PARAMETER (Nmode=1, Ndis=1)
25        REAL sigma_g(Ndis,Nmode), r_0(Ndis,Nmode), Ntot(Ndis,Nmode)
26c--Accumulation Soluble (AS)
27        DATA r_0    /0.1E-6/   !--meters
28        DATA sigma_g/1.8/
29        DATA Ntot   /1.0/
30        CHARACTER*33 chmode(Nmode)
31        DATA chmode/'BC Accumulation Soluble (AS)'/
32c
33        REAL masse,volume,surface,rho
34        PARAMETER (rho=1.550E3)  !--BC dry density kg/m3 from Cathy Liousse
35c
36c---------- RH growth parameters----------------
37c
38        INTEGER rh_int,IRH
39        PARAMETER(rh_int=12)
40        REAL RH_tab(rh_int),RH, rh_dummy
41        DATA RH_tab/0.,10.,20,30.,40.,50.,60.,70.,80.,85.,90.,95./
42        REAL rwet
43c
44c--Carrico et al. Hygroscopic growth behaviour of a carbon-dominated
45c--aerosol in Yosemite Park, Atm. Env. 39, 1393-1404, 2005
46        REAL growth(rh_int,Nmode)
47        DATA growth/
48     .  1.0,1.0,1.0,1.0,1.0,1.05,1.10,1.20,1.29,1.30,1.31,1.40/ !--AS
49c
50        REAL n_r, n_i
51c
52c-------------------------------------
53c
54        COMPLEX m          !----refractive index m=n_r-i*n_i
55        INTEGER Nmax,Nstart !--number of iterations
56        COMPLEX k2, k3, z1, z2
57        COMPLEX u1,u5,u6,u8
58        COMPLEX a(1:21000), b(1:21000)
59        COMPLEX I 
60        INTEGER n  !--loop index
61        REAL pi, nnn
62        COMPLEX nn
63        REAL Q_ext, Q_abs, Q_sca, g, omega   !--parameters for radius r
64        REAL x  !--size parameter
65        REAL r  !--radius
66        REAL sigma_sca, sigma_ext, sigma_abs
67        REAL omegatot,  gtot !--averaged parameters
68        COMPLEX ksiz2(-1:21000), psiz2(1:21000)
69        COMPLEX nu1z1(1:21010), nu1z2(1:21010)
70        COMPLEX nu3z2(0:21000)
71        REAL number, deltar
72        INTEGER bin, Nbin, k
73        PARAMETER (Nbin=700)
74c
75C---wavelengths STREAMER
76        INTEGER Nwv, NwvmaxSW
77        PARAMETER (NwvmaxSW=25)
78        REAL lambda(1:NwvmaxSW)
79        DATA lambda/0.28E-6, 0.30E-6, 0.33E-6, 0.36E-6, 0.40E-6,
80     .              0.44E-6, 0.48E-6, 0.52E-6, 0.57E-6, 0.64E-6,
81     .              0.69E-6, 0.75E-6, 0.78E-6, 0.87E-6, 1.00E-6,
82     .              1.10E-6, 1.19E-6, 1.28E-6, 1.53E-6, 1.64E-6,
83     .              2.13E-6, 2.38E-6, 2.91E-6, 3.42E-6, 4.00E-6/
84c
85        INTEGER nb, nb_lambda
86        PARAMETER (nb_lambda=5)
87        REAL lambda_ref(nb_lambda)
88        DATA lambda_ref /0.443E-6,0.550E-6,0.670E-6,0.765E-6,0.865E-6/ 
89c
90C---TOA fluxes - Streamer Cs
91        REAL weight(1:NwvmaxSW-1)
92c        DATA weight/0.839920E1, 0.231208E2, 0.322393E2, 0.465058E2,
93c     .              0.678199E2, 0.798964E2, 0.771359E2, 0.888472E2,
94c     .              0.115281E3, 0.727565E2, 0.816992E2, 0.336172E2,
95c     .              0.914603E2, 0.112706E3, 0.658840E2, 0.524470E2,
96c     .              0.391067E2, 0.883864E2, 0.276672E2, 0.681812E2,
97c     .              0.190966E2, 0.250766E2, 0.128704E2, 0.698720E1/
98C---TOA fluxes - Tad
99c        DATA weight/ 4.20, 11.56, 16.12, 23.25, 33.91, 39.95, 38.57,
100c     .              44.42, 57.64, 29.36, 47.87, 16.81, 45.74, 56.35,
101c     .              32.94, 26.22, 19.55, 44.19, 13.83, 34.09,  9.55,
102c      .              12.54,  6.44,  3.49/
103C---BOA fluxes - Tad
104        DATA weight/ 0.01,  4.05, 9.51,  15.99, 26.07, 33.10, 33.07,
105     .              39.91, 52.67, 27.89, 43.60, 13.67, 42.22, 40.12,
106     .              32.70, 14.44, 19.48, 14.23, 13.43, 16.42,  8.33,
107     .               0.95,  0.65, 2.76/
108c
109        REAL lambda_int(1:NwvmaxSW-1+nb_lambda)
110c
111        REAL final_a(1:NwvmaxSW-1+nb_lambda)
112        REAL final_g(1:NwvmaxSW-1+nb_lambda)
113        REAL final_w(1:NwvmaxSW-1+nb_lambda)
114c
115        INTEGER band, NbandSW, NbandLW
116        PARAMETER (NbandSW=6, NbandLW=5)
117c
118        REAL gcm_a(NbandSW+NbandLW)
119        REAL gcm_g(NbandSW+NbandLW)
120        REAL gcm_w(NbandSW+NbandLW)
121        REAL gcm_weight_a(NbandSW+NbandLW)
122        REAL gcm_weight_g(NbandSW+NbandLW)
123        REAL gcm_weight_w(NbandSW+NbandLW)
124c
125        REAL ss_a(NbandSW+NbandLW+nb_lambda,rh_int)
126        REAL ss_w(NbandSW+NbandLW+nb_lambda,rh_int)
127        REAL ss_g(NbandSW+NbandLW+nb_lambda,rh_int)
128c
129        INTEGER NwvmaxLW
130        PARAMETER (NwvmaxLW=100)
131        REAL Tb, Planck
132        PARAMETER (Tb=260.0)
133c
134        INTEGER wv, nb_wv, nb_wv_i
135        PARAMETER (nb_wv=200)
136        REAL wv_bc(1:nb_wv)
137        REAL index_r_bc(RH_int,Nmode,1:nb_wv)
138        REAL index_i_bc(RH_int,Nmode,1:nb_wv)
139        REAL count_n_r, count_n_i
140c
141        pi=4.*ATAN(1.)
142        I=CMPLX(0.,1.)
143c
144c------opening output files
145c
146       OPEN (unit=14, file='SEXT_bc_soluble_6bands.txt')
147       OPEN (unit=15, file='G_bc_soluble_6bands.txt')
148       OPEN (unit=16, file='SSA_bc_soluble_6bands.txt')
149c
150       OPEN (unit=24, file='SEXT_bc_insoluble_6bands.txt')
151       OPEN (unit=25, file='G_bc_insoluble_6bands.txt')
152       OPEN (unit=26, file='SSA_bc_insoluble_6bands.txt')
153c
154       OPEN (unit=34, file='SEXT_bc_soluble_5wave.txt')
155       OPEN (unit=35, file='SEXT_bc_insoluble_5wave.txt')
156       OPEN (unit=36, file='SABS_bc_soluble_5wave.txt')
157       OPEN (unit=37, file='SABS_bc_insoluble_5wave.txt')
158c
159c--initializing wavelengths
160c
161       DO Nwv=1, NwvmaxSW-1
162         lambda_int(Nwv)=( lambda(Nwv)+lambda(Nwv+1) ) /2.
163       ENDDO
164c
165       DO nb=1, nb_lambda
166         lambda_int(NwvmaxSW-1+nb)=lambda_ref(nb)
167       ENDDO
168c
169       OPEN(unit=21,file='ri_bc_AS_v2')
170c
171c---now start calculations 
172c
173       DO mode=1, Nmode
174c
175c--lecture des indices from an old file
176       DO IRH=1,rh_int     
177       DO wv=1, nb_wv
178        READ (20+mode,*) rh_dummy, wv_bc(wv),
179     .        index_r_bc(irh,mode,wv), index_i_bc(irh,mode,wv)
180       ENDDO
181       ENDDO
182c
183c--loop over RH
184c
185        DO IRH=1,rh_int     
186c
187        DO Nwv=1, NwvmaxSW-1+nb_lambda
188c
189c--sinon plus proche voisin
190        n_r=index_r_bc(irh,mode,1)
191        n_i=index_i_bc(irh,mode,1)
192        DO wv=1, nb_wv
193          IF (wv_bc(wv).LT.lambda_int(Nwv)) THEN
194            n_r=index_r_bc(irh,mode,wv)
195            n_i=index_i_bc(irh,mode,wv)
196          ENDIF 
197        ENDDO
198c
199        m=CMPLX(n_r,-n_i)
200c
201        sigma_sca=0.0
202        sigma_ext=0.0
203        sigma_abs=0.0
204        gtot=0.0
205        omegatot=0.0
206        masse = 0.0
207        volume=0.0
208        surface=0.0
209c
210        DO bin=0, Nbin !---loop on size bins
211 
212c--here r is rdry
213        r=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
214c
215        rwet=growth(IRH,mode)*r
216c
217        x=2.*pi*rwet/lambda_int(Nwv)
218        deltar=1./FLOAT(Nbin)*(log(rmax)-log(rmin))
219c
220c--rwet r*1e2 to get into cm unit
221c
222        number=0
223        DO dis=1, Ndis
224          number=number+
225     .       Ntot(dis,mode)/SQRT(2.*pi)/log(sigma_g(dis,mode))*
226     .       exp(-0.5*(log(r/(r_0(dis,mode)))/
227     .               log(sigma_g(dis,mode)))**2)
228        ENDDO
229c--dry aerosol mass
230        masse=masse  +4./3.*pi*r**3*number*deltar*rho*1.E3          !--g/m3
231        volume=volume+4./3.*pi*r**3*number*deltar
232        surface=surface+4.*pi*r**2*number*deltar
233c
234        k2=m
235        k3=CMPLX(1.0,0.0)
236
237        z2=CMPLX(x,0.0)
238        z1=m*z2
239 
240        IF (0.0.LE.x.AND.x.LE.8.) THEN
241        Nmax=INT(x+4*x**(1./3.)+1.)+2
242        ELSEIF (8..LT.x.AND.x.LT.4200.) THEN
243        Nmax=INT(x+4.05*x**(1./3.)+2.)+1
244        ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN
245        Nmax=INT(x+4*x**(1./3.)+2.)+1
246        ELSE
247        WRITE(10,*) 'x out of bound, x=', x
248        STOP
249        ENDIF
250
251        Nstart=Nmax+10
252
253C-----------loop for nu1z1, nu1z2
254
255       nu1z1(Nstart)=CMPLX(0.0,0.0)
256       nu1z2(Nstart)=CMPLX(0.0,0.0)
257       DO n=Nstart-1, 1 , -1
258       nn=CMPLX(FLOAT(n),0.0)
259       nu1z1(n)=(nn+1.)/z1 - 1./( (nn+1.)/z1 + nu1z1(n+1) ) 
260       nu1z2(n)=(nn+1.)/z2 - 1./( (nn+1.)/z2 + nu1z2(n+1) )
261       ENDDO
262
263C------------loop for nu3z2
264
265       nu3z2(0)=-I
266       DO n=1, Nmax
267       nn=CMPLX(FLOAT(n),0.0)
268       nu3z2(n)=-nn/z2 + 1./ (nn/z2 - nu3z2(n-1) ) 
269       ENDDO
270
271C-----------loop for psiz2 and ksiz2 (z2)
272       ksiz2(-1)=COS(REAL(z2))-I*SIN(REAL(z2))
273       ksiz2(0)=SIN(REAL(z2))+I*COS(REAL(z2))
274       DO n=1,Nmax
275       nn=CMPLX(FLOAT(n),0.0)
276       ksiz2(n)=(2.*nn-1.)/z2 * ksiz2(n-1) - ksiz2(n-2)
277       psiz2(n)=CMPLX(REAL(ksiz2(n)),0.0)
278       ENDDO
279
280C-----------loop for a(n) and b(n)
281   
282       DO n=1, Nmax
283        u1=k3*nu1z1(n) - k2*nu1z2(n)
284        u5=k3*nu1z1(n) - k2*nu3z2(n)
285        u6=k2*nu1z1(n) - k3*nu1z2(n)
286        u8=k2*nu1z1(n) - k3*nu3z2(n)
287        a(n)=psiz2(n)/ksiz2(n) * u1/u5
288        b(n)=psiz2(n)/ksiz2(n) * u6/u8
289       ENDDO
290
291C-----------------final loop--------------
292        Q_ext=0.0
293        Q_sca=0.0
294        g=0.0
295        DO n=Nmax-1,1,-1
296          nnn=FLOAT(n)
297          Q_ext=Q_ext+ (2.*nnn+1.) * REAL( a(n)+b(n) ) 
298          Q_sca=Q_sca+ (2.*nnn+1.) * 
299     .               REAL( a(n)*CONJG(a(n)) + b(n)*CONJG(b(n)) )
300          g=g + nnn*(nnn+2.)/(nnn+1.) * 
301     .       REAL( a(n)*CONJG(a(n+1))+b(n)*CONJG(b(n+1)) )  + 
302     .          (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n)))
303        ENDDO
304        Q_ext=2./x**2 * Q_ext
305        Q_sca=2./x**2 * Q_sca
306        Q_abs=Q_ext-Q_sca
307        IF (AIMAG(m).EQ.0.0) Q_abs=0.0
308        omega=Q_sca/Q_ext
309        g=g*4./x**2/Q_sca
310c
311        sigma_sca=sigma_sca+rwet**2*Q_sca*number*deltar
312        sigma_abs=sigma_abs+rwet**2*Q_abs*number*deltar
313        sigma_ext=sigma_ext+rwet**2*Q_ext*number*deltar
314        omegatot=omegatot+rwet**2*Q_ext*omega*number*deltar
315        gtot    =gtot+rwet**2*Q_sca*g*number*deltar 
316c
317        ENDDO   !---bin
318C------------------------------------------------------------------
319
320        sigma_sca=pi*sigma_sca
321        sigma_abs=pi*sigma_abs
322        sigma_ext=pi*sigma_ext
323        gtot=pi*gtot/sigma_sca
324        omegatot=pi*omegatot/sigma_ext
325c
326        final_g(Nwv)=gtot
327        final_w(Nwv)=omegatot
328        final_a(Nwv)=sigma_ext/masse
329c
330        ENDDO  !--loop on wavelength
331c
332c---averaging over LMDZ wavebands
333c
334        DO band=1, NbandSW
335          gcm_a(band)=0.0
336          gcm_g(band)=0.0
337          gcm_w(band)=0.0
338          gcm_weight_a(band)=0.0
339          gcm_weight_g(band)=0.0
340          gcm_weight_w(band)=0.0
341        ENDDO
342c
343c---band 1 is now in the UV, so we take the first wavelength as being representative 
344c---it doesn't matter anyway because all radiation is absorbed in the stratosphere
345       DO Nwv=1,1
346          band=1
347          gcm_a(band)=gcm_a(band)+final_a(Nwv)*weight(Nwv)
348          gcm_weight_a(band)=gcm_weight_a(band)+weight(Nwv)
349          gcm_w(band)=gcm_w(band)+
350     .                final_w(Nwv)*final_a(Nwv)*weight(Nwv)
351          gcm_weight_w(band)=gcm_weight_w(band)+
352     .                       final_a(Nwv)*weight(Nwv)
353          gcm_g(band)=gcm_g(band)+
354     .                final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
355          gcm_weight_g(band)=gcm_weight_g(band)+
356     .                       final_a(Nwv)*final_w(Nwv)*weight(Nwv)
357        ENDDO
358c
359       DO Nwv=1,NwvmaxSW-1
360c
361        IF (Nwv.LE.5) THEN          !--RRTM spectral interval 2
362          band=2
363        ELSEIF (Nwv.LE.10) THEN     !--RRTM spectral interval 3
364          band=3
365        ELSEIF (Nwv.LE.16) THEN     !--RRTM spectral interval 4
366          band=4
367        ELSEIF (Nwv.LE.21) THEN     !--RRTM spectral interval 5
368          band=5
369        ELSE                        !--RRTM spectral interval 6
370          band=6
371        ENDIF
372c
373        gcm_a(band)=gcm_a(band)+final_a(Nwv)*weight(Nwv)
374        gcm_weight_a(band)=gcm_weight_a(band)+weight(Nwv)
375        gcm_w(band)=gcm_w(band)+
376     .              final_w(Nwv)*final_a(Nwv)*weight(Nwv)
377        gcm_weight_w(band)=gcm_weight_w(band)+
378     .                     final_a(Nwv)*weight(Nwv)
379        gcm_g(band)=gcm_g(band)+
380     .              final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
381        gcm_weight_g(band)=gcm_weight_g(band)+
382     .                     final_a(Nwv)*final_w(Nwv)*weight(Nwv)
383
384        ENDDO
385c
386        DO band=1, NbandSW
387          gcm_a(band)=gcm_a(band)/gcm_weight_a(band)
388          gcm_w(band)=gcm_w(band)/gcm_weight_w(band)
389          gcm_g(band)=gcm_g(band)/gcm_weight_g(band)
390          ss_a(band,IRH)=gcm_a(band)
391          ss_w(band,IRH)=gcm_w(band)
392          ss_g(band,IRH)=gcm_g(band)
393        ENDDO
394c
395        DO nb=NbandSW+1, NbandSW+nb_lambda
396          ss_a(nb,IRH)=final_a(NwvmaxSW-1+nb-NbandSW)
397          ss_w(nb,IRH)=final_w(NwvmaxSW-1+nb-NbandSW)
398          ss_g(nb,IRH)=final_g(NwvmaxSW-1+nb-NbandSW)
399        ENDDO
400c
401        ENDDO   !--fin boucle sur RH
402c
403c--Outputs
404C
405        WRITE(14,*) '  ! '//chmode(mode)
406        DO k=1, NbandSW
407          WRITE(14,951) (ss_a(k,IRH),IRH=1,rh_int)
408        ENDDO
409        WRITE(15,*) '  ! '//chmode(mode)
410        DO k=1, NbandSW
411          WRITE(15,951) (ss_g(k,IRH),IRH=1,rh_int)
412        ENDDO
413        WRITE(16,*) '  ! '//chmode(mode)
414        DO k=1, NbandSW
415          WRITE(16,951) (ss_w(k,IRH),IRH=1,rh_int)
416        ENDDO
417c
418        WRITE(34,*) '  ! extinction '//chmode(mode)
419        DO k=NbandSW+1,NbandSW+nb_lambda
420        WRITE(34,951) (ss_a(k,IRH),IRH=1,rh_int)
421        ENDDO
422        WRITE(36,*) '  ! absorption '//chmode(mode)
423        DO k=NbandSW+1,NbandSW+nb_lambda
424        WRITE(36,951) 
425     &   ((1.0-ss_w(k,IRH))*ss_a(k,IRH),IRH=1,rh_int)
426        ENDDO
427c
428        WRITE(24,*) '  ! BC insoluble '
429        WRITE(24,954) (ss_a(k,1),k=1,NbandSW)
430        WRITE(25,*) '  ! BC insoluble '
431        WRITE(25,954) (ss_w(k,1),k=1,NbandSW)
432        WRITE(26,*) '  ! BC insoluble '
433        WRITE(26,954) (ss_g(k,1),k=1,NbandSW)
434c
435        WRITE(35,*) '  ! extinction BC insoluble '
436        WRITE(35,955) (ss_a(k,1),k=NbandSW+1,NbandSW+nb_lambda)
437        WRITE(37,*) '  ! absorption BC insoluble '
438        WRITE(37,955)
439     &   ((1.0-ss_w(k,1))*ss_a(k,1),k=NbandSW+1,NbandSW+nb_lambda)
440c
441        ENDDO !--boucle sur les modes
442
443951     FORMAT(1X,12(F6.3,','),' &')
444954     FORMAT(1X,6(F6.3,','),' &')
445955     FORMAT(1X,5(F6.3,','),' &')
446c
447101     FORMAT(F10.6,E13.6)
448c
449        CLOSE(14)
450        CLOSE(15)
451        CLOSE(16)
452c
453        CLOSE(24)
454        CLOSE(25)
455        CLOSE(26)
456c
457        CLOSE(34)
458        CLOSE(35)
459        CLOSE(36)
460        CLOSE(37)
461c
462        END
Note: See TracBrowser for help on using the repository browser.