Context Navigation

aeroptproperties.F @ 746

Last change on this file since 746 was 737, checked in by jbmadeleine, 13 years ago

Corrected aeroptproperties.F to avoid having a floating underflow in
aeroptproperties.F due to a too small minimum radius of the
interpolation grid. This can lead to wrong values of the scattering
parameters or even NaNs?.

refftabmin and refftabmax are now computed automatically.

Previously, the user had to specify them, which is dangerous for the
above mentioned reason and not flexible in the long run.

Added two dimensions (iaer and idomain) to refftab and nuefftab

so that the grid can be changed depending on the range of radiuses
contained in the lookup tables for each aerosol.

Also added two dimensions to logvratgrid and delete vratgrid, which

is only used once.

Added a message at firstcall giving the range of radiuses used for

the interpolation grid or each aerosol.

normd is now checked to make sure that it is not equal to zero (in

which case there is a division by zero at the end of the routine).
This should not occur, but we never know...

File size: 55.0 KB

Line
1	SUBROUTINE aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
2	& QVISsQREF3d,omegaVIS3d,gVIS3d,
3	& QIRsQREF3d,omegaIR3d,gIR3d,
4	& QREFvis3d,QREFir3d,
5	& omegaREFvis3d,omegaREFir3d)
6	IMPLICIT NONE
7	c =============================================================
8	c Aerosol Optical Properties
9	c
10	c Description:
11	c Compute the scattering parameters in each grid
12	c box, depending on aerosol grain sizes. Log-normal size
13	c distribution and Gauss-Legendre integration are used.
14
15	c Parameters:
16	c Don't forget to set the value of varyingnueff below; If
17	c the effective variance of the distribution for the given
18	c aerosol is considered homogeneous in the atmosphere, please
19	c set varyingnueff(iaer) to .false. Resulting computational
20	c time will be much better.
21
22	c Authors: J.-B. Madeleine, F. Forget, F. Montmessin
23	c =============================================================
24
25	#include "dimensions.h"
26	#include "dimphys.h"
27	#include "callkeys.h"
28	#include "dimradmars.h"
29	#include "yomaer.h"
30
31	c Local variables
32	c ---------------
33
34	c =============================================================
35	LOGICAL, PARAMETER :: varyingnueff(naerkind) = .false.
36	c =============================================================
37
38	c Min. and max radius of the interpolation grid (in METERS)
39	REAL, SAVE :: refftabmin(naerkind,2)
40	REAL, SAVE :: refftabmax(naerkind,2)
41	c Log of the min and max variance of the interpolation grid
42	REAL, PARAMETER :: nuefftabmin = -4.6
43	REAL, PARAMETER :: nuefftabmax = 0.
44	c Number of effective radius of the interpolation grid
45	INTEGER, PARAMETER :: refftabsize = 100
46	c Number of effective variances of the interpolation grid
47	INTEGER, PARAMETER :: nuefftabsize = 100
48	c Interpolation grid indices (reff,nueff)
49	INTEGER :: grid_i,grid_j
50	c Intermediate variable
51	REAL :: var_tmp,var3d_tmp(ngridmx,nlayermx)
52	c Bilinear interpolation factors
53	REAL :: kx,ky,k1,k2,k3,k4
54	c Size distribution parameters
55	REAL :: sizedistk1,sizedistk2
56	c Pi!
57	REAL,SAVE :: pi
58	c Variables used by the Gauss-Legendre integration:
59	INTEGER radius_id,gausind
60	REAL kint
61	REAL drad
62	INTEGER, PARAMETER :: ngau = 10
63	REAL weightgaus(ngau),radgaus(ngau)
64	SAVE weightgaus,radgaus
65	c DATA weightgaus/.2955242247,.2692667193,.2190863625,
66	c & .1494513491,.0666713443/
67	c DATA radgaus/.1488743389,.4333953941,.6794095682,
68	c & .8650633666,.9739065285/
69	DATA radgaus/0.07652652113350,0.22778585114165,
70	& 0.37370608871528,0.51086700195146,
71	& 0.63605368072468,0.74633190646476,
72	& 0.83911697181213,0.91223442826796,
73	& 0.96397192726078,0.99312859919241/
74
75	DATA weightgaus/0.15275338723120,0.14917298659407,
76	& 0.14209610937519,0.13168863843930,
77	& 0.11819453196154,0.10193011980823,
78	& 0.08327674160932,0.06267204829828,
79	& 0.04060142982019,0.01761400714091/
80
81	c Indices
82	INTEGER :: i,j,k,l,m,iaer,idomain
83	INTEGER :: ig,lg,chg
84
85	c Local saved variables
86	c ---------------------
87
88	c Radius axis of the interpolation grid
89	REAL,SAVE :: refftab(refftabsize,naerkind,2)
90	c Variance axis of the interpolation grid
91	REAL,SAVE :: nuefftab(nuefftabsize,naerkind,2)
92	c Volume ratio of the grid
93	REAL,SAVE :: logvratgrid(naerkind,2)
94	c Grid used to remember which calculation is done
95	LOGICAL,SAVE :: checkgrid(refftabsize,nuefftabsize,naerkind,2)
96	& = .false.
97	c Optical properties of the grid (VISIBLE)
98	REAL,SAVE :: qsqrefVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
99	REAL,SAVE :: qextVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
100	REAL,SAVE :: qscatVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
101	REAL,SAVE :: omegVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
102	REAL,SAVE :: gVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
103	c Optical properties of the grid (INFRARED)
104	REAL,SAVE :: qsqrefIRgrid(refftabsize,nuefftabsize,nir,naerkind)
105	REAL,SAVE :: qextIRgrid(refftabsize,nuefftabsize,nir,naerkind)
106	REAL,SAVE :: qscatIRgrid(refftabsize,nuefftabsize,nir,naerkind)
107	REAL,SAVE :: omegIRgrid(refftabsize,nuefftabsize,nir,naerkind)
108	REAL,SAVE :: gIRgrid(refftabsize,nuefftabsize,nir,naerkind)
109	c Optical properties of the grid (REFERENCE WAVELENGTHS)
110	REAL,SAVE :: qrefVISgrid(refftabsize,nuefftabsize,naerkind)
111	REAL,SAVE :: qscatrefVISgrid(refftabsize,nuefftabsize,naerkind)
112	REAL,SAVE :: qrefIRgrid(refftabsize,nuefftabsize,naerkind)
113	REAL,SAVE :: qscatrefIRgrid(refftabsize,nuefftabsize,naerkind)
114	REAL,SAVE :: omegrefVISgrid(refftabsize,nuefftabsize,naerkind)
115	REAL,SAVE :: omegrefIRgrid(refftabsize,nuefftabsize,naerkind)
116	c Firstcall
117	LOGICAL,SAVE :: firstcall = .true.
118	c Variables used by the Gauss-Legendre integration:
119	REAL,SAVE :: normd(refftabsize,nuefftabsize,naerkind,2)
120	REAL,SAVE :: dista(refftabsize,nuefftabsize,naerkind,2,ngau)
121	REAL,SAVE :: distb(refftabsize,nuefftabsize,naerkind,2,ngau)
122
123	REAL,SAVE :: radGAUSa(ngau,naerkind,2)
124	REAL,SAVE :: radGAUSb(ngau,naerkind,2)
125
126	REAL,SAVE :: qsqrefVISa(nsun,ngau,naerkind)
127	REAL,SAVE :: qrefVISa(ngau,naerkind)
128	REAL,SAVE :: qsqrefVISb(nsun,ngau,naerkind)
129	REAL,SAVE :: qrefVISb(ngau,naerkind)
130	REAL,SAVE :: omegVISa(nsun,ngau,naerkind)
131	REAL,SAVE :: omegrefVISa(ngau,naerkind)
132	REAL,SAVE :: omegVISb(nsun,ngau,naerkind)
133	REAL,SAVE :: omegrefVISb(ngau,naerkind)
134	REAL,SAVE :: gVISa(nsun,ngau,naerkind)
135	REAL,SAVE :: gVISb(nsun,ngau,naerkind)
136
137	REAL,SAVE :: qsqrefIRa(nir,ngau,naerkind)
138	REAL,SAVE :: qrefIRa(ngau,naerkind)
139	REAL,SAVE :: qsqrefIRb(nir,ngau,naerkind)
140	REAL,SAVE :: qrefIRb(ngau,naerkind)
141	REAL,SAVE :: omegIRa(nir,ngau,naerkind)
142	REAL,SAVE :: omegrefIRa(ngau,naerkind)
143	REAL,SAVE :: omegIRb(nir,ngau,naerkind)
144	REAL,SAVE :: omegrefIRb(ngau,naerkind)
145	REAL,SAVE :: gIRa(nir,ngau,naerkind)
146	REAL,SAVE :: gIRb(nir,ngau,naerkind)
147
148	REAL :: radiusm
149	REAL :: radiusr
150
151	c Inputs
152	c ------
153
154	INTEGER :: ngrid,nlayer
155	c Aerosol effective radius used for radiative transfer (meter)
156	REAL :: reffrad(ngridmx,nlayermx,naerkind)
157	c Aerosol effective variance used for radiative transfer (n.u.)
158	REAL :: nueffrad(ngridmx,nlayermx,naerkind)
159
160	c Outputs
161	c -------
162
163	REAL :: QVISsQREF3d(ngridmx,nlayermx,nsun,naerkind)
164	REAL :: omegaVIS3d(ngridmx,nlayermx,nsun,naerkind)
165	REAL :: gVIS3d(ngridmx,nlayermx,nsun,naerkind)
166
167	REAL :: QIRsQREF3d(ngridmx,nlayermx,nir,naerkind)
168	REAL :: omegaIR3d(ngridmx,nlayermx,nir,naerkind)
169	REAL :: gIR3d(ngridmx,nlayermx,nir,naerkind)
170
171	REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
172	REAL :: QREFir3d(ngridmx,nlayermx,naerkind)
173
174	REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
175	REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind)
176
177	c Tests
178	c -----
179
180	LOGICAL,SAVE :: out_qwg = .false.
181	INTEGER, PARAMETER :: out_iaer = 1
182	REAL :: out_qext(ngridmx,nlayermx)
183	REAL :: out_omeg(ngridmx,nlayermx)
184	REAL :: out_g(ngridmx,nlayermx)
185	INTEGER :: out_nchannel
186	CHARACTER*1 :: out_str
187
188	c Creating the effective radius and variance grid
189	c -----------------------------------------------
190	IF (firstcall) THEN
191	c 0.1 Pi!
192	pi = 2. * asin(1.e0)
193
194	WRITE(,) "aeroptproperties: interpolation grid"
195	DO iaer = 1, naerkind ! Loop on aerosol kind
196	DO idomain = 1, 2 ! Loop on visible or infrared channel
197
198	c 0.2 Effective radius
199	radiusm=
200	& 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain))+
201	& radiustab(iaer,idomain,1))
202	radiusr=
203	& 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain))-
204	& radiustab(iaer,idomain,1))
205	refftabmin(iaer,idomain) =
206	& radiusm-radiusr*radgaus(ngau)
207	refftabmax(iaer,idomain) =
208	& radiustab(iaer,idomain,nsize(iaer,idomain))
209
210	WRITE(,) "Scatterer: ",iaer
211	WRITE(,) "Domain: ",idomain
212	WRITE(,) "Min radius (m): ", refftabmin(iaer,idomain)
213	WRITE(,) "Max radius (m): ", refftabmax(iaer,idomain)
214
215	refftab(1,iaer,idomain) =
216	& refftabmin(iaer,idomain)
217	refftab(refftabsize,iaer,idomain) =
218	& refftabmax(iaer,idomain)
219
220	logvratgrid(iaer,idomain) =
221	& log(refftabmax(iaer,idomain)/
222	& refftabmin(iaer,idomain)) /
223	& float(refftabsize-1)*3.
224	do i = 2, refftabsize-1
225	refftab(i,iaer,idomain) = refftab(i-1,iaer,idomain)*
226	& exp(1./3.*logvratgrid(iaer,idomain))
227	enddo
228
229	c 0.3 Effective variance
230	do i = 0, nuefftabsize-1
231	nuefftab(i+1,iaer,idomain) = exp( nuefftabmin +
232	& i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
233	enddo
234
235	ENDDO
236	ENDDO
237	firstcall = .false.
238	ENDIF
239
240	DO iaer = 1, naerkind ! Loop on aerosol kind
241	IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN
242	c==================================================================
243	c If there is one single particle size, optical
244	c properties of the considered aerosol are homogeneous
245	DO lg = 1, nlayer
246	DO ig = 1, ngrid
247	DO chg = 1, nsun
248	QVISsQREF3d(ig,lg,chg,iaer)=QVISsQREF(chg,iaer,1)
249	omegaVIS3d(ig,lg,chg,iaer)=omegaVIS(chg,iaer,1)
250	gVIS3d(ig,lg,chg,iaer)=gVIS(chg,iaer,1)
251	ENDDO
252	DO chg = 1, nir
253	QIRsQREF3d(ig,lg,chg,iaer)=QIRsQREF(chg,iaer,1)
254	omegaIR3d(ig,lg,chg,iaer)=omegaIR(chg,iaer,1)
255	gIR3d(ig,lg,chg,iaer)=gIR(chg,iaer,1)
256	ENDDO
257	QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1)
258	QREFir3d(ig,lg,iaer)=QREFir(iaer,1)
259	omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1)
260	omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1)
261	ENDDO
262	ENDDO
263	ELSE ! Varying effective radius and variance
264	DO idomain = 1, 2 ! Loop on visible or infrared channel
265	c==================================================================
266
267	c 1.1 Radius middle point and range for Gauss integration
268	radiusm=
269	& 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) +
270	& radiustab(iaer,idomain,1))
271	radiusr=
272	& 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) -
273	& radiustab(iaer,idomain,1))
274
275	c 1.2 Interpolating data at the Gauss quadrature points:
276	DO gausind=1,ngau
277	drad=radiusr*radgaus(gausind)
278	radGAUSa(gausind,iaer,idomain)=radiusm-drad
279
280	radius_id=minloc(abs(radiustab(iaer,idomain,:) -
281	& (radiusm-drad)),DIM=1)
282	IF ((radiustab(iaer,idomain,radius_id) -
283	& (radiusm-drad)).GT.0) THEN
284	radius_id=radius_id-1
285	ENDIF
286	IF (radius_id.GE.nsize(iaer,idomain)) THEN
287	radius_id=nsize(iaer,idomain)-1
288	kint = 1.
289	ELSEIF (radius_id.LT.1) THEN
290	radius_id=1
291	kint = 0.
292	ELSE
293	kint = ( (radiusm-drad) -
294	& radiustab(iaer,idomain,radius_id) ) /
295	& ( radiustab(iaer,idomain,radius_id+1) -
296	& radiustab(iaer,idomain,radius_id) )
297	ENDIF
298	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
299	DO m=1,nsun
300	qsqrefVISa(m,gausind,iaer)=
301	& (1-kint)*QVISsQREF(m,iaer,radius_id) +
302	& kint*QVISsQREF(m,iaer,radius_id+1)
303	omegVISa(m,gausind,iaer)=
304	& (1-kint)*omegaVIS(m,iaer,radius_id) +
305	& kint*omegaVIS(m,iaer,radius_id+1)
306	gVISa(m,gausind,iaer)=
307	& (1-kint)*gVIS(m,iaer,radius_id) +
308	& kint*gVIS(m,iaer,radius_id+1)
309	ENDDO
310	qrefVISa(gausind,iaer)=
311	& (1-kint)*QREFvis(iaer,radius_id) +
312	& kint*QREFvis(iaer,radius_id+1)
313	omegrefVISa(gausind,iaer)=
314	& (1-kint)*omegaREFvis(iaer,radius_id) +
315	& kint*omegaREFvis(iaer,radius_id+1)
316	ELSE ! INFRARED DOMAIN ----------------------------------
317	DO m=1,nir
318	qsqrefIRa(m,gausind,iaer)=
319	& (1-kint)*QIRsQREF(m,iaer,radius_id) +
320	& kint*QIRsQREF(m,iaer,radius_id+1)
321	omegIRa(m,gausind,iaer)=
322	& (1-kint)*omegaIR(m,iaer,radius_id) +
323	& kint*omegaIR(m,iaer,radius_id+1)
324	gIRa(m,gausind,iaer)=
325	& (1-kint)*gIR(m,iaer,radius_id) +
326	& kint*gIR(m,iaer,radius_id+1)
327	ENDDO
328	qrefIRa(gausind,iaer)=
329	& (1-kint)*QREFir(iaer,radius_id) +
330	& kint*QREFir(iaer,radius_id+1)
331	omegrefIRa(gausind,iaer)=
332	& (1-kint)*omegaREFir(iaer,radius_id) +
333	& kint*omegaREFir(iaer,radius_id+1)
334	ENDIF
335	ENDDO
336
337	DO gausind=1,ngau
338	drad=radiusr*radgaus(gausind)
339	radGAUSb(gausind,iaer,idomain)=radiusm+drad
340
341	radius_id=minloc(abs(radiustab(iaer,idomain,:) -
342	& (radiusm+drad)),DIM=1)
343	IF ((radiustab(iaer,idomain,radius_id) -
344	& (radiusm+drad)).GT.0) THEN
345	radius_id=radius_id-1
346	ENDIF
347	IF (radius_id.GE.nsize(iaer,idomain)) THEN
348	radius_id=nsize(iaer,idomain)-1
349	kint = 1.
350	ELSEIF (radius_id.LT.1) THEN
351	radius_id=1
352	kint = 0.
353	ELSE
354	kint = ( (radiusm+drad) -
355	& radiustab(iaer,idomain,radius_id) ) /
356	& ( radiustab(iaer,idomain,radius_id+1) -
357	& radiustab(iaer,idomain,radius_id) )
358	ENDIF
359	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
360	DO m=1,nsun
361	qsqrefVISb(m,gausind,iaer)=
362	& (1-kint)*QVISsQREF(m,iaer,radius_id) +
363	& kint*QVISsQREF(m,iaer,radius_id+1)
364	omegVISb(m,gausind,iaer)=
365	& (1-kint)*omegaVIS(m,iaer,radius_id) +
366	& kint*omegaVIS(m,iaer,radius_id+1)
367	gVISb(m,gausind,iaer)=
368	& (1-kint)*gVIS(m,iaer,radius_id) +
369	& kint*gVIS(m,iaer,radius_id+1)
370	ENDDO
371	qrefVISb(gausind,iaer)=
372	& (1-kint)*QREFvis(iaer,radius_id) +
373	& kint*QREFvis(iaer,radius_id+1)
374	omegrefVISb(gausind,iaer)=
375	& (1-kint)*omegaREFvis(iaer,radius_id) +
376	& kint*omegaREFvis(iaer,radius_id+1)
377	ELSE ! INFRARED DOMAIN ----------------------------------
378	DO m=1,nir
379	qsqrefIRb(m,gausind,iaer)=
380	& (1-kint)*QIRsQREF(m,iaer,radius_id) +
381	& kint*QIRsQREF(m,iaer,radius_id+1)
382	omegIRb(m,gausind,iaer)=
383	& (1-kint)*omegaIR(m,iaer,radius_id) +
384	& kint*omegaIR(m,iaer,radius_id+1)
385	gIRb(m,gausind,iaer)=
386	& (1-kint)*gIR(m,iaer,radius_id) +
387	& kint*gIR(m,iaer,radius_id+1)
388	ENDDO
389	qrefIRb(gausind,iaer)=
390	& (1-kint)*QREFir(iaer,radius_id) +
391	& kint*QREFir(iaer,radius_id+1)
392	omegrefIRb(gausind,iaer)=
393	& (1-kint)*omegaREFir(iaer,radius_id) +
394	& kint*omegaREFir(iaer,radius_id+1)
395	ENDIF
396	ENDDO
397	c==================================================================
398	IF ( .NOT.varyingnueff(iaer) ) THEN ! CONSTANT NUEFF
399	c==================================================================
400	c 2. Compute the scattering parameters using linear
401	c interpolation over grain sizes and constant nueff
402	c ---------------------------------------------------
403
404	DO lg = 1,nlayer
405	DO ig = 1,ngrid
406	c 2.1 Effective radius index and kx calculation
407	var_tmp=reffrad(ig,lg,iaer)/refftabmin(iaer,idomain)
408	var_tmp=log(var_tmp)*3.
409	var_tmp=var_tmp/logvratgrid(iaer,idomain)+1.
410	grid_i=floor(var_tmp)
411	IF (grid_i.GE.refftabsize) THEN
412	c WRITE(,) 'Warning: particle size in grid box #'
413	c WRITE(,) ig,' is too large to be used by the '
414	c WRITE(,) 'radiative transfer; please extend the '
415	c WRITE(,) 'interpolation grid to larger grain sizes.'
416	grid_i=refftabsize-1
417	kx = 1.
418	ELSEIF (grid_i.LT.1) THEN
419	c WRITE(,) 'Warning: particle size in grid box #'
420	c WRITE(,) ig,' is too small to be used by the '
421	c WRITE(,) 'radiative transfer; please extend the '
422	c WRITE(,) 'interpolation grid to smaller grain sizes.'
423	grid_i=1
424	kx = 0.
425	ELSE
426	kx = ( reffrad(ig,lg,iaer)-
427	& refftab(grid_i,iaer,idomain) ) /
428	& ( refftab(grid_i+1,iaer,idomain)-
429	& refftab(grid_i,iaer,idomain) )
430	ENDIF
431	c 2.3 Integration
432	DO j=grid_i,grid_i+1
433	c 2.3.1 Check if the calculation has been done
434	IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
435	c 2.3.2 Log-normal dist., r_g and sigma_g are defined
436	c in [hansen_1974], "Light scattering in planetary
437	c atmospheres", Space Science Reviews 16 527-610.
438	c Here, sizedistk1=r_g and sizedistk2=sigma_g^2
439	sizedistk2 = log(1.+nueffrad(1,1,iaer))
440	sizedistk1 = exp(2.5*sizedistk2)
441	sizedistk1 = refftab(j,iaer,idomain) / sizedistk1
442
443	normd(j,1,iaer,idomain) = 1e-30
444	DO gausind=1,ngau
445	drad=radiusr*radgaus(gausind)
446	dista(j,1,iaer,idomain,gausind) =
447	& LOG((radiusm-drad)/sizedistk1)
448	dista(j,1,iaer,idomain,gausind) =
449	& EXP(-dista(j,1,iaer,idomain,gausind) *
450	& dista(j,1,iaer,idomain,gausind) *
451	& 0.5e0/sizedistk2)/(radiusm-drad)
452	dista(j,1,iaer,idomain,gausind) =
453	& dista(j,1,iaer,idomain,gausind) /
454	& (sqrt(2e0pisizedistk2))
455
456	distb(j,1,iaer,idomain,gausind) =
457	& LOG((radiusm+drad)/sizedistk1)
458	distb(j,1,iaer,idomain,gausind) =
459	& EXP(-distb(j,1,iaer,idomain,gausind) *
460	& distb(j,1,iaer,idomain,gausind) *
461	& 0.5e0/sizedistk2)/(radiusm+drad)
462	distb(j,1,iaer,idomain,gausind) =
463	& distb(j,1,iaer,idomain,gausind) /
464	& (sqrt(2e0pisizedistk2))
465
466	normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) +
467	& weightgaus(gausind) *
468	& (
469	& distb(j,1,iaer,idomain,gausind) * pi *
470	& radGAUSb(gausind,iaer,idomain) *
471	& radGAUSb(gausind,iaer,idomain) +
472	& dista(j,1,iaer,idomain,gausind) * pi *
473	& radGAUSa(gausind,iaer,idomain) *
474	& radGAUSa(gausind,iaer,idomain)
475	& )
476	ENDDO
477	IF (normd(j,1,iaer,idomain).EQ.1e-30) THEN
478	WRITE(,)"normd:", normd(j,1,iaer,idomain)
479	WRITE(,)"Risk of division by 0 (aeroptproperties.F)"
480	WRITE(,)"Check the size of the interpolation grid."
481	CALL ABORT
482	ENDIF
483	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
484	c 2.3.3.vis Initialization
485	qsqrefVISgrid(j,1,:,iaer)=0.
486	qextVISgrid(j,1,:,iaer)=0.
487	qscatVISgrid(j,1,:,iaer)=0.
488	omegVISgrid(j,1,:,iaer)=0.
489	gVISgrid(j,1,:,iaer)=0.
490	qrefVISgrid(j,1,iaer)=0.
491	qscatrefVISgrid(j,1,iaer)=0.
492	omegrefVISgrid(j,1,iaer)=0.
493
494	DO gausind=1,ngau
495	DO m=1,nsun
496	c Convolution:
497	qextVISgrid(j,1,m,iaer) =
498	& qextVISgrid(j,1,m,iaer) +
499	& weightgaus(gausind) *
500	& (
501	& qsqrefVISb(m,gausind,iaer) *
502	& qrefVISb(gausind,iaer) *
503	& piradGAUSb(gausind,iaer,idomain)
504	& radGAUSb(gausind,iaer,idomain) *
505	& distb(j,1,iaer,idomain,gausind) +
506	& qsqrefVISa(m,gausind,iaer) *
507	& qrefVISa(gausind,iaer) *
508	& piradGAUSa(gausind,iaer,idomain)
509	& radGAUSa(gausind,iaer,idomain) *
510	& dista(j,1,iaer,idomain,gausind)
511	& )
512	qscatVISgrid(j,1,m,iaer) =
513	& qscatVISgrid(j,1,m,iaer) +
514	& weightgaus(gausind) *
515	& (
516	& omegVISb(m,gausind,iaer) *
517	& qsqrefVISb(m,gausind,iaer) *
518	& qrefVISb(gausind,iaer) *
519	& piradGAUSb(gausind,iaer,idomain)
520	& radGAUSb(gausind,iaer,idomain) *
521	& distb(j,1,iaer,idomain,gausind) +
522	& omegVISa(m,gausind,iaer) *
523	& qsqrefVISa(m,gausind,iaer) *
524	& qrefVISa(gausind,iaer) *
525	& piradGAUSa(gausind,iaer,idomain)
526	& radGAUSa(gausind,iaer,idomain) *
527	& dista(j,1,iaer,idomain,gausind)
528	& )
529	gVISgrid(j,1,m,iaer) =
530	& gVISgrid(j,1,m,iaer) +
531	& weightgaus(gausind) *
532	& (
533	& omegVISb(m,gausind,iaer) *
534	& qsqrefVISb(m,gausind,iaer) *
535	& qrefVISb(gausind,iaer) *
536	& gVISb(m,gausind,iaer) *
537	& piradGAUSb(gausind,iaer,idomain)
538	& radGAUSb(gausind,iaer,idomain) *
539	& distb(j,1,iaer,idomain,gausind) +
540	& omegVISa(m,gausind,iaer) *
541	& qsqrefVISa(m,gausind,iaer) *
542	& qrefVISa(gausind,iaer) *
543	& gVISa(m,gausind,iaer) *
544	& piradGAUSa(gausind,iaer,idomain)
545	& radGAUSa(gausind,iaer,idomain) *
546	& dista(j,1,iaer,idomain,gausind)
547	& )
548	ENDDO
549	qrefVISgrid(j,1,iaer) =
550	& qrefVISgrid(j,1,iaer) +
551	& weightgaus(gausind) *
552	& (
553	& qrefVISb(gausind,iaer) *
554	& piradGAUSb(gausind,iaer,idomain)
555	& radGAUSb(gausind,iaer,idomain) *
556	& distb(j,1,iaer,idomain,gausind) +
557	& qrefVISa(gausind,iaer) *
558	& piradGAUSa(gausind,iaer,idomain)
559	& radGAUSa(gausind,iaer,idomain) *
560	& dista(j,1,iaer,idomain,gausind)
561	& )
562	qscatrefVISgrid(j,1,iaer) =
563	& qscatrefVISgrid(j,1,iaer) +
564	& weightgaus(gausind) *
565	& (
566	& omegrefVISb(gausind,iaer) *
567	& qrefVISb(gausind,iaer) *
568	& piradGAUSb(gausind,iaer,idomain)
569	& radGAUSb(gausind,iaer,idomain) *
570	& distb(j,1,iaer,idomain,gausind) +
571	& omegrefVISa(gausind,iaer) *
572	& qrefVISa(gausind,iaer) *
573	& piradGAUSa(gausind,iaer,idomain)
574	& radGAUSa(gausind,iaer,idomain) *
575	& dista(j,1,iaer,idomain,gausind)
576	& )
577	ENDDO
578
579	qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) /
580	& normd(j,1,iaer,idomain)
581	qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /
582	& normd(j,1,iaer,idomain)
583	omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /
584	& qrefVISgrid(j,1,iaer)
585	DO m=1,nsun
586	qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /
587	& normd(j,1,iaer,idomain)
588	qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /
589	& normd(j,1,iaer,idomain)
590	gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) /
591	& qscatVISgrid(j,1,m,iaer) /
592	& normd(j,1,iaer,idomain)
593
594	qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /
595	& qrefVISgrid(j,1,iaer)
596	omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /
597	& qextVISgrid(j,1,m,iaer)
598	ENDDO
599	ELSE ! INFRARED DOMAIN ----------
600	c 2.3.3.ir Initialization
601	qsqrefIRgrid(j,1,:,iaer)=0.
602	qextIRgrid(j,1,:,iaer)=0.
603	qscatIRgrid(j,1,:,iaer)=0.
604	omegIRgrid(j,1,:,iaer)=0.
605	gIRgrid(j,1,:,iaer)=0.
606	qrefIRgrid(j,1,iaer)=0.
607	qscatrefIRgrid(j,1,iaer)=0.
608	omegrefIRgrid(j,1,iaer)=0.
609
610	DO gausind=1,ngau
611	DO m=1,nir
612	c Convolution:
613	qextIRgrid(j,1,m,iaer) =
614	& qextIRgrid(j,1,m,iaer) +
615	& weightgaus(gausind) *
616	& (
617	& qsqrefIRb(m,gausind,iaer) *
618	& qrefVISb(gausind,iaer) *
619	& piradGAUSb(gausind,iaer,idomain)
620	& radGAUSb(gausind,iaer,idomain) *
621	& distb(j,1,iaer,idomain,gausind) +
622	& qsqrefIRa(m,gausind,iaer) *
623	& qrefVISa(gausind,iaer) *
624	& piradGAUSa(gausind,iaer,idomain)
625	& radGAUSa(gausind,iaer,idomain) *
626	& dista(j,1,iaer,idomain,gausind)
627	& )
628	qscatIRgrid(j,1,m,iaer) =
629	& qscatIRgrid(j,1,m,iaer) +
630	& weightgaus(gausind) *
631	& (
632	& omegIRb(m,gausind,iaer) *
633	& qsqrefIRb(m,gausind,iaer) *
634	& qrefVISb(gausind,iaer) *
635	& piradGAUSb(gausind,iaer,idomain)
636	& radGAUSb(gausind,iaer,idomain) *
637	& distb(j,1,iaer,idomain,gausind) +
638	& omegIRa(m,gausind,iaer) *
639	& qsqrefIRa(m,gausind,iaer) *
640	& qrefVISa(gausind,iaer) *
641	& piradGAUSa(gausind,iaer,idomain)
642	& radGAUSa(gausind,iaer,idomain) *
643	& dista(j,1,iaer,idomain,gausind)
644	& )
645	gIRgrid(j,1,m,iaer) =
646	& gIRgrid(j,1,m,iaer) +
647	& weightgaus(gausind) *
648	& (
649	& omegIRb(m,gausind,iaer) *
650	& qsqrefIRb(m,gausind,iaer) *
651	& qrefVISb(gausind,iaer) *
652	& gIRb(m,gausind,iaer) *
653	& piradGAUSb(gausind,iaer,idomain)
654	& radGAUSb(gausind,iaer,idomain) *
655	& distb(j,1,iaer,idomain,gausind) +
656	& omegIRa(m,gausind,iaer) *
657	& qsqrefIRa(m,gausind,iaer) *
658	& qrefVISa(gausind,iaer) *
659	& gIRa(m,gausind,iaer) *
660	& piradGAUSa(gausind,iaer,idomain)
661	& radGAUSa(gausind,iaer,idomain) *
662	& dista(j,1,iaer,idomain,gausind)
663	& )
664	ENDDO
665	qrefIRgrid(j,1,iaer) =
666	& qrefIRgrid(j,1,iaer) +
667	& weightgaus(gausind) *
668	& (
669	& qrefIRb(gausind,iaer) *
670	& piradGAUSb(gausind,iaer,idomain)
671	& radGAUSb(gausind,iaer,idomain) *
672	& distb(j,1,iaer,idomain,gausind) +
673	& qrefIRa(gausind,iaer) *
674	& piradGAUSa(gausind,iaer,idomain)
675	& radGAUSa(gausind,iaer,idomain) *
676	& dista(j,1,iaer,idomain,gausind)
677	& )
678	qscatrefIRgrid(j,1,iaer) =
679	& qscatrefIRgrid(j,1,iaer) +
680	& weightgaus(gausind) *
681	& (
682	& omegrefIRb(gausind,iaer) *
683	& qrefIRb(gausind,iaer) *
684	& piradGAUSb(gausind,iaer,idomain)
685	& radGAUSb(gausind,iaer,idomain) *
686	& distb(j,1,iaer,idomain,gausind) +
687	& omegrefIRa(gausind,iaer) *
688	& qrefIRa(gausind,iaer) *
689	& piradGAUSa(gausind,iaer,idomain)
690	& radGAUSa(gausind,iaer,idomain) *
691	& dista(j,1,iaer,idomain,gausind)
692	& )
693	ENDDO
694
695	qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) /
696	& normd(j,1,iaer,idomain)
697	qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /
698	& normd(j,1,iaer,idomain)
699	omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /
700	& qrefIRgrid(j,1,iaer)
701	DO m=1,nir
702	qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /
703	& normd(j,1,iaer,idomain)
704	qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /
705	& normd(j,1,iaer,idomain)
706	gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) /
707	& qscatIRgrid(j,1,m,iaer) /
708	& normd(j,1,iaer,idomain)
709
710	qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /
711	& qrefVISgrid(j,1,iaer)
712	omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /
713	& qextIRgrid(j,1,m,iaer)
714	ENDDO
715	ENDIF ! --------------------------
716	checkgrid(j,1,iaer,idomain) = .true.
717	ENDIF !checkgrid
718	ENDDO !grid_i
719	c 2.4 Linear interpolation
720	k1 = (1-kx)
721	k2 = kx
722	IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
723	DO m=1,nsun
724	QVISsQREF3d(ig,lg,m,iaer) =
725	& k1*qsqrefVISgrid(grid_i,1,m,iaer) +
726	& k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
727	omegaVIS3d(ig,lg,m,iaer) =
728	& k1*omegVISgrid(grid_i,1,m,iaer) +
729	& k2*omegVISgrid(grid_i+1,1,m,iaer)
730	gVIS3d(ig,lg,m,iaer) =
731	& k1*gVISgrid(grid_i,1,m,iaer) +
732	& k2*gVISgrid(grid_i+1,1,m,iaer)
733	ENDDO !nsun
734	QREFvis3d(ig,lg,iaer) =
735	& k1*qrefVISgrid(grid_i,1,iaer) +
736	& k2*qrefVISgrid(grid_i+1,1,iaer)
737	omegaREFvis3d(ig,lg,iaer) =
738	& k1*omegrefVISgrid(grid_i,1,iaer) +
739	& k2*omegrefVISgrid(grid_i+1,1,iaer)
740	ELSE ! INFRARED -----------------------
741	DO m=1,nir
742	QIRsQREF3d(ig,lg,m,iaer) =
743	& k1*qsqrefIRgrid(grid_i,1,m,iaer) +
744	& k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
745	omegaIR3d(ig,lg,m,iaer) =
746	& k1*omegIRgrid(grid_i,1,m,iaer) +
747	& k2*omegIRgrid(grid_i+1,1,m,iaer)
748	gIR3d(ig,lg,m,iaer) =
749	& k1*gIRgrid(grid_i,1,m,iaer) +
750	& k2*gIRgrid(grid_i+1,1,m,iaer)
751	ENDDO !nir
752	QREFir3d(ig,lg,iaer) =
753	& k1*qrefIRgrid(grid_i,1,iaer) +
754	& k2*qrefIRgrid(grid_i+1,1,iaer)
755	omegaREFir3d(ig,lg,iaer) =
756	& k1*omegrefIRgrid(grid_i,1,iaer) +
757	& k2*omegrefIRgrid(grid_i+1,1,iaer)
758	ENDIF ! --------------------------------
759	ENDDO !nlayermx
760	ENDDO !ngridmx
761	c==================================================================
762	ELSE ! VARYING NUEFF
763	c==================================================================
764	c 3. Compute the scattering parameters in each grid box
765	c using bilinear interpolation over the grain sizes
766	c and the effective variances;
767	c -----------------------------------------------------
768
769	DO lg = 1,nlayer
770	DO ig = 1,ngrid
771	c 3.1 Effective variance index and ky calculation
772	var_tmp=log(nueffrad(ig,lg,iaer))
773	grid_j=floor( (nuefftabsize-1)/(nuefftabmax-nuefftabmin)*
774	& (var_tmp-nuefftabmin)+1. )
775	IF (grid_j.GE.nuefftabsize) THEN
776	c WRITE(,) 'Warning: effective variance '
777	c WRITE(,) 'is too large to be used by the '
778	c WRITE(,) 'radiative transfer; please extend the '
779	c WRITE(,) 'interpolation grid to larger values.'
780	grid_j=nuefftabsize-1
781	ky = 1.
782	ELSEIF (grid_j.LT.1) THEN
783	c WRITE(,) 'Warning: effective variance '
784	c WRITE(,) 'is too small to be used by the '
785	c WRITE(,) 'radiative transfer; please extend the '
786	c WRITE(,) 'interpolation grid to smaller values.'
787	grid_j=1
788	ky = 0.
789	ELSE
790	ky = ( nueffrad(ig,lg,iaer)-
791	& nuefftab(grid_j,iaer,idomain) ) /
792	& ( nuefftab(grid_j+1,iaer,idomain)-
793	& nuefftab(grid_j,iaer,idomain) )
794	ENDIF
795	c 3.2 Effective radius index and kx calculation
796	var_tmp=reffrad(ig,lg,iaer)/refftabmin(iaer,idomain)
797	var_tmp=log(var_tmp)*3.
798	var_tmp=var_tmp/logvratgrid(iaer,idomain)+1.
799	grid_i=floor(var_tmp)
800	IF (grid_i.GE.refftabsize) THEN
801	c WRITE(,) 'Warning: particle size in grid box #'
802	c WRITE(,) ig,' is too large to be used by the '
803	c WRITE(,) 'radiative transfer; please extend the '
804	c WRITE(,) 'interpolation grid to larger grain sizes.'
805	grid_i=refftabsize-1
806	kx = 1.
807	ELSEIF (grid_i.LT.1) THEN
808	c WRITE(,) 'Warning: particle size in grid box #'
809	c WRITE(,) ig,' is too small to be used by the '
810	c WRITE(,) 'radiative transfer; please extend the '
811	c WRITE(,) 'interpolation grid to smaller grain sizes.'
812	grid_i=1
813	kx = 0.
814	ELSE
815	kx = ( reffrad(ig,lg,iaer)-
816	& refftab(grid_i,iaer,idomain) ) /
817	& ( refftab(grid_i+1,iaer,idomain)-
818	& refftab(grid_i,iaer,idomain) )
819	ENDIF
820	c 3.3 Integration
821	DO j=grid_i,grid_i+1
822	DO k=grid_j,grid_j+1
823	c 3.3.1 Check if the calculation has been done
824	IF (.NOT.checkgrid(j,k,iaer,idomain)) THEN
825
826	c 3.3.2 Log-normal dist., r_g and sigma_g are defined
827	c in [hansen_1974], "Light scattering in planetary
828	c atmospheres", Space Science Reviews 16 527-610.
829	c Here, sizedistk1=r_g and sizedistk2=sigma_g^2
830	sizedistk2 = log(1.+nuefftab(k,iaer,idomain))
831	sizedistk1 = exp(2.5*sizedistk2)
832	sizedistk1 = refftab(j,iaer,idomain) / sizedistk1
833
834	normd(j,k,iaer,idomain) = 1e-30
835	DO gausind=1,ngau
836	drad=radiusr*radgaus(gausind)
837
838	dista(j,k,iaer,idomain,gausind) =
839	& LOG((radiusm-drad)/sizedistk1)
840	dista(j,k,iaer,idomain,gausind) =
841	& EXP(-dista(j,k,iaer,idomain,gausind) *
842	& dista(j,k,iaer,idomain,gausind) *
843	& 0.5e0/sizedistk2)/(radiusm-drad)
844	dista(j,k,iaer,idomain,gausind) =
845	& dista(j,k,iaer,idomain,gausind) /
846	& (sqrt(2e0pisizedistk2))
847
848	distb(j,k,iaer,idomain,gausind) =
849	& LOG((radiusm+drad)/sizedistk1)
850	distb(j,k,iaer,idomain,gausind) =
851	& EXP(-distb(j,k,iaer,idomain,gausind) *
852	& distb(j,k,iaer,idomain,gausind) *
853	& 0.5e0/sizedistk2)/(radiusm+drad)
854	distb(j,k,iaer,idomain,gausind) =
855	& distb(j,k,iaer,idomain,gausind) /
856	& (sqrt(2e0pisizedistk2))
857
858	normd(j,k,iaer,idomain)=normd(j,k,iaer,idomain) +
859	& weightgaus(gausind) *
860	& (
861	& distb(j,k,iaer,idomain,gausind) * pi *
862	& radGAUSb(gausind,iaer,idomain) *
863	& radGAUSb(gausind,iaer,idomain) +
864	& dista(j,k,iaer,idomain,gausind) * pi *
865	& radGAUSa(gausind,iaer,idomain) *
866	& radGAUSa(gausind,iaer,idomain)
867	& )
868	ENDDO
869	IF (normd(j,k,iaer,idomain).EQ.1e-30) THEN
870	WRITE(,)"normd:", normd(j,k,iaer,idomain)
871	WRITE(,)"Risk of division by 0 (aeroptproperties.F)"
872	WRITE(,)"Check the size of the interpolation grid."
873	CALL ABORT
874	ENDIF
875	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
876	c 2.3.3.vis Initialization
877	qsqrefVISgrid(j,k,:,iaer)=0.
878	qextVISgrid(j,k,:,iaer)=0.
879	qscatVISgrid(j,k,:,iaer)=0.
880	omegVISgrid(j,k,:,iaer)=0.
881	gVISgrid(j,k,:,iaer)=0.
882	qrefVISgrid(j,k,iaer)=0.
883	qscatrefVISgrid(j,k,iaer)=0.
884	omegrefVISgrid(j,k,iaer)=0.
885
886	DO gausind=1,ngau
887	DO m=1,nsun
888	c Convolution:
889	qextVISgrid(j,k,m,iaer) =
890	& qextVISgrid(j,k,m,iaer) +
891	& weightgaus(gausind) *
892	& (
893	& qsqrefVISb(m,gausind,iaer) *
894	& qrefVISb(gausind,iaer) *
895	& piradGAUSb(gausind,iaer,idomain)
896	& radGAUSb(gausind,iaer,idomain) *
897	& distb(j,k,iaer,idomain,gausind) +
898	& qsqrefVISa(m,gausind,iaer) *
899	& qrefVISa(gausind,iaer) *
900	& piradGAUSa(gausind,iaer,idomain)
901	& radGAUSa(gausind,iaer,idomain) *
902	& dista(j,k,iaer,idomain,gausind)
903	& )
904	qscatVISgrid(j,k,m,iaer) =
905	& qscatVISgrid(j,k,m,iaer) +
906	& weightgaus(gausind) *
907	& (
908	& omegVISb(m,gausind,iaer) *
909	& qsqrefVISb(m,gausind,iaer) *
910	& qrefVISb(gausind,iaer) *
911	& piradGAUSb(gausind,iaer,idomain)
912	& radGAUSb(gausind,iaer,idomain) *
913	& distb(j,k,iaer,idomain,gausind) +
914	& omegVISa(m,gausind,iaer) *
915	& qsqrefVISa(m,gausind,iaer) *
916	& qrefVISa(gausind,iaer) *
917	& piradGAUSa(gausind,iaer,idomain)
918	& radGAUSa(gausind,iaer,idomain) *
919	& dista(j,k,iaer,idomain,gausind)
920	& )
921	gVISgrid(j,k,m,iaer) =
922	& gVISgrid(j,k,m,iaer) +
923	& weightgaus(gausind) *
924	& (
925	& omegVISb(m,gausind,iaer) *
926	& qsqrefVISb(m,gausind,iaer) *
927	& qrefVISb(gausind,iaer) *
928	& gVISb(m,gausind,iaer) *
929	& piradGAUSb(gausind,iaer,idomain)
930	& radGAUSb(gausind,iaer,idomain) *
931	& distb(j,k,iaer,idomain,gausind) +
932	& omegVISa(m,gausind,iaer) *
933	& qsqrefVISa(m,gausind,iaer) *
934	& qrefVISa(gausind,iaer) *
935	& gVISa(m,gausind,iaer) *
936	& piradGAUSa(gausind,iaer,idomain)
937	& radGAUSa(gausind,iaer,idomain) *
938	& dista(j,k,iaer,idomain,gausind)
939	& )
940	ENDDO
941	qrefVISgrid(j,k,iaer) =
942	& qrefVISgrid(j,k,iaer) +
943	& weightgaus(gausind) *
944	& (
945	& qrefVISb(gausind,iaer) *
946	& piradGAUSb(gausind,iaer,idomain)
947	& radGAUSb(gausind,iaer,idomain) *
948	& distb(j,k,iaer,idomain,gausind) +
949	& qrefVISa(gausind,iaer) *
950	& piradGAUSa(gausind,iaer,idomain)
951	& radGAUSa(gausind,iaer,idomain) *
952	& dista(j,k,iaer,idomain,gausind)
953	& )
954	qscatrefVISgrid(j,k,iaer) =
955	& qscatrefVISgrid(j,k,iaer) +
956	& weightgaus(gausind) *
957	& (
958	& omegrefVISb(gausind,iaer) *
959	& qrefVISb(gausind,iaer) *
960	& piradGAUSb(gausind,iaer,idomain)
961	& radGAUSb(gausind,iaer,idomain) *
962	& distb(j,k,iaer,idomain,gausind) +
963	& omegrefVISa(gausind,iaer) *
964	& qrefVISa(gausind,iaer) *
965	& piradGAUSa(gausind,iaer,idomain)
966	& radGAUSa(gausind,iaer,idomain) *
967	& dista(j,k,iaer,idomain,gausind)
968	& )
969	ENDDO
970	qrefVISgrid(j,k,iaer)=qrefVISgrid(j,k,iaer) /
971	& normd(j,k,iaer,idomain)
972	qscatrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) /
973	& normd(j,k,iaer,idomain)
974	omegrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) /
975	& qrefVISgrid(j,k,iaer)
976	DO m=1,nsun
977	qextVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) /
978	& normd(j,k,iaer,idomain)
979	qscatVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) /
980	& normd(j,k,iaer,idomain)
981	gVISgrid(j,k,m,iaer)=gVISgrid(j,k,m,iaer) /
982	& qscatVISgrid(j,k,m,iaer) /
983	& normd(j,k,iaer,idomain)
984
985	qsqrefVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) /
986	& qrefVISgrid(j,k,iaer)
987	omegVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) /
988	& qextVISgrid(j,k,m,iaer)
989	ENDDO
990	ELSE ! INFRARED DOMAIN ----------
991	c 2.3.3.ir Initialization
992	qsqrefIRgrid(j,k,:,iaer)=0.
993	qextIRgrid(j,k,:,iaer)=0.
994	qscatIRgrid(j,k,:,iaer)=0.
995	omegIRgrid(j,k,:,iaer)=0.
996	gIRgrid(j,k,:,iaer)=0.
997	qrefIRgrid(j,k,iaer)=0.
998	qscatrefIRgrid(j,k,iaer)=0.
999	omegrefIRgrid(j,k,iaer)=0.
1000
1001	DO gausind=1,ngau
1002	DO m=1,nir
1003	c Convolution:
1004	qextIRgrid(j,k,m,iaer) =
1005	& qextIRgrid(j,k,m,iaer) +
1006	& weightgaus(gausind) *
1007	& (
1008	& qsqrefIRb(m,gausind,iaer) *
1009	& qrefVISb(gausind,iaer) *
1010	& piradGAUSb(gausind,iaer,idomain)
1011	& radGAUSb(gausind,iaer,idomain) *
1012	& distb(j,k,iaer,idomain,gausind) +
1013	& qsqrefIRa(m,gausind,iaer) *
1014	& qrefVISa(gausind,iaer) *
1015	& piradGAUSa(gausind,iaer,idomain)
1016	& radGAUSa(gausind,iaer,idomain) *
1017	& dista(j,k,iaer,idomain,gausind)
1018	& )
1019	qscatIRgrid(j,k,m,iaer) =
1020	& qscatIRgrid(j,k,m,iaer) +
1021	& weightgaus(gausind) *
1022	& (
1023	& omegIRb(m,gausind,iaer) *
1024	& qsqrefIRb(m,gausind,iaer) *
1025	& qrefVISb(gausind,iaer) *
1026	& piradGAUSb(gausind,iaer,idomain)
1027	& radGAUSb(gausind,iaer,idomain) *
1028	& distb(j,k,iaer,idomain,gausind) +
1029	& omegIRa(m,gausind,iaer) *
1030	& qsqrefIRa(m,gausind,iaer) *
1031	& qrefVISa(gausind,iaer) *
1032	& piradGAUSa(gausind,iaer,idomain)
1033	& radGAUSa(gausind,iaer,idomain) *
1034	& dista(j,k,iaer,idomain,gausind)
1035	& )
1036	gIRgrid(j,k,m,iaer) =
1037	& gIRgrid(j,k,m,iaer) +
1038	& weightgaus(gausind) *
1039	& (
1040	& omegIRb(m,gausind,iaer) *
1041	& qsqrefIRb(m,gausind,iaer) *
1042	& qrefVISb(gausind,iaer) *
1043	& gIRb(m,gausind,iaer) *
1044	& piradGAUSb(gausind,iaer,idomain)
1045	& radGAUSb(gausind,iaer,idomain) *
1046	& distb(j,k,iaer,idomain,gausind) +
1047	& omegIRa(m,gausind,iaer) *
1048	& qsqrefIRa(m,gausind,iaer) *
1049	& qrefVISa(gausind,iaer) *
1050	& gIRa(m,gausind,iaer) *
1051	& piradGAUSa(gausind,iaer,idomain)
1052	& radGAUSa(gausind,iaer,idomain) *
1053	& dista(j,k,iaer,idomain,gausind)
1054	& )
1055	ENDDO
1056	qrefIRgrid(j,k,iaer) =
1057	& qrefIRgrid(j,k,iaer) +
1058	& weightgaus(gausind) *
1059	& (
1060	& qrefIRb(gausind,iaer) *
1061	& piradGAUSb(gausind,iaer,idomain)
1062	& radGAUSb(gausind,iaer,idomain) *
1063	& distb(j,k,iaer,idomain,gausind) +
1064	& qrefIRa(gausind,iaer) *
1065	& piradGAUSa(gausind,iaer,idomain)
1066	& radGAUSa(gausind,iaer,idomain) *
1067	& dista(j,k,iaer,idomain,gausind)
1068	& )
1069	qscatrefIRgrid(j,k,iaer) =
1070	& qscatrefIRgrid(j,k,iaer) +
1071	& weightgaus(gausind) *
1072	& (
1073	& omegrefIRb(gausind,iaer) *
1074	& qrefIRb(gausind,iaer) *
1075	& piradGAUSb(gausind,iaer,idomain)
1076	& radGAUSb(gausind,iaer,idomain) *
1077	& distb(j,k,iaer,idomain,gausind) +
1078	& omegrefIRa(gausind,iaer) *
1079	& qrefIRa(gausind,iaer) *
1080	& piradGAUSa(gausind,iaer,idomain)
1081	& radGAUSa(gausind,iaer,idomain) *
1082	& dista(j,k,iaer,idomain,gausind)
1083	& )
1084	ENDDO
1085	qrefIRgrid(j,k,iaer)=qrefIRgrid(j,k,iaer) /
1086	& normd(j,k,iaer,idomain)
1087	qscatrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) /
1088	& normd(j,k,iaer,idomain)
1089	omegrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) /
1090	& qrefIRgrid(j,k,iaer)
1091	DO m=1,nir
1092	qextIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) /
1093	& normd(j,k,iaer,idomain)
1094	qscatIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) /
1095	& normd(j,k,iaer,idomain)
1096	gIRgrid(j,k,m,iaer)=gIRgrid(j,k,m,iaer) /
1097	& qscatIRgrid(j,k,m,iaer) /
1098	& normd(j,k,iaer,idomain)
1099
1100	qsqrefIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) /
1101	& qrefVISgrid(j,k,iaer)
1102	omegIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) /
1103	& qextIRgrid(j,k,m,iaer)
1104	ENDDO
1105	ENDIF ! --------------------------
1106	checkgrid(j,k,iaer,idomain) = .true.
1107	ENDIF !checkgrid
1108	ENDDO !grid_j
1109	ENDDO !grid_i
1110	c 3.4 Bilinear interpolation
1111	k1 = (1-kx)*(1-ky)
1112	k2 = kx*(1-ky)
1113	k3 = kx*ky
1114	k4 = (1-kx)*ky
1115	IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
1116	DO m=1,nsun
1117	QVISsQREF3d(ig,lg,m,iaer) =
1118	& k1*qsqrefVISgrid(grid_i,grid_j,m,iaer) +
1119	& k2*qsqrefVISgrid(grid_i+1,grid_j,m,iaer) +
1120	& k3*qsqrefVISgrid(grid_i+1,grid_j+1,m,iaer) +
1121	& k4*qsqrefVISgrid(grid_i,grid_j+1,m,iaer)
1122	omegaVIS3d(ig,lg,m,iaer) =
1123	& k1*omegVISgrid(grid_i,grid_j,m,iaer) +
1124	& k2*omegVISgrid(grid_i+1,grid_j,m,iaer) +
1125	& k3*omegVISgrid(grid_i+1,grid_j+1,m,iaer) +
1126	& k4*omegVISgrid(grid_i,grid_j+1,m,iaer)
1127	gVIS3d(ig,lg,m,iaer) =
1128	& k1*gVISgrid(grid_i,grid_j,m,iaer) +
1129	& k2*gVISgrid(grid_i+1,grid_j,m,iaer) +
1130	& k3*gVISgrid(grid_i+1,grid_j+1,m,iaer) +
1131	& k4*gVISgrid(grid_i,grid_j+1,m,iaer)
1132	ENDDO !nsun
1133	QREFvis3d(ig,lg,iaer) =
1134	& k1*qrefVISgrid(grid_i,grid_j,iaer) +
1135	& k2*qrefVISgrid(grid_i+1,grid_j,iaer) +
1136	& k3*qrefVISgrid(grid_i+1,grid_j+1,iaer) +
1137	& k4*qrefVISgrid(grid_i,grid_j+1,iaer)
1138	omegaREFvis3d(ig,lg,iaer) =
1139	& k1*omegrefVISgrid(grid_i,grid_j,iaer) +
1140	& k2*omegrefVISgrid(grid_i+1,grid_j,iaer) +
1141	& k3*omegrefVISgrid(grid_i+1,grid_j+1,iaer) +
1142	& k4*omegrefVISgrid(grid_i,grid_j+1,iaer)
1143	ELSE ! INFRARED -----------------------
1144	DO m=1,nir
1145	QIRsQREF3d(ig,lg,m,iaer) =
1146	& k1*qsqrefIRgrid(grid_i,grid_j,m,iaer) +
1147	& k2*qsqrefIRgrid(grid_i+1,grid_j,m,iaer) +
1148	& k3*qsqrefIRgrid(grid_i+1,grid_j+1,m,iaer) +
1149	& k4*qsqrefIRgrid(grid_i,grid_j+1,m,iaer)
1150	omegaIR3d(ig,lg,m,iaer) =
1151	& k1*omegIRgrid(grid_i,grid_j,m,iaer) +
1152	& k2*omegIRgrid(grid_i+1,grid_j,m,iaer) +
1153	& k3*omegIRgrid(grid_i+1,grid_j+1,m,iaer) +
1154	& k4*omegIRgrid(grid_i,grid_j+1,m,iaer)
1155	gIR3d(ig,lg,m,iaer) =
1156	& k1*gIRgrid(grid_i,grid_j,m,iaer) +
1157	& k2*gIRgrid(grid_i+1,grid_j,m,iaer) +
1158	& k3*gIRgrid(grid_i+1,grid_j+1,m,iaer) +
1159	& k4*gIRgrid(grid_i,grid_j+1,m,iaer)
1160	ENDDO !nir
1161	QREFir3d(ig,lg,iaer) =
1162	& k1*qrefIRgrid(grid_i,grid_j,iaer) +
1163	& k2*qrefIRgrid(grid_i+1,grid_j,iaer) +
1164	& k3*qrefIRgrid(grid_i+1,grid_j+1,iaer) +
1165	& k4*qrefIRgrid(grid_i,grid_j+1,iaer)
1166	omegaREFir3d(ig,lg,iaer) =
1167	& k1*omegrefIRgrid(grid_i,grid_j,iaer) +
1168	& k2*omegrefIRgrid(grid_i+1,grid_j,iaer) +
1169	& k3*omegrefIRgrid(grid_i+1,grid_j+1,iaer) +
1170	& k4*omegrefIRgrid(grid_i,grid_j+1,iaer)
1171	ENDIF ! --------------------------------
1172	ENDDO !nlayermx
1173	ENDDO !ngridmx
1174
1175	ENDIF ! varyingnueff
1176	c==================================================================
1177	ENDDO ! idomain
1178
1179	ENDIF ! nsize = 1
1180
1181	ENDDO ! iaer (loop on aerosol kind)
1182
1183	c=====Radiative properties - TESTS=================================
1184	IF (out_qwg) THEN
1185	IF (ngrid.NE.1) THEN
1186	DO out_nchannel = 1, 2
1187	c -------------------------------------------------------------
1188	DO lg = 1, nlayer
1189	DO ig = 1, ngrid
1190	out_qext(ig,lg) =
1191	& QVISsQREF3d(ig,lg,out_nchannel,out_iaer)*
1192	& QREFvis3d(ig,lg,out_iaer)
1193	out_omeg(ig,lg) =
1194	& omegaVIS3d(ig,lg,out_nchannel,out_iaer)
1195	out_g(ig,lg) = gVIS3d(ig,lg,out_nchannel,out_iaer)
1196	ENDDO ! ig
1197	ENDDO ! lg
1198	write(out_str(1:1),'(i1.1)') out_nchannel
1199	call WRITEDIAGFI(ngrid,'qextvis'//out_str,"Ext.efficiency","",
1200	& 3,out_qext)
1201	call WRITEDIAGFI(ngrid,'omegvis'//out_str,"Sing.Scat.Alb.","",
1202	& 3,out_omeg)
1203	call WRITEDIAGFI(ngrid,'gvis'//out_str,"Asym.Factor","",
1204	& 3,out_g)
1205	c -------------------------------------------------------------
1206	ENDDO ! out_nchannel
1207	DO out_nchannel = 2, 4
1208	c -------------------------------------------------------------
1209	DO lg = 1, nlayer
1210	DO ig = 1, ngrid
1211	out_qext(ig,lg) =
1212	& QIRsQREF3d(ig,lg,out_nchannel,out_iaer)*
1213	& QREFir3d(ig,lg,out_iaer)
1214	out_omeg(ig,lg) =
1215	& omegaIR3d(ig,lg,out_nchannel,out_iaer)
1216	out_g(ig,lg) = gIR3d(ig,lg,out_nchannel,out_iaer)
1217	ENDDO ! ig
1218	ENDDO ! lg
1219	write(out_str(1:1),'(i1.1)') out_nchannel
1220	call WRITEDIAGFI(ngrid,'qextir'//out_str,"Ext.efficiency","",
1221	& 3,out_qext)
1222	call WRITEDIAGFI(ngrid,'omegir'//out_str,"Sing.Scat.Alb.","",
1223	& 3,out_omeg)
1224	call WRITEDIAGFI(ngrid,'gir'//out_str,"Asym.Factor","",
1225	& 3,out_g)
1226	c -------------------------------------------------------------
1227	ENDDO ! out_nchannel
1228	call WRITEDIAGFI(ngrid,"omegvisref","Sing.Scat.Alb.","",
1229	& 3,omegaREFvis3d(1,1,out_iaer))
1230	call WRITEDIAGFI(ngrid,"omegirref","Sing.Scat.Alb.","",
1231	& 3,omegaREFir3d(1,1,out_iaer))
1232	ENDIF ! ngrid.EQ.1
1233	ENDIF ! out_qwg
1234	c==================================================================
1235
1236	RETURN
1237	END

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

Download in other formats:

Original Format