subroutine growthrate(timestep,temp,pmid,ph2o,psat,
     &                      seq,rcrystal,growth)

      IMPLICIT NONE

c=======================================================================
c
c     Determination of the water ice crystal growth rate
c
c     Authors: F. Montmessin
c       Adapted for the LMD/GCM by J.-B. Madeleine (October 2011)
c     
c=======================================================================

c-----------------------------------------------------------------------
c   declarations:
c   -------------

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "tracer.h"
#include "microphys.h"

c
c   arguments:
c   ----------

      REAL timestep
      REAL temp     ! temperature in the middle of the layer (K)
      REAL pmid     ! pressure in the middle of the layer (K)
      REAL*8 ph2o   ! water vapor partial pressure (Pa)
      REAL*8 psat   ! water vapor saturation pressure (Pa) 
      REAL rcrystal ! crystal radius before condensation (m)
      REAL*8 seq    ! Equilibrium saturation ratio
      REAL*8 growth 

c   local:
c   ------

      REAL*8 k,Lv                 
      REAL*8 knudsen           ! Knudsen number (gas mean free path/particle radius)
      REAL*8 afactor,Dv,lambda       ! Intermediate computations for growth rate
      REAL*8 Rk,Rd

c-----------------------------------------------------------------------
c      Ice particle growth rate by diffusion/impegement of water molecules
c                r.dr/dt = (S-Seq) / (Seq*Rk+Rd)
c        with r the crystal radius, Rk and Rd the resistances due to 
c        latent heat release and to vapor diffusion respectively 
c----------------------------------------------------------------------- 

c     - Equilibrium saturation accounting for KeLvin Effect
c      seq=exp(2*sigh2o*mh2o/(rho_ice*rgp*t*r))
c      (already computed in improvedcloud.F)

c     - Thermal conductibility of CO2
      k  = (0.17913 * temp - 13.9789) * 4.184e-4
c     - Latent heat of h2o (J.kg-1)
      Lv = (2834.3 - 0.28 * (temp-To) - 0.004 * (temp-To)**2 ) * 1.e+3

c     - Constant to compute gas mean free path
c     l= (T/P)*a, with a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
      afactor = 0.707*rgp/(4 * pi* molco2**2  * nav)

c     - Compute Dv, water vapor diffusion coefficient
c       accounting for both kinetic and continuum regime of diffusion,
c       the nature of which depending on the Knudsen number.

      Dv = 1./3. * sqrt( 8*kbz*temp/(pi*mh2o/nav) )* kbz * temp / 
     &   ( pi * pmid * (molco2+molh2o)**2 * sqrt(1.+mh2o/mco2) )

      knudsen = temp / pmid * afactor / rcrystal
      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen)
      Dv      = Dv / (1. + lambda * knudsen)

c     - Compute Rk
      Rk = Lv**2 * rho_ice * mh2o / (k*rgp*temp**2.)
c     - Compute Rd
      Rd = rgp * temp *rho_ice / (Dv*psat*mh2o)

c     - Compute growth=rdr/dt, then r(t+1)= sqrt(r(t)**2.+2.*growth*dt)
      growth = 1. / (seq*Rk+Rd)
c       growth = (ph2o/psat-seq) / (seq*Rk+Rd)
c      rf   = sqrt( max( r**2.+2.*growth*timestep , 0. ) )
c      dr   = rf-r

      RETURN
      END