Context Navigation

aeroptproperties.F90 @ 2937

Last change on this file since 2937 was 2899, checked in by emillour, 3 years ago
Generic PCM: More code tidying: turn aeropacity, aeroptproperties, gfluxi, gfluxv, sfluxi and sfluxv into modules. EM
File size: 36.6 KB

Line
1	module aeroptproperties_mod
2
3	implicit none
4
5	contains
6
7	SUBROUTINE aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
8	QVISsQREF3d,omegaVIS3d,gVIS3d, &
9	QIRsQREF3d,omegaIR3d,gIR3d, &
10	QREFvis3d,QREFir3d)!, &
11	! omegaREFvis3d,omegaREFir3d)
12
13	use radinc_h, only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
14	use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir
15	use radcommon_h, only: qrefvis,qrefir,omegarefir !,omegarefvis
16	use radcommon_h, only: radiustab,nsize
17
18	implicit none
19
20	! =============================================================
21	! Aerosol Optical Properties
22	!
23	! Description:
24	! Compute the scattering parameters in each grid
25	! box, depending on aerosol grain sizes. Log-normal size
26	! distribution and Gauss-Legendre integration are used.
27
28	! Parameters:
29	! Don't forget to set the value of varyingnueff below; If
30	! the effective variance of the distribution for the given
31	! aerosol is considered homogeneous in the atmosphere, please
32	! set varyingnueff(iaer) to .false. Resulting computational
33	! time will be much better.
34
35	! Authors: J.-B. Madeleine, F. Forget, F. Montmessin
36	! Slightly modified and converted to F90 by R. Wordsworth (2009)
37	! Varying nueff section removed by R. Wordsworth for simplicity
38	! ==============================================================
39
40	! Local variables
41	! ---------------
42
43
44
45	! =============================================================
46	LOGICAL, PARAMETER :: varyingnueff(naerkind) = .false.
47	! =============================================================
48
49	! Min. and max radius of the interpolation grid (in METERS)
50	REAL, PARAMETER :: refftabmin = 2e-8 !2e-8
51	! REAL, PARAMETER :: refftabmax = 35e-6
52	REAL, PARAMETER :: refftabmax = 1e-3
53	! Log of the min and max variance of the interpolation grid
54	REAL, PARAMETER :: nuefftabmin = -4.6
55	REAL, PARAMETER :: nuefftabmax = 0.
56	! Number of effective radius of the interpolation grid
57	INTEGER, PARAMETER :: refftabsize = 200
58	! Number of effective variances of the interpolation grid
59	! INTEGER, PARAMETER :: nuefftabsize = 100
60	INTEGER, PARAMETER :: nuefftabsize = 1
61	! Interpolation grid indices (reff,nueff)
62	INTEGER :: grid_i,grid_j
63	! Intermediate variable
64	REAL :: var_tmp,var3d_tmp(ngrid,nlayer)
65	! Bilinear interpolation factors
66	REAL :: kx,ky,k1,k2,k3,k4
67	! Size distribution parameters
68	REAL :: sizedistk1,sizedistk2
69	! Pi!
70	REAL,SAVE :: pi
71	!$OMP THREADPRIVATE(pi)
72	! Variables used by the Gauss-Legendre integration:
73	INTEGER radius_id,gausind
74	REAL kint
75	REAL drad
76	INTEGER, PARAMETER :: ngau = 10
77	REAL weightgaus(ngau),radgaus(ngau)
78	SAVE weightgaus,radgaus
79	! DATA weightgaus/.2955242247,.2692667193,.2190863625,.1494513491,.0666713443/
80	! DATA radgaus/.1488743389,.4333953941,.6794095682,.8650633666,.9739065285/
81	DATA radgaus/0.07652652113350,0.22778585114165, &
82	0.37370608871528,0.51086700195146, &
83	0.63605368072468,0.74633190646476, &
84	0.83911697181213,0.91223442826796, &
85	0.96397192726078,0.99312859919241/
86
87	DATA weightgaus/0.15275338723120,0.14917298659407, &
88	0.14209610937519,0.13168863843930, &
89	0.11819453196154,0.10193011980823, &
90	0.08327674160932,0.06267204829828, &
91	0.04060142982019,0.01761400714091/
92	!$OMP THREADPRIVATE(radgaus,weightgaus)
93	! Indices
94	INTEGER :: i,j,k,l,m,iaer,idomain
95	INTEGER :: ig,lg,chg
96
97	! Local saved variables
98	! ---------------------
99
100	! Radius axis of the interpolation grid
101	REAL,SAVE :: refftab(refftabsize)
102	! Variance axis of the interpolation grid
103	REAL,SAVE :: nuefftab(nuefftabsize)
104	! Volume ratio of the grid
105	REAL,SAVE :: logvratgrid,vratgrid
106	! Grid used to remember which calculation is done
107	LOGICAL,SAVE :: checkgrid(refftabsize,nuefftabsize,naerkind,2) = .false.
108	!$OMP THREADPRIVATE(refftab,nuefftab,logvratgrid,vratgrid,checkgrid)
109	! Optical properties of the grid (VISIBLE)
110	REAL,SAVE :: qsqrefVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
111	REAL,SAVE :: qextVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
112	REAL,SAVE :: qscatVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
113	REAL,SAVE :: omegVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
114	REAL,SAVE :: gVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
115	!$OMP THREADPRIVATE(qsqrefVISgrid,qextVISgrid,qscatVISgrid,omegVISgrid,gVISgrid)
116	! Optical properties of the grid (INFRARED)
117	REAL,SAVE :: qsqrefIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
118	REAL,SAVE :: qextIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
119	REAL,SAVE :: qscatIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
120	REAL,SAVE :: omegIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
121	REAL,SAVE :: gIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
122	!$OMP THREADPRIVATE(qsqrefIRgrid,qextIRgrid,qscatIRgrid,omegIRgrid,gIRgrid)
123	! Optical properties of the grid (REFERENCE WAVELENGTHS)
124	REAL,SAVE :: qrefVISgrid(refftabsize,nuefftabsize,naerkind)
125	REAL,SAVE :: qscatrefVISgrid(refftabsize,nuefftabsize,naerkind)
126	REAL,SAVE :: qrefIRgrid(refftabsize,nuefftabsize,naerkind)
127	REAL,SAVE :: qscatrefIRgrid(refftabsize,nuefftabsize,naerkind)
128	REAL,SAVE :: omegrefVISgrid(refftabsize,nuefftabsize,naerkind)
129	REAL,SAVE :: omegrefIRgrid(refftabsize,nuefftabsize,naerkind)
130	!$OMP THREADPRIVATE(qrefVISgrid,qscatrefVISgrid,qrefIRgrid,qscatrefIRgrid,omegrefVISgrid,&
131	!$OMP omegrefIRgrid)
132	! Firstcall
133	LOGICAL,SAVE :: firstcall = .true.
134	!$OMP THREADPRIVATE(firstcall)
135	! Variables used by the Gauss-Legendre integration:
136	REAL,SAVE :: normd(refftabsize,nuefftabsize,naerkind,2)
137	REAL,SAVE :: dista(refftabsize,nuefftabsize,naerkind,2,ngau)
138	REAL,SAVE :: distb(refftabsize,nuefftabsize,naerkind,2,ngau)
139	!$OMP THREADPRIVATE(normd,dista,distb)
140
141	REAL,SAVE :: radGAUSa(ngau,naerkind,2)
142	REAL,SAVE :: radGAUSb(ngau,naerkind,2)
143	!$OMP THREADPRIVATE(radGAUSa,radGAUSb)
144
145	REAL,SAVE :: qsqrefVISa(L_NSPECTV,ngau,naerkind)
146	REAL,SAVE :: qrefVISa(ngau,naerkind)
147	REAL,SAVE :: qsqrefVISb(L_NSPECTV,ngau,naerkind)
148	REAL,SAVE :: qrefVISb(ngau,naerkind)
149	REAL,SAVE :: omegVISa(L_NSPECTV,ngau,naerkind)
150	REAL,SAVE :: omegrefVISa(ngau,naerkind)
151	REAL,SAVE :: omegVISb(L_NSPECTV,ngau,naerkind)
152	REAL,SAVE :: omegrefVISb(ngau,naerkind)
153	REAL,SAVE :: gVISa(L_NSPECTV,ngau,naerkind)
154	REAL,SAVE :: gVISb(L_NSPECTV,ngau,naerkind)
155	!$OMP THREADPRIVATE(qsqrefVISa,qrefVISa,qsqrefVISb,qrefVISb,omegVISa, &
156	!$OMP omegrefVISa,omegVISb,omegrefVISb,gVISa,gVISb)
157
158	REAL,SAVE :: qsqrefIRa(L_NSPECTI,ngau,naerkind)
159	REAL,SAVE :: qrefIRa(ngau,naerkind)
160	REAL,SAVE :: qsqrefIRb(L_NSPECTI,ngau,naerkind)
161	REAL,SAVE :: qrefIRb(ngau,naerkind)
162	REAL,SAVE :: omegIRa(L_NSPECTI,ngau,naerkind)
163	REAL,SAVE :: omegrefIRa(ngau,naerkind)
164	REAL,SAVE :: omegIRb(L_NSPECTI,ngau,naerkind)
165	REAL,SAVE :: omegrefIRb(ngau,naerkind)
166	REAL,SAVE :: gIRa(L_NSPECTI,ngau,naerkind)
167	REAL,SAVE :: gIRb(L_NSPECTI,ngau,naerkind)
168	!$OMP THREADPRIVATE(qsqrefIRa,qrefIRa,qsqrefIRb,qrefIRb,omegIRa,omegrefIRa,&
169	!$OMP omegIRb,omegrefIRb,gIRa,gIRb)
170
171	REAL :: radiusm
172	REAL :: radiusr
173
174	! Inputs
175	! ------
176
177	INTEGER :: ngrid,nlayer
178	! Aerosol effective radius used for radiative transfer (meter)
179	REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind)
180	! Aerosol effective variance used for radiative transfer (n.u.)
181	REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind)
182
183	! Outputs
184	! -------
185
186	REAL,INTENT(OUT) :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
187	REAL,INTENT(OUT) :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
188	REAL,INTENT(OUT) :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
189
190	REAL,INTENT(OUT) :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
191	REAL,INTENT(OUT) :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
192	REAL,INTENT(OUT) :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
193
194	REAL,INTENT(OUT) :: QREFvis3d(ngrid,nlayer,naerkind)
195	REAL,INTENT(OUT) :: QREFir3d(ngrid,nlayer,naerkind)
196
197	! REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
198	! REAL :: omegaREFir3d(ngrid,nlayer,naerkind)
199
200	DO iaer = 1, naerkind ! Loop on aerosol kind
201	IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN
202	!==================================================================
203	! If there is one single particle size, optical
204	! properties of the considered aerosol are homogeneous
205	DO lg = 1, nlayer
206	DO ig = 1, ngrid
207	DO chg = 1, L_NSPECTV
208	QVISsQREF3d(ig,lg,chg,iaer)=QVISsQREF(chg,iaer,1)
209	omegaVIS3d(ig,lg,chg,iaer)=omegaVIS(chg,iaer,1)
210	gVIS3d(ig,lg,chg,iaer)=gVIS(chg,iaer,1)
211	ENDDO
212	DO chg = 1, L_NSPECTI
213	QIRsQREF3d(ig,lg,chg,iaer)=QIRsQREF(chg,iaer,1)
214	omegaIR3d(ig,lg,chg,iaer)=omegaIR(chg,iaer,1)
215	gIR3d(ig,lg,chg,iaer)=gIR(chg,iaer,1)
216	ENDDO
217	QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1)
218	QREFir3d(ig,lg,iaer)=QREFir(iaer,1)
219	! omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1)
220	! omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1)
221	ENDDO
222	ENDDO
223
224
225	if (firstcall) then
226	print*,'Optical prop. of the aerosol are homogenous for:'
227	print*,'iaer = ',iaer
228	endif
229
230	ELSE ! Varying effective radius and variance
231	DO idomain = 1, 2 ! Loop on visible or infrared channel
232	!==================================================================
233	! 1. Creating the effective radius and variance grid
234	! --------------------------------------------------
235	IF (firstcall) THEN
236
237	! 1.1 Pi!
238	pi = 2. * asin(1.e0)
239
240	! 1.2 Effective radius
241	refftab(1) = refftabmin
242	refftab(refftabsize) = refftabmax
243
244	logvratgrid = log(refftabmax/refftabmin) / float(refftabsize-1)*3.
245	vratgrid = exp(logvratgrid)
246
247	do i = 2, refftabsize-1
248	refftab(i) = refftab(i-1)vratgrid*(1./3.)
249	enddo
250
251	! 1.3 Effective variance
252	if(nuefftabsize.eq.1)then ! addded by RDW
253	print*,'Warning: no variance range in aeroptproperties'
254	nuefftab(1)=0.2
255	else
256	do i = 0, nuefftabsize-1
257	nuefftab(i+1) = exp( nuefftabmin + i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
258	enddo
259	endif
260
261	firstcall = .false.
262	ENDIF
263
264	! 1.4 Radius middle point and range for Gauss integration
265	radiusm=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) + radiustab(iaer,idomain,1))
266	radiusr=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) - radiustab(iaer,idomain,1))
267
268	! 1.5 Interpolating data at the Gauss quadrature points:
269	DO gausind=1,ngau
270	drad=radiusr*radgaus(gausind)
271	radGAUSa(gausind,iaer,idomain)=radiusm-drad
272
273	radius_id=minloc(abs(radiustab(iaer,idomain,:) - (radiusm-drad)),DIM=1)
274	IF ((radiustab(iaer,idomain,radius_id) - (radiusm-drad)).GT.0) THEN
275	radius_id=radius_id-1
276	ENDIF
277	IF (radius_id.GE.nsize(iaer,idomain)) THEN
278	radius_id=nsize(iaer,idomain)-1
279	kint = 1.
280	ELSEIF (radius_id.LT.1) THEN
281	radius_id=1
282	kint = 0.
283	ELSE
284	kint = ( (radiusm-drad) - &
285	radiustab(iaer,idomain,radius_id) ) / &
286	( radiustab(iaer,idomain,radius_id+1) - &
287	radiustab(iaer,idomain,radius_id) )
288	ENDIF
289	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
290	DO m=1,L_NSPECTV
291	qsqrefVISa(m,gausind,iaer)= &
292	(1-kint)*QVISsQREF(m,iaer,radius_id) + &
293	kint*QVISsQREF(m,iaer,radius_id+1)
294	omegVISa(m,gausind,iaer)= &
295	(1-kint)*omegaVIS(m,iaer,radius_id) + &
296	kint*omegaVIS(m,iaer,radius_id+1)
297	gVISa(m,gausind,iaer)= &
298	(1-kint)*gVIS(m,iaer,radius_id) + &
299	kint*gVIS(m,iaer,radius_id+1)
300	ENDDO
301	qrefVISa(gausind,iaer)= &
302	(1-kint)*QREFvis(iaer,radius_id) + &
303	kint*QREFvis(iaer,radius_id+1)
304	omegrefVISa(gausind,iaer)= 0
305	! omegrefVISa(gausind,iaer)= &
306	! (1-kint)*omegaREFvis(iaer,radius_id) + &
307	! kint*omegaREFvis(iaer,radius_id+1)
308	ELSE ! INFRARED DOMAIN ----------------------------------
309	DO m=1,L_NSPECTI
310	qsqrefIRa(m,gausind,iaer)= &
311	(1-kint)*QIRsQREF(m,iaer,radius_id) + &
312	kint*QIRsQREF(m,iaer,radius_id+1)
313	omegIRa(m,gausind,iaer)= &
314	(1-kint)*omegaIR(m,iaer,radius_id) + &
315	kint*omegaIR(m,iaer,radius_id+1)
316	gIRa(m,gausind,iaer)= &
317	(1-kint)*gIR(m,iaer,radius_id) + &
318	kint*gIR(m,iaer,radius_id+1)
319	ENDDO
320	qrefIRa(gausind,iaer)= &
321	(1-kint)*QREFir(iaer,radius_id) + &
322	kint*QREFir(iaer,radius_id+1)
323	omegrefIRa(gausind,iaer)= &
324	(1-kint)*omegaREFir(iaer,radius_id) + &
325	kint*omegaREFir(iaer,radius_id+1)
326	ENDIF
327	ENDDO
328
329	DO gausind=1,ngau
330	drad=radiusr*radgaus(gausind)
331	radGAUSb(gausind,iaer,idomain)=radiusm+drad
332
333	radius_id=minloc(abs(radiustab(iaer,idomain,:) - &
334	(radiusm+drad)),DIM=1)
335	IF ((radiustab(iaer,idomain,radius_id) - &
336	(radiusm+drad)).GT.0) THEN
337	radius_id=radius_id-1
338	ENDIF
339	IF (radius_id.GE.nsize(iaer,idomain)) THEN
340	radius_id=nsize(iaer,idomain)-1
341	kint = 1.
342	ELSEIF (radius_id.LT.1) THEN
343	radius_id=1
344	kint = 0.
345	ELSE
346	kint = ( (radiusm+drad) - &
347	radiustab(iaer,idomain,radius_id) ) / &
348	( radiustab(iaer,idomain,radius_id+1) - &
349	radiustab(iaer,idomain,radius_id) )
350	ENDIF
351	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
352	DO m=1,L_NSPECTV
353	qsqrefVISb(m,gausind,iaer)= &
354	(1-kint)*QVISsQREF(m,iaer,radius_id) + &
355	kint*QVISsQREF(m,iaer,radius_id+1)
356	omegVISb(m,gausind,iaer)= &
357	(1-kint)*omegaVIS(m,iaer,radius_id) + &
358	kint*omegaVIS(m,iaer,radius_id+1)
359	gVISb(m,gausind,iaer)= &
360	(1-kint)*gVIS(m,iaer,radius_id) + &
361	kint*gVIS(m,iaer,radius_id+1)
362	ENDDO
363	qrefVISb(gausind,iaer)= &
364	(1-kint)*QREFvis(iaer,radius_id) + &
365	kint*QREFvis(iaer,radius_id+1)
366	omegrefVISb(gausind,iaer)= 0
367	! omegrefVISb(gausind,iaer)= &
368	! (1-kint)*omegaREFvis(iaer,radius_id) + &
369	! kint*omegaREFvis(iaer,radius_id+1)
370	ELSE ! INFRARED DOMAIN ----------------------------------
371	DO m=1,L_NSPECTI
372	qsqrefIRb(m,gausind,iaer)= &
373	(1-kint)*QIRsQREF(m,iaer,radius_id) + &
374	kint*QIRsQREF(m,iaer,radius_id+1)
375	omegIRb(m,gausind,iaer)= &
376	(1-kint)*omegaIR(m,iaer,radius_id) + &
377	kint*omegaIR(m,iaer,radius_id+1)
378	gIRb(m,gausind,iaer)= &
379	(1-kint)*gIR(m,iaer,radius_id) + &
380	kint*gIR(m,iaer,radius_id+1)
381	ENDDO
382	qrefIRb(gausind,iaer)= &
383	(1-kint)*QREFir(iaer,radius_id) + &
384	kint*QREFir(iaer,radius_id+1)
385	omegrefIRb(gausind,iaer)= &
386	(1-kint)*omegaREFir(iaer,radius_id) + &
387	kint*omegaREFir(iaer,radius_id+1)
388	ENDIF
389	ENDDO
390
391	!==================================================================
392	! CONSTANT NUEFF FROM HERE
393	!==================================================================
394
395	! 2. Compute the scattering parameters using linear
396	! interpolation over grain sizes and constant nueff
397	! ---------------------------------------------------
398
399	DO lg = 1,nlayer
400	DO ig = 1, ngrid
401	! 2.1 Effective radius index and kx calculation
402	var_tmp=reffrad(ig,lg,iaer)/refftabmin
403	var_tmp=log(var_tmp)*3.
404	var_tmp=var_tmp/logvratgrid+1.
405	grid_i=floor(var_tmp)
406	IF (grid_i.GE.refftabsize) THEN
407	! WRITE(,) 'Warning: particle size in grid box #'
408	! WRITE(,) ig,' is too large to be used by the '
409	! WRITE(,) 'radiative transfer; please extend the '
410	! WRITE(,) 'interpolation grid to larger grain sizes.'
411	grid_i=refftabsize-1
412	kx = 1.
413	ELSEIF (grid_i.LT.1) THEN
414	! WRITE(,) 'Warning: particle size in grid box #'
415	! WRITE(,) ig,' is too small to be used by the '
416	! WRITE(,) 'radiative transfer; please extend the '
417	! WRITE(,) 'interpolation grid to smaller grain sizes.'
418	grid_i=1
419	kx = 0.
420	ELSE
421	kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) / &
422	( refftab(grid_i+1)-refftab(grid_i) )
423	ENDIF
424	! 2.3 Integration
425	DO j=grid_i,grid_i+1
426	! 2.3.1 Check if the calculation has been done
427	IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
428	! 2.3.2 Log-normal dist., r_g and sigma_g are defined
429	! in [hansen_1974], "Light scattering in planetary
430	! atmospheres", Space Science Reviews 16 527-610.
431	! Here, sizedistk1=r_g and sizedistk2=sigma_g^2
432	sizedistk2 = log(1.+nueffrad(1,1,iaer))
433	sizedistk1 = exp(2.5*sizedistk2)
434	sizedistk1 = refftab(j) / sizedistk1
435
436	normd(j,1,iaer,idomain) = 1e-30
437	DO gausind=1,ngau
438	drad=radiusr*radgaus(gausind)
439	dista(j,1,iaer,idomain,gausind) = LOG((radiusm-drad)/sizedistk1)
440	dista(j,1,iaer,idomain,gausind) = &
441	EXP(-dista(j,1,iaer,idomain,gausind) * &
442	dista(j,1,iaer,idomain,gausind) * &
443	0.5e0/sizedistk2)/(radiusm-drad)
444	dista(j,1,iaer,idomain,gausind) = &
445	dista(j,1,iaer,idomain,gausind) / &
446	(sqrt(2e0pisizedistk2))
447
448	distb(j,1,iaer,idomain,gausind) = LOG((radiusm+drad)/sizedistk1)
449	distb(j,1,iaer,idomain,gausind) = &
450	EXP(-distb(j,1,iaer,idomain,gausind) * &
451	distb(j,1,iaer,idomain,gausind) * &
452	0.5e0/sizedistk2)/(radiusm+drad)
453	distb(j,1,iaer,idomain,gausind) = &
454	distb(j,1,iaer,idomain,gausind) / &
455	(sqrt(2e0pisizedistk2))
456
457	normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) + &
458	weightgaus(gausind) * &
459	( &
460	distb(j,1,iaer,idomain,gausind) * pi * &
461	radGAUSb(gausind,iaer,idomain) * &
462	radGAUSb(gausind,iaer,idomain) + &
463	dista(j,1,iaer,idomain,gausind) * pi * &
464	radGAUSa(gausind,iaer,idomain) * &
465	radGAUSa(gausind,iaer,idomain) &
466	)
467	ENDDO
468	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
469	! 2.3.3.vis Initialization
470	qsqrefVISgrid(j,1,:,iaer)=0.
471	qextVISgrid(j,1,:,iaer)=0.
472	qscatVISgrid(j,1,:,iaer)=0.
473	omegVISgrid(j,1,:,iaer)=0.
474	gVISgrid(j,1,:,iaer)=0.
475	qrefVISgrid(j,1,iaer)=0.
476	qscatrefVISgrid(j,1,iaer)=0.
477	omegrefVISgrid(j,1,iaer)=0.
478
479	DO gausind=1,ngau
480	DO m=1,L_NSPECTV
481	! Convolution:
482	qextVISgrid(j,1,m,iaer) = &
483	qextVISgrid(j,1,m,iaer) + &
484	weightgaus(gausind) * &
485	( &
486	qsqrefVISb(m,gausind,iaer) * &
487	qrefVISb(gausind,iaer) * &
488	piradGAUSb(gausind,iaer,idomain) &
489	radGAUSb(gausind,iaer,idomain) * &
490	distb(j,1,iaer,idomain,gausind) + &
491	qsqrefVISa(m,gausind,iaer) * &
492	qrefVISa(gausind,iaer) * &
493	piradGAUSa(gausind,iaer,idomain) &
494	radGAUSa(gausind,iaer,idomain) * &
495	dista(j,1,iaer,idomain,gausind) &
496	)
497	qscatVISgrid(j,1,m,iaer) = &
498	qscatVISgrid(j,1,m,iaer) + &
499	weightgaus(gausind) * &
500	( &
501	omegVISb(m,gausind,iaer) * &
502	qsqrefVISb(m,gausind,iaer) * &
503	qrefVISb(gausind,iaer) * &
504	piradGAUSb(gausind,iaer,idomain) &
505	radGAUSb(gausind,iaer,idomain) * &
506	distb(j,1,iaer,idomain,gausind) + &
507	omegVISa(m,gausind,iaer) * &
508	qsqrefVISa(m,gausind,iaer) * &
509	qrefVISa(gausind,iaer) * &
510	piradGAUSa(gausind,iaer,idomain) &
511	radGAUSa(gausind,iaer,idomain) * &
512	dista(j,1,iaer,idomain,gausind) &
513	)
514	gVISgrid(j,1,m,iaer) = &
515	gVISgrid(j,1,m,iaer) + &
516	weightgaus(gausind) * &
517	( &
518	omegVISb(m,gausind,iaer) * &
519	qsqrefVISb(m,gausind,iaer) * &
520	qrefVISb(gausind,iaer) * &
521	gVISb(m,gausind,iaer) * &
522	piradGAUSb(gausind,iaer,idomain) &
523	radGAUSb(gausind,iaer,idomain) * &
524	distb(j,1,iaer,idomain,gausind) + &
525	omegVISa(m,gausind,iaer) * &
526	qsqrefVISa(m,gausind,iaer) * &
527	qrefVISa(gausind,iaer) * &
528	gVISa(m,gausind,iaer) * &
529	piradGAUSa(gausind,iaer,idomain) &
530	radGAUSa(gausind,iaer,idomain) * &
531	dista(j,1,iaer,idomain,gausind) &
532	)
533	ENDDO
534	qrefVISgrid(j,1,iaer) = &
535	qrefVISgrid(j,1,iaer) + &
536	weightgaus(gausind) * &
537	( &
538	qrefVISb(gausind,iaer) * &
539	piradGAUSb(gausind,iaer,idomain) &
540	radGAUSb(gausind,iaer,idomain) * &
541	distb(j,1,iaer,idomain,gausind) + &
542	qrefVISa(gausind,iaer) * &
543	piradGAUSa(gausind,iaer,idomain) &
544	radGAUSa(gausind,iaer,idomain) * &
545	dista(j,1,iaer,idomain,gausind) &
546	)
547	qscatrefVISgrid(j,1,iaer) = &
548	qscatrefVISgrid(j,1,iaer) + &
549	weightgaus(gausind) * &
550	( &
551	omegrefVISb(gausind,iaer) * &
552	qrefVISb(gausind,iaer) * &
553	piradGAUSb(gausind,iaer,idomain) &
554	radGAUSb(gausind,iaer,idomain) * &
555	distb(j,1,iaer,idomain,gausind) + &
556	omegrefVISa(gausind,iaer) * &
557	qrefVISa(gausind,iaer) * &
558	piradGAUSa(gausind,iaer,idomain) &
559	radGAUSa(gausind,iaer,idomain) * &
560	dista(j,1,iaer,idomain,gausind) &
561	)
562	ENDDO
563
564	qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) / &
565	normd(j,1,iaer,idomain)
566	qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / &
567	normd(j,1,iaer,idomain)
568	omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / &
569	qrefVISgrid(j,1,iaer)
570	DO m=1,L_NSPECTV
571	qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / &
572	normd(j,1,iaer,idomain)
573	qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / &
574	normd(j,1,iaer,idomain)
575	gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) / &
576	qscatVISgrid(j,1,m,iaer) / &
577	normd(j,1,iaer,idomain)
578
579	qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / &
580	qrefVISgrid(j,1,iaer)
581	omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / &
582	qextVISgrid(j,1,m,iaer)
583	ENDDO
584	ELSE ! INFRARED DOMAIN ----------
585	! 2.3.3.ir Initialization
586	qsqrefIRgrid(j,1,:,iaer)=0.
587	qextIRgrid(j,1,:,iaer)=0.
588	qscatIRgrid(j,1,:,iaer)=0.
589	omegIRgrid(j,1,:,iaer)=0.
590	gIRgrid(j,1,:,iaer)=0.
591	qrefIRgrid(j,1,iaer)=0.
592	qscatrefIRgrid(j,1,iaer)=0.
593	omegrefIRgrid(j,1,iaer)=0.
594
595	DO gausind=1,ngau
596	DO m=1,L_NSPECTI
597	! Convolution:
598	qextIRgrid(j,1,m,iaer) = &
599	qextIRgrid(j,1,m,iaer) + &
600	weightgaus(gausind) * &
601	( &
602	qsqrefIRb(m,gausind,iaer) * &
603	qrefVISb(gausind,iaer) * &
604	piradGAUSb(gausind,iaer,idomain) &
605	radGAUSb(gausind,iaer,idomain) * &
606	distb(j,1,iaer,idomain,gausind) + &
607	qsqrefIRa(m,gausind,iaer) * &
608	qrefVISa(gausind,iaer) * &
609	piradGAUSa(gausind,iaer,idomain) &
610	radGAUSa(gausind,iaer,idomain) * &
611	dista(j,1,iaer,idomain,gausind) &
612	)
613	qscatIRgrid(j,1,m,iaer) = &
614	qscatIRgrid(j,1,m,iaer) + &
615	weightgaus(gausind) * &
616	( &
617	omegIRb(m,gausind,iaer) * &
618	qsqrefIRb(m,gausind,iaer) * &
619	qrefVISb(gausind,iaer) * &
620	piradGAUSb(gausind,iaer,idomain) &
621	radGAUSb(gausind,iaer,idomain) * &
622	distb(j,1,iaer,idomain,gausind) + &
623	omegIRa(m,gausind,iaer) * &
624	qsqrefIRa(m,gausind,iaer) * &
625	qrefVISa(gausind,iaer) * &
626	piradGAUSa(gausind,iaer,idomain) &
627	radGAUSa(gausind,iaer,idomain) * &
628	dista(j,1,iaer,idomain,gausind) &
629	)
630	gIRgrid(j,1,m,iaer) = &
631	gIRgrid(j,1,m,iaer) + &
632	weightgaus(gausind) * &
633	( &
634	omegIRb(m,gausind,iaer) * &
635	qsqrefIRb(m,gausind,iaer) * &
636	qrefVISb(gausind,iaer) * &
637	gIRb(m,gausind,iaer) * &
638	piradGAUSb(gausind,iaer,idomain) &
639	radGAUSb(gausind,iaer,idomain) * &
640	distb(j,1,iaer,idomain,gausind) + &
641	omegIRa(m,gausind,iaer) * &
642	qsqrefIRa(m,gausind,iaer) * &
643	qrefVISa(gausind,iaer) * &
644	gIRa(m,gausind,iaer) * &
645	piradGAUSa(gausind,iaer,idomain) &
646	radGAUSa(gausind,iaer,idomain) * &
647	dista(j,1,iaer,idomain,gausind) &
648	)
649	ENDDO
650	qrefIRgrid(j,1,iaer) = &
651	qrefIRgrid(j,1,iaer) + &
652	weightgaus(gausind) * &
653	( &
654	qrefIRb(gausind,iaer) * &
655	piradGAUSb(gausind,iaer,idomain) &
656	radGAUSb(gausind,iaer,idomain) * &
657	distb(j,1,iaer,idomain,gausind) + &
658	qrefIRa(gausind,iaer) * &
659	piradGAUSa(gausind,iaer,idomain) &
660	radGAUSa(gausind,iaer,idomain) * &
661	dista(j,1,iaer,idomain,gausind) &
662	)
663	qscatrefIRgrid(j,1,iaer) = &
664	qscatrefIRgrid(j,1,iaer) + &
665	weightgaus(gausind) * &
666	( &
667	omegrefIRb(gausind,iaer) * &
668	qrefIRb(gausind,iaer) * &
669	piradGAUSb(gausind,iaer,idomain) &
670	radGAUSb(gausind,iaer,idomain) * &
671	distb(j,1,iaer,idomain,gausind) + &
672	omegrefIRa(gausind,iaer) * &
673	qrefIRa(gausind,iaer) * &
674	piradGAUSa(gausind,iaer,idomain) &
675	radGAUSa(gausind,iaer,idomain) * &
676	dista(j,1,iaer,idomain,gausind) &
677	)
678	ENDDO
679
680	qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) / &
681	normd(j,1,iaer,idomain)
682	qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / &
683	normd(j,1,iaer,idomain)
684	omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / &
685	qrefIRgrid(j,1,iaer)
686	DO m=1,L_NSPECTI
687	qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / &
688	normd(j,1,iaer,idomain)
689	qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / &
690	normd(j,1,iaer,idomain)
691	gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) / &
692	qscatIRgrid(j,1,m,iaer) / &
693	normd(j,1,iaer,idomain)
694
695	qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / &
696	qrefVISgrid(j,1,iaer)
697	omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / &
698	qextIRgrid(j,1,m,iaer)
699	ENDDO
700	ENDIF ! --------------------------
701	checkgrid(j,1,iaer,idomain) = .true.
702	ENDIF !checkgrid
703	ENDDO !grid_i
704	! 2.4 Linear interpolation
705	k1 = (1-kx)
706	k2 = kx
707	IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
708	DO m=1,L_NSPECTV
709	QVISsQREF3d(ig,lg,m,iaer) = &
710	k1*qsqrefVISgrid(grid_i,1,m,iaer) + &
711	k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
712	omegaVIS3d(ig,lg,m,iaer) = &
713	k1*omegVISgrid(grid_i,1,m,iaer) + &
714	k2*omegVISgrid(grid_i+1,1,m,iaer)
715	gVIS3d(ig,lg,m,iaer) = &
716	k1*gVISgrid(grid_i,1,m,iaer) + &
717	k2*gVISgrid(grid_i+1,1,m,iaer)
718	ENDDO !L_NSPECTV
719	QREFvis3d(ig,lg,iaer) = &
720	k1*qrefVISgrid(grid_i,1,iaer) + &
721	k2*qrefVISgrid(grid_i+1,1,iaer)
722	! omegaREFvis3d(ig,lg,iaer) = &
723	! k1*omegrefVISgrid(grid_i,1,iaer) + &
724	! k2*omegrefVISgrid(grid_i+1,1,iaer)
725	ELSE ! INFRARED -----------------------
726	DO m=1,L_NSPECTI
727	QIRsQREF3d(ig,lg,m,iaer) = &
728	k1*qsqrefIRgrid(grid_i,1,m,iaer) + &
729	k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
730	omegaIR3d(ig,lg,m,iaer) = &
731	k1*omegIRgrid(grid_i,1,m,iaer) + &
732	k2*omegIRgrid(grid_i+1,1,m,iaer)
733	gIR3d(ig,lg,m,iaer) = &
734	k1*gIRgrid(grid_i,1,m,iaer) + &
735	k2*gIRgrid(grid_i+1,1,m,iaer)
736	ENDDO !L_NSPECTI
737	QREFir3d(ig,lg,iaer) = &
738	k1*qrefIRgrid(grid_i,1,iaer) + &
739	k2*qrefIRgrid(grid_i+1,1,iaer)
740	! omegaREFir3d(ig,lg,iaer) = &
741	! k1*omegrefIRgrid(grid_i,1,iaer) + &
742	! k2*omegrefIRgrid(grid_i+1,1,iaer)
743	ENDIF ! --------------------------------
744	ENDDO !nlayer
745	ENDDO !ngrid
746
747	!==================================================================
748
749
750
751	ENDDO ! idomain
752
753	ENDIF ! nsize = 1
754
755	ENDDO ! iaer (loop on aerosol kind)
756
757	RETURN
758	END SUBROUTINE aeroptproperties
759
760
761	end module aeroptproperties_mod

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

Download in other formats:

Original Format