Context Navigation

aeroptproperties.F @ 804

Last change on this file since 804 was 784, checked in by jbmadeleine, 12 years ago
Minor change: added optional outputs for the 1D model in aeroptproperties.F.
File size: 55.2 KB

Rev	Line
[38]	1	SUBROUTINE aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
	2	& QVISsQREF3d,omegaVIS3d,gVIS3d,
	3	& QIRsQREF3d,omegaIR3d,gIR3d,
	4	& QREFvis3d,QREFir3d,
	5	& omegaREFvis3d,omegaREFir3d)
	6	IMPLICIT NONE
	7	c =============================================================
	8	c Aerosol Optical Properties
	9	c
	10	c Description:
	11	c Compute the scattering parameters in each grid
	12	c box, depending on aerosol grain sizes. Log-normal size
	13	c distribution and Gauss-Legendre integration are used.
	14
	15	c Parameters:
	16	c Don't forget to set the value of varyingnueff below; If
	17	c the effective variance of the distribution for the given
	18	c aerosol is considered homogeneous in the atmosphere, please
	19	c set varyingnueff(iaer) to .false. Resulting computational
	20	c time will be much better.
	21
	22	c Authors: J.-B. Madeleine, F. Forget, F. Montmessin
	23	c =============================================================
	24
	25	#include "dimensions.h"
	26	#include "dimphys.h"
	27	#include "callkeys.h"
	28	#include "dimradmars.h"
	29	#include "yomaer.h"
	30
	31	c Local variables
	32	c ---------------
	33
	34	c =============================================================
	35	LOGICAL, PARAMETER :: varyingnueff(naerkind) = .false.
	36	c =============================================================
	37
	38	c Min. and max radius of the interpolation grid (in METERS)
[737]	39	REAL, SAVE :: refftabmin(naerkind,2)
	40	REAL, SAVE :: refftabmax(naerkind,2)
[38]	41	c Log of the min and max variance of the interpolation grid
	42	REAL, PARAMETER :: nuefftabmin = -4.6
	43	REAL, PARAMETER :: nuefftabmax = 0.
	44	c Number of effective radius of the interpolation grid
	45	INTEGER, PARAMETER :: refftabsize = 100
	46	c Number of effective variances of the interpolation grid
	47	INTEGER, PARAMETER :: nuefftabsize = 100
	48	c Interpolation grid indices (reff,nueff)
	49	INTEGER :: grid_i,grid_j
	50	c Intermediate variable
	51	REAL :: var_tmp,var3d_tmp(ngridmx,nlayermx)
	52	c Bilinear interpolation factors
	53	REAL :: kx,ky,k1,k2,k3,k4
	54	c Size distribution parameters
	55	REAL :: sizedistk1,sizedistk2
	56	c Pi!
	57	REAL,SAVE :: pi
	58	c Variables used by the Gauss-Legendre integration:
	59	INTEGER radius_id,gausind
	60	REAL kint
	61	REAL drad
[89]	62	INTEGER, PARAMETER :: ngau = 10
	63	REAL weightgaus(ngau),radgaus(ngau)
[38]	64	SAVE weightgaus,radgaus
[89]	65	c DATA weightgaus/.2955242247,.2692667193,.2190863625,
	66	c & .1494513491,.0666713443/
	67	c DATA radgaus/.1488743389,.4333953941,.6794095682,
	68	c & .8650633666,.9739065285/
	69	DATA radgaus/0.07652652113350,0.22778585114165,
	70	& 0.37370608871528,0.51086700195146,
	71	& 0.63605368072468,0.74633190646476,
	72	& 0.83911697181213,0.91223442826796,
	73	& 0.96397192726078,0.99312859919241/
[38]	74
[89]	75	DATA weightgaus/0.15275338723120,0.14917298659407,
	76	& 0.14209610937519,0.13168863843930,
	77	& 0.11819453196154,0.10193011980823,
	78	& 0.08327674160932,0.06267204829828,
	79	& 0.04060142982019,0.01761400714091/
	80
[38]	81	c Indices
	82	INTEGER :: i,j,k,l,m,iaer,idomain
	83	INTEGER :: ig,lg,chg
	84
	85	c Local saved variables
	86	c ---------------------
	87
	88	c Radius axis of the interpolation grid
[737]	89	REAL,SAVE :: refftab(refftabsize,naerkind,2)
[38]	90	c Variance axis of the interpolation grid
[737]	91	REAL,SAVE :: nuefftab(nuefftabsize,naerkind,2)
[38]	92	c Volume ratio of the grid
[737]	93	REAL,SAVE :: logvratgrid(naerkind,2)
[38]	94	c Grid used to remember which calculation is done
	95	LOGICAL,SAVE :: checkgrid(refftabsize,nuefftabsize,naerkind,2)
	96	& = .false.
	97	c Optical properties of the grid (VISIBLE)
	98	REAL,SAVE :: qsqrefVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
	99	REAL,SAVE :: qextVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
	100	REAL,SAVE :: qscatVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
	101	REAL,SAVE :: omegVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
	102	REAL,SAVE :: gVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
	103	c Optical properties of the grid (INFRARED)
	104	REAL,SAVE :: qsqrefIRgrid(refftabsize,nuefftabsize,nir,naerkind)
	105	REAL,SAVE :: qextIRgrid(refftabsize,nuefftabsize,nir,naerkind)
	106	REAL,SAVE :: qscatIRgrid(refftabsize,nuefftabsize,nir,naerkind)
	107	REAL,SAVE :: omegIRgrid(refftabsize,nuefftabsize,nir,naerkind)
	108	REAL,SAVE :: gIRgrid(refftabsize,nuefftabsize,nir,naerkind)
	109	c Optical properties of the grid (REFERENCE WAVELENGTHS)
	110	REAL,SAVE :: qrefVISgrid(refftabsize,nuefftabsize,naerkind)
	111	REAL,SAVE :: qscatrefVISgrid(refftabsize,nuefftabsize,naerkind)
	112	REAL,SAVE :: qrefIRgrid(refftabsize,nuefftabsize,naerkind)
	113	REAL,SAVE :: qscatrefIRgrid(refftabsize,nuefftabsize,naerkind)
	114	REAL,SAVE :: omegrefVISgrid(refftabsize,nuefftabsize,naerkind)
	115	REAL,SAVE :: omegrefIRgrid(refftabsize,nuefftabsize,naerkind)
	116	c Firstcall
	117	LOGICAL,SAVE :: firstcall = .true.
	118	c Variables used by the Gauss-Legendre integration:
	119	REAL,SAVE :: normd(refftabsize,nuefftabsize,naerkind,2)
[89]	120	REAL,SAVE :: dista(refftabsize,nuefftabsize,naerkind,2,ngau)
	121	REAL,SAVE :: distb(refftabsize,nuefftabsize,naerkind,2,ngau)
[38]	122
[89]	123	REAL,SAVE :: radGAUSa(ngau,naerkind,2)
	124	REAL,SAVE :: radGAUSb(ngau,naerkind,2)
[38]	125
[89]	126	REAL,SAVE :: qsqrefVISa(nsun,ngau,naerkind)
	127	REAL,SAVE :: qrefVISa(ngau,naerkind)
	128	REAL,SAVE :: qsqrefVISb(nsun,ngau,naerkind)
	129	REAL,SAVE :: qrefVISb(ngau,naerkind)
	130	REAL,SAVE :: omegVISa(nsun,ngau,naerkind)
	131	REAL,SAVE :: omegrefVISa(ngau,naerkind)
	132	REAL,SAVE :: omegVISb(nsun,ngau,naerkind)
	133	REAL,SAVE :: omegrefVISb(ngau,naerkind)
	134	REAL,SAVE :: gVISa(nsun,ngau,naerkind)
	135	REAL,SAVE :: gVISb(nsun,ngau,naerkind)
[38]	136
[89]	137	REAL,SAVE :: qsqrefIRa(nir,ngau,naerkind)
	138	REAL,SAVE :: qrefIRa(ngau,naerkind)
	139	REAL,SAVE :: qsqrefIRb(nir,ngau,naerkind)
	140	REAL,SAVE :: qrefIRb(ngau,naerkind)
	141	REAL,SAVE :: omegIRa(nir,ngau,naerkind)
	142	REAL,SAVE :: omegrefIRa(ngau,naerkind)
	143	REAL,SAVE :: omegIRb(nir,ngau,naerkind)
	144	REAL,SAVE :: omegrefIRb(ngau,naerkind)
	145	REAL,SAVE :: gIRa(nir,ngau,naerkind)
	146	REAL,SAVE :: gIRb(nir,ngau,naerkind)
[38]	147
	148	REAL :: radiusm
	149	REAL :: radiusr
	150
	151	c Inputs
	152	c ------
	153
	154	INTEGER :: ngrid,nlayer
	155	c Aerosol effective radius used for radiative transfer (meter)
	156	REAL :: reffrad(ngridmx,nlayermx,naerkind)
	157	c Aerosol effective variance used for radiative transfer (n.u.)
	158	REAL :: nueffrad(ngridmx,nlayermx,naerkind)
	159
	160	c Outputs
	161	c -------
	162
	163	REAL :: QVISsQREF3d(ngridmx,nlayermx,nsun,naerkind)
	164	REAL :: omegaVIS3d(ngridmx,nlayermx,nsun,naerkind)
	165	REAL :: gVIS3d(ngridmx,nlayermx,nsun,naerkind)
	166
	167	REAL :: QIRsQREF3d(ngridmx,nlayermx,nir,naerkind)
	168	REAL :: omegaIR3d(ngridmx,nlayermx,nir,naerkind)
	169	REAL :: gIR3d(ngridmx,nlayermx,nir,naerkind)
	170
	171	REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
	172	REAL :: QREFir3d(ngridmx,nlayermx,naerkind)
	173
	174	REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
	175	REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind)
	176
[89]	177	c Tests
	178	c -----
	179
	180	LOGICAL,SAVE :: out_qwg = .false.
[784]	181	INTEGER, PARAMETER :: out_iaer = 2
	182	INTEGER :: out_ndim
[89]	183	REAL :: out_qext(ngridmx,nlayermx)
	184	REAL :: out_omeg(ngridmx,nlayermx)
	185	REAL :: out_g(ngridmx,nlayermx)
	186	INTEGER :: out_nchannel
	187	CHARACTER*1 :: out_str
	188
[737]	189	c Creating the effective radius and variance grid
	190	c -----------------------------------------------
	191	IF (firstcall) THEN
	192	c 0.1 Pi!
	193	pi = 2. * asin(1.e0)
	194
	195	WRITE(,) "aeroptproperties: interpolation grid"
	196	DO iaer = 1, naerkind ! Loop on aerosol kind
	197	DO idomain = 1, 2 ! Loop on visible or infrared channel
	198
	199	c 0.2 Effective radius
	200	radiusm=
	201	& 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain))+
	202	& radiustab(iaer,idomain,1))
	203	radiusr=
	204	& 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain))-
	205	& radiustab(iaer,idomain,1))
	206	refftabmin(iaer,idomain) =
	207	& radiusm-radiusr*radgaus(ngau)
	208	refftabmax(iaer,idomain) =
	209	& radiustab(iaer,idomain,nsize(iaer,idomain))
	210
	211	WRITE(,) "Scatterer: ",iaer
	212	WRITE(,) "Domain: ",idomain
	213	WRITE(,) "Min radius (m): ", refftabmin(iaer,idomain)
	214	WRITE(,) "Max radius (m): ", refftabmax(iaer,idomain)
	215
	216	refftab(1,iaer,idomain) =
	217	& refftabmin(iaer,idomain)
	218	refftab(refftabsize,iaer,idomain) =
	219	& refftabmax(iaer,idomain)
	220
	221	logvratgrid(iaer,idomain) =
	222	& log(refftabmax(iaer,idomain)/
	223	& refftabmin(iaer,idomain)) /
	224	& float(refftabsize-1)*3.
	225	do i = 2, refftabsize-1
	226	refftab(i,iaer,idomain) = refftab(i-1,iaer,idomain)*
	227	& exp(1./3.*logvratgrid(iaer,idomain))
	228	enddo
	229
	230	c 0.3 Effective variance
	231	do i = 0, nuefftabsize-1
	232	nuefftab(i+1,iaer,idomain) = exp( nuefftabmin +
	233	& i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
	234	enddo
	235
	236	ENDDO
	237	ENDDO
	238	firstcall = .false.
	239	ENDIF
	240
[38]	241	DO iaer = 1, naerkind ! Loop on aerosol kind
	242	IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN
	243	c==================================================================
	244	c If there is one single particle size, optical
	245	c properties of the considered aerosol are homogeneous
	246	DO lg = 1, nlayer
	247	DO ig = 1, ngrid
	248	DO chg = 1, nsun
	249	QVISsQREF3d(ig,lg,chg,iaer)=QVISsQREF(chg,iaer,1)
	250	omegaVIS3d(ig,lg,chg,iaer)=omegaVIS(chg,iaer,1)
	251	gVIS3d(ig,lg,chg,iaer)=gVIS(chg,iaer,1)
	252	ENDDO
	253	DO chg = 1, nir
	254	QIRsQREF3d(ig,lg,chg,iaer)=QIRsQREF(chg,iaer,1)
	255	omegaIR3d(ig,lg,chg,iaer)=omegaIR(chg,iaer,1)
	256	gIR3d(ig,lg,chg,iaer)=gIR(chg,iaer,1)
	257	ENDDO
	258	QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1)
	259	QREFir3d(ig,lg,iaer)=QREFir(iaer,1)
	260	omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1)
	261	omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1)
	262	ENDDO
	263	ENDDO
	264	ELSE ! Varying effective radius and variance
	265	DO idomain = 1, 2 ! Loop on visible or infrared channel
	266	c==================================================================
	267
[737]	268	c 1.1 Radius middle point and range for Gauss integration
[38]	269	radiusm=
	270	& 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) +
	271	& radiustab(iaer,idomain,1))
	272	radiusr=
	273	& 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) -
	274	& radiustab(iaer,idomain,1))
	275
[737]	276	c 1.2 Interpolating data at the Gauss quadrature points:
[89]	277	DO gausind=1,ngau
[38]	278	drad=radiusr*radgaus(gausind)
	279	radGAUSa(gausind,iaer,idomain)=radiusm-drad
	280
	281	radius_id=minloc(abs(radiustab(iaer,idomain,:) -
	282	& (radiusm-drad)),DIM=1)
	283	IF ((radiustab(iaer,idomain,radius_id) -
	284	& (radiusm-drad)).GT.0) THEN
	285	radius_id=radius_id-1
	286	ENDIF
	287	IF (radius_id.GE.nsize(iaer,idomain)) THEN
	288	radius_id=nsize(iaer,idomain)-1
	289	kint = 1.
	290	ELSEIF (radius_id.LT.1) THEN
	291	radius_id=1
	292	kint = 0.
	293	ELSE
	294	kint = ( (radiusm-drad) -
	295	& radiustab(iaer,idomain,radius_id) ) /
	296	& ( radiustab(iaer,idomain,radius_id+1) -
	297	& radiustab(iaer,idomain,radius_id) )
	298	ENDIF
	299	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
	300	DO m=1,nsun
	301	qsqrefVISa(m,gausind,iaer)=
	302	& (1-kint)*QVISsQREF(m,iaer,radius_id) +
	303	& kint*QVISsQREF(m,iaer,radius_id+1)
	304	omegVISa(m,gausind,iaer)=
	305	& (1-kint)*omegaVIS(m,iaer,radius_id) +
	306	& kint*omegaVIS(m,iaer,radius_id+1)
	307	gVISa(m,gausind,iaer)=
	308	& (1-kint)*gVIS(m,iaer,radius_id) +
	309	& kint*gVIS(m,iaer,radius_id+1)
	310	ENDDO
	311	qrefVISa(gausind,iaer)=
	312	& (1-kint)*QREFvis(iaer,radius_id) +
	313	& kint*QREFvis(iaer,radius_id+1)
	314	omegrefVISa(gausind,iaer)=
	315	& (1-kint)*omegaREFvis(iaer,radius_id) +
	316	& kint*omegaREFvis(iaer,radius_id+1)
	317	ELSE ! INFRARED DOMAIN ----------------------------------
	318	DO m=1,nir
	319	qsqrefIRa(m,gausind,iaer)=
	320	& (1-kint)*QIRsQREF(m,iaer,radius_id) +
	321	& kint*QIRsQREF(m,iaer,radius_id+1)
	322	omegIRa(m,gausind,iaer)=
	323	& (1-kint)*omegaIR(m,iaer,radius_id) +
	324	& kint*omegaIR(m,iaer,radius_id+1)
	325	gIRa(m,gausind,iaer)=
	326	& (1-kint)*gIR(m,iaer,radius_id) +
	327	& kint*gIR(m,iaer,radius_id+1)
	328	ENDDO
	329	qrefIRa(gausind,iaer)=
	330	& (1-kint)*QREFir(iaer,radius_id) +
	331	& kint*QREFir(iaer,radius_id+1)
	332	omegrefIRa(gausind,iaer)=
	333	& (1-kint)*omegaREFir(iaer,radius_id) +
	334	& kint*omegaREFir(iaer,radius_id+1)
	335	ENDIF
	336	ENDDO
	337
[89]	338	DO gausind=1,ngau
[38]	339	drad=radiusr*radgaus(gausind)
	340	radGAUSb(gausind,iaer,idomain)=radiusm+drad
	341
	342	radius_id=minloc(abs(radiustab(iaer,idomain,:) -
	343	& (radiusm+drad)),DIM=1)
	344	IF ((radiustab(iaer,idomain,radius_id) -
	345	& (radiusm+drad)).GT.0) THEN
	346	radius_id=radius_id-1
	347	ENDIF
	348	IF (radius_id.GE.nsize(iaer,idomain)) THEN
	349	radius_id=nsize(iaer,idomain)-1
	350	kint = 1.
	351	ELSEIF (radius_id.LT.1) THEN
	352	radius_id=1
	353	kint = 0.
	354	ELSE
	355	kint = ( (radiusm+drad) -
	356	& radiustab(iaer,idomain,radius_id) ) /
	357	& ( radiustab(iaer,idomain,radius_id+1) -
	358	& radiustab(iaer,idomain,radius_id) )
	359	ENDIF
	360	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
	361	DO m=1,nsun
	362	qsqrefVISb(m,gausind,iaer)=
	363	& (1-kint)*QVISsQREF(m,iaer,radius_id) +
	364	& kint*QVISsQREF(m,iaer,radius_id+1)
	365	omegVISb(m,gausind,iaer)=
	366	& (1-kint)*omegaVIS(m,iaer,radius_id) +
	367	& kint*omegaVIS(m,iaer,radius_id+1)
	368	gVISb(m,gausind,iaer)=
	369	& (1-kint)*gVIS(m,iaer,radius_id) +
	370	& kint*gVIS(m,iaer,radius_id+1)
	371	ENDDO
	372	qrefVISb(gausind,iaer)=
	373	& (1-kint)*QREFvis(iaer,radius_id) +
	374	& kint*QREFvis(iaer,radius_id+1)
	375	omegrefVISb(gausind,iaer)=
	376	& (1-kint)*omegaREFvis(iaer,radius_id) +
	377	& kint*omegaREFvis(iaer,radius_id+1)
	378	ELSE ! INFRARED DOMAIN ----------------------------------
	379	DO m=1,nir
	380	qsqrefIRb(m,gausind,iaer)=
	381	& (1-kint)*QIRsQREF(m,iaer,radius_id) +
	382	& kint*QIRsQREF(m,iaer,radius_id+1)
	383	omegIRb(m,gausind,iaer)=
	384	& (1-kint)*omegaIR(m,iaer,radius_id) +
	385	& kint*omegaIR(m,iaer,radius_id+1)
	386	gIRb(m,gausind,iaer)=
	387	& (1-kint)*gIR(m,iaer,radius_id) +
	388	& kint*gIR(m,iaer,radius_id+1)
	389	ENDDO
	390	qrefIRb(gausind,iaer)=
	391	& (1-kint)*QREFir(iaer,radius_id) +
	392	& kint*QREFir(iaer,radius_id+1)
	393	omegrefIRb(gausind,iaer)=
	394	& (1-kint)*omegaREFir(iaer,radius_id) +
	395	& kint*omegaREFir(iaer,radius_id+1)
	396	ENDIF
	397	ENDDO
	398	c==================================================================
	399	IF ( .NOT.varyingnueff(iaer) ) THEN ! CONSTANT NUEFF
	400	c==================================================================
	401	c 2. Compute the scattering parameters using linear
	402	c interpolation over grain sizes and constant nueff
	403	c ---------------------------------------------------
	404
	405	DO lg = 1,nlayer
	406	DO ig = 1,ngrid
	407	c 2.1 Effective radius index and kx calculation
[737]	408	var_tmp=reffrad(ig,lg,iaer)/refftabmin(iaer,idomain)
[38]	409	var_tmp=log(var_tmp)*3.
[737]	410	var_tmp=var_tmp/logvratgrid(iaer,idomain)+1.
[38]	411	grid_i=floor(var_tmp)
	412	IF (grid_i.GE.refftabsize) THEN
	413	c WRITE(,) 'Warning: particle size in grid box #'
	414	c WRITE(,) ig,' is too large to be used by the '
	415	c WRITE(,) 'radiative transfer; please extend the '
	416	c WRITE(,) 'interpolation grid to larger grain sizes.'
	417	grid_i=refftabsize-1
	418	kx = 1.
	419	ELSEIF (grid_i.LT.1) THEN
	420	c WRITE(,) 'Warning: particle size in grid box #'
	421	c WRITE(,) ig,' is too small to be used by the '
	422	c WRITE(,) 'radiative transfer; please extend the '
	423	c WRITE(,) 'interpolation grid to smaller grain sizes.'
	424	grid_i=1
	425	kx = 0.
	426	ELSE
[737]	427	kx = ( reffrad(ig,lg,iaer)-
	428	& refftab(grid_i,iaer,idomain) ) /
	429	& ( refftab(grid_i+1,iaer,idomain)-
	430	& refftab(grid_i,iaer,idomain) )
[38]	431	ENDIF
	432	c 2.3 Integration
	433	DO j=grid_i,grid_i+1
	434	c 2.3.1 Check if the calculation has been done
	435	IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
	436	c 2.3.2 Log-normal dist., r_g and sigma_g are defined
	437	c in [hansen_1974], "Light scattering in planetary
	438	c atmospheres", Space Science Reviews 16 527-610.
	439	c Here, sizedistk1=r_g and sizedistk2=sigma_g^2
	440	sizedistk2 = log(1.+nueffrad(1,1,iaer))
	441	sizedistk1 = exp(2.5*sizedistk2)
[737]	442	sizedistk1 = refftab(j,iaer,idomain) / sizedistk1
[38]	443
	444	normd(j,1,iaer,idomain) = 1e-30
[89]	445	DO gausind=1,ngau
[38]	446	drad=radiusr*radgaus(gausind)
	447	dista(j,1,iaer,idomain,gausind) =
	448	& LOG((radiusm-drad)/sizedistk1)
	449	dista(j,1,iaer,idomain,gausind) =
	450	& EXP(-dista(j,1,iaer,idomain,gausind) *
	451	& dista(j,1,iaer,idomain,gausind) *
	452	& 0.5e0/sizedistk2)/(radiusm-drad)
	453	dista(j,1,iaer,idomain,gausind) =
	454	& dista(j,1,iaer,idomain,gausind) /
	455	& (sqrt(2e0pisizedistk2))
	456
	457	distb(j,1,iaer,idomain,gausind) =
	458	& LOG((radiusm+drad)/sizedistk1)
	459	distb(j,1,iaer,idomain,gausind) =
	460	& EXP(-distb(j,1,iaer,idomain,gausind) *
	461	& distb(j,1,iaer,idomain,gausind) *
	462	& 0.5e0/sizedistk2)/(radiusm+drad)
	463	distb(j,1,iaer,idomain,gausind) =
	464	& distb(j,1,iaer,idomain,gausind) /
	465	& (sqrt(2e0pisizedistk2))
	466
	467	normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) +
	468	& weightgaus(gausind) *
	469	& (
	470	& distb(j,1,iaer,idomain,gausind) * pi *
	471	& radGAUSb(gausind,iaer,idomain) *
	472	& radGAUSb(gausind,iaer,idomain) +
	473	& dista(j,1,iaer,idomain,gausind) * pi *
	474	& radGAUSa(gausind,iaer,idomain) *
	475	& radGAUSa(gausind,iaer,idomain)
	476	& )
	477	ENDDO
[737]	478	IF (normd(j,1,iaer,idomain).EQ.1e-30) THEN
	479	WRITE(,)"normd:", normd(j,1,iaer,idomain)
	480	WRITE(,)"Risk of division by 0 (aeroptproperties.F)"
	481	WRITE(,)"Check the size of the interpolation grid."
	482	CALL ABORT
	483	ENDIF
[38]	484	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
	485	c 2.3.3.vis Initialization
	486	qsqrefVISgrid(j,1,:,iaer)=0.
	487	qextVISgrid(j,1,:,iaer)=0.
	488	qscatVISgrid(j,1,:,iaer)=0.
	489	omegVISgrid(j,1,:,iaer)=0.
	490	gVISgrid(j,1,:,iaer)=0.
	491	qrefVISgrid(j,1,iaer)=0.
	492	qscatrefVISgrid(j,1,iaer)=0.
	493	omegrefVISgrid(j,1,iaer)=0.
	494
[89]	495	DO gausind=1,ngau
[38]	496	DO m=1,nsun
	497	c Convolution:
	498	qextVISgrid(j,1,m,iaer) =
	499	& qextVISgrid(j,1,m,iaer) +
	500	& weightgaus(gausind) *
	501	& (
	502	& qsqrefVISb(m,gausind,iaer) *
	503	& qrefVISb(gausind,iaer) *
	504	& piradGAUSb(gausind,iaer,idomain)
	505	& radGAUSb(gausind,iaer,idomain) *
	506	& distb(j,1,iaer,idomain,gausind) +
	507	& qsqrefVISa(m,gausind,iaer) *
	508	& qrefVISa(gausind,iaer) *
	509	& piradGAUSa(gausind,iaer,idomain)
	510	& radGAUSa(gausind,iaer,idomain) *
	511	& dista(j,1,iaer,idomain,gausind)
	512	& )
	513	qscatVISgrid(j,1,m,iaer) =
	514	& qscatVISgrid(j,1,m,iaer) +
	515	& weightgaus(gausind) *
	516	& (
	517	& omegVISb(m,gausind,iaer) *
	518	& qsqrefVISb(m,gausind,iaer) *
	519	& qrefVISb(gausind,iaer) *
	520	& piradGAUSb(gausind,iaer,idomain)
	521	& radGAUSb(gausind,iaer,idomain) *
	522	& distb(j,1,iaer,idomain,gausind) +
	523	& omegVISa(m,gausind,iaer) *
	524	& qsqrefVISa(m,gausind,iaer) *
	525	& qrefVISa(gausind,iaer) *
	526	& piradGAUSa(gausind,iaer,idomain)
	527	& radGAUSa(gausind,iaer,idomain) *
	528	& dista(j,1,iaer,idomain,gausind)
	529	& )
	530	gVISgrid(j,1,m,iaer) =
	531	& gVISgrid(j,1,m,iaer) +
	532	& weightgaus(gausind) *
	533	& (
	534	& omegVISb(m,gausind,iaer) *
	535	& qsqrefVISb(m,gausind,iaer) *
	536	& qrefVISb(gausind,iaer) *
	537	& gVISb(m,gausind,iaer) *
	538	& piradGAUSb(gausind,iaer,idomain)
	539	& radGAUSb(gausind,iaer,idomain) *
	540	& distb(j,1,iaer,idomain,gausind) +
	541	& omegVISa(m,gausind,iaer) *
	542	& qsqrefVISa(m,gausind,iaer) *
	543	& qrefVISa(gausind,iaer) *
	544	& gVISa(m,gausind,iaer) *
	545	& piradGAUSa(gausind,iaer,idomain)
	546	& radGAUSa(gausind,iaer,idomain) *
	547	& dista(j,1,iaer,idomain,gausind)
	548	& )
	549	ENDDO
	550	qrefVISgrid(j,1,iaer) =
	551	& qrefVISgrid(j,1,iaer) +
	552	& weightgaus(gausind) *
	553	& (
	554	& qrefVISb(gausind,iaer) *
	555	& piradGAUSb(gausind,iaer,idomain)
	556	& radGAUSb(gausind,iaer,idomain) *
	557	& distb(j,1,iaer,idomain,gausind) +
	558	& qrefVISa(gausind,iaer) *
	559	& piradGAUSa(gausind,iaer,idomain)
	560	& radGAUSa(gausind,iaer,idomain) *
	561	& dista(j,1,iaer,idomain,gausind)
	562	& )
	563	qscatrefVISgrid(j,1,iaer) =
	564	& qscatrefVISgrid(j,1,iaer) +
	565	& weightgaus(gausind) *
	566	& (
	567	& omegrefVISb(gausind,iaer) *
	568	& qrefVISb(gausind,iaer) *
	569	& piradGAUSb(gausind,iaer,idomain)
	570	& radGAUSb(gausind,iaer,idomain) *
	571	& distb(j,1,iaer,idomain,gausind) +
	572	& omegrefVISa(gausind,iaer) *
	573	& qrefVISa(gausind,iaer) *
	574	& piradGAUSa(gausind,iaer,idomain)
	575	& radGAUSa(gausind,iaer,idomain) *
	576	& dista(j,1,iaer,idomain,gausind)
	577	& )
	578	ENDDO
	579
	580	qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) /
	581	& normd(j,1,iaer,idomain)
	582	qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /
	583	& normd(j,1,iaer,idomain)
	584	omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /
	585	& qrefVISgrid(j,1,iaer)
	586	DO m=1,nsun
	587	qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /
	588	& normd(j,1,iaer,idomain)
	589	qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /
	590	& normd(j,1,iaer,idomain)
	591	gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) /
	592	& qscatVISgrid(j,1,m,iaer) /
	593	& normd(j,1,iaer,idomain)
	594
	595	qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /
	596	& qrefVISgrid(j,1,iaer)
	597	omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /
	598	& qextVISgrid(j,1,m,iaer)
	599	ENDDO
	600	ELSE ! INFRARED DOMAIN ----------
	601	c 2.3.3.ir Initialization
	602	qsqrefIRgrid(j,1,:,iaer)=0.
	603	qextIRgrid(j,1,:,iaer)=0.
	604	qscatIRgrid(j,1,:,iaer)=0.
	605	omegIRgrid(j,1,:,iaer)=0.
	606	gIRgrid(j,1,:,iaer)=0.
	607	qrefIRgrid(j,1,iaer)=0.
	608	qscatrefIRgrid(j,1,iaer)=0.
	609	omegrefIRgrid(j,1,iaer)=0.
	610
[89]	611	DO gausind=1,ngau
[38]	612	DO m=1,nir
	613	c Convolution:
	614	qextIRgrid(j,1,m,iaer) =
	615	& qextIRgrid(j,1,m,iaer) +
	616	& weightgaus(gausind) *
	617	& (
	618	& qsqrefIRb(m,gausind,iaer) *
	619	& qrefVISb(gausind,iaer) *
	620	& piradGAUSb(gausind,iaer,idomain)
	621	& radGAUSb(gausind,iaer,idomain) *
	622	& distb(j,1,iaer,idomain,gausind) +
	623	& qsqrefIRa(m,gausind,iaer) *
	624	& qrefVISa(gausind,iaer) *
	625	& piradGAUSa(gausind,iaer,idomain)
	626	& radGAUSa(gausind,iaer,idomain) *
	627	& dista(j,1,iaer,idomain,gausind)
	628	& )
	629	qscatIRgrid(j,1,m,iaer) =
	630	& qscatIRgrid(j,1,m,iaer) +
	631	& weightgaus(gausind) *
	632	& (
	633	& omegIRb(m,gausind,iaer) *
	634	& qsqrefIRb(m,gausind,iaer) *
	635	& qrefVISb(gausind,iaer) *
	636	& piradGAUSb(gausind,iaer,idomain)
	637	& radGAUSb(gausind,iaer,idomain) *
	638	& distb(j,1,iaer,idomain,gausind) +
	639	& omegIRa(m,gausind,iaer) *
	640	& qsqrefIRa(m,gausind,iaer) *
	641	& qrefVISa(gausind,iaer) *
	642	& piradGAUSa(gausind,iaer,idomain)
	643	& radGAUSa(gausind,iaer,idomain) *
	644	& dista(j,1,iaer,idomain,gausind)
	645	& )
	646	gIRgrid(j,1,m,iaer) =
	647	& gIRgrid(j,1,m,iaer) +
	648	& weightgaus(gausind) *
	649	& (
	650	& omegIRb(m,gausind,iaer) *
	651	& qsqrefIRb(m,gausind,iaer) *
	652	& qrefVISb(gausind,iaer) *
	653	& gIRb(m,gausind,iaer) *
	654	& piradGAUSb(gausind,iaer,idomain)
	655	& radGAUSb(gausind,iaer,idomain) *
	656	& distb(j,1,iaer,idomain,gausind) +
	657	& omegIRa(m,gausind,iaer) *
	658	& qsqrefIRa(m,gausind,iaer) *
	659	& qrefVISa(gausind,iaer) *
	660	& gIRa(m,gausind,iaer) *
	661	& piradGAUSa(gausind,iaer,idomain)
	662	& radGAUSa(gausind,iaer,idomain) *
	663	& dista(j,1,iaer,idomain,gausind)
	664	& )
	665	ENDDO
	666	qrefIRgrid(j,1,iaer) =
	667	& qrefIRgrid(j,1,iaer) +
	668	& weightgaus(gausind) *
	669	& (
	670	& qrefIRb(gausind,iaer) *
	671	& piradGAUSb(gausind,iaer,idomain)
	672	& radGAUSb(gausind,iaer,idomain) *
	673	& distb(j,1,iaer,idomain,gausind) +
	674	& qrefIRa(gausind,iaer) *
	675	& piradGAUSa(gausind,iaer,idomain)
	676	& radGAUSa(gausind,iaer,idomain) *
	677	& dista(j,1,iaer,idomain,gausind)
	678	& )
	679	qscatrefIRgrid(j,1,iaer) =
	680	& qscatrefIRgrid(j,1,iaer) +
	681	& weightgaus(gausind) *
	682	& (
	683	& omegrefIRb(gausind,iaer) *
	684	& qrefIRb(gausind,iaer) *
	685	& piradGAUSb(gausind,iaer,idomain)
	686	& radGAUSb(gausind,iaer,idomain) *
	687	& distb(j,1,iaer,idomain,gausind) +
	688	& omegrefIRa(gausind,iaer) *
	689	& qrefIRa(gausind,iaer) *
	690	& piradGAUSa(gausind,iaer,idomain)
	691	& radGAUSa(gausind,iaer,idomain) *
	692	& dista(j,1,iaer,idomain,gausind)
	693	& )
	694	ENDDO
	695
	696	qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) /
	697	& normd(j,1,iaer,idomain)
	698	qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /
	699	& normd(j,1,iaer,idomain)
	700	omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /
	701	& qrefIRgrid(j,1,iaer)
	702	DO m=1,nir
	703	qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /
	704	& normd(j,1,iaer,idomain)
	705	qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /
	706	& normd(j,1,iaer,idomain)
	707	gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) /
	708	& qscatIRgrid(j,1,m,iaer) /
	709	& normd(j,1,iaer,idomain)
	710
	711	qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /
	712	& qrefVISgrid(j,1,iaer)
	713	omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /
	714	& qextIRgrid(j,1,m,iaer)
	715	ENDDO
	716	ENDIF ! --------------------------
	717	checkgrid(j,1,iaer,idomain) = .true.
	718	ENDIF !checkgrid
	719	ENDDO !grid_i
	720	c 2.4 Linear interpolation
	721	k1 = (1-kx)
	722	k2 = kx
	723	IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
	724	DO m=1,nsun
	725	QVISsQREF3d(ig,lg,m,iaer) =
	726	& k1*qsqrefVISgrid(grid_i,1,m,iaer) +
	727	& k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
	728	omegaVIS3d(ig,lg,m,iaer) =
	729	& k1*omegVISgrid(grid_i,1,m,iaer) +
	730	& k2*omegVISgrid(grid_i+1,1,m,iaer)
	731	gVIS3d(ig,lg,m,iaer) =
	732	& k1*gVISgrid(grid_i,1,m,iaer) +
	733	& k2*gVISgrid(grid_i+1,1,m,iaer)
	734	ENDDO !nsun
	735	QREFvis3d(ig,lg,iaer) =
	736	& k1*qrefVISgrid(grid_i,1,iaer) +
	737	& k2*qrefVISgrid(grid_i+1,1,iaer)
	738	omegaREFvis3d(ig,lg,iaer) =
	739	& k1*omegrefVISgrid(grid_i,1,iaer) +
	740	& k2*omegrefVISgrid(grid_i+1,1,iaer)
	741	ELSE ! INFRARED -----------------------
	742	DO m=1,nir
	743	QIRsQREF3d(ig,lg,m,iaer) =
	744	& k1*qsqrefIRgrid(grid_i,1,m,iaer) +
	745	& k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
	746	omegaIR3d(ig,lg,m,iaer) =
	747	& k1*omegIRgrid(grid_i,1,m,iaer) +
	748	& k2*omegIRgrid(grid_i+1,1,m,iaer)
	749	gIR3d(ig,lg,m,iaer) =
	750	& k1*gIRgrid(grid_i,1,m,iaer) +
	751	& k2*gIRgrid(grid_i+1,1,m,iaer)
	752	ENDDO !nir
	753	QREFir3d(ig,lg,iaer) =
	754	& k1*qrefIRgrid(grid_i,1,iaer) +
	755	& k2*qrefIRgrid(grid_i+1,1,iaer)
	756	omegaREFir3d(ig,lg,iaer) =
	757	& k1*omegrefIRgrid(grid_i,1,iaer) +
	758	& k2*omegrefIRgrid(grid_i+1,1,iaer)
	759	ENDIF ! --------------------------------
	760	ENDDO !nlayermx
	761	ENDDO !ngridmx
	762	c==================================================================
	763	ELSE ! VARYING NUEFF
	764	c==================================================================
	765	c 3. Compute the scattering parameters in each grid box
	766	c using bilinear interpolation over the grain sizes
	767	c and the effective variances;
	768	c -----------------------------------------------------
	769
	770	DO lg = 1,nlayer
	771	DO ig = 1,ngrid
	772	c 3.1 Effective variance index and ky calculation
	773	var_tmp=log(nueffrad(ig,lg,iaer))
	774	grid_j=floor( (nuefftabsize-1)/(nuefftabmax-nuefftabmin)*
	775	& (var_tmp-nuefftabmin)+1. )
	776	IF (grid_j.GE.nuefftabsize) THEN
	777	c WRITE(,) 'Warning: effective variance '
	778	c WRITE(,) 'is too large to be used by the '
	779	c WRITE(,) 'radiative transfer; please extend the '
	780	c WRITE(,) 'interpolation grid to larger values.'
	781	grid_j=nuefftabsize-1
	782	ky = 1.
	783	ELSEIF (grid_j.LT.1) THEN
	784	c WRITE(,) 'Warning: effective variance '
	785	c WRITE(,) 'is too small to be used by the '
	786	c WRITE(,) 'radiative transfer; please extend the '
	787	c WRITE(,) 'interpolation grid to smaller values.'
	788	grid_j=1
	789	ky = 0.
	790	ELSE
[737]	791	ky = ( nueffrad(ig,lg,iaer)-
	792	& nuefftab(grid_j,iaer,idomain) ) /
	793	& ( nuefftab(grid_j+1,iaer,idomain)-
	794	& nuefftab(grid_j,iaer,idomain) )
[38]	795	ENDIF
	796	c 3.2 Effective radius index and kx calculation
[737]	797	var_tmp=reffrad(ig,lg,iaer)/refftabmin(iaer,idomain)
[38]	798	var_tmp=log(var_tmp)*3.
[737]	799	var_tmp=var_tmp/logvratgrid(iaer,idomain)+1.
[38]	800	grid_i=floor(var_tmp)
	801	IF (grid_i.GE.refftabsize) THEN
	802	c WRITE(,) 'Warning: particle size in grid box #'
	803	c WRITE(,) ig,' is too large to be used by the '
	804	c WRITE(,) 'radiative transfer; please extend the '
	805	c WRITE(,) 'interpolation grid to larger grain sizes.'
	806	grid_i=refftabsize-1
	807	kx = 1.
	808	ELSEIF (grid_i.LT.1) THEN
	809	c WRITE(,) 'Warning: particle size in grid box #'
	810	c WRITE(,) ig,' is too small to be used by the '
	811	c WRITE(,) 'radiative transfer; please extend the '
	812	c WRITE(,) 'interpolation grid to smaller grain sizes.'
	813	grid_i=1
	814	kx = 0.
	815	ELSE
[737]	816	kx = ( reffrad(ig,lg,iaer)-
	817	& refftab(grid_i,iaer,idomain) ) /
	818	& ( refftab(grid_i+1,iaer,idomain)-
	819	& refftab(grid_i,iaer,idomain) )
[38]	820	ENDIF
	821	c 3.3 Integration
	822	DO j=grid_i,grid_i+1
	823	DO k=grid_j,grid_j+1
	824	c 3.3.1 Check if the calculation has been done
	825	IF (.NOT.checkgrid(j,k,iaer,idomain)) THEN
	826
	827	c 3.3.2 Log-normal dist., r_g and sigma_g are defined
	828	c in [hansen_1974], "Light scattering in planetary
	829	c atmospheres", Space Science Reviews 16 527-610.
	830	c Here, sizedistk1=r_g and sizedistk2=sigma_g^2
[737]	831	sizedistk2 = log(1.+nuefftab(k,iaer,idomain))
[38]	832	sizedistk1 = exp(2.5*sizedistk2)
[737]	833	sizedistk1 = refftab(j,iaer,idomain) / sizedistk1
[38]	834
	835	normd(j,k,iaer,idomain) = 1e-30
[89]	836	DO gausind=1,ngau
[38]	837	drad=radiusr*radgaus(gausind)
	838
	839	dista(j,k,iaer,idomain,gausind) =
	840	& LOG((radiusm-drad)/sizedistk1)
	841	dista(j,k,iaer,idomain,gausind) =
	842	& EXP(-dista(j,k,iaer,idomain,gausind) *
	843	& dista(j,k,iaer,idomain,gausind) *
	844	& 0.5e0/sizedistk2)/(radiusm-drad)
	845	dista(j,k,iaer,idomain,gausind) =
	846	& dista(j,k,iaer,idomain,gausind) /
	847	& (sqrt(2e0pisizedistk2))
	848
	849	distb(j,k,iaer,idomain,gausind) =
	850	& LOG((radiusm+drad)/sizedistk1)
	851	distb(j,k,iaer,idomain,gausind) =
	852	& EXP(-distb(j,k,iaer,idomain,gausind) *
	853	& distb(j,k,iaer,idomain,gausind) *
	854	& 0.5e0/sizedistk2)/(radiusm+drad)
	855	distb(j,k,iaer,idomain,gausind) =
	856	& distb(j,k,iaer,idomain,gausind) /
	857	& (sqrt(2e0pisizedistk2))
	858
	859	normd(j,k,iaer,idomain)=normd(j,k,iaer,idomain) +
	860	& weightgaus(gausind) *
	861	& (
	862	& distb(j,k,iaer,idomain,gausind) * pi *
	863	& radGAUSb(gausind,iaer,idomain) *
	864	& radGAUSb(gausind,iaer,idomain) +
	865	& dista(j,k,iaer,idomain,gausind) * pi *
	866	& radGAUSa(gausind,iaer,idomain) *
	867	& radGAUSa(gausind,iaer,idomain)
	868	& )
	869	ENDDO
[737]	870	IF (normd(j,k,iaer,idomain).EQ.1e-30) THEN
	871	WRITE(,)"normd:", normd(j,k,iaer,idomain)
	872	WRITE(,)"Risk of division by 0 (aeroptproperties.F)"
	873	WRITE(,)"Check the size of the interpolation grid."
	874	CALL ABORT
	875	ENDIF
[38]	876	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
	877	c 2.3.3.vis Initialization
	878	qsqrefVISgrid(j,k,:,iaer)=0.
	879	qextVISgrid(j,k,:,iaer)=0.
	880	qscatVISgrid(j,k,:,iaer)=0.
	881	omegVISgrid(j,k,:,iaer)=0.
	882	gVISgrid(j,k,:,iaer)=0.
	883	qrefVISgrid(j,k,iaer)=0.
	884	qscatrefVISgrid(j,k,iaer)=0.
	885	omegrefVISgrid(j,k,iaer)=0.
	886
[89]	887	DO gausind=1,ngau
[38]	888	DO m=1,nsun
	889	c Convolution:
	890	qextVISgrid(j,k,m,iaer) =
	891	& qextVISgrid(j,k,m,iaer) +
	892	& weightgaus(gausind) *
	893	& (
	894	& qsqrefVISb(m,gausind,iaer) *
	895	& qrefVISb(gausind,iaer) *
	896	& piradGAUSb(gausind,iaer,idomain)
	897	& radGAUSb(gausind,iaer,idomain) *
	898	& distb(j,k,iaer,idomain,gausind) +
	899	& qsqrefVISa(m,gausind,iaer) *
	900	& qrefVISa(gausind,iaer) *
	901	& piradGAUSa(gausind,iaer,idomain)
	902	& radGAUSa(gausind,iaer,idomain) *
	903	& dista(j,k,iaer,idomain,gausind)
	904	& )
	905	qscatVISgrid(j,k,m,iaer) =
	906	& qscatVISgrid(j,k,m,iaer) +
	907	& weightgaus(gausind) *
	908	& (
	909	& omegVISb(m,gausind,iaer) *
	910	& qsqrefVISb(m,gausind,iaer) *
	911	& qrefVISb(gausind,iaer) *
	912	& piradGAUSb(gausind,iaer,idomain)
	913	& radGAUSb(gausind,iaer,idomain) *
	914	& distb(j,k,iaer,idomain,gausind) +
	915	& omegVISa(m,gausind,iaer) *
	916	& qsqrefVISa(m,gausind,iaer) *
	917	& qrefVISa(gausind,iaer) *
	918	& piradGAUSa(gausind,iaer,idomain)
	919	& radGAUSa(gausind,iaer,idomain) *
	920	& dista(j,k,iaer,idomain,gausind)
	921	& )
	922	gVISgrid(j,k,m,iaer) =
	923	& gVISgrid(j,k,m,iaer) +
	924	& weightgaus(gausind) *
	925	& (
	926	& omegVISb(m,gausind,iaer) *
	927	& qsqrefVISb(m,gausind,iaer) *
	928	& qrefVISb(gausind,iaer) *
	929	& gVISb(m,gausind,iaer) *
	930	& piradGAUSb(gausind,iaer,idomain)
	931	& radGAUSb(gausind,iaer,idomain) *
	932	& distb(j,k,iaer,idomain,gausind) +
	933	& omegVISa(m,gausind,iaer) *
	934	& qsqrefVISa(m,gausind,iaer) *
	935	& qrefVISa(gausind,iaer) *
	936	& gVISa(m,gausind,iaer) *
	937	& piradGAUSa(gausind,iaer,idomain)
	938	& radGAUSa(gausind,iaer,idomain) *
	939	& dista(j,k,iaer,idomain,gausind)
	940	& )
	941	ENDDO
	942	qrefVISgrid(j,k,iaer) =
	943	& qrefVISgrid(j,k,iaer) +
	944	& weightgaus(gausind) *
	945	& (
	946	& qrefVISb(gausind,iaer) *
	947	& piradGAUSb(gausind,iaer,idomain)
	948	& radGAUSb(gausind,iaer,idomain) *
	949	& distb(j,k,iaer,idomain,gausind) +
	950	& qrefVISa(gausind,iaer) *
	951	& piradGAUSa(gausind,iaer,idomain)
	952	& radGAUSa(gausind,iaer,idomain) *
	953	& dista(j,k,iaer,idomain,gausind)
	954	& )
	955	qscatrefVISgrid(j,k,iaer) =
	956	& qscatrefVISgrid(j,k,iaer) +
	957	& weightgaus(gausind) *
	958	& (
	959	& omegrefVISb(gausind,iaer) *
	960	& qrefVISb(gausind,iaer) *
	961	& piradGAUSb(gausind,iaer,idomain)
	962	& radGAUSb(gausind,iaer,idomain) *
	963	& distb(j,k,iaer,idomain,gausind) +
	964	& omegrefVISa(gausind,iaer) *
	965	& qrefVISa(gausind,iaer) *
	966	& piradGAUSa(gausind,iaer,idomain)
	967	& radGAUSa(gausind,iaer,idomain) *
	968	& dista(j,k,iaer,idomain,gausind)
	969	& )
	970	ENDDO
	971	qrefVISgrid(j,k,iaer)=qrefVISgrid(j,k,iaer) /
	972	& normd(j,k,iaer,idomain)
	973	qscatrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) /
	974	& normd(j,k,iaer,idomain)
	975	omegrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) /
	976	& qrefVISgrid(j,k,iaer)
	977	DO m=1,nsun
	978	qextVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) /
	979	& normd(j,k,iaer,idomain)
	980	qscatVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) /
	981	& normd(j,k,iaer,idomain)
	982	gVISgrid(j,k,m,iaer)=gVISgrid(j,k,m,iaer) /
	983	& qscatVISgrid(j,k,m,iaer) /
	984	& normd(j,k,iaer,idomain)
	985
	986	qsqrefVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) /
	987	& qrefVISgrid(j,k,iaer)
	988	omegVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) /
	989	& qextVISgrid(j,k,m,iaer)
	990	ENDDO
	991	ELSE ! INFRARED DOMAIN ----------
	992	c 2.3.3.ir Initialization
	993	qsqrefIRgrid(j,k,:,iaer)=0.
	994	qextIRgrid(j,k,:,iaer)=0.
	995	qscatIRgrid(j,k,:,iaer)=0.
	996	omegIRgrid(j,k,:,iaer)=0.
	997	gIRgrid(j,k,:,iaer)=0.
	998	qrefIRgrid(j,k,iaer)=0.
	999	qscatrefIRgrid(j,k,iaer)=0.
	1000	omegrefIRgrid(j,k,iaer)=0.
	1001
[89]	1002	DO gausind=1,ngau
[38]	1003	DO m=1,nir
	1004	c Convolution:
	1005	qextIRgrid(j,k,m,iaer) =
	1006	& qextIRgrid(j,k,m,iaer) +
	1007	& weightgaus(gausind) *
	1008	& (
	1009	& qsqrefIRb(m,gausind,iaer) *
	1010	& qrefVISb(gausind,iaer) *
	1011	& piradGAUSb(gausind,iaer,idomain)
	1012	& radGAUSb(gausind,iaer,idomain) *
	1013	& distb(j,k,iaer,idomain,gausind) +
	1014	& qsqrefIRa(m,gausind,iaer) *
	1015	& qrefVISa(gausind,iaer) *
	1016	& piradGAUSa(gausind,iaer,idomain)
	1017	& radGAUSa(gausind,iaer,idomain) *
	1018	& dista(j,k,iaer,idomain,gausind)
	1019	& )
	1020	qscatIRgrid(j,k,m,iaer) =
	1021	& qscatIRgrid(j,k,m,iaer) +
	1022	& weightgaus(gausind) *
	1023	& (
	1024	& omegIRb(m,gausind,iaer) *
	1025	& qsqrefIRb(m,gausind,iaer) *
	1026	& qrefVISb(gausind,iaer) *
	1027	& piradGAUSb(gausind,iaer,idomain)
	1028	& radGAUSb(gausind,iaer,idomain) *
	1029	& distb(j,k,iaer,idomain,gausind) +
	1030	& omegIRa(m,gausind,iaer) *
	1031	& qsqrefIRa(m,gausind,iaer) *
	1032	& qrefVISa(gausind,iaer) *
	1033	& piradGAUSa(gausind,iaer,idomain)
	1034	& radGAUSa(gausind,iaer,idomain) *
	1035	& dista(j,k,iaer,idomain,gausind)
	1036	& )
	1037	gIRgrid(j,k,m,iaer) =
	1038	& gIRgrid(j,k,m,iaer) +
	1039	& weightgaus(gausind) *
	1040	& (
	1041	& omegIRb(m,gausind,iaer) *
	1042	& qsqrefIRb(m,gausind,iaer) *
	1043	& qrefVISb(gausind,iaer) *
	1044	& gIRb(m,gausind,iaer) *
	1045	& piradGAUSb(gausind,iaer,idomain)
	1046	& radGAUSb(gausind,iaer,idomain) *
	1047	& distb(j,k,iaer,idomain,gausind) +
	1048	& omegIRa(m,gausind,iaer) *
	1049	& qsqrefIRa(m,gausind,iaer) *
	1050	& qrefVISa(gausind,iaer) *
	1051	& gIRa(m,gausind,iaer) *
	1052	& piradGAUSa(gausind,iaer,idomain)
	1053	& radGAUSa(gausind,iaer,idomain) *
	1054	& dista(j,k,iaer,idomain,gausind)
	1055	& )
	1056	ENDDO
	1057	qrefIRgrid(j,k,iaer) =
	1058	& qrefIRgrid(j,k,iaer) +
	1059	& weightgaus(gausind) *
	1060	& (
	1061	& qrefIRb(gausind,iaer) *
	1062	& piradGAUSb(gausind,iaer,idomain)
	1063	& radGAUSb(gausind,iaer,idomain) *
	1064	& distb(j,k,iaer,idomain,gausind) +
	1065	& qrefIRa(gausind,iaer) *
	1066	& piradGAUSa(gausind,iaer,idomain)
	1067	& radGAUSa(gausind,iaer,idomain) *
	1068	& dista(j,k,iaer,idomain,gausind)
	1069	& )
	1070	qscatrefIRgrid(j,k,iaer) =
	1071	& qscatrefIRgrid(j,k,iaer) +
	1072	& weightgaus(gausind) *
	1073	& (
	1074	& omegrefIRb(gausind,iaer) *
	1075	& qrefIRb(gausind,iaer) *
	1076	& piradGAUSb(gausind,iaer,idomain)
	1077	& radGAUSb(gausind,iaer,idomain) *
	1078	& distb(j,k,iaer,idomain,gausind) +
	1079	& omegrefIRa(gausind,iaer) *
	1080	& qrefIRa(gausind,iaer) *
	1081	& piradGAUSa(gausind,iaer,idomain)
	1082	& radGAUSa(gausind,iaer,idomain) *
	1083	& dista(j,k,iaer,idomain,gausind)
	1084	& )
	1085	ENDDO
	1086	qrefIRgrid(j,k,iaer)=qrefIRgrid(j,k,iaer) /
	1087	& normd(j,k,iaer,idomain)
	1088	qscatrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) /
	1089	& normd(j,k,iaer,idomain)
	1090	omegrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) /
	1091	& qrefIRgrid(j,k,iaer)
	1092	DO m=1,nir
	1093	qextIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) /
	1094	& normd(j,k,iaer,idomain)
	1095	qscatIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) /
	1096	& normd(j,k,iaer,idomain)
	1097	gIRgrid(j,k,m,iaer)=gIRgrid(j,k,m,iaer) /
	1098	& qscatIRgrid(j,k,m,iaer) /
	1099	& normd(j,k,iaer,idomain)
	1100
	1101	qsqrefIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) /
	1102	& qrefVISgrid(j,k,iaer)
	1103	omegIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) /
	1104	& qextIRgrid(j,k,m,iaer)
	1105	ENDDO
	1106	ENDIF ! --------------------------
	1107	checkgrid(j,k,iaer,idomain) = .true.
	1108	ENDIF !checkgrid
	1109	ENDDO !grid_j
	1110	ENDDO !grid_i
	1111	c 3.4 Bilinear interpolation
	1112	k1 = (1-kx)*(1-ky)
	1113	k2 = kx*(1-ky)
	1114	k3 = kx*ky
	1115	k4 = (1-kx)*ky
	1116	IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
	1117	DO m=1,nsun
	1118	QVISsQREF3d(ig,lg,m,iaer) =
	1119	& k1*qsqrefVISgrid(grid_i,grid_j,m,iaer) +
	1120	& k2*qsqrefVISgrid(grid_i+1,grid_j,m,iaer) +
	1121	& k3*qsqrefVISgrid(grid_i+1,grid_j+1,m,iaer) +
	1122	& k4*qsqrefVISgrid(grid_i,grid_j+1,m,iaer)
	1123	omegaVIS3d(ig,lg,m,iaer) =
	1124	& k1*omegVISgrid(grid_i,grid_j,m,iaer) +
	1125	& k2*omegVISgrid(grid_i+1,grid_j,m,iaer) +
	1126	& k3*omegVISgrid(grid_i+1,grid_j+1,m,iaer) +
	1127	& k4*omegVISgrid(grid_i,grid_j+1,m,iaer)
	1128	gVIS3d(ig,lg,m,iaer) =
	1129	& k1*gVISgrid(grid_i,grid_j,m,iaer) +
	1130	& k2*gVISgrid(grid_i+1,grid_j,m,iaer) +
	1131	& k3*gVISgrid(grid_i+1,grid_j+1,m,iaer) +
	1132	& k4*gVISgrid(grid_i,grid_j+1,m,iaer)
	1133	ENDDO !nsun
	1134	QREFvis3d(ig,lg,iaer) =
	1135	& k1*qrefVISgrid(grid_i,grid_j,iaer) +
	1136	& k2*qrefVISgrid(grid_i+1,grid_j,iaer) +
	1137	& k3*qrefVISgrid(grid_i+1,grid_j+1,iaer) +
	1138	& k4*qrefVISgrid(grid_i,grid_j+1,iaer)
	1139	omegaREFvis3d(ig,lg,iaer) =
	1140	& k1*omegrefVISgrid(grid_i,grid_j,iaer) +
	1141	& k2*omegrefVISgrid(grid_i+1,grid_j,iaer) +
	1142	& k3*omegrefVISgrid(grid_i+1,grid_j+1,iaer) +
	1143	& k4*omegrefVISgrid(grid_i,grid_j+1,iaer)
	1144	ELSE ! INFRARED -----------------------
	1145	DO m=1,nir
	1146	QIRsQREF3d(ig,lg,m,iaer) =
	1147	& k1*qsqrefIRgrid(grid_i,grid_j,m,iaer) +
	1148	& k2*qsqrefIRgrid(grid_i+1,grid_j,m,iaer) +
	1149	& k3*qsqrefIRgrid(grid_i+1,grid_j+1,m,iaer) +
	1150	& k4*qsqrefIRgrid(grid_i,grid_j+1,m,iaer)
	1151	omegaIR3d(ig,lg,m,iaer) =
	1152	& k1*omegIRgrid(grid_i,grid_j,m,iaer) +
	1153	& k2*omegIRgrid(grid_i+1,grid_j,m,iaer) +
	1154	& k3*omegIRgrid(grid_i+1,grid_j+1,m,iaer) +
	1155	& k4*omegIRgrid(grid_i,grid_j+1,m,iaer)
	1156	gIR3d(ig,lg,m,iaer) =
	1157	& k1*gIRgrid(grid_i,grid_j,m,iaer) +
	1158	& k2*gIRgrid(grid_i+1,grid_j,m,iaer) +
	1159	& k3*gIRgrid(grid_i+1,grid_j+1,m,iaer) +
	1160	& k4*gIRgrid(grid_i,grid_j+1,m,iaer)
	1161	ENDDO !nir
	1162	QREFir3d(ig,lg,iaer) =
	1163	& k1*qrefIRgrid(grid_i,grid_j,iaer) +
	1164	& k2*qrefIRgrid(grid_i+1,grid_j,iaer) +
	1165	& k3*qrefIRgrid(grid_i+1,grid_j+1,iaer) +
	1166	& k4*qrefIRgrid(grid_i,grid_j+1,iaer)
	1167	omegaREFir3d(ig,lg,iaer) =
	1168	& k1*omegrefIRgrid(grid_i,grid_j,iaer) +
	1169	& k2*omegrefIRgrid(grid_i+1,grid_j,iaer) +
	1170	& k3*omegrefIRgrid(grid_i+1,grid_j+1,iaer) +
	1171	& k4*omegrefIRgrid(grid_i,grid_j+1,iaer)
	1172	ENDIF ! --------------------------------
	1173	ENDDO !nlayermx
	1174	ENDDO !ngridmx
	1175
	1176	ENDIF ! varyingnueff
	1177	c==================================================================
	1178	ENDDO ! idomain
	1179
	1180	ENDIF ! nsize = 1
	1181
	1182	ENDDO ! iaer (loop on aerosol kind)
	1183
[89]	1184	c=====Radiative properties - TESTS=================================
	1185	IF (out_qwg) THEN
[784]	1186	c -------------------------------------------------------------
	1187	IF (ngrid.NE.1) THEN
	1188	out_ndim = 3
	1189	ELSE
	1190	out_ndim = 1
	1191	ENDIF
	1192	c -------------------------------------------------------------
[89]	1193	DO out_nchannel = 1, 2
	1194	c -------------------------------------------------------------
	1195	DO lg = 1, nlayer
	1196	DO ig = 1, ngrid
	1197	out_qext(ig,lg) =
	1198	& QVISsQREF3d(ig,lg,out_nchannel,out_iaer)*
	1199	& QREFvis3d(ig,lg,out_iaer)
	1200	out_omeg(ig,lg) =
	1201	& omegaVIS3d(ig,lg,out_nchannel,out_iaer)
	1202	out_g(ig,lg) = gVIS3d(ig,lg,out_nchannel,out_iaer)
	1203	ENDDO ! ig
	1204	ENDDO ! lg
	1205	write(out_str(1:1),'(i1.1)') out_nchannel
	1206	call WRITEDIAGFI(ngrid,'qextvis'//out_str,"Ext.efficiency","",
[784]	1207	& out_ndim,out_qext)
[89]	1208	call WRITEDIAGFI(ngrid,'omegvis'//out_str,"Sing.Scat.Alb.","",
[784]	1209	& out_ndim,out_omeg)
[89]	1210	call WRITEDIAGFI(ngrid,'gvis'//out_str,"Asym.Factor","",
[784]	1211	& out_ndim,out_g)
[89]	1212	c -------------------------------------------------------------
	1213	ENDDO ! out_nchannel
	1214	DO out_nchannel = 2, 4
	1215	c -------------------------------------------------------------
	1216	DO lg = 1, nlayer
	1217	DO ig = 1, ngrid
	1218	out_qext(ig,lg) =
	1219	& QIRsQREF3d(ig,lg,out_nchannel,out_iaer)*
	1220	& QREFir3d(ig,lg,out_iaer)
	1221	out_omeg(ig,lg) =
	1222	& omegaIR3d(ig,lg,out_nchannel,out_iaer)
	1223	out_g(ig,lg) = gIR3d(ig,lg,out_nchannel,out_iaer)
	1224	ENDDO ! ig
	1225	ENDDO ! lg
	1226	write(out_str(1:1),'(i1.1)') out_nchannel
[690]	1227	call WRITEDIAGFI(ngrid,'qextir'//out_str,"Ext.efficiency","",
[784]	1228	& out_ndim,out_qext)
[690]	1229	call WRITEDIAGFI(ngrid,'omegir'//out_str,"Sing.Scat.Alb.","",
[784]	1230	& out_ndim,out_omeg)
[690]	1231	call WRITEDIAGFI(ngrid,'gir'//out_str,"Asym.Factor","",
[784]	1232	& out_ndim,out_g)
[89]	1233	c -------------------------------------------------------------
	1234	ENDDO ! out_nchannel
	1235	call WRITEDIAGFI(ngrid,"omegvisref","Sing.Scat.Alb.","",
[784]	1236	& out_ndim,omegaREFvis3d(1,1,out_iaer))
[89]	1237	call WRITEDIAGFI(ngrid,"omegirref","Sing.Scat.Alb.","",
[784]	1238	& out_ndim,omegaREFir3d(1,1,out_iaer))
[89]	1239	ENDIF ! out_qwg
	1240	c==================================================================
	1241
[38]	1242	RETURN
	1243	END

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

Download in other formats:

Original Format