source: LMDZ5/trunk/libf/phylmd/StratAer/miecalc_aer.F90 @ 4521

Last change on this file since 4521 was 2948, checked in by oboucher, 7 years ago

Nstart is increased for the downward recurrence for more accurate results.
The calculations for the very small x (small particles and large wavelengths) has accuracy problem
File size: 22.9 KB
Line 
1SUBROUTINE MIECALC_AER(tau_strat, piz_strat, cg_strat, tau_strat_wave, tau_lw_abs_rrtm, paprs, debut)
2
3!-------Mie computations for a size distribution
4!       of homogeneous spheres.
5!
6!==========================================================
7!--Ref : Toon and Ackerman, Applied Optics, 1981
8!        Stephens, CSIRO, 1979
9! Attention : surdimensionement des tableaux
10! to be compiled with double precision option (-r8 on Sun)
11! AUTHOR: Olivier Boucher, Christoph Kleinschmitt
12!-------SIZE distribution properties----------------
13!--sigma_g : geometric standard deviation
14!--r_0     : geometric number mean radius (um)/modal radius
15!--Ntot    : total concentration in m-3
16
17  USE phys_local_var_mod, ONLY: tr_seri, mdw, alpha_bin, piz_bin, cg_bin
18  USE aerophys
19  USE aero_mod
20  USE infotrac, ONLY : nbtr, nbtr_bin, nbtr_sulgas, id_SO2_strat
21  USE dimphy
22  USE YOMCST  , ONLY : RG, RPI
23  USE mod_phys_lmdz_para, only: gather, scatter, bcast
24  USE mod_grid_phy_lmdz, ONLY : klon_glo
25  USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
26  USE print_control_mod, ONLY: prt_level, lunout
27
28  IMPLICIT NONE
29
30! Variable input
31  LOGICAL,INTENT(IN) :: debut   ! le flag de l'initialisation de la physique
32  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
33
34! Stratospheric aerosols optical properties
35  REAL, DIMENSION(klon,klev,nbands_sw_rrtm) :: tau_strat, piz_strat, cg_strat
36  REAL, DIMENSION(klon,klev,nwave_sw+nwave_lw) :: tau_strat_wave
37  REAL, DIMENSION(klon,klev,nbands_lw_rrtm) :: tau_lw_abs_rrtm
38
39!!  REAL,DIMENSION(klon_glo,klev,nbtr)     :: tr_seri_glo         ! Concentration Traceur [U/KgA] 
40
41! local variables
42  REAL Ntot
43  PARAMETER (Ntot=1.0)
44  LOGICAL, PARAMETER :: refr_ind_interpol = .TRUE. ! set interpolation of refractive index
45  REAL r_0    ! aerosol particle radius [m]
46  INTEGER bin_number, ilon, ilev
47  REAL masse,volume,surface
48  REAL rmin, rmax    !----integral bounds in  m
49
50!-------------------------------------
51
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 nnn
60  COMPLEX nn
61  REAL Q_ext, Q_abs, Q_sca, g, omega   !--parameters for radius r
62  REAL x, x_old  !--size parameter
63  REAL r, r_lower, r_upper  !--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, it
71  PARAMETER (Nbin=10)
72  LOGICAL smallx
73
74!---wavelengths STREAMER
75  INTEGER Nwv, NwvmaxSW
76  PARAMETER (NwvmaxSW=24)
77  REAL lambda(1:NwvmaxSW+1)
78  DATA lambda/0.28E-6, 0.30E-6, 0.33E-6, 0.36E-6, 0.40E-6, &
79              0.44E-6, 0.48E-6, 0.52E-6, 0.57E-6, 0.64E-6, &
80              0.69E-6, 0.75E-6, 0.78E-6, 0.87E-6, 1.00E-6, &
81              1.10E-6, 1.19E-6, 1.28E-6, 1.53E-6, 1.64E-6, &
82              2.13E-6, 2.38E-6, 2.91E-6, 3.42E-6, 4.00E-6/
83
84!---wavelengths de references
85!---be careful here the 5th wavelength is 1020 nm
86  INTEGER nb
87  REAL lambda_ref(nwave_sw+nwave_lw)
88  DATA lambda_ref /0.443E-6,0.550E-6,0.670E-6,  &
89                   0.765E-6,1.020E-6,10.E-6/
90
91!--LW
92  INTEGER NwvmaxLW
93  PARAMETER (NwvmaxLW=500)
94  REAL Tb, hh, cc, kb
95  PARAMETER (Tb=220.0, hh=6.62607e-34)
96  PARAMETER (cc=2.99792e8, kb=1.38065e-23)
97
98!---TOA fluxes - Streamer Cs
99  REAL weight(1:NwvmaxSW), weightLW
100!c        DATA weight/0.839920E1, 0.231208E2, 0.322393E2, 0.465058E2,
101!c     .              0.678199E2, 0.798964E2, 0.771359E2, 0.888472E2,
102!c     .              0.115281E3, 0.727565E2, 0.816992E2, 0.336172E2,
103!c     .              0.914603E2, 0.112706E3, 0.658840E2, 0.524470E2,
104!c     .              0.391067E2, 0.883864E2, 0.276672E2, 0.681812E2,
105!c     .              0.190966E2, 0.250766E2, 0.128704E2, 0.698720E1/
106!---TOA fluxes - Tad
107  DATA weight/ 4.20, 11.56, 16.12, 23.25, 33.91, 39.95, 38.57, &
108              44.42, 57.64, 29.36, 47.87, 16.81, 45.74, 56.35, &
109              32.94, 26.22, 19.55, 44.19, 13.83, 34.09,  9.55, &
110              12.54,  6.44,  3.49/
111!C---BOA fluxes - Tad
112!c        DATA weight/ 0.01,  4.05, 9.51,  15.99, 26.07, 33.10, 33.07,
113!c     .              39.91, 52.67, 27.89, 43.60, 13.67, 42.22, 40.12,
114!c     .              32.70, 14.44, 19.48, 14.23, 13.43, 16.42,  8.33,
115!c     .               0.95,  0.65, 2.76/
116
117  REAL lambda_int(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW), ll
118  REAL dlambda_int(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW), dl
119
120  REAL n_r(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW)
121  REAL n_i(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW)
122
123  REAL ilambda, ilambda_prev, ilambda_max, ilambda_min
124  REAL n_r_h2so4, n_i_h2so4
125  REAL n_r_h2so4_prev, n_i_h2so4_prev
126
127  REAL final_a(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW)
128  REAL final_g(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW)
129  REAL final_w(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW)
130
131  INTEGER band, bandSW, bandLW, wavenumber
132
133!---wavelengths SW RRTM
134  REAL wv_rrtm_SW(nbands_sw_rrtm+1)
135  DATA wv_rrtm_SW/  0.185E-6, 0.25E-6, 0.44E-6, 0.69E-6,  &
136                     1.19E-6, 2.38E-6, 4.00E-6/
137
138!---wavenumbers and wavelengths LW RRTM
139  REAL wn_rrtm(nbands_lw_rrtm+1), wv_rrtm(nbands_lw_rrtm+1)
140  DATA wn_rrtm/  10.,  250.,  500.,  630.,  700.,  820.,  &
141                980., 1080., 1180., 1390., 1480., 1800.,  &
142               2080., 2250., 2380., 2600., 3000./
143
144!--GCM results
145  REAL gcm_a(nbands_sw_rrtm+nbands_lw_rrtm)
146  REAL gcm_g(nbands_sw_rrtm+nbands_lw_rrtm)
147  REAL gcm_w(nbands_sw_rrtm+nbands_lw_rrtm)
148  REAL gcm_weight_a(nbands_sw_rrtm+nbands_lw_rrtm)
149  REAL gcm_weight_g(nbands_sw_rrtm+nbands_lw_rrtm)
150  REAL gcm_weight_w(nbands_sw_rrtm+nbands_lw_rrtm)
151
152  REAL ss_a(nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw)
153  REAL ss_w(nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw)
154  REAL ss_g(nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw)
155
156  INTEGER, PARAMETER :: nb_lambda_h2so4=62
157  REAL, DIMENSION (nb_lambda_h2so4,4) :: ref_ind
158  !-- fichier h2so4_0.75_300.00_hummel_1988_p_q.dat
159  ! -- wavenumber (cm-1), wavelength (um), n_r, n_i
160  DATA ref_ind /                                &
161   200.000,   50.0000,   2.01000,   6.5000E-01, &
162   250.000,   40.0000,   1.94000,   6.3000E-01, &
163   285.714,   35.0000,   1.72000,   5.2000E-01, &
164   333.333,   30.0000,   1.73000,   2.9000E-01, &
165   358.423,   27.9000,   1.78000,   2.5000E-01, &
166   400.000,   25.0000,   1.84000,   2.4000E-01, &
167   444.444,   22.5000,   1.82000,   2.9000E-01, &
168   469.484,   21.3000,   1.79000,   2.5000E-01, &
169   500.000,   20.0000,   1.81000,   2.3000E-01, &
170   540.541,   18.5000,   1.92700,   3.0200E-01, &
171   555.556,   18.0000,   1.95000,   4.1000E-01, &
172   581.395,   17.2000,   1.72400,   5.9000E-01, &
173   609.756,   16.4000,   1.52000,   4.1400E-01, &
174   666.667,   15.0000,   1.59000,   2.1100E-01, &
175   675.676,   14.8000,   1.61000,   2.0500E-01, &
176   714.286,   14.0000,   1.64000,   1.9500E-01, &
177   769.231,   13.0000,   1.69000,   1.9500E-01, &
178   800.000,   12.5000,   1.74000,   1.9800E-01, &
179   869.565,   11.5000,   1.89000,   3.7400E-01, &
180   909.091,   11.0000,   1.67000,   4.8500E-01, &
181   944.198,   10.5910,   1.72000,   3.4000E-01, &
182  1000.000,   10.0000,   1.89000,   4.5500E-01, &
183  1020.408,    9.8000,   1.91000,   6.8000E-01, &
184  1052.632,    9.5000,   1.67000,   7.5000E-01, &
185  1086.957,    9.2000,   1.60000,   5.8600E-01, &
186  1111.111,    9.0000,   1.65000,   6.3300E-01, &
187  1149.425,    8.7000,   1.53000,   7.7200E-01, &
188  1176.471,    8.5000,   1.37000,   7.5500E-01, &
189  1219.512,    8.2000,   1.20000,   6.4500E-01, &
190  1265.823,    7.9000,   1.14000,   4.8800E-01, &
191  1388.889,    7.2000,   1.21000,   1.7600E-01, &
192  1538.462,    6.5000,   1.37000,   1.2800E-01, &
193  1612.903,    6.2000,   1.42400,   1.6500E-01, &
194  1666.667,    6.0000,   1.42500,   1.9500E-01, &
195  1818.182,    5.5000,   1.33700,   1.8300E-01, &
196  2000.000,    5.0000,   1.36000,   1.2100E-01, &
197  2222.222,    4.5000,   1.38500,   1.2000E-01, &
198  2500.000,    4.0000,   1.39800,   1.2600E-01, &
199  2666.667,    3.7500,   1.39600,   1.3100E-01, &
200  2857.143,    3.5000,   1.37600,   1.5800E-01, &
201  2948.113,    3.3920,   1.35200,   1.5900E-01, &
202  3125.000,    3.2000,   1.31100,   1.3500E-01, &
203  3333.333,    3.0000,   1.29300,   9.5500E-02, &
204  3703.704,    2.7000,   1.30300,   5.7000E-03, &
205  4000.000,    2.5000,   1.34400,   3.7600E-03, &
206  4444.444,    2.2500,   1.37000,   1.8000E-03, &
207  5000.000,    2.0000,   1.38400,   1.2600E-03, &
208  5555.556,    1.8000,   1.39000,   5.5000E-04, &
209  6510.417,    1.5360,   1.40300,   1.3700E-04, &
210  7692.308,    1.3000,   1.41000,   1.0000E-05, &
211  9433.962,    1.0600,   1.42000,   1.5000E-06, &
212 11627.907,    0.8600,   1.42500,   1.7900E-07, &
213 14409.222,    0.6940,   1.42800,   1.9900E-08, &
214 15797.788,    0.6330,   1.42900,   1.4700E-08, &
215 18181.818,    0.5500,   1.43000,   1.0000E-08, &
216 19417.476,    0.5150,   1.43100,   1.0000E-08, &
217 20491.803,    0.4880,   1.43200,   1.0000E-08, &
218 25000.000,    0.4000,   1.44000,   1.0000E-08, &
219 29673.591,    0.3370,   1.45900,   1.0000E-08, &
220 33333.333,    0.3000,   1.46900,   1.0000E-08, &
221 40000.000,    0.2500,   1.48400,   1.0000E-08, &
222 50000.000,    0.2000,   1.49800,   1.0000E-08 /
223!---------------------------------------------------------
224
225  IF (debut) THEN   
226
227  !--initialising dry diameters to geometrically spaced mass/volume (see Jacobson 1994)
228      mdw(1)=mdwmin
229      IF (V_rat.LT.1.62) THEN ! compensate for dip in second bin for lower volume ratio
230        mdw(2)=mdw(1)*2.**(1./3.)
231        DO it=3, nbtr_bin
232          mdw(it)=mdw(it-1)*V_rat**(1./3.)
233        ENDDO
234      ELSE
235        DO it=2, nbtr_bin
236          mdw(it)=mdw(it-1)*V_rat**(1./3.)
237        ENDDO
238      ENDIF
239      PRINT *,'init mdw=', mdw
240
241    !--compute particle radius for a composition of 75% H2SO4 / 25% H2O at T=293K
242    DO bin_number=1, nbtr_bin
243      r_0=(dens_aer_dry/dens_aer_ref/0.75)**(1./3.)*mdw(bin_number)/2.
244    !--integral boundaries set to bin boundaries
245      rmin=r_0/sqrt(V_rat**(1./3.))
246      rmax=r_0*sqrt(V_rat**(1./3.))
247
248    !--set up SW
249      DO Nwv=1, NwvmaxSW
250        lambda_int(Nwv)=( lambda(Nwv)+lambda(Nwv+1) ) /2.
251      ENDDO
252
253      DO nb=1, nwave_sw+nwave_lw
254        lambda_int(NwvmaxSW+nb)=lambda_ref(nb)
255      ENDDO
256
257    !--set up LW
258    !--conversion wavenumber in cm-1 to wavelength in m
259      DO Nwv=1, nbands_lw_rrtm+1
260        wv_rrtm(Nwv)=10000./wn_rrtm(Nwv)*1.e-6
261      ENDDO
262
263      DO Nwv=1, NwvmaxLW
264        lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
265          exp( log(wv_rrtm(1))+float(Nwv-1)/float(NwvmaxLW-1)* &
266          (log(wv_rrtm(nbands_lw_rrtm+1))-log(wv_rrtm(1))) )
267      ENDDO
268
269!--computing the dlambdas
270      Nwv=1
271      dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
272      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)- &
273      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv+1)
274      DO Nwv=2, NwvmaxLW-1
275      dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
276      &  (lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv-1)- &
277      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv+1))/2.
278      ENDDO
279      Nwv=NwvmaxLW
280      dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
281      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv-1)- &
282      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)
283
284      IF (refr_ind_interpol) THEN
285
286        ilambda_max=ref_ind(1,2)/1.e6 !--in m
287        ilambda_min=ref_ind(nb_lambda_h2so4,2)/1.e6 !--in m
288        DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
289          IF (lambda_int(Nwv).GT.ilambda_max) THEN
290            !for lambda out of data range, take boundary values
291            n_r(Nwv)=ref_ind(1,3)
292            n_i(Nwv)=ref_ind(1,4)
293          ELSEIF (lambda_int(Nwv).LE.ilambda_min) THEN
294            n_r(Nwv)=ref_ind(nb_lambda_h2so4,3)
295            n_i(Nwv)=ref_ind(nb_lambda_h2so4,4)
296          ELSE
297            DO nb=2,nb_lambda_h2so4
298              ilambda=ref_ind(nb,2)/1.e6
299              ilambda_prev=ref_ind(nb-1,2)/1.e6
300              n_r_h2so4=ref_ind(nb,3)
301              n_r_h2so4_prev=ref_ind(nb-1,3)
302              n_i_h2so4=ref_ind(nb,4)
303              n_i_h2so4_prev=ref_ind(nb-1,4)
304              IF (lambda_int(Nwv).GT.ilambda.AND. &
305                lambda_int(Nwv).LE.ilambda_prev) THEN !--- linear interpolation
306                n_r(Nwv)=n_r_h2so4+(lambda_int(Nwv)-ilambda)/ &
307                     (ilambda_prev-ilambda)*(n_r_h2so4_prev-n_r_h2so4)
308                n_i(Nwv)=n_i_h2so4+(lambda_int(Nwv)-ilambda)/ &
309                     (ilambda_prev-ilambda)*(n_i_h2so4_prev-n_i_h2so4)
310              ENDIF
311            ENDDO
312          ENDIF
313        ENDDO
314
315      ELSE  !-- no refr_ind_interpol, closest neighbour from upper wavelength
316
317        DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
318          n_r(Nwv)=ref_ind(1,3)
319          n_i(Nwv)=ref_ind(1,4)
320          DO nb=2,nb_lambda_h2so4
321            IF (ref_ind(nb,2)/1.e6.GT.lambda_int(Nwv)) THEN !--- step function
322              n_r(Nwv)=ref_ind(nb,3)
323              n_i(Nwv)=ref_ind(nb,4)
324            ENDIF 
325          ENDDO
326        ENDDO
327      ENDIF
328
329    !---Loop on wavelengths
330      DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
331
332      m=CMPLX(n_r(Nwv),-n_i(Nwv))
333
334      I=CMPLX(0.,1.)
335
336      sigma_sca=0.0
337      sigma_ext=0.0
338      sigma_abs=0.0
339      gtot=0.0
340      omegatot=0.0
341      masse = 0.0
342      volume=0.0
343      surface=0.0
344
345      DO bin=1, Nbin !---loop on size bins
346
347      r_lower=exp(log(rmin)+FLOAT(bin-1)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
348      r_upper=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
349      deltar=r_upper-r_lower
350
351      r=sqrt(r_lower*r_upper)
352      x=2.*RPI*r/lambda_int(Nwv)
353
354!we impose a minimum value for x and extrapolate quantities for small x values
355      smallx = .FALSE.
356      IF (x.LT.0.001) THEN
357        smallx = .TRUE.
358        x_old = x
359        x = 0.001
360      ENDIF
361
362      number=Ntot*deltar/(rmax-rmin) !dN/dr constant over tracer bin
363!      masse=masse  +4./3.*RPI*(r**3)*number*deltar*ropx*1.E3  !--g/m3
364      volume=volume+4./3.*RPI*(r**3)*number*deltar
365      surface=surface+4.*RPI*r**2*number*deltar
366
367      k2=m
368      k3=CMPLX(1.0,0.0)
369
370      z2=CMPLX(x,0.0)
371      z1=m*z2
372
373      IF (0.0.LE.x.AND.x.LE.8.) THEN
374        Nmax=INT(x+4*x**(1./3.)+1.)+2
375      ELSEIF (8..LT.x.AND.x.LT.4200.) THEN
376        Nmax=INT(x+4.05*x**(1./3.)+2.)+1
377      ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN
378        Nmax=INT(x+4*x**(1./3.)+2.)+1
379      ELSE
380        PRINT *, 'x out of bound, x=', x
381        STOP
382      ENDIF
383
384      Nstart=Nmax+100
385
386    !-----------loop for nu1z1, nu1z2
387
388      nu1z1(Nstart)=CMPLX(0.0,0.0)
389      nu1z2(Nstart)=CMPLX(0.0,0.0)
390      DO n=Nstart-1, 1 , -1
391        nn=CMPLX(FLOAT(n),0.0)
392        nu1z1(n)=(nn+1.)/z1 - 1./( (nn+1.)/z1 + nu1z1(n+1) )
393        nu1z2(n)=(nn+1.)/z2 - 1./( (nn+1.)/z2 + nu1z2(n+1) )
394      ENDDO
395
396    !------------loop for nu3z2
397
398      nu3z2(0)=-I
399      DO n=1, Nmax
400        nn=CMPLX(FLOAT(n),0.0)
401        nu3z2(n)=-nn/z2 + 1./ (nn/z2 - nu3z2(n-1) )
402      ENDDO
403
404    !-----------loop for psiz2 and ksiz2 (z2)
405      ksiz2(-1)=COS(REAL(z2))-I*SIN(REAL(z2))
406      ksiz2(0)=SIN(REAL(z2))+I*COS(REAL(z2))
407      DO n=1,Nmax
408       nn=CMPLX(FLOAT(n),0.0)
409       ksiz2(n)=(2.*nn-1.)/z2 * ksiz2(n-1) - ksiz2(n-2)
410       psiz2(n)=CMPLX(REAL(ksiz2(n)),0.0)
411      ENDDO
412
413    !-----------loop for a(n) and b(n)
414
415      DO n=1, Nmax
416        u1=k3*nu1z1(n) - k2*nu1z2(n)
417        u5=k3*nu1z1(n) - k2*nu3z2(n)
418        u6=k2*nu1z1(n) - k3*nu1z2(n)
419        u8=k2*nu1z1(n) - k3*nu3z2(n)
420        a(n)=psiz2(n)/ksiz2(n) * u1/u5
421        b(n)=psiz2(n)/ksiz2(n) * u6/u8
422      ENDDO
423
424    !-----------------final loop--------------
425      Q_ext=0.0
426      Q_sca=0.0
427      g=0.0
428
429      DO n=Nmax-1,1,-1
430        nnn=FLOAT(n)
431        Q_ext=Q_ext+ (2.*nnn+1.) * REAL( a(n)+b(n) )
432        Q_sca=Q_sca+ (2.*nnn+1.) *  &
433                   REAL( a(n)*CONJG(a(n)) + b(n)*CONJG(b(n)) )
434        g=g + nnn*(nnn+2.)/(nnn+1.) *   &
435           REAL( a(n)*CONJG(a(n+1))+b(n)*CONJG(b(n+1)) )  +   &
436              (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n)))
437      ENDDO
438
439      Q_ext=2./x**2 * Q_ext
440      Q_sca=2./x**2 * Q_sca
441    !--extrapolation in case of small x values
442      IF (smallx) THEN
443        Q_ext = x_old/x * Q_ext
444        Q_sca = x_old/x * Q_sca
445      ENDIF
446
447      Q_abs=Q_ext-Q_sca
448
449      IF (AIMAG(m).EQ.0.0) Q_abs=0.0
450      omega=Q_sca/Q_ext
451
452    ! g is wrong in the smallx case (but that does not matter as long as we ignore LW scattering)
453      g=g*4./x**2/Q_sca
454
455      sigma_sca=sigma_sca+r**2*Q_sca*number
456      sigma_abs=sigma_abs+r**2*Q_abs*number
457      sigma_ext=sigma_ext+r**2*Q_ext*number
458      omegatot=omegatot+r**2*Q_ext*omega*number
459      gtot    =gtot+r**2*Q_sca*g*number
460
461      ENDDO   !---bin
462    !------------------------------------------------------------------
463
464      sigma_sca=RPI*sigma_sca
465      sigma_abs=RPI*sigma_abs
466      sigma_ext=RPI*sigma_ext
467      gtot=RPI*gtot/sigma_sca
468      omegatot=RPI*omegatot/sigma_ext
469
470      final_g(Nwv)=gtot
471      final_w(Nwv)=omegatot
472!      final_a(Nwv)=sigma_ext/masse
473      final_a(Nwv)=sigma_ext !extinction/absorption cross section per particle
474
475      ENDDO  !--loop on wavelength
476
477    !---averaging over LMDZ wavebands
478
479      DO band=1, nbands_sw_rrtm+nbands_lw_rrtm
480        gcm_a(band)=0.0
481        gcm_g(band)=0.0
482        gcm_w(band)=0.0
483        gcm_weight_a(band)=0.0
484        gcm_weight_g(band)=0.0
485        gcm_weight_w(band)=0.0
486      ENDDO
487
488    !---band 1 is now in the UV, so we take the first wavelength as being representative
489      DO Nwv=1,1
490        bandSW=1
491        gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv)
492        gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv)
493        gcm_w(bandSW)=gcm_w(bandSW)+final_w(Nwv)*final_a(Nwv)*weight(Nwv)
494        gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+final_a(Nwv)*weight(Nwv)
495        gcm_g(bandSW)=gcm_g(bandSW)+final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
496        gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+final_a(Nwv)*final_w(Nwv)*weight(Nwv)
497      ENDDO
498
499      DO Nwv=1,NwvmaxSW
500
501        IF (lambda_int(Nwv).LE.wv_rrtm_SW(3)) THEN      !--RRTM spectral interval 2
502          bandSW=2
503        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(4)) THEN  !--RRTM spectral interval 3
504          bandSW=3
505        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(5)) THEN  !--RRTM spectral interval 4
506          bandSW=4
507        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(6)) THEN  !--RRTM spectral interval 5
508          bandSW=5
509        ELSE                                            !--RRTM spectral interval 6
510          bandSW=6
511        ENDIF
512
513        gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv)
514        gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv)
515        gcm_w(bandSW)=gcm_w(bandSW)+final_w(Nwv)*final_a(Nwv)*weight(Nwv)
516        gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+final_a(Nwv)*weight(Nwv)
517        gcm_g(bandSW)=gcm_g(bandSW)+final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
518        gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+final_a(Nwv)*final_w(Nwv)*weight(Nwv)
519
520      ENDDO
521
522      DO Nwv=NwvmaxSW+nwave_sw+nwave_lw+1,NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
523        ll=lambda_int(Nwv)
524        dl=dlambda_int(Nwv)
525        weightLW=1./ll**5./(exp(hh*cc/kb/Tb/ll)-1.)*dl
526        bandLW=1  !--default value starting from the highest lambda
527        DO band=2, nbands_lw_rrtm
528          IF (ll.LT.wv_rrtm(band)) THEN   !--as long as
529            bandLW=band
530          ENDIF
531        ENDDO
532        gcm_a(nbands_sw_rrtm+bandLW)=gcm_a(nbands_sw_rrtm+bandLW)+final_a(Nwv)*   &
533             (1.-final_w(Nwv))*weightLW
534        gcm_weight_a(nbands_sw_rrtm+bandLW)=gcm_weight_a(nbands_sw_rrtm+bandLW)+weightLW
535
536        gcm_w(nbands_sw_rrtm+bandLW)=gcm_w(nbands_sw_rrtm+bandLW)+final_w(Nwv)*   &
537             final_a(Nwv)*weightLW
538        gcm_weight_w(nbands_sw_rrtm+bandLW)=gcm_weight_w(nbands_sw_rrtm+bandLW)+  &
539             final_a(Nwv)*weightLW
540
541        gcm_g(nbands_sw_rrtm+bandLW)=gcm_g(nbands_sw_rrtm+bandLW)+final_g(Nwv)*   &
542             final_a(Nwv)*final_w(Nwv)*weightLW
543        gcm_weight_g(nbands_sw_rrtm+bandLW)=gcm_weight_g(nbands_sw_rrtm+bandLW)+  &
544             final_a(Nwv)*final_w(Nwv)*weightLW
545      ENDDO
546
547      DO band=1, nbands_sw_rrtm+nbands_lw_rrtm
548        gcm_a(band)=gcm_a(band)/gcm_weight_a(band)
549        gcm_w(band)=gcm_w(band)/gcm_weight_w(band)
550        gcm_g(band)=gcm_g(band)/gcm_weight_g(band)
551        ss_a(band)=gcm_a(band)
552        ss_w(band)=gcm_w(band)
553        ss_g(band)=gcm_g(band)
554      ENDDO
555
556      DO nb=1, nwave_sw+nwave_lw
557        ss_a(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_a(NwvmaxSW+nb)
558        ss_w(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_w(NwvmaxSW+nb)
559        ss_g(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_g(NwvmaxSW+nb)
560      ENDDO
561
562      DO nb=1,nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw
563        alpha_bin(nb,bin_number)=ss_a(nb) !extinction/absorption cross section per particle
564        piz_bin(nb,bin_number)=ss_w(nb)
565        cg_bin(nb,bin_number)=ss_g(nb)
566      ENDDO
567
568    ENDDO !loop over tracer bins
569
570!!$OMP END MASTER
571!  CALL bcast(alpha_bin)
572!  CALL bcast(piz_bin)
573!  CALL bcast(cg_bin)
574!!$OMP BARRIER
575
576    !set to default values at first time step (tr_seri still zero)
577    tau_strat(:,:,:)=1.e-15
578    piz_strat(:,:,:)=1.0
579    cg_strat(:,:,:)=0.0
580    tau_lw_abs_rrtm(:,:,:)=1.e-15
581    tau_strat_wave(:,:,:)=1.e-15
582
583  ELSE  !-- not debut
584
585  !--compute optical properties of actual size distribution (from tr_seri)
586    DO ilon=1,klon
587      DO ilev=1, klev
588        DO nb=1,nbands_sw_rrtm
589          tau_strat(ilon,ilev,nb)=0.0
590          DO bin_number=1, nbtr_bin
591            tau_strat(ilon,ilev,nb)=tau_strat(ilon,ilev,nb)+alpha_bin(nb,bin_number) &
592                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
593          ENDDO
594
595          piz_strat(ilon,ilev,nb)=0.0
596          DO bin_number=1, nbtr_bin
597            piz_strat(ilon,ilev,nb)=piz_strat(ilon,ilev,nb)+piz_bin(nb,bin_number)*alpha_bin(nb,bin_number) &
598                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
599          ENDDO
600          piz_strat(ilon,ilev,nb)=piz_strat(ilon,ilev,nb)/MAX(tau_strat(ilon,ilev,nb),1.e-15)
601
602          cg_strat(ilon,ilev,nb)=0.0
603          DO bin_number=1, nbtr_bin
604            cg_strat(ilon,ilev,nb)=cg_strat(ilon,ilev,nb)+cg_bin(nb,bin_number)*piz_bin(nb,bin_number)*alpha_bin(nb,bin_number) &
605                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
606          ENDDO
607          cg_strat(ilon,ilev,nb)=cg_strat(ilon,ilev,nb)/MAX(tau_strat(ilon,ilev,nb)*piz_strat(ilon,ilev,nb),1.e-15)
608        ENDDO
609        DO nb=1,nbands_lw_rrtm
610          tau_lw_abs_rrtm(ilon,ilev,nb)=0.0
611          DO bin_number=1, nbtr_bin
612            tau_lw_abs_rrtm(ilon,ilev,nb)=tau_lw_abs_rrtm(ilon,ilev,nb)+alpha_bin(nbands_sw_rrtm+nb,bin_number) &
613                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
614          ENDDO
615        ENDDO
616        DO nb=1,nwave_sw+nwave_lw
617          tau_strat_wave(ilon,ilev,nb)=0.0
618          DO bin_number=1, nbtr_bin
619            tau_strat_wave(ilon,ilev,nb)=tau_strat_wave(ilon,ilev,nb)+alpha_bin(nbands_sw_rrtm+nbands_lw_rrtm+nb,bin_number) &
620                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
621          ENDDO
622        ENDDO
623      ENDDO
624    ENDDO
625
626  ENDIF !debut
627
628END SUBROUTINE MIECALC_AER
Note: See TracBrowser for help on using the repository browser.