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

Last change on this file since 5559 was 5559, checked in by lebasn, 14 hours ago

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