source: LMDZ6/trunk/libf/phylmd/StratAer/miecalc_aer.F90 @ 3667

Last change on this file since 3667 was 3667, checked in by oboucher, 4 years ago

Correction of a big big bug on refractive indiex and Mie calc in the StratAer? module
  • Property svn:keywords set to Id
File size: 23.1 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  !-- fichier h2so4_0.75_300.00_hummel_1988_p_q.dat
157  ! -- wavenumber (cm-1), wavelength (um), n_r, n_i
158  INTEGER, PARAMETER :: nb_lambda_h2so4=62
159  REAL, DIMENSION (nb_lambda_h2so4,4) :: ref_ind
160
161!---------------------------------------------------------
162
163  IF (debut) THEN   
164
165     ref_ind = RESHAPE( (/ &
166      200.000,   50.0000,   2.01000,   6.5000E-01, &
167      250.000,   40.0000,   1.94000,   6.3000E-01, &
168      285.714,   35.0000,   1.72000,   5.2000E-01, &
169      333.333,   30.0000,   1.73000,   2.9000E-01, &
170      358.423,   27.9000,   1.78000,   2.5000E-01, &
171      400.000,   25.0000,   1.84000,   2.4000E-01, &
172      444.444,   22.5000,   1.82000,   2.9000E-01, &
173      469.484,   21.3000,   1.79000,   2.5000E-01, &
174      500.000,   20.0000,   1.81000,   2.3000E-01, &
175      540.541,   18.5000,   1.92700,   3.0200E-01, &
176      555.556,   18.0000,   1.95000,   4.1000E-01, &
177      581.395,   17.2000,   1.72400,   5.9000E-01, &
178      609.756,   16.4000,   1.52000,   4.1400E-01, &
179      666.667,   15.0000,   1.59000,   2.1100E-01, &
180      675.676,   14.8000,   1.61000,   2.0500E-01, &
181      714.286,   14.0000,   1.64000,   1.9500E-01, &
182      769.231,   13.0000,   1.69000,   1.9500E-01, &
183      800.000,   12.5000,   1.74000,   1.9800E-01, &
184      869.565,   11.5000,   1.89000,   3.7400E-01, &
185      909.091,   11.0000,   1.67000,   4.8500E-01, &
186      944.198,   10.5910,   1.72000,   3.4000E-01, &
187     1000.000,   10.0000,   1.89000,   4.5500E-01, &
188     1020.408,    9.8000,   1.91000,   6.8000E-01, &
189     1052.632,    9.5000,   1.67000,   7.5000E-01, &
190     1086.957,    9.2000,   1.60000,   5.8600E-01, &
191     1111.111,    9.0000,   1.65000,   6.3300E-01, &
192     1149.425,    8.7000,   1.53000,   7.7200E-01, &
193     1176.471,    8.5000,   1.37000,   7.5500E-01, &
194     1219.512,    8.2000,   1.20000,   6.4500E-01, &
195     1265.823,    7.9000,   1.14000,   4.8800E-01, &
196     1388.889,    7.2000,   1.21000,   1.7600E-01, &
197     1538.462,    6.5000,   1.37000,   1.2800E-01, &
198     1612.903,    6.2000,   1.42400,   1.6500E-01, &
199     1666.667,    6.0000,   1.42500,   1.9500E-01, &
200     1818.182,    5.5000,   1.33700,   1.8300E-01, &
201     2000.000,    5.0000,   1.36000,   1.2100E-01, &
202     2222.222,    4.5000,   1.38500,   1.2000E-01, &
203     2500.000,    4.0000,   1.39800,   1.2600E-01, &
204     2666.667,    3.7500,   1.39600,   1.3100E-01, &
205     2857.143,    3.5000,   1.37600,   1.5800E-01, &
206     2948.113,    3.3920,   1.35200,   1.5900E-01, &
207     3125.000,    3.2000,   1.31100,   1.3500E-01, &
208     3333.333,    3.0000,   1.29300,   9.5500E-02, &
209     3703.704,    2.7000,   1.30300,   5.7000E-03, &
210     4000.000,    2.5000,   1.34400,   3.7600E-03, &
211     4444.444,    2.2500,   1.37000,   1.8000E-03, &
212     5000.000,    2.0000,   1.38400,   1.2600E-03, &
213     5555.556,    1.8000,   1.39000,   5.5000E-04, &
214     6510.417,    1.5360,   1.40300,   1.3700E-04, &
215     7692.308,    1.3000,   1.41000,   1.0000E-05, &
216     9433.962,    1.0600,   1.42000,   1.5000E-06, &
217    11627.907,    0.8600,   1.42500,   1.7900E-07, &
218    14409.222,    0.6940,   1.42800,   1.9900E-08, &
219    15797.788,    0.6330,   1.42900,   1.4700E-08, &
220    18181.818,    0.5500,   1.43000,   1.0000E-08, &
221    19417.476,    0.5150,   1.43100,   1.0000E-08, &
222    20491.803,    0.4880,   1.43200,   1.0000E-08, &
223    25000.000,    0.4000,   1.44000,   1.0000E-08, &
224    29673.591,    0.3370,   1.45900,   1.0000E-08, &
225    33333.333,    0.3000,   1.46900,   1.0000E-08, &
226    40000.000,    0.2500,   1.48400,   1.0000E-08, &
227    50000.000,    0.2000,   1.49800,   1.0000E-08 /), (/nb_lambda_h2so4,4/), order=(/2,1/) )
228
229  !--initialising dry diameters to geometrically spaced mass/volume (see Jacobson 1994)
230    mdw(1)=mdwmin
231    IF (V_rat.LT.1.62) THEN ! compensate for dip in second bin for lower volume ratio
232      mdw(2)=mdw(1)*2.**(1./3.)
233      DO it=3, nbtr_bin
234        mdw(it)=mdw(it-1)*V_rat**(1./3.)
235      ENDDO
236    ELSE
237      DO it=2, nbtr_bin
238        mdw(it)=mdw(it-1)*V_rat**(1./3.)
239      ENDDO
240    ENDIF
241    WRITE(lunout,*) 'init mdw=', mdw
242
243    !--compute particle radius for a composition of 75% H2SO4 / 25% H2O at T=293K
244    DO bin_number=1, nbtr_bin
245      r_0=(dens_aer_dry/dens_aer_ref/0.75)**(1./3.)*mdw(bin_number)/2.
246    !--integral boundaries set to bin boundaries
247      rmin=r_0/sqrt(V_rat**(1./3.))
248      rmax=r_0*sqrt(V_rat**(1./3.))
249
250    !--set up SW
251      DO Nwv=1, NwvmaxSW
252        lambda_int(Nwv)=( lambda(Nwv)+lambda(Nwv+1) ) /2.
253      ENDDO
254
255      DO nb=1, nwave_sw+nwave_lw
256        lambda_int(NwvmaxSW+nb)=lambda_ref(nb)
257      ENDDO
258
259    !--set up LW
260    !--conversion wavenumber in cm-1 to wavelength in m
261      DO Nwv=1, nbands_lw_rrtm+1
262        wv_rrtm(Nwv)=10000./wn_rrtm(Nwv)*1.e-6
263      ENDDO
264
265      DO Nwv=1, NwvmaxLW
266        lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
267          exp( log(wv_rrtm(1))+float(Nwv-1)/float(NwvmaxLW-1)* &
268          (log(wv_rrtm(nbands_lw_rrtm+1))-log(wv_rrtm(1))) )
269      ENDDO
270
271!--computing the dlambdas
272      Nwv=1
273      dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
274      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)- &
275      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv+1)
276      DO Nwv=2, NwvmaxLW-1
277      dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
278      &  (lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv-1)- &
279      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv+1))/2.
280      ENDDO
281      Nwv=NwvmaxLW
282      dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
283      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv-1)- &
284      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)
285
286      IF (refr_ind_interpol) THEN
287 
288        ilambda_max=ref_ind(1,2)/1.e6 !--in m
289        ilambda_min=ref_ind(nb_lambda_h2so4,2)/1.e6 !--in m
290        DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
291          IF (lambda_int(Nwv).GT.ilambda_max) THEN
292            !for lambda out of data range, take boundary values
293            n_r(Nwv)=ref_ind(1,3)
294            n_i(Nwv)=ref_ind(1,4)
295          ELSEIF (lambda_int(Nwv).LE.ilambda_min) THEN
296            n_r(Nwv)=ref_ind(nb_lambda_h2so4,3)
297            n_i(Nwv)=ref_ind(nb_lambda_h2so4,4)
298          ELSE
299            DO nb=2,nb_lambda_h2so4
300              ilambda=ref_ind(nb,2)/1.e6
301              ilambda_prev=ref_ind(nb-1,2)/1.e6
302              n_r_h2so4=ref_ind(nb,3)
303              n_r_h2so4_prev=ref_ind(nb-1,3)
304              n_i_h2so4=ref_ind(nb,4)
305              n_i_h2so4_prev=ref_ind(nb-1,4)
306              IF (lambda_int(Nwv).GT.ilambda.AND. &
307                lambda_int(Nwv).LE.ilambda_prev) THEN !--- linear interpolation
308                n_r(Nwv)=n_r_h2so4+(lambda_int(Nwv)-ilambda)/ &
309                     (ilambda_prev-ilambda)*(n_r_h2so4_prev-n_r_h2so4)
310                n_i(Nwv)=n_i_h2so4+(lambda_int(Nwv)-ilambda)/ &
311                     (ilambda_prev-ilambda)*(n_i_h2so4_prev-n_i_h2so4)
312              ENDIF
313            ENDDO
314          ENDIF
315        ENDDO
316
317      ELSE  !-- no refr_ind_interpol, closest neighbour from upper wavelength
318
319        DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
320          n_r(Nwv)=ref_ind(1,3)
321          n_i(Nwv)=ref_ind(1,4)
322          DO nb=2,nb_lambda_h2so4
323            IF (ref_ind(nb,2)/1.e6.GT.lambda_int(Nwv)) THEN !--- step function
324              n_r(Nwv)=ref_ind(nb,3)
325              n_i(Nwv)=ref_ind(nb,4)
326            ENDIF 
327          ENDDO
328        ENDDO
329      ENDIF
330
331    !---Loop on wavelengths
332      DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
333
334      m=CMPLX(n_r(Nwv),-n_i(Nwv))
335
336      I=CMPLX(0.,1.)
337
338      sigma_sca=0.0
339      sigma_ext=0.0
340      sigma_abs=0.0
341      gtot=0.0
342      omegatot=0.0
343      masse = 0.0
344      volume=0.0
345      surface=0.0
346
347      DO bin=1, Nbin !---loop on size bins
348
349      r_lower=exp(log(rmin)+FLOAT(bin-1)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
350      r_upper=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
351      deltar=r_upper-r_lower
352
353      r=sqrt(r_lower*r_upper)
354      x=2.*RPI*r/lambda_int(Nwv)
355
356!we impose a minimum value for x and extrapolate quantities for small x values
357      smallx = .FALSE.
358      IF (x.LT.0.001) THEN
359        smallx = .TRUE.
360        x_old = x
361        x = 0.001
362      ENDIF
363
364      number=Ntot*deltar/(rmax-rmin) !dN/dr constant over tracer bin
365!      masse=masse  +4./3.*RPI*(r**3)*number*deltar*ropx*1.E3  !--g/m3
366      volume=volume+4./3.*RPI*(r**3)*number*deltar
367      surface=surface+4.*RPI*r**2*number*deltar
368
369      k2=m
370      k3=CMPLX(1.0,0.0)
371
372      z2=CMPLX(x,0.0)
373      z1=m*z2
374
375      IF (0.0.LE.x.AND.x.LE.8.) THEN
376        Nmax=INT(x+4*x**(1./3.)+1.)+2
377      ELSEIF (8..LT.x.AND.x.LT.4200.) THEN
378        Nmax=INT(x+4.05*x**(1./3.)+2.)+1
379      ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN
380        Nmax=INT(x+4*x**(1./3.)+2.)+1
381      ELSE
382        WRITE(lunout,*) 'x out of bound, x=', x
383        STOP
384      ENDIF
385
386      Nstart=Nmax+100
387
388    !-----------loop for nu1z1, nu1z2
389
390      nu1z1(Nstart)=CMPLX(0.0,0.0)
391      nu1z2(Nstart)=CMPLX(0.0,0.0)
392      DO n=Nstart-1, 1 , -1
393        nn=CMPLX(FLOAT(n),0.0)
394        nu1z1(n)=(nn+1.)/z1 - 1./( (nn+1.)/z1 + nu1z1(n+1) )
395        nu1z2(n)=(nn+1.)/z2 - 1./( (nn+1.)/z2 + nu1z2(n+1) )
396      ENDDO
397
398    !------------loop for nu3z2
399
400      nu3z2(0)=-I
401      DO n=1, Nmax
402        nn=CMPLX(FLOAT(n),0.0)
403        nu3z2(n)=-nn/z2 + 1./ (nn/z2 - nu3z2(n-1) )
404      ENDDO
405
406    !-----------loop for psiz2 and ksiz2 (z2)
407      ksiz2(-1)=COS(REAL(z2))-I*SIN(REAL(z2))
408      ksiz2(0)=SIN(REAL(z2))+I*COS(REAL(z2))
409      DO n=1,Nmax
410       nn=CMPLX(FLOAT(n),0.0)
411       ksiz2(n)=(2.*nn-1.)/z2 * ksiz2(n-1) - ksiz2(n-2)
412       psiz2(n)=CMPLX(REAL(ksiz2(n)),0.0)
413      ENDDO
414
415    !-----------loop for a(n) and b(n)
416
417      DO n=1, Nmax
418        u1=k3*nu1z1(n) - k2*nu1z2(n)
419        u5=k3*nu1z1(n) - k2*nu3z2(n)
420        u6=k2*nu1z1(n) - k3*nu1z2(n)
421        u8=k2*nu1z1(n) - k3*nu3z2(n)
422        a(n)=psiz2(n)/ksiz2(n) * u1/u5
423        b(n)=psiz2(n)/ksiz2(n) * u6/u8
424      ENDDO
425
426    !-----------------final loop--------------
427      Q_ext=0.0
428      Q_sca=0.0
429      g=0.0
430
431      DO n=Nmax-1,1,-1
432        nnn=FLOAT(n)
433        Q_ext=Q_ext+ (2.*nnn+1.) * REAL( a(n)+b(n) )
434        Q_sca=Q_sca+ (2.*nnn+1.) *  &
435                   REAL( a(n)*CONJG(a(n)) + b(n)*CONJG(b(n)) )
436        g=g + nnn*(nnn+2.)/(nnn+1.) *   &
437           REAL( a(n)*CONJG(a(n+1))+b(n)*CONJG(b(n+1)) )  +   &
438              (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n)))
439      ENDDO
440
441      Q_ext=2./x**2 * Q_ext
442      Q_sca=2./x**2 * Q_sca
443    !--extrapolation in case of small x values
444      IF (smallx) THEN
445        Q_ext = x_old/x * Q_ext
446        Q_sca = x_old/x * Q_sca
447      ENDIF
448
449      Q_abs=Q_ext-Q_sca
450
451      IF (AIMAG(m).EQ.0.0) Q_abs=0.0
452      omega=Q_sca/Q_ext
453
454    ! g is wrong in the smallx case (but that does not matter as long as we ignore LW scattering)
455      g=g*4./x**2/Q_sca
456
457      sigma_sca=sigma_sca+r**2*Q_sca*number
458      sigma_abs=sigma_abs+r**2*Q_abs*number
459      sigma_ext=sigma_ext+r**2*Q_ext*number
460      omegatot=omegatot+r**2*Q_ext*omega*number
461      gtot    =gtot+r**2*Q_sca*g*number
462 
463      ENDDO   !---bin
464
465    !------------------------------------------------------------------
466
467      sigma_sca=RPI*sigma_sca
468      sigma_abs=RPI*sigma_abs
469      sigma_ext=RPI*sigma_ext
470      gtot=RPI*gtot/sigma_sca
471      omegatot=RPI*omegatot/sigma_ext
472
473      final_g(Nwv)=gtot
474      final_w(Nwv)=omegatot
475!      final_a(Nwv)=sigma_ext/masse
476      final_a(Nwv)=sigma_ext !extinction/absorption cross section per particle
477
478      ENDDO  !--loop on wavelength
479
480    !---averaging over LMDZ wavebands
481
482      DO band=1, nbands_sw_rrtm+nbands_lw_rrtm
483        gcm_a(band)=0.0
484        gcm_g(band)=0.0
485        gcm_w(band)=0.0
486        gcm_weight_a(band)=0.0
487        gcm_weight_g(band)=0.0
488        gcm_weight_w(band)=0.0
489      ENDDO
490
491    !---band 1 is now in the UV, so we take the first wavelength as being representative
492      DO Nwv=1,1
493        bandSW=1
494        gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv)
495        gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv)
496        gcm_w(bandSW)=gcm_w(bandSW)+final_w(Nwv)*final_a(Nwv)*weight(Nwv)
497        gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+final_a(Nwv)*weight(Nwv)
498        gcm_g(bandSW)=gcm_g(bandSW)+final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
499        gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+final_a(Nwv)*final_w(Nwv)*weight(Nwv)
500      ENDDO
501
502      DO Nwv=1,NwvmaxSW
503
504        IF (lambda_int(Nwv).LE.wv_rrtm_SW(3)) THEN      !--RRTM spectral interval 2
505          bandSW=2
506        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(4)) THEN  !--RRTM spectral interval 3
507          bandSW=3
508        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(5)) THEN  !--RRTM spectral interval 4
509          bandSW=4
510        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(6)) THEN  !--RRTM spectral interval 5
511          bandSW=5
512        ELSE                                            !--RRTM spectral interval 6
513          bandSW=6
514        ENDIF
515
516        gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv)
517        gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv)
518        gcm_w(bandSW)=gcm_w(bandSW)+final_w(Nwv)*final_a(Nwv)*weight(Nwv)
519        gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+final_a(Nwv)*weight(Nwv)
520        gcm_g(bandSW)=gcm_g(bandSW)+final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
521        gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+final_a(Nwv)*final_w(Nwv)*weight(Nwv)
522
523      ENDDO
524
525      DO Nwv=NwvmaxSW+nwave_sw+nwave_lw+1,NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
526        ll=lambda_int(Nwv)
527        dl=dlambda_int(Nwv)
528        weightLW=1./ll**5./(exp(hh*cc/kb/Tb/ll)-1.)*dl
529        bandLW=1  !--default value starting from the highest lambda
530        DO band=2, nbands_lw_rrtm
531          IF (ll.LT.wv_rrtm(band)) THEN   !--as long as
532            bandLW=band
533          ENDIF
534        ENDDO
535        gcm_a(nbands_sw_rrtm+bandLW)=gcm_a(nbands_sw_rrtm+bandLW)+final_a(Nwv)*   &
536             (1.-final_w(Nwv))*weightLW
537        gcm_weight_a(nbands_sw_rrtm+bandLW)=gcm_weight_a(nbands_sw_rrtm+bandLW)+weightLW
538
539        gcm_w(nbands_sw_rrtm+bandLW)=gcm_w(nbands_sw_rrtm+bandLW)+final_w(Nwv)*   &
540             final_a(Nwv)*weightLW
541        gcm_weight_w(nbands_sw_rrtm+bandLW)=gcm_weight_w(nbands_sw_rrtm+bandLW)+  &
542             final_a(Nwv)*weightLW
543
544        gcm_g(nbands_sw_rrtm+bandLW)=gcm_g(nbands_sw_rrtm+bandLW)+final_g(Nwv)*   &
545             final_a(Nwv)*final_w(Nwv)*weightLW
546        gcm_weight_g(nbands_sw_rrtm+bandLW)=gcm_weight_g(nbands_sw_rrtm+bandLW)+  &
547             final_a(Nwv)*final_w(Nwv)*weightLW
548      ENDDO
549
550      DO band=1, nbands_sw_rrtm+nbands_lw_rrtm
551        gcm_a(band)=gcm_a(band)/gcm_weight_a(band)
552        gcm_w(band)=gcm_w(band)/gcm_weight_w(band)
553        gcm_g(band)=gcm_g(band)/gcm_weight_g(band)
554        ss_a(band)=gcm_a(band)
555        ss_w(band)=gcm_w(band)
556        ss_g(band)=gcm_g(band)
557      ENDDO
558
559      DO nb=1, nwave_sw+nwave_lw
560        ss_a(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_a(NwvmaxSW+nb)
561        ss_w(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_w(NwvmaxSW+nb)
562        ss_g(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_g(NwvmaxSW+nb)
563      ENDDO
564
565      DO nb=1,nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw
566        alpha_bin(nb,bin_number)=ss_a(nb) !extinction/absorption cross section per particle
567        piz_bin(nb,bin_number)=ss_w(nb)
568        cg_bin(nb,bin_number)=ss_g(nb)
569      ENDDO
570
571    ENDDO !loop over tracer bins
572
573!!$OMP END MASTER
574!  CALL bcast(alpha_bin)
575!  CALL bcast(piz_bin)
576!  CALL bcast(cg_bin)
577!!$OMP BARRIER
578
579    !set to default values at first time step (tr_seri still zero)
580    tau_strat(:,:,:)=1.e-15
581    piz_strat(:,:,:)=1.0
582    cg_strat(:,:,:)=0.0
583    tau_lw_abs_rrtm(:,:,:)=1.e-15
584    tau_strat_wave(:,:,:)=1.e-15
585
586  ELSE  !-- not debut
587
588  !--compute optical properties of actual size distribution (from tr_seri)
589    DO ilon=1,klon
590      DO ilev=1, klev
591        DO nb=1,nbands_sw_rrtm
592          tau_strat(ilon,ilev,nb)=0.0
593          DO bin_number=1, nbtr_bin
594            tau_strat(ilon,ilev,nb)=tau_strat(ilon,ilev,nb)+alpha_bin(nb,bin_number) &
595                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
596          ENDDO
597
598          piz_strat(ilon,ilev,nb)=0.0
599          DO bin_number=1, nbtr_bin
600            piz_strat(ilon,ilev,nb)=piz_strat(ilon,ilev,nb)+piz_bin(nb,bin_number)*alpha_bin(nb,bin_number) &
601                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
602          ENDDO
603          piz_strat(ilon,ilev,nb)=piz_strat(ilon,ilev,nb)/MAX(tau_strat(ilon,ilev,nb),1.e-15)
604
605          cg_strat(ilon,ilev,nb)=0.0
606          DO bin_number=1, nbtr_bin
607            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) &
608                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
609          ENDDO
610          cg_strat(ilon,ilev,nb)=cg_strat(ilon,ilev,nb)/MAX(tau_strat(ilon,ilev,nb)*piz_strat(ilon,ilev,nb),1.e-15)
611        ENDDO
612        DO nb=1,nbands_lw_rrtm
613          tau_lw_abs_rrtm(ilon,ilev,nb)=0.0
614          DO bin_number=1, nbtr_bin
615            tau_lw_abs_rrtm(ilon,ilev,nb)=tau_lw_abs_rrtm(ilon,ilev,nb)+alpha_bin(nbands_sw_rrtm+nb,bin_number) &
616                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
617          ENDDO
618        ENDDO
619        DO nb=1,nwave_sw+nwave_lw
620          tau_strat_wave(ilon,ilev,nb)=0.0
621          DO bin_number=1, nbtr_bin
622            tau_strat_wave(ilon,ilev,nb)=tau_strat_wave(ilon,ilev,nb)+alpha_bin(nbands_sw_rrtm+nbands_lw_rrtm+nb,bin_number) &
623                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
624          ENDDO
625        ENDDO
626      ENDDO
627    ENDDO
628
629  ENDIF !debut
630
631END SUBROUTINE MIECALC_AER
Note: See TracBrowser for help on using the repository browser.