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