Context Navigation

aeroptproperties.F90 @ 3595

Last change on this file since 3595 was 3585, checked in by debatzbr, 2 weeks ago
Connecting microphysics to radiative transfer + miscellaneous cleans
File size: 40.4 KB

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

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

Download in other formats:

Original Format