Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/aeroptproperties.f90 @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago

File size: 41.3 KB

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

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

Download in other formats: