Context Navigation

miecalc_aer.F90 @ 3634

Last change on this file since 3634 was 3526, checked in by Laurent Fairhead, 5 years ago
Olivier's modifications + addition of Id keyword
Property svn:keywords set to `Id`
File size: 22.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format