Context Navigation

aeroptproperties.F90 @ 3562

Last change on this file since 3562 was 2972, checked in by emillour, 20 months ago
Generic PCM: Make number of scatterers fully dynamic (i.e. set in callphys.def and no longer by compilation option "-s #"). One should now specify naerkind = # in callphys.def (default is 0). EM
File size: 39.2 KB

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

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

Download in other formats:

Original Format