Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

aeroptproperties.F90

Last change on this file was 3966, checked in by debatzbr, 5 weeks ago
Pluto PCM: Add sanity check in the radiatif transfer (usefull for debug mode). BBT
File size: 41.8 KB

Line
1	module aeroptproperties_mod
2
3	implicit none
4
5	contains
6
7	SUBROUTINE aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
8	QVISsQREF3d,omegaVIS3d,gVIS3d, &
9	QIRsQREF3d,omegaIR3d,gIR3d, &
10	QREFvis3d,QREFir3d)!, &
11	! omegaREFvis3d,omegaREFir3d)
12
13	use radinc_h, only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
14	use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir
15	use radcommon_h, only: qrefvis,qrefir,omegarefir !,omegarefvis
16	use radcommon_h, only: radiustab,nsize
17
18	implicit none
19
20	! =============================================================
21	! Aerosol Optical Properties
22	!
23	! Description:
24	! Compute the scattering parameters in each grid
25	! box, depending on aerosol grain sizes. Log-normal size
26	! distribution and Gauss-Legendre integration are used.
27
28	! Parameters:
29	! Don't forget to set the value of varyingnueff below; If
30	! the effective variance of the distribution for the given
31	! aerosol is considered homogeneous in the atmosphere, please
32	! set varyingnueff(iaer) to .false. Resulting computational
33	! time will be much better.
34
35	! Authors: J.-B. Madeleine, F. Forget, F. Montmessin
36	! Slightly modified and converted to F90 by R. Wordsworth (2009)
37	! Varying nueff section removed by R. Wordsworth for simplicity
38	! ==============================================================
39
40	! INPUTS
41	! ------
42	INTEGER :: ngrid,nlayer
43	! Aerosol effective radius used for radiative transfer (meter)
44	REAL :: reffrad(ngrid,nlayer,naerkind)
45	! Aerosol effective variance used for radiative transfer (n.u.)
46	REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind)
47
48	! OUTPUTS
49	! -------
50	REAL,INTENT(OUT) :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
51	REAL,INTENT(OUT) :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
52	REAL,INTENT(OUT) :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
53
54	REAL,INTENT(OUT) :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
55	REAL,INTENT(OUT) :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
56	REAL,INTENT(OUT) :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
57
58	REAL,INTENT(OUT) :: QREFvis3d(ngrid,nlayer,naerkind)
59	REAL,INTENT(OUT) :: QREFir3d(ngrid,nlayer,naerkind)
60
61	! REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
62	! REAL :: omegaREFir3d(ngrid,nlayer,naerkind)
63
64	! =============================================================
65	! LOGICAL, PARAMETER :: varyingnueff(naerkind) = .false. ! not used!
66	! =============================================================
67
68	! Min. and max radius of the interpolation grid (in METERS)
69	REAL, PARAMETER :: refftabmin = 4e-9 !2e-8
70	! REAL, PARAMETER :: refftabmax = 35e-6
71	REAL, PARAMETER :: refftabmax = 1e-3
72	! Log of the min and max variance of the interpolation grid
73	REAL, PARAMETER :: nuefftabmin = -4.6
74	REAL, PARAMETER :: nuefftabmax = 0.
75	! Number of effective radius of the interpolation grid
76	INTEGER, PARAMETER :: refftabsize = 200
77	! Number of effective variances of the interpolation grid
78	! INTEGER, PARAMETER :: nuefftabsize = 100
79	INTEGER, PARAMETER :: nuefftabsize = 1
80	! Interpolation grid indices (reff,nueff)
81	INTEGER :: grid_i,grid_j
82	! Intermediate variable
83	REAL :: var_tmp,var3d_tmp(ngrid,nlayer)
84	! Bilinear interpolation factors
85	REAL :: kx,ky,k1,k2,k3,k4
86	! Size distribution parameters
87	REAL :: sizedistk1,sizedistk2
88	! Pi!
89	REAL,SAVE :: pi
90	!$OMP THREADPRIVATE(pi)
91	! Variables used by the Gauss-Legendre integration:
92	INTEGER radius_id,gausind
93	REAL kint
94	REAL drad
95	INTEGER, PARAMETER :: ngau = 10
96	REAL weightgaus(ngau),radgaus(ngau)
97	SAVE weightgaus,radgaus
98	! DATA weightgaus/.2955242247,.2692667193,.2190863625,.1494513491,.0666713443/
99	! DATA radgaus/.1488743389,.4333953941,.6794095682,.8650633666,.9739065285/
100	DATA radgaus/0.07652652113350,0.22778585114165, &
101	0.37370608871528,0.51086700195146, &
102	0.63605368072468,0.74633190646476, &
103	0.83911697181213,0.91223442826796, &
104	0.96397192726078,0.99312859919241/
105
106	DATA weightgaus/0.15275338723120,0.14917298659407, &
107	0.14209610937519,0.13168863843930, &
108	0.11819453196154,0.10193011980823, &
109	0.08327674160932,0.06267204829828, &
110	0.04060142982019,0.01761400714091/
111	!$OMP THREADPRIVATE(radgaus,weightgaus)
112	! Indices
113	INTEGER :: i,j,k,l,m,iaer,idomain
114	INTEGER :: ig,lg,chg
115
116
117	! Local saved variables
118	! ---------------------
119	! Radius axis of the interpolation grid
120	REAL,SAVE :: refftab(refftabsize)
121	! Variance axis of the interpolation grid
122	REAL,SAVE :: nuefftab(nuefftabsize)
123	! Volume ratio of the grid
124	REAL,SAVE :: logvratgrid,vratgrid
125	! Grid used to remember which calculation is done
126	LOGICAL,SAVE,ALLOCATABLE :: checkgrid(:,:,:,:)
127	!$OMP THREADPRIVATE(refftab,nuefftab,logvratgrid,vratgrid,checkgrid)
128	! Optical properties of the grid (VISIBLE)
129	REAL,SAVE,ALLOCATABLE :: qsqrefVISgrid(:,:,:,:)
130	REAL,SAVE,ALLOCATABLE :: qextVISgrid(:,:,:,:)
131	REAL,SAVE,ALLOCATABLE :: qscatVISgrid(:,:,:,:)
132	REAL,SAVE,ALLOCATABLE :: omegVISgrid(:,:,:,:)
133	REAL,SAVE,ALLOCATABLE :: gVISgrid(:,:,:,:)
134	!$OMP THREADPRIVATE(qsqrefVISgrid,qextVISgrid,qscatVISgrid,omegVISgrid,gVISgrid)
135	! Optical properties of the grid (INFRARED)
136	REAL,SAVE,ALLOCATABLE :: qsqrefIRgrid(:,:,:,:)
137	REAL,SAVE,ALLOCATABLE :: qextIRgrid(:,:,:,:)
138	REAL,SAVE,ALLOCATABLE :: qscatIRgrid(:,:,:,:)
139	REAL,SAVE,ALLOCATABLE :: omegIRgrid(:,:,:,:)
140	REAL,SAVE,ALLOCATABLE :: gIRgrid(:,:,:,:)
141	!$OMP THREADPRIVATE(qsqrefIRgrid,qextIRgrid,qscatIRgrid,omegIRgrid,gIRgrid)
142	! Optical properties of the grid (REFERENCE WAVELENGTHS)
143	REAL,SAVE,ALLOCATABLE :: qrefVISgrid(:,:,:)
144	REAL,SAVE,ALLOCATABLE :: qscatrefVISgrid(:,:,:)
145	REAL,SAVE,ALLOCATABLE :: qrefIRgrid(:,:,:)
146	REAL,SAVE,ALLOCATABLE :: qscatrefIRgrid(:,:,:)
147	REAL,SAVE,ALLOCATABLE :: omegrefVISgrid(:,:,:)
148	REAL,SAVE,ALLOCATABLE :: omegrefIRgrid(:,:,:)
149	!$OMP THREADPRIVATE(qrefVISgrid,qscatrefVISgrid,qrefIRgrid,qscatrefIRgrid,omegrefVISgrid,&
150	!$OMP omegrefIRgrid)
151	! Firstcall
152	LOGICAL,SAVE :: firstcall = .true.
153	LOGICAL,SAVE :: first_allocate=.true.
154	!$OMP THREADPRIVATE(firstcall,first_allocate)
155	! Variables used by the Gauss-Legendre integration:
156	REAL,SAVE,ALLOCATABLE :: normd(:,:,:,:)
157	REAL,SAVE,ALLOCATABLE :: dista(:,:,:,:,:)
158	REAL,SAVE,ALLOCATABLE :: distb(:,:,:,:,:)
159	!$OMP THREADPRIVATE(normd,dista,distb)
160
161	REAL,SAVE,ALLOCATABLE :: radGAUSa(:,:,:)
162	REAL,SAVE,ALLOCATABLE :: radGAUSb(:,:,:)
163	!$OMP THREADPRIVATE(radGAUSa,radGAUSb)
164
165	REAL,SAVE,ALLOCATABLE :: qsqrefVISa(:,:,:)
166	REAL,SAVE,ALLOCATABLE :: qrefVISa(:,:)
167	REAL,SAVE,ALLOCATABLE :: qsqrefVISb(:,:,:)
168	REAL,SAVE,ALLOCATABLE :: qrefVISb(:,:)
169	!$OMP THREADPRIVATE(qsqrefVISa,qrefVISa,qsqrefVISb,qrefVISb)
170	REAL,SAVE,ALLOCATABLE :: omegVISa(:,:,:)
171	REAL,SAVE,ALLOCATABLE :: omegrefVISa(:,:)
172	REAL,SAVE,ALLOCATABLE :: omegVISb(:,:,:)
173	REAL,SAVE,ALLOCATABLE :: omegrefVISb(:,:)
174	REAL,SAVE,ALLOCATABLE :: gVISa(:,:,:)
175	REAL,SAVE,ALLOCATABLE :: gVISb(:,:,:)
176	!$OMP THREADPRIVATE(omegVISa,omegrefVISa,omegVISb,omegrefVISb,gVISa,gVISb)
177
178	REAL,SAVE,ALLOCATABLE :: qsqrefIRa(:,:,:)
179	REAL,SAVE,ALLOCATABLE :: qrefIRa(:,:)
180	REAL,SAVE,ALLOCATABLE :: qsqrefIRb(:,:,:)
181	REAL,SAVE,ALLOCATABLE :: qrefIRb(:,:)
182	!$OMP THREADPRIVATE(qsqrefIRa,qrefIRa,qsqrefIRb,qrefIRb)
183	REAL,SAVE,ALLOCATABLE :: omegIRa(:,:,:)
184	REAL,SAVE,ALLOCATABLE :: omegrefIRa(:,:)
185	REAL,SAVE,ALLOCATABLE :: omegIRb(:,:,:)
186	REAL,SAVE,ALLOCATABLE :: omegrefIRb(:,:)
187	REAL,SAVE,ALLOCATABLE :: gIRa(:,:,:)
188	REAL,SAVE,ALLOCATABLE :: gIRb(:,:,:)
189	!$OMP THREADPRIVATE(omegIRa,omegrefIRa,omegIRb,omegrefIRb,gIRa,gIRb)
190
191	! Local variables:
192	!-----------------
193	REAL :: radiusm
194	REAL :: radiusr
195
196	REAL :: minrad ! minimal radius in table .dat radiustab (outside 0)
197	REAL :: maxrad
198
199	! 0. Allocate local saved arrays at firstcall
200	! --------------------------------------------------
201	IF (first_allocate) THEN
202	! Grid used to remember computations already done at previous calls
203	ALLOCATE(checkgrid(refftabsize,nuefftabsize,naerkind,2))
204	checkgrid(:,:,:,:)=.false.
205	! Optical properties of the grid (VISIBLE)
206	ALLOCATE(qsqrefVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind))
207	ALLOCATE(qextVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind))
208	ALLOCATE(qscatVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind))
209	ALLOCATE(omegVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind))
210	ALLOCATE(gVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind))
211	! Optical properties of the grid (INFRARED)
212	ALLOCATE(qsqrefIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind))
213	ALLOCATE(qextIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind))
214	ALLOCATE(qscatIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind))
215	ALLOCATE(omegIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind))
216	ALLOCATE(gIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind))
217	! Optical properties of the grid (REFERENCE WAVELENGTHS)
218	ALLOCATE(qrefVISgrid(refftabsize,nuefftabsize,naerkind))
219	ALLOCATE(qscatrefVISgrid(refftabsize,nuefftabsize,naerkind))
220	ALLOCATE(qrefIRgrid(refftabsize,nuefftabsize,naerkind))
221	ALLOCATE(qscatrefIRgrid(refftabsize,nuefftabsize,naerkind))
222	ALLOCATE(omegrefVISgrid(refftabsize,nuefftabsize,naerkind))
223	ALLOCATE(omegrefIRgrid(refftabsize,nuefftabsize,naerkind))
224	! Variables used by the Gauss-Legendre integration:
225	ALLOCATE(normd(refftabsize,nuefftabsize,naerkind,2))
226	ALLOCATE(dista(refftabsize,nuefftabsize,naerkind,2,ngau))
227	ALLOCATE(distb(refftabsize,nuefftabsize,naerkind,2,ngau))
228	ALLOCATE(radGAUSa(ngau,naerkind,2))
229	ALLOCATE(radGAUSb(ngau,naerkind,2))
230	!
231	ALLOCATE(qsqrefVISa(L_NSPECTV,ngau,naerkind))
232	ALLOCATE(qrefVISa(ngau,naerkind))
233	ALLOCATE(qsqrefVISb(L_NSPECTV,ngau,naerkind))
234	ALLOCATE(qrefVISb(ngau,naerkind))
235	ALLOCATE(omegVISa(L_NSPECTV,ngau,naerkind))
236	ALLOCATE(omegrefVISa(ngau,naerkind))
237	ALLOCATE(omegVISb(L_NSPECTV,ngau,naerkind))
238	ALLOCATE(omegrefVISb(ngau,naerkind))
239	ALLOCATE(gVISa(L_NSPECTV,ngau,naerkind))
240	ALLOCATE(gVISb(L_NSPECTV,ngau,naerkind))
241	!
242	ALLOCATE(qsqrefIRa(L_NSPECTI,ngau,naerkind))
243	ALLOCATE(qrefIRa(ngau,naerkind))
244	ALLOCATE(qsqrefIRb(L_NSPECTI,ngau,naerkind))
245	ALLOCATE(qrefIRb(ngau,naerkind))
246
247	ALLOCATE(omegIRa(L_NSPECTI,ngau,naerkind))
248	ALLOCATE(omegrefIRa(ngau,naerkind))
249	ALLOCATE(omegIRb(L_NSPECTI,ngau,naerkind))
250	ALLOCATE(omegrefIRb(ngau,naerkind))
251	ALLOCATE(gIRa(L_NSPECTI,ngau,naerkind))
252	ALLOCATE(gIRb(L_NSPECTI,ngau,naerkind))
253
254	first_allocate=.false.
255	ENDIF ! of IF (first_allocate)
256
257
258	DO iaer = 1, naerkind ! Loop on aerosol kind
259	!==================================================================
260	! If there is one single particle size, optical
261	! properties of the considered aerosol are homogeneous
262	IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN
263
264	DO lg = 1, nlayer
265	DO ig = 1, ngrid
266	DO chg = 1, L_NSPECTV
267	QVISsQREF3d(ig,lg,chg,iaer)=QVISsQREF(chg,iaer,1)
268	omegaVIS3d(ig,lg,chg,iaer)=omegaVIS(chg,iaer,1)
269	gVIS3d(ig,lg,chg,iaer)=gVIS(chg,iaer,1)
270	ENDDO
271	DO chg = 1, L_NSPECTI
272	QIRsQREF3d(ig,lg,chg,iaer)=QIRsQREF(chg,iaer,1)
273	omegaIR3d(ig,lg,chg,iaer)=omegaIR(chg,iaer,1)
274	gIR3d(ig,lg,chg,iaer)=gIR(chg,iaer,1)
275	ENDDO
276	QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1)
277	QREFir3d(ig,lg,iaer)=QREFir(iaer,1)
278	!omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1)
279	!omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1)
280	ENDDO
281	ENDDO
282
283	if (firstcall) then
284	print*,'Optical prop. of the aerosol are homogenous for:'
285	print*,'iaer = ',iaer
286	endif
287
288	!==================================================================
289	! Various particule size
290	! Varying effective radius and variance
291	ELSE
292
293	DO idomain = 1, 2 ! Loop on visible or infrared channel
294	! 1. Creating the effective radius and variance grid
295	! --------------------------------------------------
296	IF (firstcall) THEN
297
298	! 1.1 Pi
299	!-------
300	pi = 2. * asin(1.e0)
301
302	! 1.2 Effective radius
303	!---------------------
304	refftab(1) = refftabmin
305	refftab(refftabsize) = refftabmax
306
307	logvratgrid = log(refftabmax/refftabmin) / float(refftabsize-1)*3.
308	vratgrid = exp(logvratgrid)
309
310	do i = 2, refftabsize-1
311	refftab(i) = refftab(i-1)vratgrid*(1./3.)
312	enddo
313
314	! 1.3 Effective variance
315	!-----------------------
316	if(nuefftabsize.eq.1)then ! addded by RDW
317	print*,'Warning: no variance range in aeroptproperties'
318	nuefftab(1)=0.2
319	else
320	do i = 0, nuefftabsize-1
321	nuefftab(i+1) = exp( nuefftabmin + i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
322	enddo
323	endif
324
325	! check that radiustab (table.dat) cover the reffrad used for haze : iaer=1
326	minrad=min(MINVAL(radiustab(1,1,1:nsize(1,1))),MINVAL(radiustab(1,2,1:nsize(1,2))))
327	maxrad=min(MAXVAL(radiustab(1,1,1:nsize(1,1))),MAXVAL(radiustab(1,2,1:nsize(1,2))))
328	IF ((MINVAL(reffrad).LT.minrad).OR.(MAXVAL(reffrad).GT.maxrad)) then
329	WRITE(,) 'Warning: particle size in grid box #'
330	WRITE(,) i,' is larger than optical properties. '
331	WRITE(,) 'reffrad=',MINVAL(reffrad),'-',MAXVAL(reffrad)
332	WRITE(,) 'radiustab=',minrad,'-',maxrad
333
334	! ensure reffrad is within bounds of radiustab
335	WHERE(reffrad.LT.minrad)
336	reffrad=minrad
337	ENDWHERE
338	WHERE(reffrad.GT.maxrad)
339	reffrad=maxrad
340	ENDWHERE
341	WRITE(,) 'Truncated reffrad within radiustab bounds:'
342	WRITE(,) 'New reffrad=',MINVAL(reffrad),'-',MAXVAL(reffrad)
343	ENDIF
344
345	firstcall = .false.
346	ENDIF ! of IF (firstcall)
347
348	! 1.4 Radius middle point and range for Gauss integration
349	!--------------------------------------------------------
350	radiusm=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) + radiustab(iaer,idomain,1))
351	radiusr=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) - radiustab(iaer,idomain,1))
352
353	! 1.5 Interpolating data at the Gauss quadrature points:
354	!-------------------------------------------------------
355	DO gausind=1,ngau
356	drad=radiusr*radgaus(gausind)
357	radGAUSa(gausind,iaer,idomain)=radiusm-drad
358
359	radius_id=minloc(abs(radiustab(iaer,idomain,:) - (radiusm-drad)),DIM=1)
360	IF ((radiustab(iaer,idomain,radius_id) - (radiusm-drad)).GT.0) THEN
361	radius_id=radius_id-1
362	ENDIF
363	IF (radius_id.GE.nsize(iaer,idomain)) THEN
364	radius_id=nsize(iaer,idomain)-1
365	kint = 1.
366	ELSEIF (radius_id.LT.1) THEN
367	radius_id=1
368	kint = 0.
369	ELSE
370	kint = ( (radiusm-drad) - &
371	radiustab(iaer,idomain,radius_id) ) / &
372	( radiustab(iaer,idomain,radius_id+1) - &
373	radiustab(iaer,idomain,radius_id) )
374	ENDIF
375	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
376	DO m=1,L_NSPECTV
377	qsqrefVISa(m,gausind,iaer)= &
378	(1-kint)*QVISsQREF(m,iaer,radius_id) + &
379	kint*QVISsQREF(m,iaer,radius_id+1)
380	omegVISa(m,gausind,iaer)= &
381	(1-kint)*omegaVIS(m,iaer,radius_id) + &
382	kint*omegaVIS(m,iaer,radius_id+1)
383	gVISa(m,gausind,iaer)= &
384	(1-kint)*gVIS(m,iaer,radius_id) + &
385	kint*gVIS(m,iaer,radius_id+1)
386	ENDDO
387	qrefVISa(gausind,iaer)= &
388	(1-kint)*QREFvis(iaer,radius_id) + &
389	kint*QREFvis(iaer,radius_id+1)
390	omegrefVISa(gausind,iaer)= 0
391	! omegrefVISa(gausind,iaer)= &
392	! (1-kint)*omegaREFvis(iaer,radius_id) + &
393	! kint*omegaREFvis(iaer,radius_id+1)
394	ELSE ! INFRARED DOMAIN ----------------------------------
395	DO m=1,L_NSPECTI
396	qsqrefIRa(m,gausind,iaer)= &
397	(1-kint)*QIRsQREF(m,iaer,radius_id) + &
398	kint*QIRsQREF(m,iaer,radius_id+1)
399	omegIRa(m,gausind,iaer)= &
400	(1-kint)*omegaIR(m,iaer,radius_id) + &
401	kint*omegaIR(m,iaer,radius_id+1)
402	gIRa(m,gausind,iaer)= &
403	(1-kint)*gIR(m,iaer,radius_id) + &
404	kint*gIR(m,iaer,radius_id+1)
405	ENDDO
406	qrefIRa(gausind,iaer)= &
407	(1-kint)*QREFir(iaer,radius_id) + &
408	kint*QREFir(iaer,radius_id+1)
409	omegrefIRa(gausind,iaer)= &
410	(1-kint)*omegaREFir(iaer,radius_id) + &
411	kint*omegaREFir(iaer,radius_id+1)
412	ENDIF
413	ENDDO
414
415	DO gausind=1,ngau
416	drad=radiusr*radgaus(gausind)
417	radGAUSb(gausind,iaer,idomain)=radiusm+drad
418
419	radius_id=minloc(abs(radiustab(iaer,idomain,:) - &
420	(radiusm+drad)),DIM=1)
421	IF ((radiustab(iaer,idomain,radius_id) - &
422	(radiusm+drad)).GT.0) THEN
423	radius_id=radius_id-1
424	ENDIF
425	IF (radius_id.GE.nsize(iaer,idomain)) THEN
426	radius_id=nsize(iaer,idomain)-1
427	kint = 1.
428	ELSEIF (radius_id.LT.1) THEN
429	radius_id=1
430	kint = 0.
431	ELSE
432	kint = ( (radiusm+drad) - &
433	radiustab(iaer,idomain,radius_id) ) / &
434	( radiustab(iaer,idomain,radius_id+1) - &
435	radiustab(iaer,idomain,radius_id) )
436	ENDIF
437	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
438	DO m=1,L_NSPECTV
439	qsqrefVISb(m,gausind,iaer)= &
440	(1-kint)*QVISsQREF(m,iaer,radius_id) + &
441	kint*QVISsQREF(m,iaer,radius_id+1)
442	omegVISb(m,gausind,iaer)= &
443	(1-kint)*omegaVIS(m,iaer,radius_id) + &
444	kint*omegaVIS(m,iaer,radius_id+1)
445	gVISb(m,gausind,iaer)= &
446	(1-kint)*gVIS(m,iaer,radius_id) + &
447	kint*gVIS(m,iaer,radius_id+1)
448	ENDDO
449	qrefVISb(gausind,iaer)= &
450	(1-kint)*QREFvis(iaer,radius_id) + &
451	kint*QREFvis(iaer,radius_id+1)
452	omegrefVISb(gausind,iaer)= 0
453	! omegrefVISb(gausind,iaer)= &
454	! (1-kint)*omegaREFvis(iaer,radius_id) + &
455	! kint*omegaREFvis(iaer,radius_id+1)
456	ELSE ! INFRARED DOMAIN ----------------------------------
457	DO m=1,L_NSPECTI
458	qsqrefIRb(m,gausind,iaer)= &
459	(1-kint)*QIRsQREF(m,iaer,radius_id) + &
460	kint*QIRsQREF(m,iaer,radius_id+1)
461	omegIRb(m,gausind,iaer)= &
462	(1-kint)*omegaIR(m,iaer,radius_id) + &
463	kint*omegaIR(m,iaer,radius_id+1)
464	gIRb(m,gausind,iaer)= &
465	(1-kint)*gIR(m,iaer,radius_id) + &
466	kint*gIR(m,iaer,radius_id+1)
467	ENDDO
468	qrefIRb(gausind,iaer)= &
469	(1-kint)*QREFir(iaer,radius_id) + &
470	kint*QREFir(iaer,radius_id+1)
471	omegrefIRb(gausind,iaer)= &
472	(1-kint)*omegaREFir(iaer,radius_id) + &
473	kint*omegaREFir(iaer,radius_id+1)
474	ENDIF
475	ENDDO
476
477	!==================================================================
478	! CONSTANT NUEFF FROM HERE
479	!==================================================================
480
481	! 2. Compute the scattering parameters using linear
482	! interpolation over grain sizes and constant nueff
483	! ---------------------------------------------------
484
485	DO lg = 1,nlayer
486	DO ig = 1, ngrid
487	! 2.1 Effective radius index and kx calculation
488	var_tmp=reffrad(ig,lg,iaer)/refftabmin
489	if (var_tmp.GT.0.) then ! for debug
490	var_tmp=log(var_tmp)*3.
491	else
492	var_tmp=0.
493	endif
494	var_tmp=var_tmp/logvratgrid+1.
495	grid_i=floor(var_tmp)
496	IF (grid_i.GE.refftabsize) THEN
497	! WRITE(,) 'Warning: particle size in grid box #'
498	! WRITE(,) ig,' is too large to be used by the '
499	! WRITE(,) 'radiative transfer; please extend the '
500	! WRITE(,) 'interpolation grid to larger grain sizes.'
501	grid_i=refftabsize-1
502	kx = 1.
503	ELSEIF (grid_i.LT.1) THEN
504	! WRITE(,) 'Warning: particle size in grid box #'
505	! WRITE(,) ig,' is too small to be used by the '
506	! WRITE(,) 'radiative transfer; please extend the '
507	! WRITE(,) 'interpolation grid to smaller grain sizes.'
508	grid_i=1
509	kx = 0.
510	ELSE
511	kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) / &
512	( refftab(grid_i+1)-refftab(grid_i) )
513	ENDIF
514	! 2.3 Integration
515	DO j=grid_i,grid_i+1
516	! 2.3.1 Check if the calculation has been done
517	IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
518	! 2.3.2 Log-normal dist., r_g and sigma_g are defined
519	! in [hansen_1974], "Light scattering in planetary
520	! atmospheres", Space Science Reviews 16 527-610.
521	! Here, sizedistk1=r_g and sizedistk2=sigma_g^2
522	sizedistk2 = log(1.+nueffrad(1,1,iaer))
523	sizedistk1 = exp(2.5*sizedistk2)
524	sizedistk1 = refftab(j) / sizedistk1
525
526	normd(j,1,iaer,idomain) = 1e-30
527	DO gausind=1,ngau
528	drad=radiusr*radgaus(gausind)
529	dista(j,1,iaer,idomain,gausind) = LOG((radiusm-drad)/sizedistk1)
530	dista(j,1,iaer,idomain,gausind) = &
531	EXP(-dista(j,1,iaer,idomain,gausind) * &
532	dista(j,1,iaer,idomain,gausind) * &
533	0.5e0/sizedistk2)/(radiusm-drad)
534	dista(j,1,iaer,idomain,gausind) = &
535	dista(j,1,iaer,idomain,gausind) / &
536	(sqrt(2e0pisizedistk2))
537
538	distb(j,1,iaer,idomain,gausind) = LOG((radiusm+drad)/sizedistk1)
539	distb(j,1,iaer,idomain,gausind) = &
540	EXP(-distb(j,1,iaer,idomain,gausind) * &
541	distb(j,1,iaer,idomain,gausind) * &
542	0.5e0/sizedistk2)/(radiusm+drad)
543	distb(j,1,iaer,idomain,gausind) = &
544	distb(j,1,iaer,idomain,gausind) / &
545	(sqrt(2e0pisizedistk2))
546
547	normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) + &
548	weightgaus(gausind) * &
549	( &
550	distb(j,1,iaer,idomain,gausind) * pi * &
551	radGAUSb(gausind,iaer,idomain) * &
552	radGAUSb(gausind,iaer,idomain) + &
553	dista(j,1,iaer,idomain,gausind) * pi * &
554	radGAUSa(gausind,iaer,idomain) * &
555	radGAUSa(gausind,iaer,idomain) &
556	)
557	ENDDO
558	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
559	! 2.3.3.vis Initialization
560	qsqrefVISgrid(j,1,:,iaer)=0.
561	qextVISgrid(j,1,:,iaer)=0.
562	qscatVISgrid(j,1,:,iaer)=0.
563	omegVISgrid(j,1,:,iaer)=0.
564	gVISgrid(j,1,:,iaer)=0.
565	qrefVISgrid(j,1,iaer)=0.
566	qscatrefVISgrid(j,1,iaer)=0.
567	omegrefVISgrid(j,1,iaer)=0.
568
569	DO gausind=1,ngau
570	DO m=1,L_NSPECTV
571	! Convolution:
572	qextVISgrid(j,1,m,iaer) = &
573	qextVISgrid(j,1,m,iaer) + &
574	weightgaus(gausind) * &
575	( &
576	qsqrefVISb(m,gausind,iaer) * &
577	qrefVISb(gausind,iaer) * &
578	piradGAUSb(gausind,iaer,idomain) &
579	radGAUSb(gausind,iaer,idomain) * &
580	distb(j,1,iaer,idomain,gausind) + &
581	qsqrefVISa(m,gausind,iaer) * &
582	qrefVISa(gausind,iaer) * &
583	piradGAUSa(gausind,iaer,idomain) &
584	radGAUSa(gausind,iaer,idomain) * &
585	dista(j,1,iaer,idomain,gausind) &
586	)
587	qscatVISgrid(j,1,m,iaer) = &
588	qscatVISgrid(j,1,m,iaer) + &
589	weightgaus(gausind) * &
590	( &
591	omegVISb(m,gausind,iaer) * &
592	qsqrefVISb(m,gausind,iaer) * &
593	qrefVISb(gausind,iaer) * &
594	piradGAUSb(gausind,iaer,idomain) &
595	radGAUSb(gausind,iaer,idomain) * &
596	distb(j,1,iaer,idomain,gausind) + &
597	omegVISa(m,gausind,iaer) * &
598	qsqrefVISa(m,gausind,iaer) * &
599	qrefVISa(gausind,iaer) * &
600	piradGAUSa(gausind,iaer,idomain) &
601	radGAUSa(gausind,iaer,idomain) * &
602	dista(j,1,iaer,idomain,gausind) &
603	)
604	gVISgrid(j,1,m,iaer) = &
605	gVISgrid(j,1,m,iaer) + &
606	weightgaus(gausind) * &
607	( &
608	omegVISb(m,gausind,iaer) * &
609	qsqrefVISb(m,gausind,iaer) * &
610	qrefVISb(gausind,iaer) * &
611	gVISb(m,gausind,iaer) * &
612	piradGAUSb(gausind,iaer,idomain) &
613	radGAUSb(gausind,iaer,idomain) * &
614	distb(j,1,iaer,idomain,gausind) + &
615	omegVISa(m,gausind,iaer) * &
616	qsqrefVISa(m,gausind,iaer) * &
617	qrefVISa(gausind,iaer) * &
618	gVISa(m,gausind,iaer) * &
619	piradGAUSa(gausind,iaer,idomain) &
620	radGAUSa(gausind,iaer,idomain) * &
621	dista(j,1,iaer,idomain,gausind) &
622	)
623	ENDDO
624	qrefVISgrid(j,1,iaer) = &
625	qrefVISgrid(j,1,iaer) + &
626	weightgaus(gausind) * &
627	( &
628	qrefVISb(gausind,iaer) * &
629	piradGAUSb(gausind,iaer,idomain) &
630	radGAUSb(gausind,iaer,idomain) * &
631	distb(j,1,iaer,idomain,gausind) + &
632	qrefVISa(gausind,iaer) * &
633	piradGAUSa(gausind,iaer,idomain) &
634	radGAUSa(gausind,iaer,idomain) * &
635	dista(j,1,iaer,idomain,gausind) &
636	)
637	qscatrefVISgrid(j,1,iaer) = &
638	qscatrefVISgrid(j,1,iaer) + &
639	weightgaus(gausind) * &
640	( &
641	omegrefVISb(gausind,iaer) * &
642	qrefVISb(gausind,iaer) * &
643	piradGAUSb(gausind,iaer,idomain) &
644	radGAUSb(gausind,iaer,idomain) * &
645	distb(j,1,iaer,idomain,gausind) + &
646	omegrefVISa(gausind,iaer) * &
647	qrefVISa(gausind,iaer) * &
648	piradGAUSa(gausind,iaer,idomain) &
649	radGAUSa(gausind,iaer,idomain) * &
650	dista(j,1,iaer,idomain,gausind) &
651	)
652	ENDDO
653
654	qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) / &
655	normd(j,1,iaer,idomain)
656	qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / &
657	normd(j,1,iaer,idomain)
658	if (qrefVISgrid(j,1,iaer).gt.0.d0) then
659	omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / &
660	qrefVISgrid(j,1,iaer)
661	else
662	omegrefVISgrid(j,1,iaer)=0.
663	endif
664	DO m=1,L_NSPECTV
665	qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / &
666	normd(j,1,iaer,idomain)
667	qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / &
668	normd(j,1,iaer,idomain)
669	if (qscatVISgrid(j,1,m,iaer).gt.0.d0) then
670	gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) / &
671	qscatVISgrid(j,1,m,iaer) / &
672	normd(j,1,iaer,idomain)
673	else
674	gVISgrid(j,1,m,iaer)=0.
675	endif
676
677	if (qrefVISgrid(j,1,iaer).gt.0.d0) then
678	qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /&
679	qrefVISgrid(j,1,iaer)
680	else
681	qsqrefVISgrid(j,1,m,iaer)=0.
682	endif
683	if (qextVISgrid(j,1,m,iaer).gt.0.d0) then
684	omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / &
685	qextVISgrid(j,1,m,iaer)
686	else
687	omegVISgrid(j,1,m,iaer)=0.
688	endif
689	ENDDO
690	ELSE ! INFRARED DOMAIN ----------
691	! 2.3.3.ir Initialization
692	qsqrefIRgrid(j,1,:,iaer)=0.
693	qextIRgrid(j,1,:,iaer)=0.
694	qscatIRgrid(j,1,:,iaer)=0.
695	omegIRgrid(j,1,:,iaer)=0.
696	gIRgrid(j,1,:,iaer)=0.
697	qrefIRgrid(j,1,iaer)=0.
698	qscatrefIRgrid(j,1,iaer)=0.
699	omegrefIRgrid(j,1,iaer)=0.
700
701	DO gausind=1,ngau
702	DO m=1,L_NSPECTI
703	! Convolution:
704	qextIRgrid(j,1,m,iaer) = &
705	qextIRgrid(j,1,m,iaer) + &
706	weightgaus(gausind) * &
707	( &
708	qsqrefIRb(m,gausind,iaer) * &
709	qrefVISb(gausind,iaer) * &
710	piradGAUSb(gausind,iaer,idomain) &
711	radGAUSb(gausind,iaer,idomain) * &
712	distb(j,1,iaer,idomain,gausind) + &
713	qsqrefIRa(m,gausind,iaer) * &
714	qrefVISa(gausind,iaer) * &
715	piradGAUSa(gausind,iaer,idomain) &
716	radGAUSa(gausind,iaer,idomain) * &
717	dista(j,1,iaer,idomain,gausind) &
718	)
719	qscatIRgrid(j,1,m,iaer) = &
720	qscatIRgrid(j,1,m,iaer) + &
721	weightgaus(gausind) * &
722	( &
723	omegIRb(m,gausind,iaer) * &
724	qsqrefIRb(m,gausind,iaer) * &
725	qrefVISb(gausind,iaer) * &
726	piradGAUSb(gausind,iaer,idomain) &
727	radGAUSb(gausind,iaer,idomain) * &
728	distb(j,1,iaer,idomain,gausind) + &
729	omegIRa(m,gausind,iaer) * &
730	qsqrefIRa(m,gausind,iaer) * &
731	qrefVISa(gausind,iaer) * &
732	piradGAUSa(gausind,iaer,idomain) &
733	radGAUSa(gausind,iaer,idomain) * &
734	dista(j,1,iaer,idomain,gausind) &
735	)
736	gIRgrid(j,1,m,iaer) = &
737	gIRgrid(j,1,m,iaer) + &
738	weightgaus(gausind) * &
739	( &
740	omegIRb(m,gausind,iaer) * &
741	qsqrefIRb(m,gausind,iaer) * &
742	qrefVISb(gausind,iaer) * &
743	gIRb(m,gausind,iaer) * &
744	piradGAUSb(gausind,iaer,idomain) &
745	radGAUSb(gausind,iaer,idomain) * &
746	distb(j,1,iaer,idomain,gausind) + &
747	omegIRa(m,gausind,iaer) * &
748	qsqrefIRa(m,gausind,iaer) * &
749	qrefVISa(gausind,iaer) * &
750	gIRa(m,gausind,iaer) * &
751	piradGAUSa(gausind,iaer,idomain) &
752	radGAUSa(gausind,iaer,idomain) * &
753	dista(j,1,iaer,idomain,gausind) &
754	)
755	ENDDO
756	qrefIRgrid(j,1,iaer) = &
757	qrefIRgrid(j,1,iaer) + &
758	weightgaus(gausind) * &
759	( &
760	qrefIRb(gausind,iaer) * &
761	piradGAUSb(gausind,iaer,idomain) &
762	radGAUSb(gausind,iaer,idomain) * &
763	distb(j,1,iaer,idomain,gausind) + &
764	qrefIRa(gausind,iaer) * &
765	piradGAUSa(gausind,iaer,idomain) &
766	radGAUSa(gausind,iaer,idomain) * &
767	dista(j,1,iaer,idomain,gausind) &
768	)
769	qscatrefIRgrid(j,1,iaer) = &
770	qscatrefIRgrid(j,1,iaer) + &
771	weightgaus(gausind) * &
772	( &
773	omegrefIRb(gausind,iaer) * &
774	qrefIRb(gausind,iaer) * &
775	piradGAUSb(gausind,iaer,idomain) &
776	radGAUSb(gausind,iaer,idomain) * &
777	distb(j,1,iaer,idomain,gausind) + &
778	omegrefIRa(gausind,iaer) * &
779	qrefIRa(gausind,iaer) * &
780	piradGAUSa(gausind,iaer,idomain) &
781	radGAUSa(gausind,iaer,idomain) * &
782	dista(j,1,iaer,idomain,gausind) &
783	)
784	ENDDO
785
786	qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) / &
787	normd(j,1,iaer,idomain)
788	qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / &
789	normd(j,1,iaer,idomain)
790	if (qrefIRgrid(j,1,iaer).gt.0.d0) then
791	omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / &
792	qrefIRgrid(j,1,iaer)
793	else
794	omegrefIRgrid(j,1,iaer)=0.
795	endif
796	DO m=1,L_NSPECTI
797	qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / &
798	normd(j,1,iaer,idomain)
799	qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / &
800	normd(j,1,iaer,idomain)
801
802	if (qscatIRgrid(j,1,m,iaer).gt.0.d0) then
803	gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) / &
804	qscatIRgrid(j,1,m,iaer) / &
805	normd(j,1,iaer,idomain)
806	else
807	gIRgrid(j,1,m,iaer)=0.
808	endif
809
810	if (qrefIRgrid(j,1,iaer).gt.0.d0) then
811	qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /&
812	qrefVISgrid(j,1,iaer)
813	else
814	qsqrefIRgrid(j,1,m,iaer)=0.
815	endif
816	if (qextIRgrid(j,1,m,iaer).gt.0.d0) then
817	omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / &
818	qextIRgrid(j,1,m,iaer)
819	else
820	omegIRgrid(j,1,m,iaer)=0.
821	endif
822	ENDDO
823	ENDIF ! --------------------------
824	checkgrid(j,1,iaer,idomain) = .true.
825	ENDIF !checkgrid
826	ENDDO !grid_i
827	! 2.4 Linear interpolation
828	k1 = (1-kx)
829	k2 = kx
830	IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
831	DO m=1,L_NSPECTV
832	QVISsQREF3d(ig,lg,m,iaer) = &
833	k1*qsqrefVISgrid(grid_i,1,m,iaer) + &
834	k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
835	omegaVIS3d(ig,lg,m,iaer) = &
836	k1*omegVISgrid(grid_i,1,m,iaer) + &
837	k2*omegVISgrid(grid_i+1,1,m,iaer)
838	gVIS3d(ig,lg,m,iaer) = &
839	k1*gVISgrid(grid_i,1,m,iaer) + &
840	k2*gVISgrid(grid_i+1,1,m,iaer)
841	ENDDO !L_NSPECTV
842	QREFvis3d(ig,lg,iaer) = &
843	k1*qrefVISgrid(grid_i,1,iaer) + &
844	k2*qrefVISgrid(grid_i+1,1,iaer)
845	! omegaREFvis3d(ig,lg,iaer) = &
846	! k1*omegrefVISgrid(grid_i,1,iaer) + &
847	! k2*omegrefVISgrid(grid_i+1,1,iaer)
848	ELSE ! INFRARED -----------------------
849	DO m=1,L_NSPECTI
850	QIRsQREF3d(ig,lg,m,iaer) = &
851	k1*qsqrefIRgrid(grid_i,1,m,iaer) + &
852	k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
853	omegaIR3d(ig,lg,m,iaer) = &
854	k1*omegIRgrid(grid_i,1,m,iaer) + &
855	k2*omegIRgrid(grid_i+1,1,m,iaer)
856	gIR3d(ig,lg,m,iaer) = &
857	k1*gIRgrid(grid_i,1,m,iaer) + &
858	k2*gIRgrid(grid_i+1,1,m,iaer)
859	ENDDO !L_NSPECTI
860	QREFir3d(ig,lg,iaer) = &
861	k1*qrefIRgrid(grid_i,1,iaer) + &
862	k2*qrefIRgrid(grid_i+1,1,iaer)
863	! omegaREFir3d(ig,lg,iaer) = &
864	! k1*omegrefIRgrid(grid_i,1,iaer) + &
865	! k2*omegrefIRgrid(grid_i+1,1,iaer)
866	ENDIF ! --------------------------------
867	ENDDO !nlayer
868	ENDDO !ngrid
869
870	!==================================================================
871
872	ENDDO ! idomain
873
874	ENDIF ! nsize = 1
875
876	ENDDO ! iaer (loop on aerosol kind)
877
878	END SUBROUTINE aeroptproperties
879
880
881	end module aeroptproperties_mod

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

Download in other formats:

Original Format