Context Navigation

aeroptproperties.F @ 2156

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

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

Download in other formats:

Original Format