Context Navigation

aeroptproperties.F @ 89

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

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

Download in other formats:

Original Format