Context Navigation

aeroptproperties.F @ 784

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

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

Download in other formats:

Original Format