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