Context Navigation

aeroptproperties.F90 @ 3506

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

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

Download in other formats:

Original Format