Context Navigation

aeroptproperties.F @ 1047

Last change on this file since 1047 was 1047, checked in by emillour, 11 years ago

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

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

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

Download in other formats:

Original Format