Context Navigation

aeroptproperties.F90 @ 1834

Last change on this file since 1834 was 1529, checked in by emillour, 9 years ago

Generic GCM: Towards a cleaner separation between physics and dynamics

Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
Removed module "comhdiff_mod.F90", as it is only used by module surf_heat_transp_mod.F90, moved module variables there.
Addedin "surf_heat_transp_mod" local versions of some arrays and routines (from dyn3d) required to compute gradient, divergence, etc. on the global dynamics grid. As before, the slab ocean only works in serial.

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

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

Download in other formats:

Original Format