Context Navigation

aeroptproperties.F @ 1292

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

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

Download in other formats:

Original Format