Context Navigation

aeroptproperties.F90 @ 2613

Last change on this file since 2613 was 2033, checked in by emillour, 6 years ago

Generic GCM:

Cleanup around aeroptproperties.90, remove arrays that were not initialized but still used, probably for nothing since removing them leads to no change in results. Possibly to revisit and further clean up later.

File size: 36.5 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
[2033]	9	use radcommon_h, only: qrefvis,qrefir,omegarefir !,omegarefvis
[253]	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
[2033]	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)
[2033]	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)
[2033]	298	omegrefVISa(gausind,iaer)= 0
	299	! omegrefVISa(gausind,iaer)= &
	300	! (1-kint)*omegaREFvis(iaer,radius_id) + &
	301	! kint*omegaREFvis(iaer,radius_id+1)
[253]	302	ELSE ! INFRARED DOMAIN ----------------------------------
	303	DO m=1,L_NSPECTI
	304	qsqrefIRa(m,gausind,iaer)= &
	305	(1-kint)*QIRsQREF(m,iaer,radius_id) + &
	306	kint*QIRsQREF(m,iaer,radius_id+1)
	307	omegIRa(m,gausind,iaer)= &
	308	(1-kint)*omegaIR(m,iaer,radius_id) + &
	309	kint*omegaIR(m,iaer,radius_id+1)
	310	gIRa(m,gausind,iaer)= &
	311	(1-kint)*gIR(m,iaer,radius_id) + &
	312	kint*gIR(m,iaer,radius_id+1)
	313	ENDDO
	314	qrefIRa(gausind,iaer)= &
	315	(1-kint)*QREFir(iaer,radius_id) + &
	316	kint*QREFir(iaer,radius_id+1)
	317	omegrefIRa(gausind,iaer)= &
	318	(1-kint)*omegaREFir(iaer,radius_id) + &
	319	kint*omegaREFir(iaer,radius_id+1)
	320	ENDIF
	321	ENDDO
	322
	323	DO gausind=1,ngau
	324	drad=radiusr*radgaus(gausind)
	325	radGAUSb(gausind,iaer,idomain)=radiusm+drad
	326
	327	radius_id=minloc(abs(radiustab(iaer,idomain,:) - &
	328	(radiusm+drad)),DIM=1)
	329	IF ((radiustab(iaer,idomain,radius_id) - &
	330	(radiusm+drad)).GT.0) THEN
	331	radius_id=radius_id-1
	332	ENDIF
	333	IF (radius_id.GE.nsize(iaer,idomain)) THEN
	334	radius_id=nsize(iaer,idomain)-1
	335	kint = 1.
	336	ELSEIF (radius_id.LT.1) THEN
	337	radius_id=1
	338	kint = 0.
	339	ELSE
	340	kint = ( (radiusm+drad) - &
	341	radiustab(iaer,idomain,radius_id) ) / &
	342	( radiustab(iaer,idomain,radius_id+1) - &
	343	radiustab(iaer,idomain,radius_id) )
	344	ENDIF
	345	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
	346	DO m=1,L_NSPECTV
	347	qsqrefVISb(m,gausind,iaer)= &
	348	(1-kint)*QVISsQREF(m,iaer,radius_id) + &
	349	kint*QVISsQREF(m,iaer,radius_id+1)
	350	omegVISb(m,gausind,iaer)= &
	351	(1-kint)*omegaVIS(m,iaer,radius_id) + &
	352	kint*omegaVIS(m,iaer,radius_id+1)
	353	gVISb(m,gausind,iaer)= &
	354	(1-kint)*gVIS(m,iaer,radius_id) + &
	355	kint*gVIS(m,iaer,radius_id+1)
	356	ENDDO
	357	qrefVISb(gausind,iaer)= &
	358	(1-kint)*QREFvis(iaer,radius_id) + &
	359	kint*QREFvis(iaer,radius_id+1)
[2033]	360	omegrefVISb(gausind,iaer)= 0
	361	! omegrefVISb(gausind,iaer)= &
	362	! (1-kint)*omegaREFvis(iaer,radius_id) + &
	363	! kint*omegaREFvis(iaer,radius_id+1)
[253]	364	ELSE ! INFRARED DOMAIN ----------------------------------
	365	DO m=1,L_NSPECTI
	366	qsqrefIRb(m,gausind,iaer)= &
	367	(1-kint)*QIRsQREF(m,iaer,radius_id) + &
	368	kint*QIRsQREF(m,iaer,radius_id+1)
	369	omegIRb(m,gausind,iaer)= &
	370	(1-kint)*omegaIR(m,iaer,radius_id) + &
	371	kint*omegaIR(m,iaer,radius_id+1)
	372	gIRb(m,gausind,iaer)= &
	373	(1-kint)*gIR(m,iaer,radius_id) + &
	374	kint*gIR(m,iaer,radius_id+1)
	375	ENDDO
	376	qrefIRb(gausind,iaer)= &
	377	(1-kint)*QREFir(iaer,radius_id) + &
	378	kint*QREFir(iaer,radius_id+1)
	379	omegrefIRb(gausind,iaer)= &
	380	(1-kint)*omegaREFir(iaer,radius_id) + &
	381	kint*omegaREFir(iaer,radius_id+1)
	382	ENDIF
	383	ENDDO
	384
	385	!==================================================================
	386	! CONSTANT NUEFF FROM HERE
	387	!==================================================================
	388
[135]	389	! 2. Compute the scattering parameters using linear
	390	! interpolation over grain sizes and constant nueff
	391	! ---------------------------------------------------
	392
	393	DO lg = 1,nlayer
[787]	394	DO ig = 1, ngrid
[135]	395	! 2.1 Effective radius index and kx calculation
	396	var_tmp=reffrad(ig,lg,iaer)/refftabmin
	397	var_tmp=log(var_tmp)*3.
	398	var_tmp=var_tmp/logvratgrid+1.
	399	grid_i=floor(var_tmp)
	400	IF (grid_i.GE.refftabsize) THEN
[253]	401	! WRITE(,) 'Warning: particle size in grid box #'
	402	! WRITE(,) ig,' is too large to be used by the '
	403	! WRITE(,) 'radiative transfer; please extend the '
	404	! WRITE(,) 'interpolation grid to larger grain sizes.'
[135]	405	grid_i=refftabsize-1
	406	kx = 1.
	407	ELSEIF (grid_i.LT.1) THEN
[253]	408	! WRITE(,) 'Warning: particle size in grid box #'
	409	! WRITE(,) ig,' is too small to be used by the '
	410	! WRITE(,) 'radiative transfer; please extend the '
	411	! WRITE(,) 'interpolation grid to smaller grain sizes.'
[135]	412	grid_i=1
	413	kx = 0.
	414	ELSE
[253]	415	kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) / &
	416	( refftab(grid_i+1)-refftab(grid_i) )
[135]	417	ENDIF
	418	! 2.3 Integration
	419	DO j=grid_i,grid_i+1
[253]	420	! 2.3.1 Check if the calculation has been done
[135]	421	IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
[253]	422	! 2.3.2 Log-normal dist., r_g and sigma_g are defined
	423	! in [hansen_1974], "Light scattering in planetary
	424	! atmospheres", Space Science Reviews 16 527-610.
	425	! Here, sizedistk1=r_g and sizedistk2=sigma_g^2
	426	sizedistk2 = log(1.+nueffrad(1,1,iaer))
	427	sizedistk1 = exp(2.5*sizedistk2)
	428	sizedistk1 = refftab(j) / sizedistk1
	429
	430	normd(j,1,iaer,idomain) = 1e-30
	431	DO gausind=1,ngau
	432	drad=radiusr*radgaus(gausind)
	433	dista(j,1,iaer,idomain,gausind) = LOG((radiusm-drad)/sizedistk1)
	434	dista(j,1,iaer,idomain,gausind) = &
	435	EXP(-dista(j,1,iaer,idomain,gausind) * &
	436	dista(j,1,iaer,idomain,gausind) * &
	437	0.5e0/sizedistk2)/(radiusm-drad)
	438	dista(j,1,iaer,idomain,gausind) = &
	439	dista(j,1,iaer,idomain,gausind) / &
	440	(sqrt(2e0pisizedistk2))
	441
	442	distb(j,1,iaer,idomain,gausind) = LOG((radiusm+drad)/sizedistk1)
	443	distb(j,1,iaer,idomain,gausind) = &
	444	EXP(-distb(j,1,iaer,idomain,gausind) * &
	445	distb(j,1,iaer,idomain,gausind) * &
	446	0.5e0/sizedistk2)/(radiusm+drad)
	447	distb(j,1,iaer,idomain,gausind) = &
	448	distb(j,1,iaer,idomain,gausind) / &
	449	(sqrt(2e0pisizedistk2))
	450
	451	normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) + &
	452	weightgaus(gausind) * &
	453	( &
	454	distb(j,1,iaer,idomain,gausind) * pi * &
	455	radGAUSb(gausind,iaer,idomain) * &
	456	radGAUSb(gausind,iaer,idomain) + &
	457	dista(j,1,iaer,idomain,gausind) * pi * &
	458	radGAUSa(gausind,iaer,idomain) * &
	459	radGAUSa(gausind,iaer,idomain) &
	460	)
	461	ENDDO
[135]	462	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
	463	! 2.3.3.vis Initialization
[253]	464	qsqrefVISgrid(j,1,:,iaer)=0.
	465	qextVISgrid(j,1,:,iaer)=0.
	466	qscatVISgrid(j,1,:,iaer)=0.
[135]	467	omegVISgrid(j,1,:,iaer)=0.
	468	gVISgrid(j,1,:,iaer)=0.
[253]	469	qrefVISgrid(j,1,iaer)=0.
	470	qscatrefVISgrid(j,1,iaer)=0.
	471	omegrefVISgrid(j,1,iaer)=0.
	472
	473	DO gausind=1,ngau
	474	DO m=1,L_NSPECTV
	475	! Convolution:
	476	qextVISgrid(j,1,m,iaer) = &
	477	qextVISgrid(j,1,m,iaer) + &
	478	weightgaus(gausind) * &
	479	( &
	480	qsqrefVISb(m,gausind,iaer) * &
	481	qrefVISb(gausind,iaer) * &
	482	piradGAUSb(gausind,iaer,idomain) &
	483	radGAUSb(gausind,iaer,idomain) * &
	484	distb(j,1,iaer,idomain,gausind) + &
	485	qsqrefVISa(m,gausind,iaer) * &
	486	qrefVISa(gausind,iaer) * &
	487	piradGAUSa(gausind,iaer,idomain) &
	488	radGAUSa(gausind,iaer,idomain) * &
	489	dista(j,1,iaer,idomain,gausind) &
	490	)
	491	qscatVISgrid(j,1,m,iaer) = &
	492	qscatVISgrid(j,1,m,iaer) + &
	493	weightgaus(gausind) * &
	494	( &
	495	omegVISb(m,gausind,iaer) * &
	496	qsqrefVISb(m,gausind,iaer) * &
	497	qrefVISb(gausind,iaer) * &
	498	piradGAUSb(gausind,iaer,idomain) &
	499	radGAUSb(gausind,iaer,idomain) * &
	500	distb(j,1,iaer,idomain,gausind) + &
	501	omegVISa(m,gausind,iaer) * &
	502	qsqrefVISa(m,gausind,iaer) * &
	503	qrefVISa(gausind,iaer) * &
	504	piradGAUSa(gausind,iaer,idomain) &
	505	radGAUSa(gausind,iaer,idomain) * &
	506	dista(j,1,iaer,idomain,gausind) &
	507	)
	508	gVISgrid(j,1,m,iaer) = &
	509	gVISgrid(j,1,m,iaer) + &
	510	weightgaus(gausind) * &
	511	( &
	512	omegVISb(m,gausind,iaer) * &
	513	qsqrefVISb(m,gausind,iaer) * &
	514	qrefVISb(gausind,iaer) * &
	515	gVISb(m,gausind,iaer) * &
	516	piradGAUSb(gausind,iaer,idomain) &
	517	radGAUSb(gausind,iaer,idomain) * &
	518	distb(j,1,iaer,idomain,gausind) + &
	519	omegVISa(m,gausind,iaer) * &
	520	qsqrefVISa(m,gausind,iaer) * &
	521	qrefVISa(gausind,iaer) * &
	522	gVISa(m,gausind,iaer) * &
	523	piradGAUSa(gausind,iaer,idomain) &
	524	radGAUSa(gausind,iaer,idomain) * &
	525	dista(j,1,iaer,idomain,gausind) &
	526	)
	527	ENDDO
	528	qrefVISgrid(j,1,iaer) = &
	529	qrefVISgrid(j,1,iaer) + &
	530	weightgaus(gausind) * &
	531	( &
	532	qrefVISb(gausind,iaer) * &
	533	piradGAUSb(gausind,iaer,idomain) &
	534	radGAUSb(gausind,iaer,idomain) * &
	535	distb(j,1,iaer,idomain,gausind) + &
	536	qrefVISa(gausind,iaer) * &
	537	piradGAUSa(gausind,iaer,idomain) &
	538	radGAUSa(gausind,iaer,idomain) * &
	539	dista(j,1,iaer,idomain,gausind) &
	540	)
	541	qscatrefVISgrid(j,1,iaer) = &
	542	qscatrefVISgrid(j,1,iaer) + &
	543	weightgaus(gausind) * &
	544	( &
	545	omegrefVISb(gausind,iaer) * &
	546	qrefVISb(gausind,iaer) * &
	547	piradGAUSb(gausind,iaer,idomain) &
	548	radGAUSb(gausind,iaer,idomain) * &
	549	distb(j,1,iaer,idomain,gausind) + &
	550	omegrefVISa(gausind,iaer) * &
	551	qrefVISa(gausind,iaer) * &
	552	piradGAUSa(gausind,iaer,idomain) &
	553	radGAUSa(gausind,iaer,idomain) * &
	554	dista(j,1,iaer,idomain,gausind) &
	555	)
[135]	556	ENDDO
[253]	557
	558	qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) / &
	559	normd(j,1,iaer,idomain)
	560	qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / &
	561	normd(j,1,iaer,idomain)
	562	omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / &
	563	qrefVISgrid(j,1,iaer)
	564	DO m=1,L_NSPECTV
	565	qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / &
	566	normd(j,1,iaer,idomain)
	567	qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / &
	568	normd(j,1,iaer,idomain)
	569	gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) / &
	570	qscatVISgrid(j,1,m,iaer) / &
	571	normd(j,1,iaer,idomain)
	572
	573	qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / &
	574	qrefVISgrid(j,1,iaer)
	575	omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / &
	576	qextVISgrid(j,1,m,iaer)
	577	ENDDO
[135]	578	ELSE ! INFRARED DOMAIN ----------
	579	! 2.3.3.ir Initialization
[253]	580	qsqrefIRgrid(j,1,:,iaer)=0.
	581	qextIRgrid(j,1,:,iaer)=0.
	582	qscatIRgrid(j,1,:,iaer)=0.
[135]	583	omegIRgrid(j,1,:,iaer)=0.
	584	gIRgrid(j,1,:,iaer)=0.
[253]	585	qrefIRgrid(j,1,iaer)=0.
	586	qscatrefIRgrid(j,1,iaer)=0.
	587	omegrefIRgrid(j,1,iaer)=0.
	588
	589	DO gausind=1,ngau
	590	DO m=1,L_NSPECTI
	591	! Convolution:
	592	qextIRgrid(j,1,m,iaer) = &
	593	qextIRgrid(j,1,m,iaer) + &
	594	weightgaus(gausind) * &
	595	( &
	596	qsqrefIRb(m,gausind,iaer) * &
	597	qrefVISb(gausind,iaer) * &
	598	piradGAUSb(gausind,iaer,idomain) &
	599	radGAUSb(gausind,iaer,idomain) * &
	600	distb(j,1,iaer,idomain,gausind) + &
	601	qsqrefIRa(m,gausind,iaer) * &
	602	qrefVISa(gausind,iaer) * &
	603	piradGAUSa(gausind,iaer,idomain) &
	604	radGAUSa(gausind,iaer,idomain) * &
	605	dista(j,1,iaer,idomain,gausind) &
	606	)
	607	qscatIRgrid(j,1,m,iaer) = &
	608	qscatIRgrid(j,1,m,iaer) + &
	609	weightgaus(gausind) * &
	610	( &
	611	omegIRb(m,gausind,iaer) * &
	612	qsqrefIRb(m,gausind,iaer) * &
	613	qrefVISb(gausind,iaer) * &
	614	piradGAUSb(gausind,iaer,idomain) &
	615	radGAUSb(gausind,iaer,idomain) * &
	616	distb(j,1,iaer,idomain,gausind) + &
	617	omegIRa(m,gausind,iaer) * &
	618	qsqrefIRa(m,gausind,iaer) * &
	619	qrefVISa(gausind,iaer) * &
	620	piradGAUSa(gausind,iaer,idomain) &
	621	radGAUSa(gausind,iaer,idomain) * &
	622	dista(j,1,iaer,idomain,gausind) &
	623	)
	624	gIRgrid(j,1,m,iaer) = &
	625	gIRgrid(j,1,m,iaer) + &
	626	weightgaus(gausind) * &
	627	( &
	628	omegIRb(m,gausind,iaer) * &
	629	qsqrefIRb(m,gausind,iaer) * &
	630	qrefVISb(gausind,iaer) * &
	631	gIRb(m,gausind,iaer) * &
	632	piradGAUSb(gausind,iaer,idomain) &
	633	radGAUSb(gausind,iaer,idomain) * &
	634	distb(j,1,iaer,idomain,gausind) + &
	635	omegIRa(m,gausind,iaer) * &
	636	qsqrefIRa(m,gausind,iaer) * &
	637	qrefVISa(gausind,iaer) * &
	638	gIRa(m,gausind,iaer) * &
	639	piradGAUSa(gausind,iaer,idomain) &
	640	radGAUSa(gausind,iaer,idomain) * &
	641	dista(j,1,iaer,idomain,gausind) &
	642	)
	643	ENDDO
	644	qrefIRgrid(j,1,iaer) = &
	645	qrefIRgrid(j,1,iaer) + &
	646	weightgaus(gausind) * &
	647	( &
	648	qrefIRb(gausind,iaer) * &
	649	piradGAUSb(gausind,iaer,idomain) &
	650	radGAUSb(gausind,iaer,idomain) * &
	651	distb(j,1,iaer,idomain,gausind) + &
	652	qrefIRa(gausind,iaer) * &
	653	piradGAUSa(gausind,iaer,idomain) &
	654	radGAUSa(gausind,iaer,idomain) * &
	655	dista(j,1,iaer,idomain,gausind) &
	656	)
	657	qscatrefIRgrid(j,1,iaer) = &
	658	qscatrefIRgrid(j,1,iaer) + &
	659	weightgaus(gausind) * &
	660	( &
	661	omegrefIRb(gausind,iaer) * &
	662	qrefIRb(gausind,iaer) * &
	663	piradGAUSb(gausind,iaer,idomain) &
	664	radGAUSb(gausind,iaer,idomain) * &
	665	distb(j,1,iaer,idomain,gausind) + &
	666	omegrefIRa(gausind,iaer) * &
	667	qrefIRa(gausind,iaer) * &
	668	piradGAUSa(gausind,iaer,idomain) &
	669	radGAUSa(gausind,iaer,idomain) * &
	670	dista(j,1,iaer,idomain,gausind) &
	671	)
[135]	672	ENDDO
[253]	673
	674	qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) / &
	675	normd(j,1,iaer,idomain)
	676	qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / &
	677	normd(j,1,iaer,idomain)
	678	omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / &
	679	qrefIRgrid(j,1,iaer)
	680	DO m=1,L_NSPECTI
	681	qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / &
	682	normd(j,1,iaer,idomain)
	683	qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / &
	684	normd(j,1,iaer,idomain)
	685	gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) / &
	686	qscatIRgrid(j,1,m,iaer) / &
	687	normd(j,1,iaer,idomain)
	688
	689	qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / &
	690	qrefVISgrid(j,1,iaer)
	691	omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / &
	692	qextIRgrid(j,1,m,iaer)
	693	ENDDO
[135]	694	ENDIF ! --------------------------
	695	checkgrid(j,1,iaer,idomain) = .true.
	696	ENDIF !checkgrid
	697	ENDDO !grid_i
	698	! 2.4 Linear interpolation
	699	k1 = (1-kx)
	700	k2 = kx
	701	IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
	702	DO m=1,L_NSPECTV
[253]	703	QVISsQREF3d(ig,lg,m,iaer) = &
	704	k1*qsqrefVISgrid(grid_i,1,m,iaer) + &
	705	k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
	706	omegaVIS3d(ig,lg,m,iaer) = &
	707	k1*omegVISgrid(grid_i,1,m,iaer) + &
	708	k2*omegVISgrid(grid_i+1,1,m,iaer)
	709	gVIS3d(ig,lg,m,iaer) = &
	710	k1*gVISgrid(grid_i,1,m,iaer) + &
	711	k2*gVISgrid(grid_i+1,1,m,iaer)
[135]	712	ENDDO !L_NSPECTV
[253]	713	QREFvis3d(ig,lg,iaer) = &
	714	k1*qrefVISgrid(grid_i,1,iaer) + &
	715	k2*qrefVISgrid(grid_i+1,1,iaer)
[2033]	716	! omegaREFvis3d(ig,lg,iaer) = &
	717	! k1*omegrefVISgrid(grid_i,1,iaer) + &
	718	! k2*omegrefVISgrid(grid_i+1,1,iaer)
[135]	719	ELSE ! INFRARED -----------------------
	720	DO m=1,L_NSPECTI
[253]	721	QIRsQREF3d(ig,lg,m,iaer) = &
	722	k1*qsqrefIRgrid(grid_i,1,m,iaer) + &
	723	k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
	724	omegaIR3d(ig,lg,m,iaer) = &
	725	k1*omegIRgrid(grid_i,1,m,iaer) + &
	726	k2*omegIRgrid(grid_i+1,1,m,iaer)
	727	gIR3d(ig,lg,m,iaer) = &
	728	k1*gIRgrid(grid_i,1,m,iaer) + &
	729	k2*gIRgrid(grid_i+1,1,m,iaer)
[135]	730	ENDDO !L_NSPECTI
[253]	731	QREFir3d(ig,lg,iaer) = &
	732	k1*qrefIRgrid(grid_i,1,iaer) + &
	733	k2*qrefIRgrid(grid_i+1,1,iaer)
[2033]	734	! omegaREFir3d(ig,lg,iaer) = &
	735	! k1*omegrefIRgrid(grid_i,1,iaer) + &
	736	! k2*omegrefIRgrid(grid_i+1,1,iaer)
[135]	737	ENDIF ! --------------------------------
[1308]	738	ENDDO !nlayer
[787]	739	ENDDO !ngrid
[135]	740
[253]	741	!==================================================================
	742
	743
	744
[135]	745	ENDDO ! idomain
	746
[253]	747	ENDIF ! nsize = 1
[135]	748
	749	ENDDO ! iaer (loop on aerosol kind)
	750
[253]	751	RETURN
	752	END SUBROUTINE aeroptproperties
[135]	753
[253]	754
	755
	756

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

Download in other formats:

Original Format