source: LMDZ6/trunk/libf/phylmd/StratAer/miecalc_aer.f90 @ 5354

Last change on this file since 5354 was 5268, checked in by abarral, 7 weeks ago

.f90 <-> .F90 depending on cpp key use
  • Property svn:keywords set to Id
File size: 22.7 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, ONLY: dens_aer_dry, dens_aer_ref, V_rat
19  USE aero_mod
20  USE infotrac_phy, ONLY : nbtr, nbtr_bin, nbtr_sulgas, id_SO2_strat
21  USE dimphy
22  USE yomcst_mod_h  , 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    !--compute particle radius for a composition of 75% H2SO4 / 25% H2O at T=293K
230    DO bin_number=1, nbtr_bin
231      r_0=(dens_aer_dry/dens_aer_ref/0.75)**(1./3.)*mdw(bin_number)/2.
232    !--integral boundaries set to bin boundaries
233      rmin=r_0/sqrt(V_rat**(1./3.))
234      rmax=r_0*sqrt(V_rat**(1./3.))
235
236    !--set up SW
237      DO Nwv=1, NwvmaxSW
238        lambda_int(Nwv)=( lambda(Nwv)+lambda(Nwv+1) ) /2.
239      ENDDO
240
241      DO nb=1, nwave_sw+nwave_lw
242        lambda_int(NwvmaxSW+nb)=lambda_ref(nb)
243      ENDDO
244
245    !--set up LW
246    !--conversion wavenumber in cm-1 to wavelength in m
247      DO Nwv=1, nbands_lw_rrtm+1
248        wv_rrtm(Nwv)=10000./wn_rrtm(Nwv)*1.e-6
249      ENDDO
250
251      DO Nwv=1, NwvmaxLW
252        lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
253          exp( log(wv_rrtm(1))+float(Nwv-1)/float(NwvmaxLW-1)* &
254          (log(wv_rrtm(nbands_lw_rrtm+1))-log(wv_rrtm(1))) )
255      ENDDO
256
257!--computing the dlambdas
258      Nwv=1
259      dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
260      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)- &
261      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv+1)
262      DO Nwv=2, NwvmaxLW-1
263      dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
264      &  (lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv-1)- &
265      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv+1))/2.
266      ENDDO
267      Nwv=NwvmaxLW
268      dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= &
269      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv-1)- &
270      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)
271
272      IF (refr_ind_interpol) THEN
273 
274        ilambda_max=ref_ind(1,2)/1.e6 !--in m
275        ilambda_min=ref_ind(nb_lambda_h2so4,2)/1.e6 !--in m
276        DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
277          IF (lambda_int(Nwv).GT.ilambda_max) THEN
278            !for lambda out of data range, take boundary values
279            n_r(Nwv)=ref_ind(1,3)
280            n_i(Nwv)=ref_ind(1,4)
281          ELSEIF (lambda_int(Nwv).LE.ilambda_min) THEN
282            n_r(Nwv)=ref_ind(nb_lambda_h2so4,3)
283            n_i(Nwv)=ref_ind(nb_lambda_h2so4,4)
284          ELSE
285            DO nb=2,nb_lambda_h2so4
286              ilambda=ref_ind(nb,2)/1.e6
287              ilambda_prev=ref_ind(nb-1,2)/1.e6
288              n_r_h2so4=ref_ind(nb,3)
289              n_r_h2so4_prev=ref_ind(nb-1,3)
290              n_i_h2so4=ref_ind(nb,4)
291              n_i_h2so4_prev=ref_ind(nb-1,4)
292              IF (lambda_int(Nwv).GT.ilambda.AND. &
293                lambda_int(Nwv).LE.ilambda_prev) THEN !--- linear interpolation
294                n_r(Nwv)=n_r_h2so4+(lambda_int(Nwv)-ilambda)/ &
295                     (ilambda_prev-ilambda)*(n_r_h2so4_prev-n_r_h2so4)
296                n_i(Nwv)=n_i_h2so4+(lambda_int(Nwv)-ilambda)/ &
297                     (ilambda_prev-ilambda)*(n_i_h2so4_prev-n_i_h2so4)
298              ENDIF
299            ENDDO
300          ENDIF
301        ENDDO
302
303      ELSE  !-- no refr_ind_interpol, closest neighbour from upper wavelength
304
305        DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
306          n_r(Nwv)=ref_ind(1,3)
307          n_i(Nwv)=ref_ind(1,4)
308          DO nb=2,nb_lambda_h2so4
309            IF (ref_ind(nb,2)/1.e6.GT.lambda_int(Nwv)) THEN !--- step function
310              n_r(Nwv)=ref_ind(nb,3)
311              n_i(Nwv)=ref_ind(nb,4)
312            ENDIF 
313          ENDDO
314        ENDDO
315      ENDIF
316
317    !---Loop on wavelengths
318      DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
319
320      m=CMPLX(n_r(Nwv),-n_i(Nwv))
321
322      I=CMPLX(0.,1.)
323
324      sigma_sca=0.0
325      sigma_ext=0.0
326      sigma_abs=0.0
327      gtot=0.0
328      omegatot=0.0
329      masse = 0.0
330      volume=0.0
331      surface=0.0
332
333      DO bin=1, Nbin !---loop on size bins
334
335      r_lower=exp(log(rmin)+FLOAT(bin-1)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
336      r_upper=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
337      deltar=r_upper-r_lower
338
339      r=sqrt(r_lower*r_upper)
340      x=2.*RPI*r/lambda_int(Nwv)
341
342!we impose a minimum value for x and extrapolate quantities for small x values
343      smallx = .FALSE.
344      IF (x.LT.0.001) THEN
345        smallx = .TRUE.
346        x_old = x
347        x = 0.001
348      ENDIF
349
350      number=Ntot*deltar/(rmax-rmin) !dN/dr constant over tracer bin
351!      masse=masse  +4./3.*RPI*(r**3)*number*deltar*ropx*1.E3  !--g/m3
352      volume=volume+4./3.*RPI*(r**3)*number*deltar
353      surface=surface+4.*RPI*r**2*number*deltar
354
355      k2=m
356      k3=CMPLX(1.0,0.0)
357
358      z2=CMPLX(x,0.0)
359      z1=m*z2
360
361      IF (0.0.LE.x.AND.x.LE.8.) THEN
362        Nmax=INT(x+4*x**(1./3.)+1.)+2
363      ELSEIF (8..LT.x.AND.x.LT.4200.) THEN
364        Nmax=INT(x+4.05*x**(1./3.)+2.)+1
365      ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN
366        Nmax=INT(x+4*x**(1./3.)+2.)+1
367      ELSE
368        WRITE(lunout,*) 'x out of bound, x=', x
369        STOP
370      ENDIF
371
372      Nstart=Nmax+100
373
374    !-----------loop for nu1z1, nu1z2
375
376      nu1z1(Nstart)=CMPLX(0.0,0.0)
377      nu1z2(Nstart)=CMPLX(0.0,0.0)
378      DO n=Nstart-1, 1 , -1
379        nn=CMPLX(FLOAT(n),0.0)
380        nu1z1(n)=(nn+1.)/z1 - 1./( (nn+1.)/z1 + nu1z1(n+1) )
381        nu1z2(n)=(nn+1.)/z2 - 1./( (nn+1.)/z2 + nu1z2(n+1) )
382      ENDDO
383
384    !------------loop for nu3z2
385
386      nu3z2(0)=-I
387      DO n=1, Nmax
388        nn=CMPLX(FLOAT(n),0.0)
389        nu3z2(n)=-nn/z2 + 1./ (nn/z2 - nu3z2(n-1) )
390      ENDDO
391
392    !-----------loop for psiz2 and ksiz2 (z2)
393      ksiz2(-1)=COS(REAL(z2))-I*SIN(REAL(z2))
394      ksiz2(0)=SIN(REAL(z2))+I*COS(REAL(z2))
395      DO n=1,Nmax
396       nn=CMPLX(FLOAT(n),0.0)
397       ksiz2(n)=(2.*nn-1.)/z2 * ksiz2(n-1) - ksiz2(n-2)
398       psiz2(n)=CMPLX(REAL(ksiz2(n)),0.0)
399      ENDDO
400
401    !-----------loop for a(n) and b(n)
402
403      DO n=1, Nmax
404        u1=k3*nu1z1(n) - k2*nu1z2(n)
405        u5=k3*nu1z1(n) - k2*nu3z2(n)
406        u6=k2*nu1z1(n) - k3*nu1z2(n)
407        u8=k2*nu1z1(n) - k3*nu3z2(n)
408        a(n)=psiz2(n)/ksiz2(n) * u1/u5
409        b(n)=psiz2(n)/ksiz2(n) * u6/u8
410      ENDDO
411
412    !-----------------final loop--------------
413      Q_ext=0.0
414      Q_sca=0.0
415      g=0.0
416
417      DO n=Nmax-1,1,-1
418        nnn=FLOAT(n)
419        Q_ext=Q_ext+ (2.*nnn+1.) * REAL( a(n)+b(n) )
420        Q_sca=Q_sca+ (2.*nnn+1.) *  &
421                   REAL( a(n)*CONJG(a(n)) + b(n)*CONJG(b(n)) )
422        g=g + nnn*(nnn+2.)/(nnn+1.) *   &
423           REAL( a(n)*CONJG(a(n+1))+b(n)*CONJG(b(n+1)) )  +   &
424              (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n)))
425      ENDDO
426
427      Q_ext=2./x**2 * Q_ext
428      Q_sca=2./x**2 * Q_sca
429    !--extrapolation in case of small x values
430      IF (smallx) THEN
431        Q_ext = x_old/x * Q_ext
432        Q_sca = x_old/x * Q_sca
433      ENDIF
434
435      Q_abs=Q_ext-Q_sca
436
437      IF (AIMAG(m).EQ.0.0) Q_abs=0.0
438      omega=Q_sca/Q_ext
439
440    ! g is wrong in the smallx case (but that does not matter as long as we ignore LW scattering)
441      g=g*4./x**2/Q_sca
442
443      sigma_sca=sigma_sca+r**2*Q_sca*number
444      sigma_abs=sigma_abs+r**2*Q_abs*number
445      sigma_ext=sigma_ext+r**2*Q_ext*number
446      omegatot=omegatot+r**2*Q_ext*omega*number
447      gtot    =gtot+r**2*Q_sca*g*number
448 
449      ENDDO   !---bin
450
451    !------------------------------------------------------------------
452
453      sigma_sca=RPI*sigma_sca
454      sigma_abs=RPI*sigma_abs
455      sigma_ext=RPI*sigma_ext
456      gtot=RPI*gtot/sigma_sca
457      omegatot=RPI*omegatot/sigma_ext
458
459      final_g(Nwv)=gtot
460      final_w(Nwv)=omegatot
461!      final_a(Nwv)=sigma_ext/masse
462      final_a(Nwv)=sigma_ext !extinction/absorption cross section per particle
463
464      ENDDO  !--loop on wavelength
465
466    !---averaging over LMDZ wavebands
467
468      DO band=1, nbands_sw_rrtm+nbands_lw_rrtm
469        gcm_a(band)=0.0
470        gcm_g(band)=0.0
471        gcm_w(band)=0.0
472        gcm_weight_a(band)=0.0
473        gcm_weight_g(band)=0.0
474        gcm_weight_w(band)=0.0
475      ENDDO
476
477    !---band 1 is now in the UV, so we take the first wavelength as being representative
478      DO Nwv=1,1
479        bandSW=1
480        gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv)
481        gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv)
482        gcm_w(bandSW)=gcm_w(bandSW)+final_w(Nwv)*final_a(Nwv)*weight(Nwv)
483        gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+final_a(Nwv)*weight(Nwv)
484        gcm_g(bandSW)=gcm_g(bandSW)+final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
485        gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+final_a(Nwv)*final_w(Nwv)*weight(Nwv)
486      ENDDO
487
488      DO Nwv=1,NwvmaxSW
489
490        IF (lambda_int(Nwv).LE.wv_rrtm_SW(3)) THEN      !--RRTM spectral interval 2
491          bandSW=2
492        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(4)) THEN  !--RRTM spectral interval 3
493          bandSW=3
494        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(5)) THEN  !--RRTM spectral interval 4
495          bandSW=4
496        ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(6)) THEN  !--RRTM spectral interval 5
497          bandSW=5
498        ELSE                                            !--RRTM spectral interval 6
499          bandSW=6
500        ENDIF
501
502        gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv)
503        gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv)
504        gcm_w(bandSW)=gcm_w(bandSW)+final_w(Nwv)*final_a(Nwv)*weight(Nwv)
505        gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+final_a(Nwv)*weight(Nwv)
506        gcm_g(bandSW)=gcm_g(bandSW)+final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
507        gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+final_a(Nwv)*final_w(Nwv)*weight(Nwv)
508
509      ENDDO
510
511      DO Nwv=NwvmaxSW+nwave_sw+nwave_lw+1,NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
512        ll=lambda_int(Nwv)
513        dl=dlambda_int(Nwv)
514        weightLW=1./ll**5./(exp(hh*cc/kb/Tb/ll)-1.)*dl
515        bandLW=1  !--default value starting from the highest lambda
516        DO band=2, nbands_lw_rrtm
517          IF (ll.LT.wv_rrtm(band)) THEN   !--as long as
518            bandLW=band
519          ENDIF
520        ENDDO
521        gcm_a(nbands_sw_rrtm+bandLW)=gcm_a(nbands_sw_rrtm+bandLW)+final_a(Nwv)*   &
522             (1.-final_w(Nwv))*weightLW
523        gcm_weight_a(nbands_sw_rrtm+bandLW)=gcm_weight_a(nbands_sw_rrtm+bandLW)+weightLW
524
525        gcm_w(nbands_sw_rrtm+bandLW)=gcm_w(nbands_sw_rrtm+bandLW)+final_w(Nwv)*   &
526             final_a(Nwv)*weightLW
527        gcm_weight_w(nbands_sw_rrtm+bandLW)=gcm_weight_w(nbands_sw_rrtm+bandLW)+  &
528             final_a(Nwv)*weightLW
529
530        gcm_g(nbands_sw_rrtm+bandLW)=gcm_g(nbands_sw_rrtm+bandLW)+final_g(Nwv)*   &
531             final_a(Nwv)*final_w(Nwv)*weightLW
532        gcm_weight_g(nbands_sw_rrtm+bandLW)=gcm_weight_g(nbands_sw_rrtm+bandLW)+  &
533             final_a(Nwv)*final_w(Nwv)*weightLW
534      ENDDO
535
536      DO band=1, nbands_sw_rrtm+nbands_lw_rrtm
537        gcm_a(band)=gcm_a(band)/gcm_weight_a(band)
538        gcm_w(band)=gcm_w(band)/gcm_weight_w(band)
539        gcm_g(band)=gcm_g(band)/gcm_weight_g(band)
540        ss_a(band)=gcm_a(band)
541        ss_w(band)=gcm_w(band)
542        ss_g(band)=gcm_g(band)
543      ENDDO
544
545      DO nb=1, nwave_sw+nwave_lw
546        ss_a(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_a(NwvmaxSW+nb)
547        ss_w(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_w(NwvmaxSW+nb)
548        ss_g(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_g(NwvmaxSW+nb)
549      ENDDO
550
551      DO nb=1,nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw
552        alpha_bin(nb,bin_number)=ss_a(nb) !extinction/absorption cross section per particle
553        piz_bin(nb,bin_number)=ss_w(nb)
554        cg_bin(nb,bin_number)=ss_g(nb)
555      ENDDO
556
557    ENDDO !loop over tracer bins
558
559!!$OMP END MASTER
560!  CALL bcast(alpha_bin)
561!  CALL bcast(piz_bin)
562!  CALL bcast(cg_bin)
563!!$OMP BARRIER
564
565    !set to default values at first time step (tr_seri still zero)
566    tau_strat(:,:,:)=1.e-15
567    piz_strat(:,:,:)=1.0
568    cg_strat(:,:,:)=0.0
569    tau_lw_abs_rrtm(:,:,:)=1.e-15
570    tau_strat_wave(:,:,:)=1.e-15
571
572  ELSE  !-- not debut
573
574  !--compute optical properties of actual size distribution (from tr_seri)
575    DO ilon=1,klon
576      DO ilev=1, klev
577        DO nb=1,nbands_sw_rrtm
578          tau_strat(ilon,ilev,nb)=0.0
579          DO bin_number=1, nbtr_bin
580            tau_strat(ilon,ilev,nb)=tau_strat(ilon,ilev,nb)+alpha_bin(nb,bin_number) &
581                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
582          ENDDO
583
584          piz_strat(ilon,ilev,nb)=0.0
585          DO bin_number=1, nbtr_bin
586            piz_strat(ilon,ilev,nb)=piz_strat(ilon,ilev,nb)+piz_bin(nb,bin_number)*alpha_bin(nb,bin_number) &
587                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
588          ENDDO
589          piz_strat(ilon,ilev,nb)=piz_strat(ilon,ilev,nb)/MAX(tau_strat(ilon,ilev,nb),1.e-15)
590
591          cg_strat(ilon,ilev,nb)=0.0
592          DO bin_number=1, nbtr_bin
593            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) &
594                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
595          ENDDO
596          cg_strat(ilon,ilev,nb)=cg_strat(ilon,ilev,nb)/MAX(tau_strat(ilon,ilev,nb)*piz_strat(ilon,ilev,nb),1.e-15)
597        ENDDO
598        DO nb=1,nbands_lw_rrtm
599          tau_lw_abs_rrtm(ilon,ilev,nb)=0.0
600          DO bin_number=1, nbtr_bin
601            tau_lw_abs_rrtm(ilon,ilev,nb)=tau_lw_abs_rrtm(ilon,ilev,nb)+alpha_bin(nbands_sw_rrtm+nb,bin_number) &
602                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
603          ENDDO
604        ENDDO
605        DO nb=1,nwave_sw+nwave_lw
606          tau_strat_wave(ilon,ilev,nb)=0.0
607          DO bin_number=1, nbtr_bin
608            tau_strat_wave(ilon,ilev,nb)=tau_strat_wave(ilon,ilev,nb)+alpha_bin(nbands_sw_rrtm+nbands_lw_rrtm+nb,bin_number) &
609                                *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG
610          ENDDO
611        ENDDO
612      ENDDO
613    ENDDO
614
615  ENDIF !debut
616
617END SUBROUTINE MIECALC_AER
Note: See TracBrowser for help on using the repository browser.