      SUBROUTINE aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
     &                            QVISsQREF3d,omegaVIS3d,gVIS3d,
     &                            QIRsQREF3d,omegaIR3d,gIR3d,
     &                            QREFvis3d,QREFir3d,
     &                            omegaREFvis3d,omegaREFir3d)
      IMPLICIT NONE
c     =============================================================
c     Aerosol Optical Properties
c
c     Description:
c       Compute the scattering parameters in each grid
c       box, depending on aerosol grain sizes. Log-normal size
c       distribution and Gauss-Legendre integration are used.

c     Parameters:
c       Don't forget to set the value of varyingnueff below; If
c       the effective variance of the distribution for the given
c       aerosol is considered homogeneous in the atmosphere, please
c       set varyingnueff(iaer) to .false. Resulting computational
c       time will be much better.

c     Authors: J.-B. Madeleine, F. Forget, F. Montmessin
c     =============================================================

#include "dimensions.h"
#include "dimphys.h"
#include "callkeys.h"
#include "dimradmars.h"
#include "yomaer.h"

c     Local variables 
c     ---------------

c     =============================================================
      LOGICAL, PARAMETER :: varyingnueff(naerkind) = .false.
c     =============================================================

c     Min. and max radius of the interpolation grid (in METERS)
      REAL, SAVE :: refftabmin(naerkind,2)
      REAL, SAVE :: refftabmax(naerkind,2)
c     Log of the min and max variance of the interpolation grid
      REAL, PARAMETER :: nuefftabmin = -4.6
      REAL, PARAMETER :: nuefftabmax = 0.
c     Number of effective radius of the interpolation grid
      INTEGER, PARAMETER :: refftabsize = 100
c     Number of effective variances of the interpolation grid
      INTEGER, PARAMETER :: nuefftabsize = 100
c     Interpolation grid indices (reff,nueff)
      INTEGER :: grid_i,grid_j
c     Intermediate variable
      REAL :: var_tmp,var3d_tmp(ngridmx,nlayermx)
c     Bilinear interpolation factors
      REAL :: kx,ky,k1,k2,k3,k4
c     Size distribution parameters
      REAL :: sizedistk1,sizedistk2
c     Pi!
      REAL,SAVE :: pi
c     Variables used by the Gauss-Legendre integration:
      INTEGER radius_id,gausind
      REAL kint
      REAL drad
      INTEGER, PARAMETER :: ngau = 10
      REAL weightgaus(ngau),radgaus(ngau)
      SAVE weightgaus,radgaus
c     DATA weightgaus/.2955242247,.2692667193,.2190863625,
c    &                .1494513491,.0666713443/
c     DATA radgaus/.1488743389,.4333953941,.6794095682,
c    &             .8650633666,.9739065285/
      DATA    radgaus/0.07652652113350,0.22778585114165,
     &                0.37370608871528,0.51086700195146,
     &                0.63605368072468,0.74633190646476,
     &                0.83911697181213,0.91223442826796,
     &                0.96397192726078,0.99312859919241/

      DATA weightgaus/0.15275338723120,0.14917298659407,
     &                0.14209610937519,0.13168863843930,
     &                0.11819453196154,0.10193011980823,
     &                0.08327674160932,0.06267204829828,
     &                0.04060142982019,0.01761400714091/

c     Indices
      INTEGER :: i,j,k,l,m,iaer,idomain
      INTEGER :: ig,lg,chg

c     Local saved variables
c     ---------------------

c     Radius axis of the interpolation grid
      REAL,SAVE :: refftab(refftabsize,naerkind,2)
c     Variance axis of the interpolation grid
      REAL,SAVE :: nuefftab(nuefftabsize,naerkind,2)
c     Volume ratio of the grid
      REAL,SAVE :: logvratgrid(naerkind,2)
c     Grid used to remember which calculation is done
      LOGICAL,SAVE :: checkgrid(refftabsize,nuefftabsize,naerkind,2)
     &                = .false.
c     Optical properties of the grid (VISIBLE)
      REAL,SAVE :: qsqrefVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
      REAL,SAVE :: qextVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
      REAL,SAVE :: qscatVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
      REAL,SAVE :: omegVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
      REAL,SAVE :: gVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
c     Optical properties of the grid (INFRARED)
      REAL,SAVE :: qsqrefIRgrid(refftabsize,nuefftabsize,nir,naerkind)
      REAL,SAVE :: qextIRgrid(refftabsize,nuefftabsize,nir,naerkind)
      REAL,SAVE :: qscatIRgrid(refftabsize,nuefftabsize,nir,naerkind)
      REAL,SAVE :: omegIRgrid(refftabsize,nuefftabsize,nir,naerkind)
      REAL,SAVE :: gIRgrid(refftabsize,nuefftabsize,nir,naerkind)
c     Optical properties of the grid (REFERENCE WAVELENGTHS)
      REAL,SAVE :: qrefVISgrid(refftabsize,nuefftabsize,naerkind)
      REAL,SAVE :: qscatrefVISgrid(refftabsize,nuefftabsize,naerkind)
      REAL,SAVE :: qrefIRgrid(refftabsize,nuefftabsize,naerkind)
      REAL,SAVE :: qscatrefIRgrid(refftabsize,nuefftabsize,naerkind)
      REAL,SAVE :: omegrefVISgrid(refftabsize,nuefftabsize,naerkind)
      REAL,SAVE :: omegrefIRgrid(refftabsize,nuefftabsize,naerkind)
c     Firstcall
      LOGICAL,SAVE :: firstcall = .true.
c     Variables used by the Gauss-Legendre integration:
      REAL,SAVE :: normd(refftabsize,nuefftabsize,naerkind,2)
      REAL,SAVE :: dista(refftabsize,nuefftabsize,naerkind,2,ngau)
      REAL,SAVE :: distb(refftabsize,nuefftabsize,naerkind,2,ngau)

      REAL,SAVE :: radGAUSa(ngau,naerkind,2)
      REAL,SAVE :: radGAUSb(ngau,naerkind,2)

      REAL,SAVE :: qsqrefVISa(nsun,ngau,naerkind)
      REAL,SAVE :: qrefVISa(ngau,naerkind)
      REAL,SAVE :: qsqrefVISb(nsun,ngau,naerkind)
      REAL,SAVE :: qrefVISb(ngau,naerkind)
      REAL,SAVE :: omegVISa(nsun,ngau,naerkind)
      REAL,SAVE :: omegrefVISa(ngau,naerkind)
      REAL,SAVE :: omegVISb(nsun,ngau,naerkind)
      REAL,SAVE :: omegrefVISb(ngau,naerkind)
      REAL,SAVE :: gVISa(nsun,ngau,naerkind)
      REAL,SAVE :: gVISb(nsun,ngau,naerkind)

      REAL,SAVE :: qsqrefIRa(nir,ngau,naerkind)
      REAL,SAVE :: qrefIRa(ngau,naerkind)
      REAL,SAVE :: qsqrefIRb(nir,ngau,naerkind)
      REAL,SAVE :: qrefIRb(ngau,naerkind)
      REAL,SAVE :: omegIRa(nir,ngau,naerkind)
      REAL,SAVE :: omegrefIRa(ngau,naerkind)
      REAL,SAVE :: omegIRb(nir,ngau,naerkind)
      REAL,SAVE :: omegrefIRb(ngau,naerkind)
      REAL,SAVE :: gIRa(nir,ngau,naerkind)
      REAL,SAVE :: gIRb(nir,ngau,naerkind)

      REAL :: radiusm
      REAL :: radiusr

c     Inputs
c     ------

      INTEGER :: ngrid,nlayer
c     Aerosol effective radius used for radiative transfer (meter)
      REAL :: reffrad(ngridmx,nlayermx,naerkind)
c     Aerosol effective variance used for radiative transfer (n.u.)
      REAL :: nueffrad(ngridmx,nlayermx,naerkind)

c     Outputs
c     -------

      REAL :: QVISsQREF3d(ngridmx,nlayermx,nsun,naerkind)
      REAL :: omegaVIS3d(ngridmx,nlayermx,nsun,naerkind)
      REAL :: gVIS3d(ngridmx,nlayermx,nsun,naerkind)

      REAL :: QIRsQREF3d(ngridmx,nlayermx,nir,naerkind)
      REAL :: omegaIR3d(ngridmx,nlayermx,nir,naerkind)
      REAL :: gIR3d(ngridmx,nlayermx,nir,naerkind)

      REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
      REAL :: QREFir3d(ngridmx,nlayermx,naerkind)

      REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
      REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind)

c     Tests
c     -----

      LOGICAL,SAVE :: out_qwg = .false.
      INTEGER, PARAMETER :: out_iaer = 2
      INTEGER :: out_ndim
      REAL :: out_qext(ngridmx,nlayermx)
      REAL :: out_omeg(ngridmx,nlayermx)
      REAL :: out_g(ngridmx,nlayermx)
      INTEGER :: out_nchannel
      CHARACTER*1 :: out_str

c     Creating the effective radius and variance grid
c     -----------------------------------------------
      IF (firstcall) THEN
c       0.1 Pi!
        pi = 2. * asin(1.e0)

        WRITE(*,*) "aeroptproperties: interpolation grid"
        DO iaer = 1, naerkind ! Loop on aerosol kind
          DO idomain = 1, 2 ! Loop on visible or infrared channel

c           0.2 Effective radius
            radiusm=
     &             0.5*(radiustab(iaer,idomain,nsize(iaer,idomain))+
     &                  radiustab(iaer,idomain,1))
            radiusr=
     &             0.5*(radiustab(iaer,idomain,nsize(iaer,idomain))-
     &                  radiustab(iaer,idomain,1))
            refftabmin(iaer,idomain) = 
     &        radiusm-radiusr*radgaus(ngau)
            refftabmax(iaer,idomain) = 
     &        radiustab(iaer,idomain,nsize(iaer,idomain))

            WRITE(*,*) "Scatterer: ",iaer
            WRITE(*,*) "Domain: ",idomain
            WRITE(*,*) "Min radius (m): ", refftabmin(iaer,idomain)
            WRITE(*,*) "Max radius (m): ", refftabmax(iaer,idomain)

            refftab(1,iaer,idomain) = 
     &        refftabmin(iaer,idomain)
            refftab(refftabsize,iaer,idomain) = 
     &        refftabmax(iaer,idomain)

            logvratgrid(iaer,idomain) = 
     &                    log(refftabmax(iaer,idomain)/
     &                      refftabmin(iaer,idomain)) /
     &                    float(refftabsize-1)*3.
            do i = 2, refftabsize-1
              refftab(i,iaer,idomain) = refftab(i-1,iaer,idomain)*
     &                             exp(1./3.*logvratgrid(iaer,idomain))
            enddo

c           0.3 Effective variance
            do i = 0, nuefftabsize-1
              nuefftab(i+1,iaer,idomain) = exp( nuefftabmin +
     &                 i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
            enddo

          ENDDO
        ENDDO
        firstcall = .false.
      ENDIF

      DO iaer = 1, naerkind ! Loop on aerosol kind
        IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN
c==================================================================
c       If there is one single particle size, optical
c         properties of the considered aerosol are homogeneous
          DO lg = 1, nlayer
            DO ig = 1, ngrid
              DO chg = 1, nsun
                QVISsQREF3d(ig,lg,chg,iaer)=QVISsQREF(chg,iaer,1)
                omegaVIS3d(ig,lg,chg,iaer)=omegaVIS(chg,iaer,1)
                gVIS3d(ig,lg,chg,iaer)=gVIS(chg,iaer,1)
              ENDDO
              DO chg = 1, nir
                QIRsQREF3d(ig,lg,chg,iaer)=QIRsQREF(chg,iaer,1)
                omegaIR3d(ig,lg,chg,iaer)=omegaIR(chg,iaer,1)
                gIR3d(ig,lg,chg,iaer)=gIR(chg,iaer,1)
              ENDDO
              QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1)
              QREFir3d(ig,lg,iaer)=QREFir(iaer,1)
              omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1)
              omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1)
            ENDDO
          ENDDO
        ELSE ! Varying effective radius and variance
      DO idomain = 1, 2 ! Loop on visible or infrared channel
c==================================================================

c       1.1 Radius middle point and range for Gauss integration
        radiusm=
     &         0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) +
     &              radiustab(iaer,idomain,1))
        radiusr=
     &         0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) -
     &              radiustab(iaer,idomain,1))

c       1.2 Interpolating data at the Gauss quadrature points:
        DO gausind=1,ngau
          drad=radiusr*radgaus(gausind)
          radGAUSa(gausind,iaer,idomain)=radiusm-drad

          radius_id=minloc(abs(radiustab(iaer,idomain,:) -
     &                         (radiusm-drad)),DIM=1)
          IF ((radiustab(iaer,idomain,radius_id) - 
     &         (radiusm-drad)).GT.0) THEN
            radius_id=radius_id-1
          ENDIF
          IF (radius_id.GE.nsize(iaer,idomain)) THEN
            radius_id=nsize(iaer,idomain)-1
            kint = 1.
          ELSEIF (radius_id.LT.1) THEN
            radius_id=1
            kint = 0.
          ELSE
          kint = ( (radiusm-drad) - 
     &             radiustab(iaer,idomain,radius_id) ) /
     &           ( radiustab(iaer,idomain,radius_id+1) - 
     &             radiustab(iaer,idomain,radius_id) )
          ENDIF
          IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
            DO m=1,nsun
            qsqrefVISa(m,gausind,iaer)=
     &              (1-kint)*QVISsQREF(m,iaer,radius_id) + 
     &              kint*QVISsQREF(m,iaer,radius_id+1)
            omegVISa(m,gausind,iaer)=
     &              (1-kint)*omegaVIS(m,iaer,radius_id) + 
     &              kint*omegaVIS(m,iaer,radius_id+1)
            gVISa(m,gausind,iaer)=
     &              (1-kint)*gVIS(m,iaer,radius_id) + 
     &              kint*gVIS(m,iaer,radius_id+1)
            ENDDO
            qrefVISa(gausind,iaer)=
     &              (1-kint)*QREFvis(iaer,radius_id) + 
     &              kint*QREFvis(iaer,radius_id+1)
            omegrefVISa(gausind,iaer)=
     &              (1-kint)*omegaREFvis(iaer,radius_id) + 
     &              kint*omegaREFvis(iaer,radius_id+1)
          ELSE ! INFRARED DOMAIN ----------------------------------
            DO m=1,nir
            qsqrefIRa(m,gausind,iaer)=
     &              (1-kint)*QIRsQREF(m,iaer,radius_id) + 
     &              kint*QIRsQREF(m,iaer,radius_id+1)
            omegIRa(m,gausind,iaer)=
     &              (1-kint)*omegaIR(m,iaer,radius_id) + 
     &              kint*omegaIR(m,iaer,radius_id+1)
            gIRa(m,gausind,iaer)=
     &              (1-kint)*gIR(m,iaer,radius_id) + 
     &              kint*gIR(m,iaer,radius_id+1)
            ENDDO
            qrefIRa(gausind,iaer)=
     &              (1-kint)*QREFir(iaer,radius_id) + 
     &              kint*QREFir(iaer,radius_id+1)
            omegrefIRa(gausind,iaer)=
     &              (1-kint)*omegaREFir(iaer,radius_id) + 
     &              kint*omegaREFir(iaer,radius_id+1)
          ENDIF
        ENDDO

        DO gausind=1,ngau
          drad=radiusr*radgaus(gausind)
          radGAUSb(gausind,iaer,idomain)=radiusm+drad

          radius_id=minloc(abs(radiustab(iaer,idomain,:) - 
     &                         (radiusm+drad)),DIM=1)
          IF ((radiustab(iaer,idomain,radius_id) - 
     &         (radiusm+drad)).GT.0) THEN
            radius_id=radius_id-1
          ENDIF
          IF (radius_id.GE.nsize(iaer,idomain)) THEN
            radius_id=nsize(iaer,idomain)-1
            kint = 1.
          ELSEIF (radius_id.LT.1) THEN
            radius_id=1
            kint = 0.
          ELSE
            kint = ( (radiusm+drad) -
     &               radiustab(iaer,idomain,radius_id) ) /
     &             ( radiustab(iaer,idomain,radius_id+1) - 
     &               radiustab(iaer,idomain,radius_id) )
          ENDIF
          IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
            DO m=1,nsun
            qsqrefVISb(m,gausind,iaer)=
     &              (1-kint)*QVISsQREF(m,iaer,radius_id) + 
     &              kint*QVISsQREF(m,iaer,radius_id+1)
            omegVISb(m,gausind,iaer)=
     &              (1-kint)*omegaVIS(m,iaer,radius_id) + 
     &              kint*omegaVIS(m,iaer,radius_id+1)
            gVISb(m,gausind,iaer)=
     &              (1-kint)*gVIS(m,iaer,radius_id) + 
     &              kint*gVIS(m,iaer,radius_id+1)
            ENDDO
            qrefVISb(gausind,iaer)=
     &              (1-kint)*QREFvis(iaer,radius_id) + 
     &              kint*QREFvis(iaer,radius_id+1)
            omegrefVISb(gausind,iaer)=
     &              (1-kint)*omegaREFvis(iaer,radius_id) + 
     &              kint*omegaREFvis(iaer,radius_id+1)
          ELSE ! INFRARED DOMAIN ----------------------------------
            DO m=1,nir
            qsqrefIRb(m,gausind,iaer)=
     &              (1-kint)*QIRsQREF(m,iaer,radius_id) + 
     &              kint*QIRsQREF(m,iaer,radius_id+1)
            omegIRb(m,gausind,iaer)=
     &              (1-kint)*omegaIR(m,iaer,radius_id) + 
     &              kint*omegaIR(m,iaer,radius_id+1)
            gIRb(m,gausind,iaer)=
     &              (1-kint)*gIR(m,iaer,radius_id) + 
     &              kint*gIR(m,iaer,radius_id+1)
            ENDDO
            qrefIRb(gausind,iaer)=
     &              (1-kint)*QREFir(iaer,radius_id) + 
     &              kint*QREFir(iaer,radius_id+1)
            omegrefIRb(gausind,iaer)=
     &              (1-kint)*omegaREFir(iaer,radius_id) + 
     &              kint*omegaREFir(iaer,radius_id+1)
          ENDIF
        ENDDO
c==================================================================
      IF ( .NOT.varyingnueff(iaer) ) THEN          ! CONSTANT NUEFF
c==================================================================
c     2. Compute the scattering parameters using linear
c       interpolation over grain sizes and constant nueff
c     ---------------------------------------------------

      DO lg = 1,nlayer
        DO ig = 1,ngrid
c         2.1 Effective radius index and kx calculation
          var_tmp=reffrad(ig,lg,iaer)/refftabmin(iaer,idomain)
          var_tmp=log(var_tmp)*3.
          var_tmp=var_tmp/logvratgrid(iaer,idomain)+1.
          grid_i=floor(var_tmp)
          IF (grid_i.GE.refftabsize) THEN
c           WRITE(*,*) 'Warning: particle size in grid box #'
c           WRITE(*,*) ig,' is too large to be used by the '
c           WRITE(*,*) 'radiative transfer; please extend the '
c           WRITE(*,*) 'interpolation grid to larger grain sizes.'
            grid_i=refftabsize-1
            kx = 1.
          ELSEIF (grid_i.LT.1) THEN
c           WRITE(*,*) 'Warning: particle size in grid box #'
c           WRITE(*,*) ig,' is too small to be used by the '
c           WRITE(*,*) 'radiative transfer; please extend the '
c           WRITE(*,*) 'interpolation grid to smaller grain sizes.'
            grid_i=1
            kx = 0.
          ELSE
            kx = ( reffrad(ig,lg,iaer)-
     &             refftab(grid_i,iaer,idomain) ) /
     &           ( refftab(grid_i+1,iaer,idomain)-
     &             refftab(grid_i,iaer,idomain) )
          ENDIF
c         2.3 Integration
          DO j=grid_i,grid_i+1
c             2.3.1 Check if the calculation has been done
              IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
c               2.3.2 Log-normal dist., r_g and sigma_g are defined
c                 in [hansen_1974], "Light scattering in planetary
c                 atmospheres", Space Science Reviews 16 527-610.
c                 Here, sizedistk1=r_g and sizedistk2=sigma_g^2
                sizedistk2 = log(1.+nueffrad(1,1,iaer))
                sizedistk1 = exp(2.5*sizedistk2)
                sizedistk1 = refftab(j,iaer,idomain) / sizedistk1

                normd(j,1,iaer,idomain) = 1e-30
                DO gausind=1,ngau
                  drad=radiusr*radgaus(gausind)
                  dista(j,1,iaer,idomain,gausind) =
     &              LOG((radiusm-drad)/sizedistk1)
                  dista(j,1,iaer,idomain,gausind) =
     &              EXP(-dista(j,1,iaer,idomain,gausind) * 
     &              dista(j,1,iaer,idomain,gausind) *
     &              0.5e0/sizedistk2)/(radiusm-drad)
                  dista(j,1,iaer,idomain,gausind) =
     &              dista(j,1,iaer,idomain,gausind) / 
     &              (sqrt(2e0*pi*sizedistk2))

                  distb(j,1,iaer,idomain,gausind) = 
     &              LOG((radiusm+drad)/sizedistk1)
                  distb(j,1,iaer,idomain,gausind) = 
     &              EXP(-distb(j,1,iaer,idomain,gausind) * 
     &              distb(j,1,iaer,idomain,gausind) *
     &              0.5e0/sizedistk2)/(radiusm+drad)
                  distb(j,1,iaer,idomain,gausind) = 
     &              distb(j,1,iaer,idomain,gausind) / 
     &              (sqrt(2e0*pi*sizedistk2))

                  normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) +
     &              weightgaus(gausind) *
     &              (
     &              distb(j,1,iaer,idomain,gausind) * pi * 
     &              radGAUSb(gausind,iaer,idomain) * 
     &              radGAUSb(gausind,iaer,idomain) +
     &              dista(j,1,iaer,idomain,gausind) * pi * 
     &              radGAUSa(gausind,iaer,idomain) * 
     &              radGAUSa(gausind,iaer,idomain)
     &              )
                ENDDO
                IF (normd(j,1,iaer,idomain).EQ.1e-30) THEN
                  WRITE(*,*)"normd:", normd(j,1,iaer,idomain)
                  WRITE(*,*)"Risk of division by 0 (aeroptproperties.F)"
                  WRITE(*,*)"Check the size of the interpolation grid."
                  CALL ABORT
                ENDIF
                IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
c                 2.3.3.vis Initialization
                  qsqrefVISgrid(j,1,:,iaer)=0.
                  qextVISgrid(j,1,:,iaer)=0.
                  qscatVISgrid(j,1,:,iaer)=0.
                  omegVISgrid(j,1,:,iaer)=0.
                  gVISgrid(j,1,:,iaer)=0.
                  qrefVISgrid(j,1,iaer)=0.
                  qscatrefVISgrid(j,1,iaer)=0.
                  omegrefVISgrid(j,1,iaer)=0.

                  DO gausind=1,ngau
                    DO m=1,nsun
c                     Convolution:
                      qextVISgrid(j,1,m,iaer) = 
     &                  qextVISgrid(j,1,m,iaer) + 
     &                  weightgaus(gausind) *
     &                  (
     &                  qsqrefVISb(m,gausind,iaer) * 
     &                  qrefVISb(gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) * 
     &                  radGAUSb(gausind,iaer,idomain) * 
     &                  distb(j,1,iaer,idomain,gausind) +
     &                  qsqrefVISa(m,gausind,iaer) * 
     &                  qrefVISa(gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) * 
     &                  radGAUSa(gausind,iaer,idomain) * 
     &                  dista(j,1,iaer,idomain,gausind)
     &                  )
                      qscatVISgrid(j,1,m,iaer) = 
     &                  qscatVISgrid(j,1,m,iaer) + 
     &                  weightgaus(gausind) *
     &                  (
     &                  omegVISb(m,gausind,iaer) *
     &                  qsqrefVISb(m,gausind,iaer) *
     &                  qrefVISb(gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) *
     &                  radGAUSb(gausind,iaer,idomain) *
     &                  distb(j,1,iaer,idomain,gausind) +
     &                  omegVISa(m,gausind,iaer) *
     &                  qsqrefVISa(m,gausind,iaer) *
     &                  qrefVISa(gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) *
     &                  radGAUSa(gausind,iaer,idomain) *
     &                  dista(j,1,iaer,idomain,gausind)
     &                  )
                      gVISgrid(j,1,m,iaer) =
     &                  gVISgrid(j,1,m,iaer) +
     &                  weightgaus(gausind) *
     &                  (
     &                  omegVISb(m,gausind,iaer) *
     &                  qsqrefVISb(m,gausind,iaer) *
     &                  qrefVISb(gausind,iaer) *
     &                  gVISb(m,gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) *
     &                  radGAUSb(gausind,iaer,idomain) *
     &                  distb(j,1,iaer,idomain,gausind) +
     &                  omegVISa(m,gausind,iaer) *
     &                  qsqrefVISa(m,gausind,iaer) *
     &                  qrefVISa(gausind,iaer) *
     &                  gVISa(m,gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) *
     &                  radGAUSa(gausind,iaer,idomain) *
     &                  dista(j,1,iaer,idomain,gausind)
     &                  )
                    ENDDO
                    qrefVISgrid(j,1,iaer) = 
     &                qrefVISgrid(j,1,iaer) + 
     &                weightgaus(gausind) *
     &                (
     &                qrefVISb(gausind,iaer) *
     &                pi*radGAUSb(gausind,iaer,idomain) * 
     &                radGAUSb(gausind,iaer,idomain) * 
     &                distb(j,1,iaer,idomain,gausind) +
     &                qrefVISa(gausind,iaer) *
     &                pi*radGAUSa(gausind,iaer,idomain) * 
     &                radGAUSa(gausind,iaer,idomain) * 
     &                dista(j,1,iaer,idomain,gausind)
     &                )
                    qscatrefVISgrid(j,1,iaer) = 
     &                qscatrefVISgrid(j,1,iaer) + 
     &                weightgaus(gausind) *
     &                (
     &                omegrefVISb(gausind,iaer) *
     &                qrefVISb(gausind,iaer) *
     &                pi*radGAUSb(gausind,iaer,idomain) * 
     &                radGAUSb(gausind,iaer,idomain) * 
     &                distb(j,1,iaer,idomain,gausind) +
     &                omegrefVISa(gausind,iaer) *
     &                qrefVISa(gausind,iaer) *
     &                pi*radGAUSa(gausind,iaer,idomain) * 
     &                radGAUSa(gausind,iaer,idomain) * 
     &                dista(j,1,iaer,idomain,gausind)
     &                )
                  ENDDO

                  qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) /
     &                          normd(j,1,iaer,idomain)
                  qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /
     &                          normd(j,1,iaer,idomain)
                  omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /
     &                         qrefVISgrid(j,1,iaer)
                  DO m=1,nsun
                    qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /
     &                          normd(j,1,iaer,idomain)
                    qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /
     &                          normd(j,1,iaer,idomain)
                    gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) / 
     &                          qscatVISgrid(j,1,m,iaer) /
     &                          normd(j,1,iaer,idomain)

                    qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /
     &                          qrefVISgrid(j,1,iaer)
                    omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / 
     &                          qextVISgrid(j,1,m,iaer)
                  ENDDO
                ELSE                   ! INFRARED DOMAIN ----------
c                 2.3.3.ir Initialization
                  qsqrefIRgrid(j,1,:,iaer)=0.
                  qextIRgrid(j,1,:,iaer)=0.
                  qscatIRgrid(j,1,:,iaer)=0.
                  omegIRgrid(j,1,:,iaer)=0.
                  gIRgrid(j,1,:,iaer)=0.
                  qrefIRgrid(j,1,iaer)=0.
                  qscatrefIRgrid(j,1,iaer)=0.
                  omegrefIRgrid(j,1,iaer)=0.

                  DO gausind=1,ngau
                    DO m=1,nir
c                     Convolution:
                      qextIRgrid(j,1,m,iaer) = 
     &                  qextIRgrid(j,1,m,iaer) + 
     &                  weightgaus(gausind) *
     &                  (
     &                  qsqrefIRb(m,gausind,iaer) * 
     &                  qrefVISb(gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) * 
     &                  radGAUSb(gausind,iaer,idomain) * 
     &                  distb(j,1,iaer,idomain,gausind) +
     &                  qsqrefIRa(m,gausind,iaer) * 
     &                  qrefVISa(gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) * 
     &                  radGAUSa(gausind,iaer,idomain) * 
     &                  dista(j,1,iaer,idomain,gausind)
     &                  )
                      qscatIRgrid(j,1,m,iaer) = 
     &                  qscatIRgrid(j,1,m,iaer) + 
     &                  weightgaus(gausind) *
     &                  (
     &                  omegIRb(m,gausind,iaer) *
     &                  qsqrefIRb(m,gausind,iaer) *
     &                  qrefVISb(gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) *
     &                  radGAUSb(gausind,iaer,idomain) *
     &                  distb(j,1,iaer,idomain,gausind) +
     &                  omegIRa(m,gausind,iaer) *
     &                  qsqrefIRa(m,gausind,iaer) *
     &                  qrefVISa(gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) *
     &                  radGAUSa(gausind,iaer,idomain) *
     &                  dista(j,1,iaer,idomain,gausind)
     &                  )
                      gIRgrid(j,1,m,iaer) =
     &                  gIRgrid(j,1,m,iaer) +
     &                  weightgaus(gausind) *
     &                  (
     &                  omegIRb(m,gausind,iaer) *
     &                  qsqrefIRb(m,gausind,iaer) *
     &                  qrefVISb(gausind,iaer) *
     &                  gIRb(m,gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) *
     &                  radGAUSb(gausind,iaer,idomain) *
     &                  distb(j,1,iaer,idomain,gausind) +
     &                  omegIRa(m,gausind,iaer) *
     &                  qsqrefIRa(m,gausind,iaer) *
     &                  qrefVISa(gausind,iaer) *
     &                  gIRa(m,gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) *
     &                  radGAUSa(gausind,iaer,idomain) *
     &                  dista(j,1,iaer,idomain,gausind)
     &                  )
                    ENDDO
                    qrefIRgrid(j,1,iaer) = 
     &                qrefIRgrid(j,1,iaer) + 
     &                weightgaus(gausind) *
     &                (
     &                qrefIRb(gausind,iaer) *
     &                pi*radGAUSb(gausind,iaer,idomain) * 
     &                radGAUSb(gausind,iaer,idomain) * 
     &                distb(j,1,iaer,idomain,gausind) +
     &                qrefIRa(gausind,iaer) *
     &                pi*radGAUSa(gausind,iaer,idomain) * 
     &                radGAUSa(gausind,iaer,idomain) * 
     &                dista(j,1,iaer,idomain,gausind)
     &                )
                    qscatrefIRgrid(j,1,iaer) = 
     &                qscatrefIRgrid(j,1,iaer) + 
     &                weightgaus(gausind) *
     &                (
     &                omegrefIRb(gausind,iaer) *
     &                qrefIRb(gausind,iaer) *
     &                pi*radGAUSb(gausind,iaer,idomain) * 
     &                radGAUSb(gausind,iaer,idomain) * 
     &                distb(j,1,iaer,idomain,gausind) +
     &                omegrefIRa(gausind,iaer) *
     &                qrefIRa(gausind,iaer) *
     &                pi*radGAUSa(gausind,iaer,idomain) * 
     &                radGAUSa(gausind,iaer,idomain) * 
     &                dista(j,1,iaer,idomain,gausind)
     &                )
                  ENDDO

                  qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) /
     &                          normd(j,1,iaer,idomain)
                  qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /
     &                          normd(j,1,iaer,idomain)
                  omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /
     &                         qrefIRgrid(j,1,iaer)
                  DO m=1,nir
                    qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /
     &                          normd(j,1,iaer,idomain)
                    qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /
     &                          normd(j,1,iaer,idomain)
                    gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) / 
     &                          qscatIRgrid(j,1,m,iaer) /
     &                          normd(j,1,iaer,idomain)

                    qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /
     &                          qrefVISgrid(j,1,iaer)
                    omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / 
     &                          qextIRgrid(j,1,m,iaer)
                  ENDDO
                ENDIF                  ! --------------------------
                checkgrid(j,1,iaer,idomain) = .true.
              ENDIF !checkgrid
          ENDDO !grid_i
c         2.4 Linear interpolation
          k1 = (1-kx)
          k2 = kx
          IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
          DO m=1,nsun
            QVISsQREF3d(ig,lg,m,iaer) =
     &                  k1*qsqrefVISgrid(grid_i,1,m,iaer) +
     &                  k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
            omegaVIS3d(ig,lg,m,iaer) =
     &                  k1*omegVISgrid(grid_i,1,m,iaer) +
     &                  k2*omegVISgrid(grid_i+1,1,m,iaer)
            gVIS3d(ig,lg,m,iaer) =
     &                  k1*gVISgrid(grid_i,1,m,iaer) +
     &                  k2*gVISgrid(grid_i+1,1,m,iaer)
          ENDDO !nsun
          QREFvis3d(ig,lg,iaer) =
     &                  k1*qrefVISgrid(grid_i,1,iaer) +
     &                  k2*qrefVISgrid(grid_i+1,1,iaer)
          omegaREFvis3d(ig,lg,iaer) =
     &                  k1*omegrefVISgrid(grid_i,1,iaer) +
     &                  k2*omegrefVISgrid(grid_i+1,1,iaer)
          ELSE                   ! INFRARED -----------------------
          DO m=1,nir
            QIRsQREF3d(ig,lg,m,iaer) =
     &                  k1*qsqrefIRgrid(grid_i,1,m,iaer) +
     &                  k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
            omegaIR3d(ig,lg,m,iaer) =
     &                  k1*omegIRgrid(grid_i,1,m,iaer) +
     &                  k2*omegIRgrid(grid_i+1,1,m,iaer)
            gIR3d(ig,lg,m,iaer) =
     &                  k1*gIRgrid(grid_i,1,m,iaer) +
     &                  k2*gIRgrid(grid_i+1,1,m,iaer)
          ENDDO !nir
          QREFir3d(ig,lg,iaer) =
     &                  k1*qrefIRgrid(grid_i,1,iaer) +
     &                  k2*qrefIRgrid(grid_i+1,1,iaer)
          omegaREFir3d(ig,lg,iaer) =
     &                  k1*omegrefIRgrid(grid_i,1,iaer) +
     &                  k2*omegrefIRgrid(grid_i+1,1,iaer)
          ENDIF                  ! --------------------------------
        ENDDO !nlayermx
      ENDDO !ngridmx
c==================================================================
      ELSE                                          ! VARYING NUEFF
c==================================================================
c     3. Compute the scattering parameters in each grid box
c       using bilinear interpolation over the grain sizes
c       and the effective variances;
c     -----------------------------------------------------

      DO lg = 1,nlayer
        DO ig = 1,ngrid
c         3.1 Effective variance index and ky calculation
          var_tmp=log(nueffrad(ig,lg,iaer))
          grid_j=floor( (nuefftabsize-1)/(nuefftabmax-nuefftabmin)*
     &       (var_tmp-nuefftabmin)+1. )
          IF (grid_j.GE.nuefftabsize) THEN
c           WRITE(*,*) 'Warning: effective variance '
c           WRITE(*,*) 'is too large to be used by the '
c           WRITE(*,*) 'radiative transfer; please extend the '
c           WRITE(*,*) 'interpolation grid to larger values.'
            grid_j=nuefftabsize-1
            ky = 1.
          ELSEIF (grid_j.LT.1) THEN
c           WRITE(*,*) 'Warning: effective variance '
c           WRITE(*,*) 'is too small to be used by the '
c           WRITE(*,*) 'radiative transfer; please extend the '
c           WRITE(*,*) 'interpolation grid to smaller values.'
            grid_j=1
            ky = 0.
          ELSE
            ky = ( nueffrad(ig,lg,iaer)-
     &             nuefftab(grid_j,iaer,idomain) ) /
     &           ( nuefftab(grid_j+1,iaer,idomain)-
     &             nuefftab(grid_j,iaer,idomain) )
          ENDIF
c         3.2 Effective radius index and kx calculation
          var_tmp=reffrad(ig,lg,iaer)/refftabmin(iaer,idomain)
          var_tmp=log(var_tmp)*3.
          var_tmp=var_tmp/logvratgrid(iaer,idomain)+1.
          grid_i=floor(var_tmp)
          IF (grid_i.GE.refftabsize) THEN
c           WRITE(*,*) 'Warning: particle size in grid box #'
c           WRITE(*,*) ig,' is too large to be used by the '
c           WRITE(*,*) 'radiative transfer; please extend the '
c           WRITE(*,*) 'interpolation grid to larger grain sizes.'
            grid_i=refftabsize-1
            kx = 1.
          ELSEIF (grid_i.LT.1) THEN
c           WRITE(*,*) 'Warning: particle size in grid box #'
c           WRITE(*,*) ig,' is too small to be used by the '
c           WRITE(*,*) 'radiative transfer; please extend the '
c           WRITE(*,*) 'interpolation grid to smaller grain sizes.'
            grid_i=1
            kx = 0.
          ELSE
            kx = ( reffrad(ig,lg,iaer)-
     &             refftab(grid_i,iaer,idomain) ) /
     &           ( refftab(grid_i+1,iaer,idomain)-
     &             refftab(grid_i,iaer,idomain) )
          ENDIF
c         3.3 Integration
          DO j=grid_i,grid_i+1
            DO k=grid_j,grid_j+1
c             3.3.1 Check if the calculation has been done
              IF (.NOT.checkgrid(j,k,iaer,idomain)) THEN

c               3.3.2 Log-normal dist., r_g and sigma_g are defined
c                 in [hansen_1974], "Light scattering in planetary
c                 atmospheres", Space Science Reviews 16 527-610.
c                 Here, sizedistk1=r_g and sizedistk2=sigma_g^2
                sizedistk2 = log(1.+nuefftab(k,iaer,idomain))
                sizedistk1 = exp(2.5*sizedistk2)
                sizedistk1 = refftab(j,iaer,idomain) / sizedistk1

                normd(j,k,iaer,idomain) = 1e-30
                DO gausind=1,ngau
                  drad=radiusr*radgaus(gausind)

                  dista(j,k,iaer,idomain,gausind) =
     &              LOG((radiusm-drad)/sizedistk1)
                  dista(j,k,iaer,idomain,gausind) =
     &              EXP(-dista(j,k,iaer,idomain,gausind) *
     &              dista(j,k,iaer,idomain,gausind) *
     &              0.5e0/sizedistk2)/(radiusm-drad)
                  dista(j,k,iaer,idomain,gausind) =
     &              dista(j,k,iaer,idomain,gausind) /
     &              (sqrt(2e0*pi*sizedistk2))

                  distb(j,k,iaer,idomain,gausind) =
     &              LOG((radiusm+drad)/sizedistk1)
                  distb(j,k,iaer,idomain,gausind) =
     &              EXP(-distb(j,k,iaer,idomain,gausind) *
     &              distb(j,k,iaer,idomain,gausind) *
     &              0.5e0/sizedistk2)/(radiusm+drad)
                  distb(j,k,iaer,idomain,gausind) =
     &              distb(j,k,iaer,idomain,gausind) /
     &              (sqrt(2e0*pi*sizedistk2))

                  normd(j,k,iaer,idomain)=normd(j,k,iaer,idomain) +
     &              weightgaus(gausind) *
     &              (
     &              distb(j,k,iaer,idomain,gausind) * pi *
     &              radGAUSb(gausind,iaer,idomain) *
     &              radGAUSb(gausind,iaer,idomain) +
     &              dista(j,k,iaer,idomain,gausind) * pi *
     &              radGAUSa(gausind,iaer,idomain) *
     &              radGAUSa(gausind,iaer,idomain)
     &              )
                ENDDO
                IF (normd(j,k,iaer,idomain).EQ.1e-30) THEN
                  WRITE(*,*)"normd:", normd(j,k,iaer,idomain)
                  WRITE(*,*)"Risk of division by 0 (aeroptproperties.F)"
                  WRITE(*,*)"Check the size of the interpolation grid."
                  CALL ABORT
                ENDIF
                IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
c                 2.3.3.vis Initialization
                  qsqrefVISgrid(j,k,:,iaer)=0.
                  qextVISgrid(j,k,:,iaer)=0.
                  qscatVISgrid(j,k,:,iaer)=0.
                  omegVISgrid(j,k,:,iaer)=0.
                  gVISgrid(j,k,:,iaer)=0.
                  qrefVISgrid(j,k,iaer)=0.
                  qscatrefVISgrid(j,k,iaer)=0.
                  omegrefVISgrid(j,k,iaer)=0.

                  DO gausind=1,ngau
                    DO m=1,nsun
c                     Convolution:
                      qextVISgrid(j,k,m,iaer) =
     &                  qextVISgrid(j,k,m,iaer) +
     &                  weightgaus(gausind) *
     &                  (
     &                  qsqrefVISb(m,gausind,iaer) *
     &                  qrefVISb(gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) *
     &                  radGAUSb(gausind,iaer,idomain) *
     &                  distb(j,k,iaer,idomain,gausind) +
     &                  qsqrefVISa(m,gausind,iaer) *
     &                  qrefVISa(gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) *
     &                  radGAUSa(gausind,iaer,idomain) *
     &                  dista(j,k,iaer,idomain,gausind)
     &                  )
                      qscatVISgrid(j,k,m,iaer) =
     &                  qscatVISgrid(j,k,m,iaer) +
     &                  weightgaus(gausind) *
     &                  (
     &                  omegVISb(m,gausind,iaer) *
     &                  qsqrefVISb(m,gausind,iaer) *
     &                  qrefVISb(gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) *
     &                  radGAUSb(gausind,iaer,idomain) *
     &                  distb(j,k,iaer,idomain,gausind) +
     &                  omegVISa(m,gausind,iaer) *
     &                  qsqrefVISa(m,gausind,iaer) *
     &                  qrefVISa(gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) *
     &                  radGAUSa(gausind,iaer,idomain) *
     &                  dista(j,k,iaer,idomain,gausind)
     &                  )
                      gVISgrid(j,k,m,iaer) =
     &                  gVISgrid(j,k,m,iaer) +
     &                  weightgaus(gausind) *
     &                  (
     &                  omegVISb(m,gausind,iaer) *
     &                  qsqrefVISb(m,gausind,iaer) *
     &                  qrefVISb(gausind,iaer) *
     &                  gVISb(m,gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) *
     &                  radGAUSb(gausind,iaer,idomain) *
     &                  distb(j,k,iaer,idomain,gausind) +
     &                  omegVISa(m,gausind,iaer) *
     &                  qsqrefVISa(m,gausind,iaer) *
     &                  qrefVISa(gausind,iaer) *
     &                  gVISa(m,gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) *
     &                  radGAUSa(gausind,iaer,idomain) *
     &                  dista(j,k,iaer,idomain,gausind)
     &                  )
                    ENDDO
                    qrefVISgrid(j,k,iaer) =
     &                qrefVISgrid(j,k,iaer) +
     &                weightgaus(gausind) *
     &                (
     &                qrefVISb(gausind,iaer) *
     &                pi*radGAUSb(gausind,iaer,idomain) *
     &                radGAUSb(gausind,iaer,idomain) *
     &                distb(j,k,iaer,idomain,gausind) +
     &                qrefVISa(gausind,iaer) *
     &                pi*radGAUSa(gausind,iaer,idomain) *
     &                radGAUSa(gausind,iaer,idomain) *
     &                dista(j,k,iaer,idomain,gausind)
     &                )
                    qscatrefVISgrid(j,k,iaer) =
     &                qscatrefVISgrid(j,k,iaer) +
     &                weightgaus(gausind) *
     &                (
     &                omegrefVISb(gausind,iaer) *
     &                qrefVISb(gausind,iaer) *
     &                pi*radGAUSb(gausind,iaer,idomain) *
     &                radGAUSb(gausind,iaer,idomain) *
     &                distb(j,k,iaer,idomain,gausind) +
     &                omegrefVISa(gausind,iaer) *
     &                qrefVISa(gausind,iaer) *
     &                pi*radGAUSa(gausind,iaer,idomain) *
     &                radGAUSa(gausind,iaer,idomain) *
     &                dista(j,k,iaer,idomain,gausind)
     &                )
                  ENDDO
                  qrefVISgrid(j,k,iaer)=qrefVISgrid(j,k,iaer) /
     &                          normd(j,k,iaer,idomain)
                  qscatrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) /
     &                          normd(j,k,iaer,idomain)
                  omegrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) /
     &                         qrefVISgrid(j,k,iaer)
                  DO m=1,nsun
                    qextVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) /
     &                          normd(j,k,iaer,idomain)
                    qscatVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) /
     &                          normd(j,k,iaer,idomain)
                    gVISgrid(j,k,m,iaer)=gVISgrid(j,k,m,iaer) /
     &                          qscatVISgrid(j,k,m,iaer) /
     &                          normd(j,k,iaer,idomain)

                    qsqrefVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) /
     &                          qrefVISgrid(j,k,iaer)
                    omegVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) /
     &                          qextVISgrid(j,k,m,iaer)
                  ENDDO
                ELSE                   ! INFRARED DOMAIN ----------
c                 2.3.3.ir Initialization
                  qsqrefIRgrid(j,k,:,iaer)=0.
                  qextIRgrid(j,k,:,iaer)=0.
                  qscatIRgrid(j,k,:,iaer)=0.
                  omegIRgrid(j,k,:,iaer)=0.
                  gIRgrid(j,k,:,iaer)=0.
                  qrefIRgrid(j,k,iaer)=0.
                  qscatrefIRgrid(j,k,iaer)=0.
                  omegrefIRgrid(j,k,iaer)=0.

                  DO gausind=1,ngau
                    DO m=1,nir
c                     Convolution:
                      qextIRgrid(j,k,m,iaer) =
     &                  qextIRgrid(j,k,m,iaer) +
     &                  weightgaus(gausind) *
     &                  (
     &                  qsqrefIRb(m,gausind,iaer) *
     &                  qrefVISb(gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) *
     &                  radGAUSb(gausind,iaer,idomain) *
     &                  distb(j,k,iaer,idomain,gausind) +
     &                  qsqrefIRa(m,gausind,iaer) *
     &                  qrefVISa(gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) *
     &                  radGAUSa(gausind,iaer,idomain) *
     &                  dista(j,k,iaer,idomain,gausind)
     &                  )
                      qscatIRgrid(j,k,m,iaer) =
     &                  qscatIRgrid(j,k,m,iaer) +
     &                  weightgaus(gausind) *
     &                  (
     &                  omegIRb(m,gausind,iaer) *
     &                  qsqrefIRb(m,gausind,iaer) *
     &                  qrefVISb(gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) *
     &                  radGAUSb(gausind,iaer,idomain) *
     &                  distb(j,k,iaer,idomain,gausind) +
     &                  omegIRa(m,gausind,iaer) *
     &                  qsqrefIRa(m,gausind,iaer) *
     &                  qrefVISa(gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) *
     &                  radGAUSa(gausind,iaer,idomain) *
     &                  dista(j,k,iaer,idomain,gausind)
     &                  )
                      gIRgrid(j,k,m,iaer) =
     &                  gIRgrid(j,k,m,iaer) +
     &                  weightgaus(gausind) *
     &                  (
     &                  omegIRb(m,gausind,iaer) *
     &                  qsqrefIRb(m,gausind,iaer) *
     &                  qrefVISb(gausind,iaer) *
     &                  gIRb(m,gausind,iaer) *
     &                  pi*radGAUSb(gausind,iaer,idomain) *
     &                  radGAUSb(gausind,iaer,idomain) *
     &                  distb(j,k,iaer,idomain,gausind) +
     &                  omegIRa(m,gausind,iaer) *
     &                  qsqrefIRa(m,gausind,iaer) *
     &                  qrefVISa(gausind,iaer) *
     &                  gIRa(m,gausind,iaer) *
     &                  pi*radGAUSa(gausind,iaer,idomain) *
     &                  radGAUSa(gausind,iaer,idomain) *
     &                  dista(j,k,iaer,idomain,gausind)
     &                  )
                    ENDDO
                    qrefIRgrid(j,k,iaer) =
     &                qrefIRgrid(j,k,iaer) +
     &                weightgaus(gausind) *
     &                (
     &                qrefIRb(gausind,iaer) *
     &                pi*radGAUSb(gausind,iaer,idomain) *
     &                radGAUSb(gausind,iaer,idomain) *
     &                distb(j,k,iaer,idomain,gausind) +
     &                qrefIRa(gausind,iaer) *
     &                pi*radGAUSa(gausind,iaer,idomain) *
     &                radGAUSa(gausind,iaer,idomain) *
     &                dista(j,k,iaer,idomain,gausind)
     &                )
                    qscatrefIRgrid(j,k,iaer) =
     &                qscatrefIRgrid(j,k,iaer) +
     &                weightgaus(gausind) *
     &                (
     &                omegrefIRb(gausind,iaer) *
     &                qrefIRb(gausind,iaer) *
     &                pi*radGAUSb(gausind,iaer,idomain) *
     &                radGAUSb(gausind,iaer,idomain) *
     &                distb(j,k,iaer,idomain,gausind) +
     &                omegrefIRa(gausind,iaer) *
     &                qrefIRa(gausind,iaer) *
     &                pi*radGAUSa(gausind,iaer,idomain) *
     &                radGAUSa(gausind,iaer,idomain) *
     &                dista(j,k,iaer,idomain,gausind)
     &                )
                  ENDDO
                  qrefIRgrid(j,k,iaer)=qrefIRgrid(j,k,iaer) /
     &                          normd(j,k,iaer,idomain)
                  qscatrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) /
     &                          normd(j,k,iaer,idomain)
                  omegrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) /
     &                         qrefIRgrid(j,k,iaer)
                  DO m=1,nir
                    qextIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) /
     &                          normd(j,k,iaer,idomain)
                    qscatIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) /
     &                          normd(j,k,iaer,idomain)
                    gIRgrid(j,k,m,iaer)=gIRgrid(j,k,m,iaer) /
     &                          qscatIRgrid(j,k,m,iaer) /
     &                          normd(j,k,iaer,idomain)

                    qsqrefIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) /
     &                          qrefVISgrid(j,k,iaer)
                    omegIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) /
     &                          qextIRgrid(j,k,m,iaer)
                  ENDDO
                ENDIF                  ! --------------------------
                checkgrid(j,k,iaer,idomain) = .true.
              ENDIF !checkgrid
            ENDDO !grid_j
          ENDDO !grid_i
c         3.4 Bilinear interpolation
          k1 = (1-kx)*(1-ky)
          k2 = kx*(1-ky)
          k3 = kx*ky
          k4 = (1-kx)*ky
          IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
          DO m=1,nsun
            QVISsQREF3d(ig,lg,m,iaer) = 
     &                  k1*qsqrefVISgrid(grid_i,grid_j,m,iaer) +
     &                  k2*qsqrefVISgrid(grid_i+1,grid_j,m,iaer) +
     &                  k3*qsqrefVISgrid(grid_i+1,grid_j+1,m,iaer) +
     &                  k4*qsqrefVISgrid(grid_i,grid_j+1,m,iaer)
            omegaVIS3d(ig,lg,m,iaer) = 
     &                  k1*omegVISgrid(grid_i,grid_j,m,iaer) +
     &                  k2*omegVISgrid(grid_i+1,grid_j,m,iaer) +
     &                  k3*omegVISgrid(grid_i+1,grid_j+1,m,iaer) +
     &                  k4*omegVISgrid(grid_i,grid_j+1,m,iaer)
            gVIS3d(ig,lg,m,iaer) = 
     &                  k1*gVISgrid(grid_i,grid_j,m,iaer) +
     &                  k2*gVISgrid(grid_i+1,grid_j,m,iaer) +
     &                  k3*gVISgrid(grid_i+1,grid_j+1,m,iaer) +
     &                  k4*gVISgrid(grid_i,grid_j+1,m,iaer)
          ENDDO !nsun
          QREFvis3d(ig,lg,iaer) = 
     &                  k1*qrefVISgrid(grid_i,grid_j,iaer) +
     &                  k2*qrefVISgrid(grid_i+1,grid_j,iaer) +
     &                  k3*qrefVISgrid(grid_i+1,grid_j+1,iaer) +
     &                  k4*qrefVISgrid(grid_i,grid_j+1,iaer)
          omegaREFvis3d(ig,lg,iaer) = 
     &                  k1*omegrefVISgrid(grid_i,grid_j,iaer) +
     &                  k2*omegrefVISgrid(grid_i+1,grid_j,iaer) +
     &                  k3*omegrefVISgrid(grid_i+1,grid_j+1,iaer) +
     &                  k4*omegrefVISgrid(grid_i,grid_j+1,iaer)
          ELSE                   ! INFRARED -----------------------
          DO m=1,nir
            QIRsQREF3d(ig,lg,m,iaer) = 
     &                  k1*qsqrefIRgrid(grid_i,grid_j,m,iaer) + 
     &                  k2*qsqrefIRgrid(grid_i+1,grid_j,m,iaer) +
     &                  k3*qsqrefIRgrid(grid_i+1,grid_j+1,m,iaer) +
     &                  k4*qsqrefIRgrid(grid_i,grid_j+1,m,iaer)
            omegaIR3d(ig,lg,m,iaer) = 
     &                  k1*omegIRgrid(grid_i,grid_j,m,iaer) + 
     &                  k2*omegIRgrid(grid_i+1,grid_j,m,iaer) +
     &                  k3*omegIRgrid(grid_i+1,grid_j+1,m,iaer) +
     &                  k4*omegIRgrid(grid_i,grid_j+1,m,iaer)
            gIR3d(ig,lg,m,iaer) =
     &                  k1*gIRgrid(grid_i,grid_j,m,iaer) + 
     &                  k2*gIRgrid(grid_i+1,grid_j,m,iaer) +
     &                  k3*gIRgrid(grid_i+1,grid_j+1,m,iaer) +
     &                  k4*gIRgrid(grid_i,grid_j+1,m,iaer)
          ENDDO !nir
          QREFir3d(ig,lg,iaer) =
     &                  k1*qrefIRgrid(grid_i,grid_j,iaer) +
     &                  k2*qrefIRgrid(grid_i+1,grid_j,iaer) +
     &                  k3*qrefIRgrid(grid_i+1,grid_j+1,iaer) +
     &                  k4*qrefIRgrid(grid_i,grid_j+1,iaer)
          omegaREFir3d(ig,lg,iaer) =
     &                  k1*omegrefIRgrid(grid_i,grid_j,iaer) +
     &                  k2*omegrefIRgrid(grid_i+1,grid_j,iaer) +
     &                  k3*omegrefIRgrid(grid_i+1,grid_j+1,iaer) +
     &                  k4*omegrefIRgrid(grid_i,grid_j+1,iaer)
          ENDIF                  ! --------------------------------
        ENDDO !nlayermx
      ENDDO !ngridmx

      ENDIF ! varyingnueff
c==================================================================
      ENDDO ! idomain

      ENDIF ! nsize = 1

      ENDDO ! iaer (loop on aerosol kind)

c=====Radiative properties - TESTS=================================
      IF (out_qwg) THEN
c     -------------------------------------------------------------
        IF (ngrid.NE.1) THEN
          out_ndim = 3
        ELSE
          out_ndim = 1
        ENDIF
c     -------------------------------------------------------------
        DO out_nchannel = 1, 2
c     -------------------------------------------------------------
          DO lg = 1, nlayer
            DO ig = 1, ngrid
              out_qext(ig,lg) = 
     &           QVISsQREF3d(ig,lg,out_nchannel,out_iaer)*
     &           QREFvis3d(ig,lg,out_iaer)
              out_omeg(ig,lg) =
     &           omegaVIS3d(ig,lg,out_nchannel,out_iaer)
              out_g(ig,lg)    = gVIS3d(ig,lg,out_nchannel,out_iaer)
            ENDDO ! ig
          ENDDO ! lg
          write(out_str(1:1),'(i1.1)') out_nchannel
          call WRITEDIAGFI(ngrid,'qextvis'//out_str,"Ext.efficiency","",
     &                     out_ndim,out_qext)
          call WRITEDIAGFI(ngrid,'omegvis'//out_str,"Sing.Scat.Alb.","",
     &                     out_ndim,out_omeg)
          call WRITEDIAGFI(ngrid,'gvis'//out_str,"Asym.Factor","",
     &                     out_ndim,out_g)
c     -------------------------------------------------------------
        ENDDO ! out_nchannel 
        DO out_nchannel = 2, 4
c     -------------------------------------------------------------
          DO lg = 1, nlayer
            DO ig = 1, ngrid
              out_qext(ig,lg) =
     &          QIRsQREF3d(ig,lg,out_nchannel,out_iaer)*
     &          QREFir3d(ig,lg,out_iaer)
              out_omeg(ig,lg) =
     &          omegaIR3d(ig,lg,out_nchannel,out_iaer)
              out_g(ig,lg)    = gIR3d(ig,lg,out_nchannel,out_iaer)
            ENDDO ! ig
          ENDDO ! lg
          write(out_str(1:1),'(i1.1)') out_nchannel
       call WRITEDIAGFI(ngrid,'qextir'//out_str,"Ext.efficiency","",
     &                     out_ndim,out_qext)
       call WRITEDIAGFI(ngrid,'omegir'//out_str,"Sing.Scat.Alb.","",
     &                     out_ndim,out_omeg)
       call WRITEDIAGFI(ngrid,'gir'//out_str,"Asym.Factor","",
     &                     out_ndim,out_g)
c     -------------------------------------------------------------
        ENDDO ! out_nchannel 
        call WRITEDIAGFI(ngrid,"omegvisref","Sing.Scat.Alb.","",
     &                   out_ndim,omegaREFvis3d(1,1,out_iaer))
        call WRITEDIAGFI(ngrid,"omegirref","Sing.Scat.Alb.","",
     &                   out_ndim,omegaREFir3d(1,1,out_iaer))
      ENDIF ! out_qwg
c==================================================================

      RETURN
      END
