Context Navigation

miecalc_aer.F90 @ 4996

Last change on this file since 4996 was 3668, checked in by oboucher, 5 years ago
Correction of a big big bug on ref_ind array initialisation and Mie calc in StratAer?
File size: 23.1 KB

Line
1	SUBROUTINE 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_0sqrt(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.RPIr/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)numberdeltarropx*1.E3 !--g/m3
366	volume=volume+4./3.RPI(r*3)number*deltar
367	surface=surface+4.RPIr*2number*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+4x*(1./3.)+1.)+2
377	ELSEIF (8..LT.x.AND.x.LT.4200.) THEN
378	Nmax=INT(x+4.05x*(1./3.)+2.)+1
379	ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN
380	Nmax=INT(x+4x*(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=k3nu1z1(n) - k2nu1z2(n)
419	u5=k3nu1z1(n) - k2nu3z2(n)
420	u6=k2nu1z1(n) - k3nu1z2(n)
421	u8=k2nu1z1(n) - k3nu3z2(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=g4./x*2/Q_sca
456
457	sigma_sca=sigma_sca+r*2Q_sca*number
458	sigma_abs=sigma_abs+r*2Q_abs*number
459	sigma_ext=sigma_ext+r*2Q_ext*number
460	omegatot=omegatot+r*2Q_extomeganumber
461	gtot =gtot+r*2Q_scagnumber
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(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