Context Navigation

aeroptproperties.F @ 735

Last change on this file since 735 was 690, checked in by acolaitis, 13 years ago
Code re-organization in diverse parts of the GCM code. These are NOT cosmetic changes, but are needed for compilation of the Mesoscale model in NESTED configuration
File size: 52.9 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)
[89]	39	REAL, PARAMETER :: refftabmin = 1e-8 !5e-8
[38]	40	REAL, PARAMETER :: refftabmax = 35e-6
	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
	89	REAL,SAVE :: refftab(refftabsize)
	90	c Variance axis of the interpolation grid
	91	REAL,SAVE :: nuefftab(nuefftabsize)
	92	c Volume ratio of the grid
	93	REAL,SAVE :: logvratgrid,vratgrid
	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.
	181	INTEGER, PARAMETER :: out_iaer = 1
	182	REAL :: out_qext(ngridmx,nlayermx)
	183	REAL :: out_omeg(ngridmx,nlayermx)
	184	REAL :: out_g(ngridmx,nlayermx)
	185	INTEGER :: out_nchannel
	186	CHARACTER*1 :: out_str
	187
[38]	188	DO iaer = 1, naerkind ! Loop on aerosol kind
	189	IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN
	190	c==================================================================
	191	c If there is one single particle size, optical
	192	c properties of the considered aerosol are homogeneous
	193	DO lg = 1, nlayer
	194	DO ig = 1, ngrid
	195	DO chg = 1, nsun
	196	QVISsQREF3d(ig,lg,chg,iaer)=QVISsQREF(chg,iaer,1)
	197	omegaVIS3d(ig,lg,chg,iaer)=omegaVIS(chg,iaer,1)
	198	gVIS3d(ig,lg,chg,iaer)=gVIS(chg,iaer,1)
	199	ENDDO
	200	DO chg = 1, nir
	201	QIRsQREF3d(ig,lg,chg,iaer)=QIRsQREF(chg,iaer,1)
	202	omegaIR3d(ig,lg,chg,iaer)=omegaIR(chg,iaer,1)
	203	gIR3d(ig,lg,chg,iaer)=gIR(chg,iaer,1)
	204	ENDDO
	205	QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1)
	206	QREFir3d(ig,lg,iaer)=QREFir(iaer,1)
	207	omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1)
	208	omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1)
	209	ENDDO
	210	ENDDO
	211	ELSE ! Varying effective radius and variance
	212	DO idomain = 1, 2 ! Loop on visible or infrared channel
	213	c==================================================================
	214	c 1. Creating the effective radius and variance grid
	215	c --------------------------------------------------
	216	IF (firstcall) THEN
	217
	218	c 1.1 Pi!
	219	pi = 2. * asin(1.e0)
	220
	221	c 1.2 Effective radius
	222	refftab(1) = refftabmin
	223	refftab(refftabsize) = refftabmax
	224
	225	logvratgrid = log(refftabmax/refftabmin) /
	226	& float(refftabsize-1)*3.
	227	vratgrid = exp(logvratgrid)
	228
	229	do i = 2, refftabsize-1
	230	refftab(i) = refftab(i-1)vratgrid*(1./3.)
	231	enddo
	232
	233	c 1.3 Effective variance
	234	do i = 0, nuefftabsize-1
	235	nuefftab(i+1) = exp( nuefftabmin +
	236	& i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
	237	enddo
	238	firstcall = .false.
	239	ENDIF
	240
	241	c 1.4 Radius middle point and range for Gauss integration
	242	radiusm=
	243	& 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) +
	244	& radiustab(iaer,idomain,1))
	245	radiusr=
	246	& 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) -
	247	& radiustab(iaer,idomain,1))
	248
	249	c 1.5 Interpolating data at the Gauss quadrature points:
[89]	250	DO gausind=1,ngau
[38]	251	drad=radiusr*radgaus(gausind)
	252	radGAUSa(gausind,iaer,idomain)=radiusm-drad
	253
	254	radius_id=minloc(abs(radiustab(iaer,idomain,:) -
	255	& (radiusm-drad)),DIM=1)
	256	IF ((radiustab(iaer,idomain,radius_id) -
	257	& (radiusm-drad)).GT.0) THEN
	258	radius_id=radius_id-1
	259	ENDIF
	260	IF (radius_id.GE.nsize(iaer,idomain)) THEN
	261	radius_id=nsize(iaer,idomain)-1
	262	kint = 1.
	263	ELSEIF (radius_id.LT.1) THEN
	264	radius_id=1
	265	kint = 0.
	266	ELSE
	267	kint = ( (radiusm-drad) -
	268	& radiustab(iaer,idomain,radius_id) ) /
	269	& ( radiustab(iaer,idomain,radius_id+1) -
	270	& radiustab(iaer,idomain,radius_id) )
	271	ENDIF
	272	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
	273	DO m=1,nsun
	274	qsqrefVISa(m,gausind,iaer)=
	275	& (1-kint)*QVISsQREF(m,iaer,radius_id) +
	276	& kint*QVISsQREF(m,iaer,radius_id+1)
	277	omegVISa(m,gausind,iaer)=
	278	& (1-kint)*omegaVIS(m,iaer,radius_id) +
	279	& kint*omegaVIS(m,iaer,radius_id+1)
	280	gVISa(m,gausind,iaer)=
	281	& (1-kint)*gVIS(m,iaer,radius_id) +
	282	& kint*gVIS(m,iaer,radius_id+1)
	283	ENDDO
	284	qrefVISa(gausind,iaer)=
	285	& (1-kint)*QREFvis(iaer,radius_id) +
	286	& kint*QREFvis(iaer,radius_id+1)
	287	omegrefVISa(gausind,iaer)=
	288	& (1-kint)*omegaREFvis(iaer,radius_id) +
	289	& kint*omegaREFvis(iaer,radius_id+1)
	290	ELSE ! INFRARED DOMAIN ----------------------------------
	291	DO m=1,nir
	292	qsqrefIRa(m,gausind,iaer)=
	293	& (1-kint)*QIRsQREF(m,iaer,radius_id) +
	294	& kint*QIRsQREF(m,iaer,radius_id+1)
	295	omegIRa(m,gausind,iaer)=
	296	& (1-kint)*omegaIR(m,iaer,radius_id) +
	297	& kint*omegaIR(m,iaer,radius_id+1)
	298	gIRa(m,gausind,iaer)=
	299	& (1-kint)*gIR(m,iaer,radius_id) +
	300	& kint*gIR(m,iaer,radius_id+1)
	301	ENDDO
	302	qrefIRa(gausind,iaer)=
	303	& (1-kint)*QREFir(iaer,radius_id) +
	304	& kint*QREFir(iaer,radius_id+1)
	305	omegrefIRa(gausind,iaer)=
	306	& (1-kint)*omegaREFir(iaer,radius_id) +
	307	& kint*omegaREFir(iaer,radius_id+1)
	308	ENDIF
	309	ENDDO
	310
[89]	311	DO gausind=1,ngau
[38]	312	drad=radiusr*radgaus(gausind)
	313	radGAUSb(gausind,iaer,idomain)=radiusm+drad
	314
	315	radius_id=minloc(abs(radiustab(iaer,idomain,:) -
	316	& (radiusm+drad)),DIM=1)
	317	IF ((radiustab(iaer,idomain,radius_id) -
	318	& (radiusm+drad)).GT.0) THEN
	319	radius_id=radius_id-1
	320	ENDIF
	321	IF (radius_id.GE.nsize(iaer,idomain)) THEN
	322	radius_id=nsize(iaer,idomain)-1
	323	kint = 1.
	324	ELSEIF (radius_id.LT.1) THEN
	325	radius_id=1
	326	kint = 0.
	327	ELSE
	328	kint = ( (radiusm+drad) -
	329	& radiustab(iaer,idomain,radius_id) ) /
	330	& ( radiustab(iaer,idomain,radius_id+1) -
	331	& radiustab(iaer,idomain,radius_id) )
	332	ENDIF
	333	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
	334	DO m=1,nsun
	335	qsqrefVISb(m,gausind,iaer)=
	336	& (1-kint)*QVISsQREF(m,iaer,radius_id) +
	337	& kint*QVISsQREF(m,iaer,radius_id+1)
	338	omegVISb(m,gausind,iaer)=
	339	& (1-kint)*omegaVIS(m,iaer,radius_id) +
	340	& kint*omegaVIS(m,iaer,radius_id+1)
	341	gVISb(m,gausind,iaer)=
	342	& (1-kint)*gVIS(m,iaer,radius_id) +
	343	& kint*gVIS(m,iaer,radius_id+1)
	344	ENDDO
	345	qrefVISb(gausind,iaer)=
	346	& (1-kint)*QREFvis(iaer,radius_id) +
	347	& kint*QREFvis(iaer,radius_id+1)
	348	omegrefVISb(gausind,iaer)=
	349	& (1-kint)*omegaREFvis(iaer,radius_id) +
	350	& kint*omegaREFvis(iaer,radius_id+1)
	351	ELSE ! INFRARED DOMAIN ----------------------------------
	352	DO m=1,nir
	353	qsqrefIRb(m,gausind,iaer)=
	354	& (1-kint)*QIRsQREF(m,iaer,radius_id) +
	355	& kint*QIRsQREF(m,iaer,radius_id+1)
	356	omegIRb(m,gausind,iaer)=
	357	& (1-kint)*omegaIR(m,iaer,radius_id) +
	358	& kint*omegaIR(m,iaer,radius_id+1)
	359	gIRb(m,gausind,iaer)=
	360	& (1-kint)*gIR(m,iaer,radius_id) +
	361	& kint*gIR(m,iaer,radius_id+1)
	362	ENDDO
	363	qrefIRb(gausind,iaer)=
	364	& (1-kint)*QREFir(iaer,radius_id) +
	365	& kint*QREFir(iaer,radius_id+1)
	366	omegrefIRb(gausind,iaer)=
	367	& (1-kint)*omegaREFir(iaer,radius_id) +
	368	& kint*omegaREFir(iaer,radius_id+1)
	369	ENDIF
	370	ENDDO
	371	c==================================================================
	372	IF ( .NOT.varyingnueff(iaer) ) THEN ! CONSTANT NUEFF
	373	c==================================================================
	374	c 2. Compute the scattering parameters using linear
	375	c interpolation over grain sizes and constant nueff
	376	c ---------------------------------------------------
	377
	378	DO lg = 1,nlayer
	379	DO ig = 1,ngrid
	380	c 2.1 Effective radius index and kx calculation
	381	var_tmp=reffrad(ig,lg,iaer)/refftabmin
	382	var_tmp=log(var_tmp)*3.
	383	var_tmp=var_tmp/logvratgrid+1.
	384	grid_i=floor(var_tmp)
	385	IF (grid_i.GE.refftabsize) THEN
	386	c WRITE(,) 'Warning: particle size in grid box #'
	387	c WRITE(,) ig,' is too large to be used by the '
	388	c WRITE(,) 'radiative transfer; please extend the '
	389	c WRITE(,) 'interpolation grid to larger grain sizes.'
	390	grid_i=refftabsize-1
	391	kx = 1.
	392	ELSEIF (grid_i.LT.1) THEN
	393	c WRITE(,) 'Warning: particle size in grid box #'
	394	c WRITE(,) ig,' is too small to be used by the '
	395	c WRITE(,) 'radiative transfer; please extend the '
	396	c WRITE(,) 'interpolation grid to smaller grain sizes.'
	397	grid_i=1
	398	kx = 0.
	399	ELSE
	400	kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) /
	401	& ( refftab(grid_i+1)-refftab(grid_i) )
	402	ENDIF
	403	c 2.3 Integration
	404	DO j=grid_i,grid_i+1
	405	c 2.3.1 Check if the calculation has been done
	406	IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
	407	c 2.3.2 Log-normal dist., r_g and sigma_g are defined
	408	c in [hansen_1974], "Light scattering in planetary
	409	c atmospheres", Space Science Reviews 16 527-610.
	410	c Here, sizedistk1=r_g and sizedistk2=sigma_g^2
	411	sizedistk2 = log(1.+nueffrad(1,1,iaer))
	412	sizedistk1 = exp(2.5*sizedistk2)
	413	sizedistk1 = refftab(j) / sizedistk1
	414
	415	normd(j,1,iaer,idomain) = 1e-30
[89]	416	DO gausind=1,ngau
[38]	417	drad=radiusr*radgaus(gausind)
	418	dista(j,1,iaer,idomain,gausind) =
	419	& LOG((radiusm-drad)/sizedistk1)
	420	dista(j,1,iaer,idomain,gausind) =
	421	& EXP(-dista(j,1,iaer,idomain,gausind) *
	422	& dista(j,1,iaer,idomain,gausind) *
	423	& 0.5e0/sizedistk2)/(radiusm-drad)
	424	dista(j,1,iaer,idomain,gausind) =
	425	& dista(j,1,iaer,idomain,gausind) /
	426	& (sqrt(2e0pisizedistk2))
	427
	428	distb(j,1,iaer,idomain,gausind) =
	429	& LOG((radiusm+drad)/sizedistk1)
	430	distb(j,1,iaer,idomain,gausind) =
	431	& EXP(-distb(j,1,iaer,idomain,gausind) *
	432	& distb(j,1,iaer,idomain,gausind) *
	433	& 0.5e0/sizedistk2)/(radiusm+drad)
	434	distb(j,1,iaer,idomain,gausind) =
	435	& distb(j,1,iaer,idomain,gausind) /
	436	& (sqrt(2e0pisizedistk2))
	437
	438	normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) +
	439	& weightgaus(gausind) *
	440	& (
	441	& distb(j,1,iaer,idomain,gausind) * pi *
	442	& radGAUSb(gausind,iaer,idomain) *
	443	& radGAUSb(gausind,iaer,idomain) +
	444	& dista(j,1,iaer,idomain,gausind) * pi *
	445	& radGAUSa(gausind,iaer,idomain) *
	446	& radGAUSa(gausind,iaer,idomain)
	447	& )
	448	ENDDO
	449	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
	450	c 2.3.3.vis Initialization
	451	qsqrefVISgrid(j,1,:,iaer)=0.
	452	qextVISgrid(j,1,:,iaer)=0.
	453	qscatVISgrid(j,1,:,iaer)=0.
	454	omegVISgrid(j,1,:,iaer)=0.
	455	gVISgrid(j,1,:,iaer)=0.
	456	qrefVISgrid(j,1,iaer)=0.
	457	qscatrefVISgrid(j,1,iaer)=0.
	458	omegrefVISgrid(j,1,iaer)=0.
	459
[89]	460	DO gausind=1,ngau
[38]	461	DO m=1,nsun
	462	c Convolution:
	463	qextVISgrid(j,1,m,iaer) =
	464	& qextVISgrid(j,1,m,iaer) +
	465	& weightgaus(gausind) *
	466	& (
	467	& qsqrefVISb(m,gausind,iaer) *
	468	& qrefVISb(gausind,iaer) *
	469	& piradGAUSb(gausind,iaer,idomain)
	470	& radGAUSb(gausind,iaer,idomain) *
	471	& distb(j,1,iaer,idomain,gausind) +
	472	& qsqrefVISa(m,gausind,iaer) *
	473	& qrefVISa(gausind,iaer) *
	474	& piradGAUSa(gausind,iaer,idomain)
	475	& radGAUSa(gausind,iaer,idomain) *
	476	& dista(j,1,iaer,idomain,gausind)
	477	& )
	478	qscatVISgrid(j,1,m,iaer) =
	479	& qscatVISgrid(j,1,m,iaer) +
	480	& weightgaus(gausind) *
	481	& (
	482	& omegVISb(m,gausind,iaer) *
	483	& qsqrefVISb(m,gausind,iaer) *
	484	& qrefVISb(gausind,iaer) *
	485	& piradGAUSb(gausind,iaer,idomain)
	486	& radGAUSb(gausind,iaer,idomain) *
	487	& distb(j,1,iaer,idomain,gausind) +
	488	& omegVISa(m,gausind,iaer) *
	489	& qsqrefVISa(m,gausind,iaer) *
	490	& qrefVISa(gausind,iaer) *
	491	& piradGAUSa(gausind,iaer,idomain)
	492	& radGAUSa(gausind,iaer,idomain) *
	493	& dista(j,1,iaer,idomain,gausind)
	494	& )
	495	gVISgrid(j,1,m,iaer) =
	496	& gVISgrid(j,1,m,iaer) +
	497	& weightgaus(gausind) *
	498	& (
	499	& omegVISb(m,gausind,iaer) *
	500	& qsqrefVISb(m,gausind,iaer) *
	501	& qrefVISb(gausind,iaer) *
	502	& gVISb(m,gausind,iaer) *
	503	& piradGAUSb(gausind,iaer,idomain)
	504	& radGAUSb(gausind,iaer,idomain) *
	505	& distb(j,1,iaer,idomain,gausind) +
	506	& omegVISa(m,gausind,iaer) *
	507	& qsqrefVISa(m,gausind,iaer) *
	508	& qrefVISa(gausind,iaer) *
	509	& gVISa(m,gausind,iaer) *
	510	& piradGAUSa(gausind,iaer,idomain)
	511	& radGAUSa(gausind,iaer,idomain) *
	512	& dista(j,1,iaer,idomain,gausind)
	513	& )
	514	ENDDO
	515	qrefVISgrid(j,1,iaer) =
	516	& qrefVISgrid(j,1,iaer) +
	517	& weightgaus(gausind) *
	518	& (
	519	& qrefVISb(gausind,iaer) *
	520	& piradGAUSb(gausind,iaer,idomain)
	521	& radGAUSb(gausind,iaer,idomain) *
	522	& distb(j,1,iaer,idomain,gausind) +
	523	& qrefVISa(gausind,iaer) *
	524	& piradGAUSa(gausind,iaer,idomain)
	525	& radGAUSa(gausind,iaer,idomain) *
	526	& dista(j,1,iaer,idomain,gausind)
	527	& )
	528	qscatrefVISgrid(j,1,iaer) =
	529	& qscatrefVISgrid(j,1,iaer) +
	530	& weightgaus(gausind) *
	531	& (
	532	& omegrefVISb(gausind,iaer) *
	533	& qrefVISb(gausind,iaer) *
	534	& piradGAUSb(gausind,iaer,idomain)
	535	& radGAUSb(gausind,iaer,idomain) *
	536	& distb(j,1,iaer,idomain,gausind) +
	537	& omegrefVISa(gausind,iaer) *
	538	& qrefVISa(gausind,iaer) *
	539	& piradGAUSa(gausind,iaer,idomain)
	540	& radGAUSa(gausind,iaer,idomain) *
	541	& dista(j,1,iaer,idomain,gausind)
	542	& )
	543	ENDDO
	544
	545	qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) /
	546	& normd(j,1,iaer,idomain)
	547	qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /
	548	& normd(j,1,iaer,idomain)
	549	omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /
	550	& qrefVISgrid(j,1,iaer)
	551	DO m=1,nsun
	552	qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /
	553	& normd(j,1,iaer,idomain)
	554	qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /
	555	& normd(j,1,iaer,idomain)
	556	gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) /
	557	& qscatVISgrid(j,1,m,iaer) /
	558	& normd(j,1,iaer,idomain)
	559
	560	qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /
	561	& qrefVISgrid(j,1,iaer)
	562	omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /
	563	& qextVISgrid(j,1,m,iaer)
	564	ENDDO
	565	ELSE ! INFRARED DOMAIN ----------
	566	c 2.3.3.ir Initialization
	567	qsqrefIRgrid(j,1,:,iaer)=0.
	568	qextIRgrid(j,1,:,iaer)=0.
	569	qscatIRgrid(j,1,:,iaer)=0.
	570	omegIRgrid(j,1,:,iaer)=0.
	571	gIRgrid(j,1,:,iaer)=0.
	572	qrefIRgrid(j,1,iaer)=0.
	573	qscatrefIRgrid(j,1,iaer)=0.
	574	omegrefIRgrid(j,1,iaer)=0.
	575
[89]	576	DO gausind=1,ngau
[38]	577	DO m=1,nir
	578	c Convolution:
	579	qextIRgrid(j,1,m,iaer) =
	580	& qextIRgrid(j,1,m,iaer) +
	581	& weightgaus(gausind) *
	582	& (
	583	& qsqrefIRb(m,gausind,iaer) *
	584	& qrefVISb(gausind,iaer) *
	585	& piradGAUSb(gausind,iaer,idomain)
	586	& radGAUSb(gausind,iaer,idomain) *
	587	& distb(j,1,iaer,idomain,gausind) +
	588	& qsqrefIRa(m,gausind,iaer) *
	589	& qrefVISa(gausind,iaer) *
	590	& piradGAUSa(gausind,iaer,idomain)
	591	& radGAUSa(gausind,iaer,idomain) *
	592	& dista(j,1,iaer,idomain,gausind)
	593	& )
	594	qscatIRgrid(j,1,m,iaer) =
	595	& qscatIRgrid(j,1,m,iaer) +
	596	& weightgaus(gausind) *
	597	& (
	598	& omegIRb(m,gausind,iaer) *
	599	& qsqrefIRb(m,gausind,iaer) *
	600	& qrefVISb(gausind,iaer) *
	601	& piradGAUSb(gausind,iaer,idomain)
	602	& radGAUSb(gausind,iaer,idomain) *
	603	& distb(j,1,iaer,idomain,gausind) +
	604	& omegIRa(m,gausind,iaer) *
	605	& qsqrefIRa(m,gausind,iaer) *
	606	& qrefVISa(gausind,iaer) *
	607	& piradGAUSa(gausind,iaer,idomain)
	608	& radGAUSa(gausind,iaer,idomain) *
	609	& dista(j,1,iaer,idomain,gausind)
	610	& )
	611	gIRgrid(j,1,m,iaer) =
	612	& gIRgrid(j,1,m,iaer) +
	613	& weightgaus(gausind) *
	614	& (
	615	& omegIRb(m,gausind,iaer) *
	616	& qsqrefIRb(m,gausind,iaer) *
	617	& qrefVISb(gausind,iaer) *
	618	& gIRb(m,gausind,iaer) *
	619	& piradGAUSb(gausind,iaer,idomain)
	620	& radGAUSb(gausind,iaer,idomain) *
	621	& distb(j,1,iaer,idomain,gausind) +
	622	& omegIRa(m,gausind,iaer) *
	623	& qsqrefIRa(m,gausind,iaer) *
	624	& qrefVISa(gausind,iaer) *
	625	& gIRa(m,gausind,iaer) *
	626	& piradGAUSa(gausind,iaer,idomain)
	627	& radGAUSa(gausind,iaer,idomain) *
	628	& dista(j,1,iaer,idomain,gausind)
	629	& )
	630	ENDDO
	631	qrefIRgrid(j,1,iaer) =
	632	& qrefIRgrid(j,1,iaer) +
	633	& weightgaus(gausind) *
	634	& (
	635	& qrefIRb(gausind,iaer) *
	636	& piradGAUSb(gausind,iaer,idomain)
	637	& radGAUSb(gausind,iaer,idomain) *
	638	& distb(j,1,iaer,idomain,gausind) +
	639	& qrefIRa(gausind,iaer) *
	640	& piradGAUSa(gausind,iaer,idomain)
	641	& radGAUSa(gausind,iaer,idomain) *
	642	& dista(j,1,iaer,idomain,gausind)
	643	& )
	644	qscatrefIRgrid(j,1,iaer) =
	645	& qscatrefIRgrid(j,1,iaer) +
	646	& weightgaus(gausind) *
	647	& (
	648	& omegrefIRb(gausind,iaer) *
	649	& qrefIRb(gausind,iaer) *
	650	& piradGAUSb(gausind,iaer,idomain)
	651	& radGAUSb(gausind,iaer,idomain) *
	652	& distb(j,1,iaer,idomain,gausind) +
	653	& omegrefIRa(gausind,iaer) *
	654	& qrefIRa(gausind,iaer) *
	655	& piradGAUSa(gausind,iaer,idomain)
	656	& radGAUSa(gausind,iaer,idomain) *
	657	& dista(j,1,iaer,idomain,gausind)
	658	& )
	659	ENDDO
	660
	661	qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) /
	662	& normd(j,1,iaer,idomain)
	663	qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /
	664	& normd(j,1,iaer,idomain)
	665	omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /
	666	& qrefIRgrid(j,1,iaer)
	667	DO m=1,nir
	668	qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /
	669	& normd(j,1,iaer,idomain)
	670	qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /
	671	& normd(j,1,iaer,idomain)
	672	gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) /
	673	& qscatIRgrid(j,1,m,iaer) /
	674	& normd(j,1,iaer,idomain)
	675
	676	qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /
	677	& qrefVISgrid(j,1,iaer)
	678	omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /
	679	& qextIRgrid(j,1,m,iaer)
	680	ENDDO
	681	ENDIF ! --------------------------
	682	checkgrid(j,1,iaer,idomain) = .true.
	683	ENDIF !checkgrid
	684	ENDDO !grid_i
	685	c 2.4 Linear interpolation
	686	k1 = (1-kx)
	687	k2 = kx
	688	IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
	689	DO m=1,nsun
	690	QVISsQREF3d(ig,lg,m,iaer) =
	691	& k1*qsqrefVISgrid(grid_i,1,m,iaer) +
	692	& k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
	693	omegaVIS3d(ig,lg,m,iaer) =
	694	& k1*omegVISgrid(grid_i,1,m,iaer) +
	695	& k2*omegVISgrid(grid_i+1,1,m,iaer)
	696	gVIS3d(ig,lg,m,iaer) =
	697	& k1*gVISgrid(grid_i,1,m,iaer) +
	698	& k2*gVISgrid(grid_i+1,1,m,iaer)
	699	ENDDO !nsun
	700	QREFvis3d(ig,lg,iaer) =
	701	& k1*qrefVISgrid(grid_i,1,iaer) +
	702	& k2*qrefVISgrid(grid_i+1,1,iaer)
	703	omegaREFvis3d(ig,lg,iaer) =
	704	& k1*omegrefVISgrid(grid_i,1,iaer) +
	705	& k2*omegrefVISgrid(grid_i+1,1,iaer)
	706	ELSE ! INFRARED -----------------------
	707	DO m=1,nir
	708	QIRsQREF3d(ig,lg,m,iaer) =
	709	& k1*qsqrefIRgrid(grid_i,1,m,iaer) +
	710	& k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
	711	omegaIR3d(ig,lg,m,iaer) =
	712	& k1*omegIRgrid(grid_i,1,m,iaer) +
	713	& k2*omegIRgrid(grid_i+1,1,m,iaer)
	714	gIR3d(ig,lg,m,iaer) =
	715	& k1*gIRgrid(grid_i,1,m,iaer) +
	716	& k2*gIRgrid(grid_i+1,1,m,iaer)
	717	ENDDO !nir
	718	QREFir3d(ig,lg,iaer) =
	719	& k1*qrefIRgrid(grid_i,1,iaer) +
	720	& k2*qrefIRgrid(grid_i+1,1,iaer)
	721	omegaREFir3d(ig,lg,iaer) =
	722	& k1*omegrefIRgrid(grid_i,1,iaer) +
	723	& k2*omegrefIRgrid(grid_i+1,1,iaer)
	724	ENDIF ! --------------------------------
	725	ENDDO !nlayermx
	726	ENDDO !ngridmx
	727	c==================================================================
	728	ELSE ! VARYING NUEFF
	729	c==================================================================
	730	c 3. Compute the scattering parameters in each grid box
	731	c using bilinear interpolation over the grain sizes
	732	c and the effective variances;
	733	c -----------------------------------------------------
	734
	735	DO lg = 1,nlayer
	736	DO ig = 1,ngrid
	737	c 3.1 Effective variance index and ky calculation
	738	var_tmp=log(nueffrad(ig,lg,iaer))
	739	grid_j=floor( (nuefftabsize-1)/(nuefftabmax-nuefftabmin)*
	740	& (var_tmp-nuefftabmin)+1. )
	741	IF (grid_j.GE.nuefftabsize) THEN
	742	c WRITE(,) 'Warning: effective variance '
	743	c WRITE(,) 'is too large to be used by the '
	744	c WRITE(,) 'radiative transfer; please extend the '
	745	c WRITE(,) 'interpolation grid to larger values.'
	746	grid_j=nuefftabsize-1
	747	ky = 1.
	748	ELSEIF (grid_j.LT.1) THEN
	749	c WRITE(,) 'Warning: effective variance '
	750	c WRITE(,) 'is too small to be used by the '
	751	c WRITE(,) 'radiative transfer; please extend the '
	752	c WRITE(,) 'interpolation grid to smaller values.'
	753	grid_j=1
	754	ky = 0.
	755	ELSE
	756	ky = ( nueffrad(ig,lg,iaer)-nuefftab(grid_j) ) /
	757	& ( nuefftab(grid_j+1)-nuefftab(grid_j) )
	758	ENDIF
	759	c 3.2 Effective radius index and kx calculation
	760	var_tmp=reffrad(ig,lg,iaer)/refftabmin
	761	var_tmp=log(var_tmp)*3.
	762	var_tmp=var_tmp/logvratgrid+1.
	763	grid_i=floor(var_tmp)
	764	IF (grid_i.GE.refftabsize) THEN
	765	c WRITE(,) 'Warning: particle size in grid box #'
	766	c WRITE(,) ig,' is too large to be used by the '
	767	c WRITE(,) 'radiative transfer; please extend the '
	768	c WRITE(,) 'interpolation grid to larger grain sizes.'
	769	grid_i=refftabsize-1
	770	kx = 1.
	771	ELSEIF (grid_i.LT.1) THEN
	772	c WRITE(,) 'Warning: particle size in grid box #'
	773	c WRITE(,) ig,' is too small to be used by the '
	774	c WRITE(,) 'radiative transfer; please extend the '
	775	c WRITE(,) 'interpolation grid to smaller grain sizes.'
	776	grid_i=1
	777	kx = 0.
	778	ELSE
	779	kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) /
	780	& ( refftab(grid_i+1)-refftab(grid_i) )
	781	ENDIF
	782	c 3.3 Integration
	783	DO j=grid_i,grid_i+1
	784	DO k=grid_j,grid_j+1
	785	c 3.3.1 Check if the calculation has been done
	786	IF (.NOT.checkgrid(j,k,iaer,idomain)) THEN
	787
	788	c 3.3.2 Log-normal dist., r_g and sigma_g are defined
	789	c in [hansen_1974], "Light scattering in planetary
	790	c atmospheres", Space Science Reviews 16 527-610.
	791	c Here, sizedistk1=r_g and sizedistk2=sigma_g^2
	792	sizedistk2 = log(1.+nuefftab(k))
	793	sizedistk1 = exp(2.5*sizedistk2)
	794	sizedistk1 = refftab(j) / sizedistk1
	795
	796	normd(j,k,iaer,idomain) = 1e-30
[89]	797	DO gausind=1,ngau
[38]	798	drad=radiusr*radgaus(gausind)
	799
	800	dista(j,k,iaer,idomain,gausind) =
	801	& LOG((radiusm-drad)/sizedistk1)
	802	dista(j,k,iaer,idomain,gausind) =
	803	& EXP(-dista(j,k,iaer,idomain,gausind) *
	804	& dista(j,k,iaer,idomain,gausind) *
	805	& 0.5e0/sizedistk2)/(radiusm-drad)
	806	dista(j,k,iaer,idomain,gausind) =
	807	& dista(j,k,iaer,idomain,gausind) /
	808	& (sqrt(2e0pisizedistk2))
	809
	810	distb(j,k,iaer,idomain,gausind) =
	811	& LOG((radiusm+drad)/sizedistk1)
	812	distb(j,k,iaer,idomain,gausind) =
	813	& EXP(-distb(j,k,iaer,idomain,gausind) *
	814	& distb(j,k,iaer,idomain,gausind) *
	815	& 0.5e0/sizedistk2)/(radiusm+drad)
	816	distb(j,k,iaer,idomain,gausind) =
	817	& distb(j,k,iaer,idomain,gausind) /
	818	& (sqrt(2e0pisizedistk2))
	819
	820	normd(j,k,iaer,idomain)=normd(j,k,iaer,idomain) +
	821	& weightgaus(gausind) *
	822	& (
	823	& distb(j,k,iaer,idomain,gausind) * pi *
	824	& radGAUSb(gausind,iaer,idomain) *
	825	& radGAUSb(gausind,iaer,idomain) +
	826	& dista(j,k,iaer,idomain,gausind) * pi *
	827	& radGAUSa(gausind,iaer,idomain) *
	828	& radGAUSa(gausind,iaer,idomain)
	829	& )
	830	ENDDO
	831	IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
	832	c 2.3.3.vis Initialization
	833	qsqrefVISgrid(j,k,:,iaer)=0.
	834	qextVISgrid(j,k,:,iaer)=0.
	835	qscatVISgrid(j,k,:,iaer)=0.
	836	omegVISgrid(j,k,:,iaer)=0.
	837	gVISgrid(j,k,:,iaer)=0.
	838	qrefVISgrid(j,k,iaer)=0.
	839	qscatrefVISgrid(j,k,iaer)=0.
	840	omegrefVISgrid(j,k,iaer)=0.
	841
[89]	842	DO gausind=1,ngau
[38]	843	DO m=1,nsun
	844	c Convolution:
	845	qextVISgrid(j,k,m,iaer) =
	846	& qextVISgrid(j,k,m,iaer) +
	847	& weightgaus(gausind) *
	848	& (
	849	& qsqrefVISb(m,gausind,iaer) *
	850	& qrefVISb(gausind,iaer) *
	851	& piradGAUSb(gausind,iaer,idomain)
	852	& radGAUSb(gausind,iaer,idomain) *
	853	& distb(j,k,iaer,idomain,gausind) +
	854	& qsqrefVISa(m,gausind,iaer) *
	855	& qrefVISa(gausind,iaer) *
	856	& piradGAUSa(gausind,iaer,idomain)
	857	& radGAUSa(gausind,iaer,idomain) *
	858	& dista(j,k,iaer,idomain,gausind)
	859	& )
	860	qscatVISgrid(j,k,m,iaer) =
	861	& qscatVISgrid(j,k,m,iaer) +
	862	& weightgaus(gausind) *
	863	& (
	864	& omegVISb(m,gausind,iaer) *
	865	& qsqrefVISb(m,gausind,iaer) *
	866	& qrefVISb(gausind,iaer) *
	867	& piradGAUSb(gausind,iaer,idomain)
	868	& radGAUSb(gausind,iaer,idomain) *
	869	& distb(j,k,iaer,idomain,gausind) +
	870	& omegVISa(m,gausind,iaer) *
	871	& qsqrefVISa(m,gausind,iaer) *
	872	& qrefVISa(gausind,iaer) *
	873	& piradGAUSa(gausind,iaer,idomain)
	874	& radGAUSa(gausind,iaer,idomain) *
	875	& dista(j,k,iaer,idomain,gausind)
	876	& )
	877	gVISgrid(j,k,m,iaer) =
	878	& gVISgrid(j,k,m,iaer) +
	879	& weightgaus(gausind) *
	880	& (
	881	& omegVISb(m,gausind,iaer) *
	882	& qsqrefVISb(m,gausind,iaer) *
	883	& qrefVISb(gausind,iaer) *
	884	& gVISb(m,gausind,iaer) *
	885	& piradGAUSb(gausind,iaer,idomain)
	886	& radGAUSb(gausind,iaer,idomain) *
	887	& distb(j,k,iaer,idomain,gausind) +
	888	& omegVISa(m,gausind,iaer) *
	889	& qsqrefVISa(m,gausind,iaer) *
	890	& qrefVISa(gausind,iaer) *
	891	& gVISa(m,gausind,iaer) *
	892	& piradGAUSa(gausind,iaer,idomain)
	893	& radGAUSa(gausind,iaer,idomain) *
	894	& dista(j,k,iaer,idomain,gausind)
	895	& )
	896	ENDDO
	897	qrefVISgrid(j,k,iaer) =
	898	& qrefVISgrid(j,k,iaer) +
	899	& weightgaus(gausind) *
	900	& (
	901	& qrefVISb(gausind,iaer) *
	902	& piradGAUSb(gausind,iaer,idomain)
	903	& radGAUSb(gausind,iaer,idomain) *
	904	& distb(j,k,iaer,idomain,gausind) +
	905	& qrefVISa(gausind,iaer) *
	906	& piradGAUSa(gausind,iaer,idomain)
	907	& radGAUSa(gausind,iaer,idomain) *
	908	& dista(j,k,iaer,idomain,gausind)
	909	& )
	910	qscatrefVISgrid(j,k,iaer) =
	911	& qscatrefVISgrid(j,k,iaer) +
	912	& weightgaus(gausind) *
	913	& (
	914	& omegrefVISb(gausind,iaer) *
	915	& qrefVISb(gausind,iaer) *
	916	& piradGAUSb(gausind,iaer,idomain)
	917	& radGAUSb(gausind,iaer,idomain) *
	918	& distb(j,k,iaer,idomain,gausind) +
	919	& omegrefVISa(gausind,iaer) *
	920	& qrefVISa(gausind,iaer) *
	921	& piradGAUSa(gausind,iaer,idomain)
	922	& radGAUSa(gausind,iaer,idomain) *
	923	& dista(j,k,iaer,idomain,gausind)
	924	& )
	925	ENDDO
	926	qrefVISgrid(j,k,iaer)=qrefVISgrid(j,k,iaer) /
	927	& normd(j,k,iaer,idomain)
	928	qscatrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) /
	929	& normd(j,k,iaer,idomain)
	930	omegrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) /
	931	& qrefVISgrid(j,k,iaer)
	932	DO m=1,nsun
	933	qextVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) /
	934	& normd(j,k,iaer,idomain)
	935	qscatVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) /
	936	& normd(j,k,iaer,idomain)
	937	gVISgrid(j,k,m,iaer)=gVISgrid(j,k,m,iaer) /
	938	& qscatVISgrid(j,k,m,iaer) /
	939	& normd(j,k,iaer,idomain)
	940
	941	qsqrefVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) /
	942	& qrefVISgrid(j,k,iaer)
	943	omegVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) /
	944	& qextVISgrid(j,k,m,iaer)
	945	ENDDO
	946	ELSE ! INFRARED DOMAIN ----------
	947	c 2.3.3.ir Initialization
	948	qsqrefIRgrid(j,k,:,iaer)=0.
	949	qextIRgrid(j,k,:,iaer)=0.
	950	qscatIRgrid(j,k,:,iaer)=0.
	951	omegIRgrid(j,k,:,iaer)=0.
	952	gIRgrid(j,k,:,iaer)=0.
	953	qrefIRgrid(j,k,iaer)=0.
	954	qscatrefIRgrid(j,k,iaer)=0.
	955	omegrefIRgrid(j,k,iaer)=0.
	956
[89]	957	DO gausind=1,ngau
[38]	958	DO m=1,nir
	959	c Convolution:
	960	qextIRgrid(j,k,m,iaer) =
	961	& qextIRgrid(j,k,m,iaer) +
	962	& weightgaus(gausind) *
	963	& (
	964	& qsqrefIRb(m,gausind,iaer) *
	965	& qrefVISb(gausind,iaer) *
	966	& piradGAUSb(gausind,iaer,idomain)
	967	& radGAUSb(gausind,iaer,idomain) *
	968	& distb(j,k,iaer,idomain,gausind) +
	969	& qsqrefIRa(m,gausind,iaer) *
	970	& qrefVISa(gausind,iaer) *
	971	& piradGAUSa(gausind,iaer,idomain)
	972	& radGAUSa(gausind,iaer,idomain) *
	973	& dista(j,k,iaer,idomain,gausind)
	974	& )
	975	qscatIRgrid(j,k,m,iaer) =
	976	& qscatIRgrid(j,k,m,iaer) +
	977	& weightgaus(gausind) *
	978	& (
	979	& omegIRb(m,gausind,iaer) *
	980	& qsqrefIRb(m,gausind,iaer) *
	981	& qrefVISb(gausind,iaer) *
	982	& piradGAUSb(gausind,iaer,idomain)
	983	& radGAUSb(gausind,iaer,idomain) *
	984	& distb(j,k,iaer,idomain,gausind) +
	985	& omegIRa(m,gausind,iaer) *
	986	& qsqrefIRa(m,gausind,iaer) *
	987	& qrefVISa(gausind,iaer) *
	988	& piradGAUSa(gausind,iaer,idomain)
	989	& radGAUSa(gausind,iaer,idomain) *
	990	& dista(j,k,iaer,idomain,gausind)
	991	& )
	992	gIRgrid(j,k,m,iaer) =
	993	& gIRgrid(j,k,m,iaer) +
	994	& weightgaus(gausind) *
	995	& (
	996	& omegIRb(m,gausind,iaer) *
	997	& qsqrefIRb(m,gausind,iaer) *
	998	& qrefVISb(gausind,iaer) *
	999	& gIRb(m,gausind,iaer) *
	1000	& piradGAUSb(gausind,iaer,idomain)
	1001	& radGAUSb(gausind,iaer,idomain) *
	1002	& distb(j,k,iaer,idomain,gausind) +
	1003	& omegIRa(m,gausind,iaer) *
	1004	& qsqrefIRa(m,gausind,iaer) *
	1005	& qrefVISa(gausind,iaer) *
	1006	& gIRa(m,gausind,iaer) *
	1007	& piradGAUSa(gausind,iaer,idomain)
	1008	& radGAUSa(gausind,iaer,idomain) *
	1009	& dista(j,k,iaer,idomain,gausind)
	1010	& )
	1011	ENDDO
	1012	qrefIRgrid(j,k,iaer) =
	1013	& qrefIRgrid(j,k,iaer) +
	1014	& weightgaus(gausind) *
	1015	& (
	1016	& qrefIRb(gausind,iaer) *
	1017	& piradGAUSb(gausind,iaer,idomain)
	1018	& radGAUSb(gausind,iaer,idomain) *
	1019	& distb(j,k,iaer,idomain,gausind) +
	1020	& qrefIRa(gausind,iaer) *
	1021	& piradGAUSa(gausind,iaer,idomain)
	1022	& radGAUSa(gausind,iaer,idomain) *
	1023	& dista(j,k,iaer,idomain,gausind)
	1024	& )
	1025	qscatrefIRgrid(j,k,iaer) =
	1026	& qscatrefIRgrid(j,k,iaer) +
	1027	& weightgaus(gausind) *
	1028	& (
	1029	& omegrefIRb(gausind,iaer) *
	1030	& qrefIRb(gausind,iaer) *
	1031	& piradGAUSb(gausind,iaer,idomain)
	1032	& radGAUSb(gausind,iaer,idomain) *
	1033	& distb(j,k,iaer,idomain,gausind) +
	1034	& omegrefIRa(gausind,iaer) *
	1035	& qrefIRa(gausind,iaer) *
	1036	& piradGAUSa(gausind,iaer,idomain)
	1037	& radGAUSa(gausind,iaer,idomain) *
	1038	& dista(j,k,iaer,idomain,gausind)
	1039	& )
	1040	ENDDO
	1041	qrefIRgrid(j,k,iaer)=qrefIRgrid(j,k,iaer) /
	1042	& normd(j,k,iaer,idomain)
	1043	qscatrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) /
	1044	& normd(j,k,iaer,idomain)
	1045	omegrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) /
	1046	& qrefIRgrid(j,k,iaer)
	1047	DO m=1,nir
	1048	qextIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) /
	1049	& normd(j,k,iaer,idomain)
	1050	qscatIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) /
	1051	& normd(j,k,iaer,idomain)
	1052	gIRgrid(j,k,m,iaer)=gIRgrid(j,k,m,iaer) /
	1053	& qscatIRgrid(j,k,m,iaer) /
	1054	& normd(j,k,iaer,idomain)
	1055
	1056	qsqrefIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) /
	1057	& qrefVISgrid(j,k,iaer)
	1058	omegIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) /
	1059	& qextIRgrid(j,k,m,iaer)
	1060	ENDDO
	1061	ENDIF ! --------------------------
	1062	checkgrid(j,k,iaer,idomain) = .true.
	1063	ENDIF !checkgrid
	1064	ENDDO !grid_j
	1065	ENDDO !grid_i
	1066	c 3.4 Bilinear interpolation
	1067	k1 = (1-kx)*(1-ky)
	1068	k2 = kx*(1-ky)
	1069	k3 = kx*ky
	1070	k4 = (1-kx)*ky
	1071	IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
	1072	DO m=1,nsun
	1073	QVISsQREF3d(ig,lg,m,iaer) =
	1074	& k1*qsqrefVISgrid(grid_i,grid_j,m,iaer) +
	1075	& k2*qsqrefVISgrid(grid_i+1,grid_j,m,iaer) +
	1076	& k3*qsqrefVISgrid(grid_i+1,grid_j+1,m,iaer) +
	1077	& k4*qsqrefVISgrid(grid_i,grid_j+1,m,iaer)
	1078	omegaVIS3d(ig,lg,m,iaer) =
	1079	& k1*omegVISgrid(grid_i,grid_j,m,iaer) +
	1080	& k2*omegVISgrid(grid_i+1,grid_j,m,iaer) +
	1081	& k3*omegVISgrid(grid_i+1,grid_j+1,m,iaer) +
	1082	& k4*omegVISgrid(grid_i,grid_j+1,m,iaer)
	1083	gVIS3d(ig,lg,m,iaer) =
	1084	& k1*gVISgrid(grid_i,grid_j,m,iaer) +
	1085	& k2*gVISgrid(grid_i+1,grid_j,m,iaer) +
	1086	& k3*gVISgrid(grid_i+1,grid_j+1,m,iaer) +
	1087	& k4*gVISgrid(grid_i,grid_j+1,m,iaer)
	1088	ENDDO !nsun
	1089	QREFvis3d(ig,lg,iaer) =
	1090	& k1*qrefVISgrid(grid_i,grid_j,iaer) +
	1091	& k2*qrefVISgrid(grid_i+1,grid_j,iaer) +
	1092	& k3*qrefVISgrid(grid_i+1,grid_j+1,iaer) +
	1093	& k4*qrefVISgrid(grid_i,grid_j+1,iaer)
	1094	omegaREFvis3d(ig,lg,iaer) =
	1095	& k1*omegrefVISgrid(grid_i,grid_j,iaer) +
	1096	& k2*omegrefVISgrid(grid_i+1,grid_j,iaer) +
	1097	& k3*omegrefVISgrid(grid_i+1,grid_j+1,iaer) +
	1098	& k4*omegrefVISgrid(grid_i,grid_j+1,iaer)
	1099	ELSE ! INFRARED -----------------------
	1100	DO m=1,nir
	1101	QIRsQREF3d(ig,lg,m,iaer) =
	1102	& k1*qsqrefIRgrid(grid_i,grid_j,m,iaer) +
	1103	& k2*qsqrefIRgrid(grid_i+1,grid_j,m,iaer) +
	1104	& k3*qsqrefIRgrid(grid_i+1,grid_j+1,m,iaer) +
	1105	& k4*qsqrefIRgrid(grid_i,grid_j+1,m,iaer)
	1106	omegaIR3d(ig,lg,m,iaer) =
	1107	& k1*omegIRgrid(grid_i,grid_j,m,iaer) +
	1108	& k2*omegIRgrid(grid_i+1,grid_j,m,iaer) +
	1109	& k3*omegIRgrid(grid_i+1,grid_j+1,m,iaer) +
	1110	& k4*omegIRgrid(grid_i,grid_j+1,m,iaer)
	1111	gIR3d(ig,lg,m,iaer) =
	1112	& k1*gIRgrid(grid_i,grid_j,m,iaer) +
	1113	& k2*gIRgrid(grid_i+1,grid_j,m,iaer) +
	1114	& k3*gIRgrid(grid_i+1,grid_j+1,m,iaer) +
	1115	& k4*gIRgrid(grid_i,grid_j+1,m,iaer)
	1116	ENDDO !nir
	1117	QREFir3d(ig,lg,iaer) =
	1118	& k1*qrefIRgrid(grid_i,grid_j,iaer) +
	1119	& k2*qrefIRgrid(grid_i+1,grid_j,iaer) +
	1120	& k3*qrefIRgrid(grid_i+1,grid_j+1,iaer) +
	1121	& k4*qrefIRgrid(grid_i,grid_j+1,iaer)
	1122	omegaREFir3d(ig,lg,iaer) =
	1123	& k1*omegrefIRgrid(grid_i,grid_j,iaer) +
	1124	& k2*omegrefIRgrid(grid_i+1,grid_j,iaer) +
	1125	& k3*omegrefIRgrid(grid_i+1,grid_j+1,iaer) +
	1126	& k4*omegrefIRgrid(grid_i,grid_j+1,iaer)
	1127	ENDIF ! --------------------------------
	1128	ENDDO !nlayermx
	1129	ENDDO !ngridmx
	1130
	1131	ENDIF ! varyingnueff
	1132	c==================================================================
	1133	ENDDO ! idomain
	1134
	1135	ENDIF ! nsize = 1
	1136
	1137	ENDDO ! iaer (loop on aerosol kind)
	1138
[89]	1139	c=====Radiative properties - TESTS=================================
	1140	IF (out_qwg) THEN
	1141	IF (ngrid.NE.1) THEN
	1142	DO out_nchannel = 1, 2
	1143	c -------------------------------------------------------------
	1144	DO lg = 1, nlayer
	1145	DO ig = 1, ngrid
	1146	out_qext(ig,lg) =
	1147	& QVISsQREF3d(ig,lg,out_nchannel,out_iaer)*
	1148	& QREFvis3d(ig,lg,out_iaer)
	1149	out_omeg(ig,lg) =
	1150	& omegaVIS3d(ig,lg,out_nchannel,out_iaer)
	1151	out_g(ig,lg) = gVIS3d(ig,lg,out_nchannel,out_iaer)
	1152	ENDDO ! ig
	1153	ENDDO ! lg
	1154	write(out_str(1:1),'(i1.1)') out_nchannel
	1155	call WRITEDIAGFI(ngrid,'qextvis'//out_str,"Ext.efficiency","",
	1156	& 3,out_qext)
	1157	call WRITEDIAGFI(ngrid,'omegvis'//out_str,"Sing.Scat.Alb.","",
	1158	& 3,out_omeg)
	1159	call WRITEDIAGFI(ngrid,'gvis'//out_str,"Asym.Factor","",
	1160	& 3,out_g)
	1161	c -------------------------------------------------------------
	1162	ENDDO ! out_nchannel
	1163	DO out_nchannel = 2, 4
	1164	c -------------------------------------------------------------
	1165	DO lg = 1, nlayer
	1166	DO ig = 1, ngrid
	1167	out_qext(ig,lg) =
	1168	& QIRsQREF3d(ig,lg,out_nchannel,out_iaer)*
	1169	& QREFir3d(ig,lg,out_iaer)
	1170	out_omeg(ig,lg) =
	1171	& omegaIR3d(ig,lg,out_nchannel,out_iaer)
	1172	out_g(ig,lg) = gIR3d(ig,lg,out_nchannel,out_iaer)
	1173	ENDDO ! ig
	1174	ENDDO ! lg
	1175	write(out_str(1:1),'(i1.1)') out_nchannel
[690]	1176	call WRITEDIAGFI(ngrid,'qextir'//out_str,"Ext.efficiency","",
[89]	1177	& 3,out_qext)
[690]	1178	call WRITEDIAGFI(ngrid,'omegir'//out_str,"Sing.Scat.Alb.","",
[89]	1179	& 3,out_omeg)
[690]	1180	call WRITEDIAGFI(ngrid,'gir'//out_str,"Asym.Factor","",
[89]	1181	& 3,out_g)
	1182	c -------------------------------------------------------------
	1183	ENDDO ! out_nchannel
	1184	call WRITEDIAGFI(ngrid,"omegvisref","Sing.Scat.Alb.","",
	1185	& 3,omegaREFvis3d(1,1,out_iaer))
	1186	call WRITEDIAGFI(ngrid,"omegirref","Sing.Scat.Alb.","",
	1187	& 3,omegaREFir3d(1,1,out_iaer))
	1188	ENDIF ! ngrid.EQ.1
	1189	ENDIF ! out_qwg
	1190	c==================================================================
	1191
[38]	1192	RETURN
	1193	END

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

Download in other formats:

Original Format