Context Navigation

aeroptproperties.F90 @ 1351

Last change on this file since 1351 was 1315, checked in by milmd, 10 years ago
LMDZ.GENERIC. OpenMP directives added in generic physic. When running in pure OpenMP or hybrid OpenMP/MPI, may have some bugs with condense_cloud and wstats routines.
File size: 36.5 KB

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

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

Download in other formats:

Original Format