Context Navigation

aeroptproperties.F @ 1266

Last change on this file since 1266 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

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

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

Download in other formats:

Original Format