source: trunk/LMDZ.MARS/libf/phymars/updatereffrad.F @ 627

      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)

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

c      LOGICAL firstcall
c      DATA firstcall/.true./
c      SAVE firstcall

      REAL CBRT
      EXTERNAL CBRT

      REAL,SAVE :: nueffdust(ngridmx,nlayermx) ! Dust effective variance

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

c==================================================================

c      IF (firstcall) THEN
c       At firstcall, rdust and rice are not known; therefore
c         they need to be computed below.

c      Correction TN 17/04: rdust and rice must be updated at all steps, 
c      otherwise it is a possible source of bugs

c       1.1 Dust particles
c       ------------------
        IF (doubleq.AND.active) THEN
          DO l=1,nlayer
            DO ig=1, ngrid
              rdust(ig,l) = 
     &          CBRT(r3n_q*pq(ig,l,igcm_dust_mass)/
     &          max(pq(ig,l,igcm_dust_number),0.01))
              rdust(ig,l)=min(max(rdust(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
              rdust(ig,l) = 0.8E-6
              nueffdust(ig,l) = 0.3
            ENDDO
          ENDDO
        ENDIF
c       1.2 Water-ice particles
c       -----------------------
        IF (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
        ENDIF ! of if (water.AND.activice)
        
c        firstcall = .false.
c      ENDIF ! of if firstcall

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) = rdust(ig,l) * 
     &          (1.e0 + nueffdust(ig,l))**2.5
              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) = rdust(ig,l) * ref_r0
              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
c             About reffice, do not confuse the mass mean radius
c             (rayon moyen massique) and the number median radius
c             (or geometric mean radius, rayon moyen géométrique).
c             rice is a mass mean radius, whereas rdust
c             is a geometric mean radius:
c             number median rad = mass mean rad x exp(-1.5 sigma0^2)
c             (Montmessin et al. 2004 paragraph 30). Therefore:
              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