Context Navigation

aeroptproperties.F90 @ 3436

Last change on this file since 3436 was 2972, checked in by emillour, 18 months ago
Generic PCM: Make number of scatterers fully dynamic (i.e. set in callphys.def and no longer by compilation option "-s #"). One should now specify naerkind = # in callphys.def (default is 0). EM
File size: 39.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,INTENT(IN) :: 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	! 0. Allocate local saved arrays at firstcall
202	! --------------------------------------------------
203	IF (first_allocate) THEN
204	! Grid used to remember computations already done at previous calls
205	ALLOCATE(checkgrid(refftabsize,nuefftabsize,naerkind,2))
206	checkgrid(:,:,:,:)=.false.
207	! Optical properties of the grid (VISIBLE)
208	ALLOCATE(qsqrefVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind))
209	ALLOCATE(qextVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind))
210	ALLOCATE(qscatVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind))
211	ALLOCATE(omegVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind))
212	ALLOCATE(gVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind))
213	! Optical properties of the grid (INFRARED)
214	ALLOCATE(qsqrefIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind))
215	ALLOCATE(qextIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind))
216	ALLOCATE(qscatIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind))
217	ALLOCATE(omegIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind))
218	ALLOCATE(gIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind))
219	! Optical properties of the grid (REFERENCE WAVELENGTHS)
220	ALLOCATE(qrefVISgrid(refftabsize,nuefftabsize,naerkind))
221	ALLOCATE(qscatrefVISgrid(refftabsize,nuefftabsize,naerkind))
222	ALLOCATE(qrefIRgrid(refftabsize,nuefftabsize,naerkind))
223	ALLOCATE(qscatrefIRgrid(refftabsize,nuefftabsize,naerkind))
224	ALLOCATE(omegrefVISgrid(refftabsize,nuefftabsize,naerkind))
225	ALLOCATE(omegrefIRgrid(refftabsize,nuefftabsize,naerkind))
226	! Variables used by the Gauss-Legendre integration:
227	ALLOCATE(normd(refftabsize,nuefftabsize,naerkind,2))
228	ALLOCATE(dista(refftabsize,nuefftabsize,naerkind,2,ngau))
229	ALLOCATE(distb(refftabsize,nuefftabsize,naerkind,2,ngau))
230	ALLOCATE(radGAUSa(ngau,naerkind,2))
231	ALLOCATE(radGAUSb(ngau,naerkind,2))
232	!
233	ALLOCATE(qsqrefVISa(L_NSPECTV,ngau,naerkind))
234	ALLOCATE(qrefVISa(ngau,naerkind))
235	ALLOCATE(qsqrefVISb(L_NSPECTV,ngau,naerkind))
236	ALLOCATE(qrefVISb(ngau,naerkind))
237	ALLOCATE(omegVISa(L_NSPECTV,ngau,naerkind))
238	ALLOCATE(omegrefVISa(ngau,naerkind))
239	ALLOCATE(omegVISb(L_NSPECTV,ngau,naerkind))
240	ALLOCATE(omegrefVISb(ngau,naerkind))
241	ALLOCATE(gVISa(L_NSPECTV,ngau,naerkind))
242	ALLOCATE(gVISb(L_NSPECTV,ngau,naerkind))
243	!
244	ALLOCATE(qsqrefIRa(L_NSPECTI,ngau,naerkind))
245	ALLOCATE(qrefIRa(ngau,naerkind))
246	ALLOCATE(qsqrefIRb(L_NSPECTI,ngau,naerkind))
247	ALLOCATE(qrefIRb(ngau,naerkind))
248
249	ALLOCATE(omegIRa(L_NSPECTI,ngau,naerkind))
250	ALLOCATE(omegrefIRa(ngau,naerkind))
251	ALLOCATE(omegIRb(L_NSPECTI,ngau,naerkind))
252	ALLOCATE(omegrefIRb(ngau,naerkind))
253	ALLOCATE(gIRa(L_NSPECTI,ngau,naerkind))
254	ALLOCATE(gIRb(L_NSPECTI,ngau,naerkind))
255
256	first_allocate=.false.
257	ENDIF ! of IF (first_allocate)
258
259	DO iaer = 1, naerkind ! Loop on aerosol kind
260	IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN
261	!==================================================================
262	! If there is one single particle size, optical
263	! properties of the considered aerosol are homogeneous
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
284	if (firstcall) then
285	print*,'Optical prop. of the aerosol are homogenous for:'
286	print*,'iaer = ',iaer
287	endif
288
289	ELSE ! Varying effective radius and variance
290	DO idomain = 1, 2 ! Loop on visible or infrared channel
291	!==================================================================
292	! 1. Creating the effective radius and variance grid
293	! --------------------------------------------------
294	IF (firstcall) THEN
295
296	! 1.1 Pi!
297	pi = 2. * asin(1.e0)
298
299	! 1.2 Effective radius
300	refftab(1) = refftabmin
301	refftab(refftabsize) = refftabmax
302
303	logvratgrid = log(refftabmax/refftabmin) / float(refftabsize-1)*3.
304	vratgrid = exp(logvratgrid)
305
306	do i = 2, refftabsize-1
307	refftab(i) = refftab(i-1)vratgrid*(1./3.)
308	enddo
309
310	! 1.3 Effective variance
311	if(nuefftabsize.eq.1)then ! addded by RDW
312	print*,'Warning: no variance range in aeroptproperties'
313	nuefftab(1)=0.2
314	else
315	do i = 0, nuefftabsize-1
316	nuefftab(i+1) = exp( nuefftabmin + i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
317	enddo
318	endif
319
320	firstcall = .false.
321	ENDIF ! of IF (firstcall)
322
323	! 1.4 Radius middle point and range for Gauss integration
324	radiusm=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) + radiustab(iaer,idomain,1))
325	radiusr=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) - radiustab(iaer,idomain,1))
326
327	! 1.5 Interpolating data at the Gauss quadrature points:
328	DO gausind=1,ngau
329	drad=radiusr*radgaus(gausind)
330	radGAUSa(gausind,iaer,idomain)=radiusm-drad
331
332	radius_id=minloc(abs(radiustab(iaer,idomain,:) - (radiusm-drad)),DIM=1)
333	IF ((radiustab(iaer,idomain,radius_id) - (radiusm-drad)).GT.0) THEN
334	radius_id=radius_id-1
335	ENDIF
336	IF (radius_id.GE.nsize(iaer,idomain)) THEN
337	radius_id=nsize(iaer,idomain)-1
338	kint = 1.
339	ELSEIF (radius_id.LT.1) THEN
340	radius_id=1
341	kint = 0.
342	ELSE
343	kint = ( (radiusm-drad) - &
344	radiustab(iaer,idomain,radius_id) ) / &
345	( radiustab(iaer,idomain,radius_id+1) - &
346	radiustab(iaer,idomain,radius_id) )
347	ENDIF
348	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
349	DO m=1,L_NSPECTV
350	qsqrefVISa(m,gausind,iaer)= &
351	(1-kint)*QVISsQREF(m,iaer,radius_id) + &
352	kint*QVISsQREF(m,iaer,radius_id+1)
353	omegVISa(m,gausind,iaer)= &
354	(1-kint)*omegaVIS(m,iaer,radius_id) + &
355	kint*omegaVIS(m,iaer,radius_id+1)
356	gVISa(m,gausind,iaer)= &
357	(1-kint)*gVIS(m,iaer,radius_id) + &
358	kint*gVIS(m,iaer,radius_id+1)
359	ENDDO
360	qrefVISa(gausind,iaer)= &
361	(1-kint)*QREFvis(iaer,radius_id) + &
362	kint*QREFvis(iaer,radius_id+1)
363	omegrefVISa(gausind,iaer)= 0
364	! omegrefVISa(gausind,iaer)= &
365	! (1-kint)*omegaREFvis(iaer,radius_id) + &
366	! kint*omegaREFvis(iaer,radius_id+1)
367	ELSE ! INFRARED DOMAIN ----------------------------------
368	DO m=1,L_NSPECTI
369	qsqrefIRa(m,gausind,iaer)= &
370	(1-kint)*QIRsQREF(m,iaer,radius_id) + &
371	kint*QIRsQREF(m,iaer,radius_id+1)
372	omegIRa(m,gausind,iaer)= &
373	(1-kint)*omegaIR(m,iaer,radius_id) + &
374	kint*omegaIR(m,iaer,radius_id+1)
375	gIRa(m,gausind,iaer)= &
376	(1-kint)*gIR(m,iaer,radius_id) + &
377	kint*gIR(m,iaer,radius_id+1)
378	ENDDO
379	qrefIRa(gausind,iaer)= &
380	(1-kint)*QREFir(iaer,radius_id) + &
381	kint*QREFir(iaer,radius_id+1)
382	omegrefIRa(gausind,iaer)= &
383	(1-kint)*omegaREFir(iaer,radius_id) + &
384	kint*omegaREFir(iaer,radius_id+1)
385	ENDIF
386	ENDDO
387
388	DO gausind=1,ngau
389	drad=radiusr*radgaus(gausind)
390	radGAUSb(gausind,iaer,idomain)=radiusm+drad
391
392	radius_id=minloc(abs(radiustab(iaer,idomain,:) - &
393	(radiusm+drad)),DIM=1)
394	IF ((radiustab(iaer,idomain,radius_id) - &
395	(radiusm+drad)).GT.0) THEN
396	radius_id=radius_id-1
397	ENDIF
398	IF (radius_id.GE.nsize(iaer,idomain)) THEN
399	radius_id=nsize(iaer,idomain)-1
400	kint = 1.
401	ELSEIF (radius_id.LT.1) THEN
402	radius_id=1
403	kint = 0.
404	ELSE
405	kint = ( (radiusm+drad) - &
406	radiustab(iaer,idomain,radius_id) ) / &
407	( radiustab(iaer,idomain,radius_id+1) - &
408	radiustab(iaer,idomain,radius_id) )
409	ENDIF
410	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
411	DO m=1,L_NSPECTV
412	qsqrefVISb(m,gausind,iaer)= &
413	(1-kint)*QVISsQREF(m,iaer,radius_id) + &
414	kint*QVISsQREF(m,iaer,radius_id+1)
415	omegVISb(m,gausind,iaer)= &
416	(1-kint)*omegaVIS(m,iaer,radius_id) + &
417	kint*omegaVIS(m,iaer,radius_id+1)
418	gVISb(m,gausind,iaer)= &
419	(1-kint)*gVIS(m,iaer,radius_id) + &
420	kint*gVIS(m,iaer,radius_id+1)
421	ENDDO
422	qrefVISb(gausind,iaer)= &
423	(1-kint)*QREFvis(iaer,radius_id) + &
424	kint*QREFvis(iaer,radius_id+1)
425	omegrefVISb(gausind,iaer)= 0
426	! omegrefVISb(gausind,iaer)= &
427	! (1-kint)*omegaREFvis(iaer,radius_id) + &
428	! kint*omegaREFvis(iaer,radius_id+1)
429	ELSE ! INFRARED DOMAIN ----------------------------------
430	DO m=1,L_NSPECTI
431	qsqrefIRb(m,gausind,iaer)= &
432	(1-kint)*QIRsQREF(m,iaer,radius_id) + &
433	kint*QIRsQREF(m,iaer,radius_id+1)
434	omegIRb(m,gausind,iaer)= &
435	(1-kint)*omegaIR(m,iaer,radius_id) + &
436	kint*omegaIR(m,iaer,radius_id+1)
437	gIRb(m,gausind,iaer)= &
438	(1-kint)*gIR(m,iaer,radius_id) + &
439	kint*gIR(m,iaer,radius_id+1)
440	ENDDO
441	qrefIRb(gausind,iaer)= &
442	(1-kint)*QREFir(iaer,radius_id) + &
443	kint*QREFir(iaer,radius_id+1)
444	omegrefIRb(gausind,iaer)= &
445	(1-kint)*omegaREFir(iaer,radius_id) + &
446	kint*omegaREFir(iaer,radius_id+1)
447	ENDIF
448	ENDDO
449
450	!==================================================================
451	! CONSTANT NUEFF FROM HERE
452	!==================================================================
453
454	! 2. Compute the scattering parameters using linear
455	! interpolation over grain sizes and constant nueff
456	! ---------------------------------------------------
457
458	DO lg = 1,nlayer
459	DO ig = 1, ngrid
460	! 2.1 Effective radius index and kx calculation
461	var_tmp=reffrad(ig,lg,iaer)/refftabmin
462	var_tmp=log(var_tmp)*3.
463	var_tmp=var_tmp/logvratgrid+1.
464	grid_i=floor(var_tmp)
465	IF (grid_i.GE.refftabsize) THEN
466	! WRITE(,) 'Warning: particle size in grid box #'
467	! WRITE(,) ig,' is too large to be used by the '
468	! WRITE(,) 'radiative transfer; please extend the '
469	! WRITE(,) 'interpolation grid to larger grain sizes.'
470	grid_i=refftabsize-1
471	kx = 1.
472	ELSEIF (grid_i.LT.1) THEN
473	! WRITE(,) 'Warning: particle size in grid box #'
474	! WRITE(,) ig,' is too small to be used by the '
475	! WRITE(,) 'radiative transfer; please extend the '
476	! WRITE(,) 'interpolation grid to smaller grain sizes.'
477	grid_i=1
478	kx = 0.
479	ELSE
480	kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) / &
481	( refftab(grid_i+1)-refftab(grid_i) )
482	ENDIF
483	! 2.3 Integration
484	DO j=grid_i,grid_i+1
485	! 2.3.1 Check if the calculation has been done
486	IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
487	! 2.3.2 Log-normal dist., r_g and sigma_g are defined
488	! in [hansen_1974], "Light scattering in planetary
489	! atmospheres", Space Science Reviews 16 527-610.
490	! Here, sizedistk1=r_g and sizedistk2=sigma_g^2
491	sizedistk2 = log(1.+nueffrad(1,1,iaer))
492	sizedistk1 = exp(2.5*sizedistk2)
493	sizedistk1 = refftab(j) / sizedistk1
494
495	normd(j,1,iaer,idomain) = 1e-30
496	DO gausind=1,ngau
497	drad=radiusr*radgaus(gausind)
498	dista(j,1,iaer,idomain,gausind) = LOG((radiusm-drad)/sizedistk1)
499	dista(j,1,iaer,idomain,gausind) = &
500	EXP(-dista(j,1,iaer,idomain,gausind) * &
501	dista(j,1,iaer,idomain,gausind) * &
502	0.5e0/sizedistk2)/(radiusm-drad)
503	dista(j,1,iaer,idomain,gausind) = &
504	dista(j,1,iaer,idomain,gausind) / &
505	(sqrt(2e0pisizedistk2))
506
507	distb(j,1,iaer,idomain,gausind) = LOG((radiusm+drad)/sizedistk1)
508	distb(j,1,iaer,idomain,gausind) = &
509	EXP(-distb(j,1,iaer,idomain,gausind) * &
510	distb(j,1,iaer,idomain,gausind) * &
511	0.5e0/sizedistk2)/(radiusm+drad)
512	distb(j,1,iaer,idomain,gausind) = &
513	distb(j,1,iaer,idomain,gausind) / &
514	(sqrt(2e0pisizedistk2))
515
516	normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) + &
517	weightgaus(gausind) * &
518	( &
519	distb(j,1,iaer,idomain,gausind) * pi * &
520	radGAUSb(gausind,iaer,idomain) * &
521	radGAUSb(gausind,iaer,idomain) + &
522	dista(j,1,iaer,idomain,gausind) * pi * &
523	radGAUSa(gausind,iaer,idomain) * &
524	radGAUSa(gausind,iaer,idomain) &
525	)
526	ENDDO
527	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
528	! 2.3.3.vis Initialization
529	qsqrefVISgrid(j,1,:,iaer)=0.
530	qextVISgrid(j,1,:,iaer)=0.
531	qscatVISgrid(j,1,:,iaer)=0.
532	omegVISgrid(j,1,:,iaer)=0.
533	gVISgrid(j,1,:,iaer)=0.
534	qrefVISgrid(j,1,iaer)=0.
535	qscatrefVISgrid(j,1,iaer)=0.
536	omegrefVISgrid(j,1,iaer)=0.
537
538	DO gausind=1,ngau
539	DO m=1,L_NSPECTV
540	! Convolution:
541	qextVISgrid(j,1,m,iaer) = &
542	qextVISgrid(j,1,m,iaer) + &
543	weightgaus(gausind) * &
544	( &
545	qsqrefVISb(m,gausind,iaer) * &
546	qrefVISb(gausind,iaer) * &
547	piradGAUSb(gausind,iaer,idomain) &
548	radGAUSb(gausind,iaer,idomain) * &
549	distb(j,1,iaer,idomain,gausind) + &
550	qsqrefVISa(m,gausind,iaer) * &
551	qrefVISa(gausind,iaer) * &
552	piradGAUSa(gausind,iaer,idomain) &
553	radGAUSa(gausind,iaer,idomain) * &
554	dista(j,1,iaer,idomain,gausind) &
555	)
556	qscatVISgrid(j,1,m,iaer) = &
557	qscatVISgrid(j,1,m,iaer) + &
558	weightgaus(gausind) * &
559	( &
560	omegVISb(m,gausind,iaer) * &
561	qsqrefVISb(m,gausind,iaer) * &
562	qrefVISb(gausind,iaer) * &
563	piradGAUSb(gausind,iaer,idomain) &
564	radGAUSb(gausind,iaer,idomain) * &
565	distb(j,1,iaer,idomain,gausind) + &
566	omegVISa(m,gausind,iaer) * &
567	qsqrefVISa(m,gausind,iaer) * &
568	qrefVISa(gausind,iaer) * &
569	piradGAUSa(gausind,iaer,idomain) &
570	radGAUSa(gausind,iaer,idomain) * &
571	dista(j,1,iaer,idomain,gausind) &
572	)
573	gVISgrid(j,1,m,iaer) = &
574	gVISgrid(j,1,m,iaer) + &
575	weightgaus(gausind) * &
576	( &
577	omegVISb(m,gausind,iaer) * &
578	qsqrefVISb(m,gausind,iaer) * &
579	qrefVISb(gausind,iaer) * &
580	gVISb(m,gausind,iaer) * &
581	piradGAUSb(gausind,iaer,idomain) &
582	radGAUSb(gausind,iaer,idomain) * &
583	distb(j,1,iaer,idomain,gausind) + &
584	omegVISa(m,gausind,iaer) * &
585	qsqrefVISa(m,gausind,iaer) * &
586	qrefVISa(gausind,iaer) * &
587	gVISa(m,gausind,iaer) * &
588	piradGAUSa(gausind,iaer,idomain) &
589	radGAUSa(gausind,iaer,idomain) * &
590	dista(j,1,iaer,idomain,gausind) &
591	)
592	ENDDO
593	qrefVISgrid(j,1,iaer) = &
594	qrefVISgrid(j,1,iaer) + &
595	weightgaus(gausind) * &
596	( &
597	qrefVISb(gausind,iaer) * &
598	piradGAUSb(gausind,iaer,idomain) &
599	radGAUSb(gausind,iaer,idomain) * &
600	distb(j,1,iaer,idomain,gausind) + &
601	qrefVISa(gausind,iaer) * &
602	piradGAUSa(gausind,iaer,idomain) &
603	radGAUSa(gausind,iaer,idomain) * &
604	dista(j,1,iaer,idomain,gausind) &
605	)
606	qscatrefVISgrid(j,1,iaer) = &
607	qscatrefVISgrid(j,1,iaer) + &
608	weightgaus(gausind) * &
609	( &
610	omegrefVISb(gausind,iaer) * &
611	qrefVISb(gausind,iaer) * &
612	piradGAUSb(gausind,iaer,idomain) &
613	radGAUSb(gausind,iaer,idomain) * &
614	distb(j,1,iaer,idomain,gausind) + &
615	omegrefVISa(gausind,iaer) * &
616	qrefVISa(gausind,iaer) * &
617	piradGAUSa(gausind,iaer,idomain) &
618	radGAUSa(gausind,iaer,idomain) * &
619	dista(j,1,iaer,idomain,gausind) &
620	)
621	ENDDO
622
623	qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) / &
624	normd(j,1,iaer,idomain)
625	qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / &
626	normd(j,1,iaer,idomain)
627	omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / &
628	qrefVISgrid(j,1,iaer)
629	DO m=1,L_NSPECTV
630	qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / &
631	normd(j,1,iaer,idomain)
632	qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / &
633	normd(j,1,iaer,idomain)
634	gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) / &
635	qscatVISgrid(j,1,m,iaer) / &
636	normd(j,1,iaer,idomain)
637
638	qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / &
639	qrefVISgrid(j,1,iaer)
640	omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / &
641	qextVISgrid(j,1,m,iaer)
642	ENDDO
643	ELSE ! INFRARED DOMAIN ----------
644	! 2.3.3.ir Initialization
645	qsqrefIRgrid(j,1,:,iaer)=0.
646	qextIRgrid(j,1,:,iaer)=0.
647	qscatIRgrid(j,1,:,iaer)=0.
648	omegIRgrid(j,1,:,iaer)=0.
649	gIRgrid(j,1,:,iaer)=0.
650	qrefIRgrid(j,1,iaer)=0.
651	qscatrefIRgrid(j,1,iaer)=0.
652	omegrefIRgrid(j,1,iaer)=0.
653
654	DO gausind=1,ngau
655	DO m=1,L_NSPECTI
656	! Convolution:
657	qextIRgrid(j,1,m,iaer) = &
658	qextIRgrid(j,1,m,iaer) + &
659	weightgaus(gausind) * &
660	( &
661	qsqrefIRb(m,gausind,iaer) * &
662	qrefVISb(gausind,iaer) * &
663	piradGAUSb(gausind,iaer,idomain) &
664	radGAUSb(gausind,iaer,idomain) * &
665	distb(j,1,iaer,idomain,gausind) + &
666	qsqrefIRa(m,gausind,iaer) * &
667	qrefVISa(gausind,iaer) * &
668	piradGAUSa(gausind,iaer,idomain) &
669	radGAUSa(gausind,iaer,idomain) * &
670	dista(j,1,iaer,idomain,gausind) &
671	)
672	qscatIRgrid(j,1,m,iaer) = &
673	qscatIRgrid(j,1,m,iaer) + &
674	weightgaus(gausind) * &
675	( &
676	omegIRb(m,gausind,iaer) * &
677	qsqrefIRb(m,gausind,iaer) * &
678	qrefVISb(gausind,iaer) * &
679	piradGAUSb(gausind,iaer,idomain) &
680	radGAUSb(gausind,iaer,idomain) * &
681	distb(j,1,iaer,idomain,gausind) + &
682	omegIRa(m,gausind,iaer) * &
683	qsqrefIRa(m,gausind,iaer) * &
684	qrefVISa(gausind,iaer) * &
685	piradGAUSa(gausind,iaer,idomain) &
686	radGAUSa(gausind,iaer,idomain) * &
687	dista(j,1,iaer,idomain,gausind) &
688	)
689	gIRgrid(j,1,m,iaer) = &
690	gIRgrid(j,1,m,iaer) + &
691	weightgaus(gausind) * &
692	( &
693	omegIRb(m,gausind,iaer) * &
694	qsqrefIRb(m,gausind,iaer) * &
695	qrefVISb(gausind,iaer) * &
696	gIRb(m,gausind,iaer) * &
697	piradGAUSb(gausind,iaer,idomain) &
698	radGAUSb(gausind,iaer,idomain) * &
699	distb(j,1,iaer,idomain,gausind) + &
700	omegIRa(m,gausind,iaer) * &
701	qsqrefIRa(m,gausind,iaer) * &
702	qrefVISa(gausind,iaer) * &
703	gIRa(m,gausind,iaer) * &
704	piradGAUSa(gausind,iaer,idomain) &
705	radGAUSa(gausind,iaer,idomain) * &
706	dista(j,1,iaer,idomain,gausind) &
707	)
708	ENDDO
709	qrefIRgrid(j,1,iaer) = &
710	qrefIRgrid(j,1,iaer) + &
711	weightgaus(gausind) * &
712	( &
713	qrefIRb(gausind,iaer) * &
714	piradGAUSb(gausind,iaer,idomain) &
715	radGAUSb(gausind,iaer,idomain) * &
716	distb(j,1,iaer,idomain,gausind) + &
717	qrefIRa(gausind,iaer) * &
718	piradGAUSa(gausind,iaer,idomain) &
719	radGAUSa(gausind,iaer,idomain) * &
720	dista(j,1,iaer,idomain,gausind) &
721	)
722	qscatrefIRgrid(j,1,iaer) = &
723	qscatrefIRgrid(j,1,iaer) + &
724	weightgaus(gausind) * &
725	( &
726	omegrefIRb(gausind,iaer) * &
727	qrefIRb(gausind,iaer) * &
728	piradGAUSb(gausind,iaer,idomain) &
729	radGAUSb(gausind,iaer,idomain) * &
730	distb(j,1,iaer,idomain,gausind) + &
731	omegrefIRa(gausind,iaer) * &
732	qrefIRa(gausind,iaer) * &
733	piradGAUSa(gausind,iaer,idomain) &
734	radGAUSa(gausind,iaer,idomain) * &
735	dista(j,1,iaer,idomain,gausind) &
736	)
737	ENDDO
738
739	qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) / &
740	normd(j,1,iaer,idomain)
741	qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / &
742	normd(j,1,iaer,idomain)
743	omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / &
744	qrefIRgrid(j,1,iaer)
745	DO m=1,L_NSPECTI
746	qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / &
747	normd(j,1,iaer,idomain)
748	qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / &
749	normd(j,1,iaer,idomain)
750	gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) / &
751	qscatIRgrid(j,1,m,iaer) / &
752	normd(j,1,iaer,idomain)
753
754	qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / &
755	qrefVISgrid(j,1,iaer)
756	omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / &
757	qextIRgrid(j,1,m,iaer)
758	ENDDO
759	ENDIF ! --------------------------
760	checkgrid(j,1,iaer,idomain) = .true.
761	ENDIF !checkgrid
762	ENDDO !grid_i
763	! 2.4 Linear interpolation
764	k1 = (1-kx)
765	k2 = kx
766	IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
767	DO m=1,L_NSPECTV
768	QVISsQREF3d(ig,lg,m,iaer) = &
769	k1*qsqrefVISgrid(grid_i,1,m,iaer) + &
770	k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
771	omegaVIS3d(ig,lg,m,iaer) = &
772	k1*omegVISgrid(grid_i,1,m,iaer) + &
773	k2*omegVISgrid(grid_i+1,1,m,iaer)
774	gVIS3d(ig,lg,m,iaer) = &
775	k1*gVISgrid(grid_i,1,m,iaer) + &
776	k2*gVISgrid(grid_i+1,1,m,iaer)
777	ENDDO !L_NSPECTV
778	QREFvis3d(ig,lg,iaer) = &
779	k1*qrefVISgrid(grid_i,1,iaer) + &
780	k2*qrefVISgrid(grid_i+1,1,iaer)
781	! omegaREFvis3d(ig,lg,iaer) = &
782	! k1*omegrefVISgrid(grid_i,1,iaer) + &
783	! k2*omegrefVISgrid(grid_i+1,1,iaer)
784	ELSE ! INFRARED -----------------------
785	DO m=1,L_NSPECTI
786	QIRsQREF3d(ig,lg,m,iaer) = &
787	k1*qsqrefIRgrid(grid_i,1,m,iaer) + &
788	k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
789	omegaIR3d(ig,lg,m,iaer) = &
790	k1*omegIRgrid(grid_i,1,m,iaer) + &
791	k2*omegIRgrid(grid_i+1,1,m,iaer)
792	gIR3d(ig,lg,m,iaer) = &
793	k1*gIRgrid(grid_i,1,m,iaer) + &
794	k2*gIRgrid(grid_i+1,1,m,iaer)
795	ENDDO !L_NSPECTI
796	QREFir3d(ig,lg,iaer) = &
797	k1*qrefIRgrid(grid_i,1,iaer) + &
798	k2*qrefIRgrid(grid_i+1,1,iaer)
799	! omegaREFir3d(ig,lg,iaer) = &
800	! k1*omegrefIRgrid(grid_i,1,iaer) + &
801	! k2*omegrefIRgrid(grid_i+1,1,iaer)
802	ENDIF ! --------------------------------
803	ENDDO !nlayer
804	ENDDO !ngrid
805
806	!==================================================================
807
808
809
810	ENDDO ! idomain
811
812	ENDIF ! nsize = 1
813
814	ENDDO ! iaer (loop on aerosol kind)
815
816	END SUBROUTINE aeroptproperties
817
818
819	end module aeroptproperties_mod

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

Download in other formats:

Original Format