Context Navigation

aeroptproperties.F @ 2541

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

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

Download in other formats:

Original Format