source: TOOLS/CMIP6_FORCING/AER_OPTICS/mie_6bands_nitrate_v4.f @ 4049

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

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

File size: 14.7 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 nitrate
22c--NDis=1 distribution per mode
23        INTEGER Nmode, Ndis, mode, dis
24        PARAMETER (Nmode=2, Ndis=1)
25        REAL sigma_g(Ndis,Nmode), r_0(Ndis,Nmode), Ntot(Ndis,Nmode)
26c--Accumulation Soluble (AS)
27        DATA r_0    /0.433E-6,0.1E-6/   !--meters
28        DATA sigma_g/2.0,1.8/
29        DATA Ntot   /1.0,1.0/
30        CHARACTER*33 chmode(Nmode)
31        DATA chmode/'Nitrate Coarse Soluble (CS)',
32     .              'Nitrate Accumulation Soluble (AS)'/
33c
34        REAL masse,volume,surface,rho
35        PARAMETER (rho=1.725E3)  !--dry density nitrate kg/m3
36c
37c---------- RH growth parameters----------------
38c
39        INTEGER rh_int,IRH
40        PARAMETER(rh_int=12)
41        REAL RH_tab(rh_int),RH, rh_dummy
42        DATA RH_tab/0.,10.,20,30.,40.,50.,60.,70.,80.,85.,90.,95./
43        REAL rwet
44c
45        REAL growth(rh_int,Nmode)
46        DATA growth/
47     .       1.0,1.0,1.0,1.05,1.10,1.15,1.21,1.33,1.50,1.65,1.85,1.90, !--CS
48     .       1.0,1.0,1.0,1.05,1.10,1.15,1.21,1.33,1.50,1.65,1.85,1.90/ !--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_nitrate(1:nb_wv)
137        REAL index_r_nitrate(RH_int,Nmode,1:nb_wv)
138        REAL index_i_nitrate(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_nitrate_soluble_6bands.txt')
147       OPEN (unit=15, file='G_nitrate_soluble_6bands.txt')
148       OPEN (unit=16, file='SSA_nitrate_soluble_6bands.txt')
149c
150       OPEN (unit=34, file='SEXT_nitrate_soluble_5wave.txt')
151       OPEN (unit=35, file='SEXT_nitrate_insoluble_5wave.txt')
152       OPEN (unit=36, file='SABS_nitrate_soluble_5wave.txt')
153       OPEN (unit=37, file='SABS_nitrate_insoluble_5wave.txt')
154c
155       OPEN (unit=24, file='SEXT_nitrate_insoluble_6bands.txt')
156       OPEN (unit=25, file='G_nitrate_insoluble_6bands.txt')
157       OPEN (unit=26, file='SSA_nitrate_insoluble_6bands.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
169c
170c---now start calculations 
171c
172       DO mode=1, Nmode
173c
174c--lecture des indices from an old file
175       OPEN(unit=20,file='ri_nitrate_AS_v2')  !--same data for AS and CS
176       DO IRH=1,rh_int     
177       DO wv=1, nb_wv
178        READ (20,*) rh_dummy, wv_nitrate(wv),
179     .        index_r_nitrate(irh,mode,wv), index_i_nitrate(irh,mode,wv)
180       ENDDO
181       ENDDO
182       CLOSE(20)
183c
184c--loop over RH
185c
186        DO IRH=1,rh_int     
187c
188        DO Nwv=1, NwvmaxSW-1+nb_lambda
189c
190c--sinon plus proche voisin
191        n_r=index_r_nitrate(irh,mode,1)
192        n_i=index_i_nitrate(irh,mode,1)
193        DO wv=1, nb_wv
194          IF (wv_nitrate(wv).LT.lambda_int(Nwv)) THEN
195            n_r=index_r_nitrate(irh,mode,wv)
196            n_i=index_i_nitrate(irh,mode,wv)
197          ENDIF 
198        ENDDO
199c
200        m=CMPLX(n_r,-n_i)
201c
202        sigma_sca=0.0
203        sigma_ext=0.0
204        sigma_abs=0.0
205        gtot=0.0
206        omegatot=0.0
207        masse = 0.0
208        volume=0.0
209        surface=0.0
210c
211        DO bin=0, Nbin !---loop on size bins
212 
213c--here r is rdry
214        r=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
215c
216        rwet=growth(IRH,mode)*r
217c
218        x=2.*pi*rwet/lambda_int(Nwv)
219        deltar=1./FLOAT(Nbin)*(log(rmax)-log(rmin))
220c
221c--rwet r*1e2 to get into cm unit
222c
223        number=0
224        DO dis=1, Ndis
225          number=number+
226     .       Ntot(dis,mode)/SQRT(2.*pi)/log(sigma_g(dis,mode))*
227     .       exp(-0.5*(log(r/(r_0(dis,mode)))/
228     .               log(sigma_g(dis,mode)))**2)
229        ENDDO
230c--dry aerosol mass
231        masse=masse  +4./3.*pi*r**3*number*deltar*rho*1.E3          !--g/m3
232        volume=volume+4./3.*pi*r**3*number*deltar
233        surface=surface+4.*pi*r**2*number*deltar
234c
235        k2=m
236        k3=CMPLX(1.0,0.0)
237
238        z2=CMPLX(x,0.0)
239        z1=m*z2
240 
241        IF (0.0.LE.x.AND.x.LE.8.) THEN
242        Nmax=INT(x+4*x**(1./3.)+1.)+2
243        ELSEIF (8..LT.x.AND.x.LT.4200.) THEN
244        Nmax=INT(x+4.05*x**(1./3.)+2.)+1
245        ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN
246        Nmax=INT(x+4*x**(1./3.)+2.)+1
247        ELSE
248        WRITE(10,*) 'x out of bound, x=', x
249        STOP
250        ENDIF
251
252        Nstart=Nmax+10
253
254C-----------loop for nu1z1, nu1z2
255
256       nu1z1(Nstart)=CMPLX(0.0,0.0)
257       nu1z2(Nstart)=CMPLX(0.0,0.0)
258       DO n=Nstart-1, 1 , -1
259       nn=CMPLX(FLOAT(n),0.0)
260       nu1z1(n)=(nn+1.)/z1 - 1./( (nn+1.)/z1 + nu1z1(n+1) ) 
261       nu1z2(n)=(nn+1.)/z2 - 1./( (nn+1.)/z2 + nu1z2(n+1) )
262       ENDDO
263
264C------------loop for nu3z2
265
266       nu3z2(0)=-I
267       DO n=1, Nmax
268       nn=CMPLX(FLOAT(n),0.0)
269       nu3z2(n)=-nn/z2 + 1./ (nn/z2 - nu3z2(n-1) ) 
270       ENDDO
271
272C-----------loop for psiz2 and ksiz2 (z2)
273       ksiz2(-1)=COS(REAL(z2))-I*SIN(REAL(z2))
274       ksiz2(0)=SIN(REAL(z2))+I*COS(REAL(z2))
275       DO n=1,Nmax
276       nn=CMPLX(FLOAT(n),0.0)
277       ksiz2(n)=(2.*nn-1.)/z2 * ksiz2(n-1) - ksiz2(n-2)
278       psiz2(n)=CMPLX(REAL(ksiz2(n)),0.0)
279       ENDDO
280
281C-----------loop for a(n) and b(n)
282   
283       DO n=1, Nmax
284        u1=k3*nu1z1(n) - k2*nu1z2(n)
285        u5=k3*nu1z1(n) - k2*nu3z2(n)
286        u6=k2*nu1z1(n) - k3*nu1z2(n)
287        u8=k2*nu1z1(n) - k3*nu3z2(n)
288        a(n)=psiz2(n)/ksiz2(n) * u1/u5
289        b(n)=psiz2(n)/ksiz2(n) * u6/u8
290       ENDDO
291
292C-----------------final loop--------------
293        Q_ext=0.0
294        Q_sca=0.0
295        g=0.0
296        DO n=Nmax-1,1,-1
297          nnn=FLOAT(n)
298          Q_ext=Q_ext+ (2.*nnn+1.) * REAL( a(n)+b(n) ) 
299          Q_sca=Q_sca+ (2.*nnn+1.) * 
300     .               REAL( a(n)*CONJG(a(n)) + b(n)*CONJG(b(n)) )
301          g=g + nnn*(nnn+2.)/(nnn+1.) * 
302     .       REAL( a(n)*CONJG(a(n+1))+b(n)*CONJG(b(n+1)) )  + 
303     .          (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n)))
304        ENDDO
305        Q_ext=2./x**2 * Q_ext
306        Q_sca=2./x**2 * Q_sca
307        Q_abs=Q_ext-Q_sca
308        IF (AIMAG(m).EQ.0.0) Q_abs=0.0
309        omega=Q_sca/Q_ext
310        g=g*4./x**2/Q_sca
311c
312        sigma_sca=sigma_sca+rwet**2*Q_sca*number*deltar
313        sigma_abs=sigma_abs+rwet**2*Q_abs*number*deltar
314        sigma_ext=sigma_ext+rwet**2*Q_ext*number*deltar
315        omegatot=omegatot+rwet**2*Q_ext*omega*number*deltar
316        gtot    =gtot+rwet**2*Q_sca*g*number*deltar 
317c
318        ENDDO   !---bin
319C------------------------------------------------------------------
320
321        sigma_sca=pi*sigma_sca
322        sigma_abs=pi*sigma_abs
323        sigma_ext=pi*sigma_ext
324        gtot=pi*gtot/sigma_sca
325        omegatot=pi*omegatot/sigma_ext
326c
327        final_g(Nwv)=gtot
328        final_w(Nwv)=omegatot
329        final_a(Nwv)=sigma_ext/masse
330c
331        ENDDO  !--loop on wavelength
332c
333c---averaging over LMDZ wavebands
334c
335        DO band=1, NbandSW
336          gcm_a(band)=0.0
337          gcm_g(band)=0.0
338          gcm_w(band)=0.0
339          gcm_weight_a(band)=0.0
340          gcm_weight_g(band)=0.0
341          gcm_weight_w(band)=0.0
342        ENDDO
343c
344c---band 1 is now in the UV, so we take the first wavelength as being representative 
345c---it doesn't matter anyway because all radiation is absorbed in the stratosphere
346       DO Nwv=1,1
347          band=1
348          gcm_a(band)=gcm_a(band)+final_a(Nwv)*weight(Nwv)
349          gcm_weight_a(band)=gcm_weight_a(band)+weight(Nwv)
350          gcm_w(band)=gcm_w(band)+
351     .                final_w(Nwv)*final_a(Nwv)*weight(Nwv)
352          gcm_weight_w(band)=gcm_weight_w(band)+
353     .                       final_a(Nwv)*weight(Nwv)
354          gcm_g(band)=gcm_g(band)+
355     .                final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
356          gcm_weight_g(band)=gcm_weight_g(band)+
357     .                       final_a(Nwv)*final_w(Nwv)*weight(Nwv)
358        ENDDO
359c
360       DO Nwv=1,NwvmaxSW-1
361c
362        IF (Nwv.LE.5) THEN          !--RRTM spectral interval 2
363          band=2
364        ELSEIF (Nwv.LE.10) THEN     !--RRTM spectral interval 3
365          band=3
366        ELSEIF (Nwv.LE.16) THEN     !--RRTM spectral interval 4
367          band=4
368        ELSEIF (Nwv.LE.21) THEN     !--RRTM spectral interval 5
369          band=5
370        ELSE                        !--RRTM spectral interval 6
371          band=6
372        ENDIF
373c
374        gcm_a(band)=gcm_a(band)+final_a(Nwv)*weight(Nwv)
375        gcm_weight_a(band)=gcm_weight_a(band)+weight(Nwv)
376        gcm_w(band)=gcm_w(band)+
377     .              final_w(Nwv)*final_a(Nwv)*weight(Nwv)
378        gcm_weight_w(band)=gcm_weight_w(band)+
379     .                     final_a(Nwv)*weight(Nwv)
380        gcm_g(band)=gcm_g(band)+
381     .              final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
382        gcm_weight_g(band)=gcm_weight_g(band)+
383     .                     final_a(Nwv)*final_w(Nwv)*weight(Nwv)
384
385        ENDDO
386c
387        DO band=1, NbandSW
388          gcm_a(band)=gcm_a(band)/gcm_weight_a(band)
389          gcm_w(band)=gcm_w(band)/gcm_weight_w(band)
390          gcm_g(band)=gcm_g(band)/gcm_weight_g(band)
391          ss_a(band,IRH)=gcm_a(band)
392          ss_w(band,IRH)=gcm_w(band)
393          ss_g(band,IRH)=gcm_g(band)
394        ENDDO
395c
396        DO nb=NbandSW+1, NbandSW+nb_lambda
397          ss_a(nb,IRH)=final_a(NwvmaxSW-1+nb-NbandSW)
398          ss_w(nb,IRH)=final_w(NwvmaxSW-1+nb-NbandSW)
399          ss_g(nb,IRH)=final_g(NwvmaxSW-1+nb-NbandSW)
400        ENDDO
401c
402        ENDDO   !--fin boucle sur RH
403c
404c--Outputs soluble
405C
406        WRITE(14,*) '  ! '//chmode(mode)
407        DO k=1, NbandSW
408          WRITE(14,951) (ss_a(k,IRH),IRH=1,rh_int)
409        ENDDO
410        WRITE(15,*) '  ! '//chmode(mode)
411        DO k=1, NbandSW
412          WRITE(15,951) (ss_g(k,IRH),IRH=1,rh_int)
413        ENDDO
414        WRITE(16,*) '  ! '//chmode(mode)
415        DO k=1, NbandSW
416          WRITE(16,951) (ss_w(k,IRH),IRH=1,rh_int)
417        ENDDO
418c
419        WRITE(34,*) '  ! extinction '//chmode(mode)
420        DO k=NbandSW+1,NbandSW+nb_lambda
421        WRITE(34,951) (ss_a(k,IRH),IRH=1,rh_int)
422        ENDDO
423        WRITE(36,*) '  ! absorption '//chmode(mode)
424        DO k=NbandSW+1,NbandSW+nb_lambda
425        WRITE(36,951) 
426     .    ((1.0-ss_w(k,IRH))*ss_a(k,IRH),IRH=1,rh_int)
427        ENDDO
428c
429c--Outputs insoluble
430C 
431        IF (mode.EQ.1) THEN
432c
433        WRITE(24,*) '  ! Nitrate insoluble'
434        WRITE(24,953) (ss_a(k,1),k=1,NbandSW)
435        WRITE(25,*) '  ! Nitrate insoluble'
436        WRITE(25,953) (ss_g(k,1),k=1,NbandSW)
437        WRITE(26,*) '  ! Nitrate insoluble'
438        WRITE(26,953) (ss_w(k,1),k=1,NbandSW)
439c
440        WRITE(35,*) '  ! extinction Nitrate insoluble'
441        WRITE(35,954) (ss_a(k,1),k=NbandSW+1,NbandSW+nb_lambda)
442        WRITE(37,*) '  ! absorption Nitrate insoluble'
443        WRITE(37,954)
444     .     ((1.0-ss_w(k,1))*ss_a(k,1),k=NbandSW+1,NbandSW+nb_lambda)
445c
446        ENDIF
447c
448        ENDDO !--boucle sur les modes
449
450951     FORMAT(1X,12(F6.3,','),' &')
451953     FORMAT(1X,5(F6.3,','),F6.3,' /')
452954     FORMAT(1X,4(F6.3,','),F6.3,' /')
453c
454101     FORMAT(F10.6,E13.6)
455c
456        CLOSE(14)
457        CLOSE(15)
458        CLOSE(16)
459c
460        CLOSE(24)
461        CLOSE(25)
462        CLOSE(26)
463c
464        CLOSE(34)
465        CLOSE(35)
466        CLOSE(36)
467        CLOSE(37)
468c
469        END
Note: See TracBrowser for help on using the repository browser.