Context Navigation

aeroptproperties.F @ 3325

Last change on this file since 3325 was 2932, checked in by romain.vande, 20 months ago

Mars PCM:
Add a new routine to write output called write_output.
It needs to be imported (for example : use write_output_mod, only: write_output)
Then, it requires only 4 arguments : the name of the variable, its title, its units and the variable itself.
It detects the dimension of the variable and decide to output either in diagfi or diagsoil.
It is also compatible with XIOS (xml file needs to be adapted)
Writediagfi and writediagsoil routines are still available in case.
RV

File size: 58.9 KB

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

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

Download in other formats:

Original Format