Context Navigation

aeroptproperties.F90 @ 1600

Last change on this file since 1600 was 1529, checked in by emillour, 9 years ago

Generic GCM: Towards a cleaner separation between physics and dynamics

Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
Removed module "comhdiff_mod.F90", as it is only used by module surf_heat_transp_mod.F90, moved module variables there.
Addedin "surf_heat_transp_mod" local versions of some arrays and routines (from dyn3d) required to compute gradient, divergence, etc. on the global dynamics grid. As before, the slab ocean only works in serial.

File size: 36.4 KB

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

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

Download in other formats:

Original Format