SUBROUTINE updatereffrad(ngrid,nlayer,
     &                rdust,rice,nuice,
     &                reffrad,nueffrad,
     &                pq)

       IMPLICIT NONE
c=======================================================================
c   subject:
c   --------
c   Subroutine designed to update the aerosol size distribution used by
c     the radiative transfer scheme. This size distribution is assumed
c     to be a log-normal distribution, with effective radius "reffrad" and
c     variance "nueffrad".
c   At firstcall, "rice" and "nuice" are not known, because
c     the H2O ice microphysical scheme is called after the radiative
c     transfer in physiq.F. That's why we assess the size of the 
c     water-ice particles at firstcall (see part 1.2 below).
c
c   author:   
c   ------
c   J.-B. Madeleine (2009-2010)
c
c=======================================================================
c
c    Declarations :
c    -------------
c
#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "dimradmars.h"
#include "tracer.h"
#include "aerkind.h"
#include "yomaer.h"

c-----------------------------------------------------------------------
c     Inputs:
c     ------

      INTEGER ngrid,nlayer
c     Ice geometric mean radius (m)
      REAL :: rice(ngridmx,nlayermx)
c     Estimated effective variance of the size distribution (n.u.)
      REAL :: nuice(ngridmx,nlayermx)
c     Tracer mass mixing ratio (kg/kg)
      REAL pq(ngrid,nlayer,nqmx)
      real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
      real reffdust(ngridmx,nlayermx) ! Dust effective radius (m)
      real nueffdust(ngridmx,nlayermx) ! Dust effective variance

c     Outputs:
c     -------

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     Local variables:
c     ---------------

      INTEGER :: ig,l          ! 3D grid indices
      INTEGER :: iaer          ! Aerosol index

c     Number of cloud condensation nuclei near the surface
c     (only used at firstcall). This value is taken from 
c     Montmessin et al. 2004 JGR 109 E10004 p5 (2E6 part m-3), and
c     converted to part kg-1 using a typical atmospheric density.

      REAL, PARAMETER :: ccn0 = 1.3E8

      LOGICAL firstcall
      DATA firstcall/.true./
      SAVE firstcall

      REAL CBRT
      EXTERNAL CBRT

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

c==================================================================
c 1. Radius used in the physical subroutines
c   (but not in the rad. transfer)
c==================================================================

c     1.1 Dust particles
c     ------------------

      IF (doubleq.AND.active) THEN
        DO l=1,nlayer
          DO ig=1, ngrid
            reffdust(ig,l) = ref_r0 *
     &        CBRT(r3n_q*pq(ig,l,igcm_dust_mass)/
     &        max(pq(ig,l,igcm_dust_number),0.01))
            reffdust(ig,l)=min(max(reffdust(ig,l),1.e-10),500.e-6)
            nueffdust(ig,l) = exp(varian**2.)-1.
          ENDDO
        ENDDO
      ELSE
        DO l=1,nlayer
          DO ig=1, ngrid
            reffdust(ig,l) = 1.5E-6
            nueffdust(ig,l) = 0.3
          ENDDO
        ENDDO
      ENDIF

      DO l=1,nlayer
        DO ig=1, ngrid
c         Geometric mean radius = Effective radius / (1+nueff)^5/2
          rdust(ig,l) = reffdust(ig,l)/(1.+nueffdust(ig,l))**2.5
        ENDDO
      ENDDO

c     1.2 Water-ice particles
c     -----------------------

      IF (firstcall.AND.water.AND.activice) THEN
        DO l=1,nlayer
          DO ig=1,ngrid
            rice(ig,l) = max( CBRT(
     &        (pq(ig,l,igcm_h2o_ice)/rho_ice +
     &        ccn0*(4./3.)*pi*rdust(ig,l)**3.) /
     &        (ccn0*4./3.*pi)),rdust(ig,l) )
            nuice(ig,l) = nuice_ref
          ENDDO
        ENDDO
        firstcall = .false.
      ENDIF

c==================================================================
c 2. Radius used in the radiative transfer code (reffrad)
c==================================================================

      DO iaer = 1, naerkind ! Loop on aerosol kind
        aerkind: SELECT CASE (name_iaer(iaer))
c==================================================================
        CASE("dust_conrath") aerkind         ! Typical dust profile
c==================================================================
          DO l=1,nlayer
            DO ig=1,ngrid
              reffrad(ig,l,iaer) = reffdust(ig,l)
              nueffrad(ig,l,iaer) = nueffdust(ig,l)
            ENDDO
          ENDDO
c==================================================================
        CASE("dust_doubleq") aerkind! Two-moment scheme for dust
c==================================================================
          DO l=1,nlayer
            DO ig=1,ngrid
              reffrad(ig,l,iaer) = reffdust(ig,l)
              nueffrad(ig,l,iaer) = nueffdust(ig,l)
            ENDDO
          ENDDO
c==================================================================
        CASE("dust_submicron") aerkind   ! Small dust population
c==================================================================
          DO l=1,nlayer
            DO ig=1,ngrid
              reffrad(ig,l,iaer)=radius(igcm_dust_submicron)
              nueffrad(ig,l,iaer)=0.03
            ENDDO
          ENDDO     
c==================================================================
        CASE("h2o_ice") aerkind             ! Water ice crystals
c==================================================================
          DO l=1,nlayer
            DO ig=1,ngrid
              reffrad(ig,l,iaer)=rice(ig,l)*(1.+nuice_ref)
              nueffrad(ig,l,iaer)=nuice_ref
            ENDDO
          ENDDO
c==================================================================
        END SELECT aerkind
      ENDDO ! iaer (loop on aerosol kind)

      RETURN
      END