Context Navigation

aeroptproperties.F @ 3026

Last change on this file since 3026 was 2932, checked in by romain.vande, 21 months ago

Mars PCM:
Add a new routine to write output called write_output.
It needs to be imported (for example : use write_output_mod, only: write_output)
Then, it requires only 4 arguments : the name of the variable, its title, its units and the variable itself.
It detects the dimension of the variable and decide to output either in diagfi or diagsoil.
It is also compatible with XIOS (xml file needs to be adapted)
Writediagfi and writediagsoil routines are still available in case.
RV

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

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

Download in other formats:

Original Format