Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

aeroptproperties.F90

Last change on this file was 227, checked in by milmd, 10 years ago
Last LMDZ version (1315) with OpenMP directives and other stuff
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.

Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/aeroptproperties.F90

Download in other formats: