Context Navigation

aeroptproperties.F90 @ 1145

Last change on this file since 1145 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
RUN RCM1D : perfect match in normal mode.
RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

Additional arguments to some subroutines (ngrid and nq)
F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this

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

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

Download in other formats:

Original Format