Context Navigation

aeroptproperties.F90 @ 3436

Last change on this file since 3436 was 3353, checked in by afalco, 6 months ago
Pluto PCM: Added zrecast & old sedim ; Choose haze file ; AF
File size: 38.2 KB

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

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

Download in other formats:

Original Format