Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2592)
+++ /trunk/LMDZ.MARS/README	(revision 2593)
@@ -3539,2 +3539,5 @@
 == 06/12/2021 == CM
 MPCO2: remove old file Meteo_flux_Plane.dat, change file name Meteo_flux_Plane_v2.dat to Meteo_flux_Plane.dat
+
+== 10/12/2021 == MV
+hdo_surfex_mod.F, vdifc_mod.F: addings to account for the effect of kinetics in the fractionation by condensation of HDO at the surface
Index: /trunk/LMDZ.MARS/libf/phymars/hdo_surfex_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/hdo_surfex_mod.F	(revision 2592)
+++ /trunk/LMDZ.MARS/libf/phymars/hdo_surfex_mod.F	(revision 2593)
@@ -6,6 +6,6 @@
 
       subroutine hdo_surfex(ngrid,nlay,nq,ptimestep,
-     &                      zt,zq,pqsurf,
-     &                      old_h2o_vap,pdqsdif,dwatercap_dif,
+     &                      zt,pplay,zq,pqsurf,
+     &                      old_h2o_vap,qsat,pdqsdif,dwatercap_dif,
      &                      hdoflux)
 
@@ -15,4 +15,5 @@
       use surfdat_h, only: watercaptag
       use geometry_mod, only: longitude_deg,latitude_deg
+      use comcstfi_h, only: pi
 
       implicit none
@@ -23,5 +24,5 @@
 c------------------------------------------------------------------
       include "callkeys.h"
-
+      include "microphys.h"
 c------------------------------------------------------------------
 c     Arguments:
@@ -33,8 +34,10 @@
       REAL, INTENT(IN) :: ptimestep             ! pas de temps physique (s)
       REAL, INTENT(IN) :: zt(ngrid,nlay)       ! local value of temperature
+      REAL, INTENT(IN) :: pplay(ngrid,nlay)     ! pression au milieu des couches (Pa) 
       REAL, INTENT(IN) :: zq(ngrid,nlay,nq)    ! local value of tracers
       REAL, INTENT(IN) :: pqsurf(ngrid,nq)
       REAL, INTENT(IN) :: old_h2o_vap(ngrid)     ! traceur d'eau avant
                                            !traitement de l'eau (kg/kg)
+      REAL, INTENT(IN) :: qsat(ngrid) ! saturation mixing ratio
       REAL, INTENT(IN) :: dwatercap_dif(ngrid)  ! trend related to permanent ice
       REAL, INTENT(INOUT) :: pdqsdif(ngrid,nq)    ! tendance towards surface 
@@ -47,10 +50,13 @@
 c     Local variables:
 
-      REAL alpha_c(ngrid)  ! fractionation factor
+      REAL alpha(ngrid)    ! equilibrium fractionation factor
+      REAL alpha_c(ngrid)  ! real fractionation factor
       REAL extrasublim ! sublimation in excess of surface ice
       REAL tmpratio(ngrid)   ! D/H ratio in flux to surf
       REAL h2oflux(ngrid)       ! value of vapour flux of H2O
                                       ! same sign as pdqsdif
-
+      REAL*8 satu(ngrid)          ! Water vapor saturation ratio over ice
+      REAL zt1(ngrid),pplay1(ngrid)
+      REAL Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
       INTEGER ig,l
 
@@ -59,6 +65,9 @@
 c-----------------------------------------------------------------------
 c    Calculation of the fluxes for HDO
-        
-        alpha_c(1:ngrid)=0.
+        !! Calculation of the saturation ratio in the layer above the surface
+        satu(1:ngrid)=old_h2o_vap(1:ngrid) / qsat(1:ngrid)
+        !! Initialisation of the fractionation coefficient
+        alpha(1:ngrid)=1.        
+        alpha_c(1:ngrid)=1.
 
         DO ig=1,ngrid
@@ -105,6 +114,20 @@
 
                if (hdofrac) then !do we use fractionation?
-c               alpha_c(ig) = exp(16288./zt(ig,1)**2.-9.34e-2)
-                alpha_c = exp(13525./zt(ig,1)**2.-5.59e-2) !Lamb
+                !! Calculation of the H2O vapor diffusion coefficient
+                Dv = 1./3. * sqrt( 8*kbz*zt(ig,1)/(pi*mh2o/nav) )
+     &             * kbz * zt(ig,1) / 
+     &            ( pi * pplay(ig,1) * (molco2+molh2o)*(molco2+molh2o) 
+     &             * sqrt(1.+mh2o/mco2) )
+                !! Calculation of the HDO vapor diffusion coefficient
+                Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,1)/(pi*mhdo/nav) )
+     &             * kbz * zt(ig,1) / 
+     &            ( pi * pplay(ig,1) * (molco2+molhdo)*(molco2+molhdo) 
+     &             * sqrt(1.+mhdo/mco2) )
+                !! Calculation of the "equilibrium" fractionation coefficient
+c               alpha(ig) = exp(16288./zt(ig,1)**2.-9.34e-2)
+                alpha(ig) = exp(13525./zt(ig,1)**2.-5.59e-2) !Lamb
+                !! Calculation of the 'real' fractionnation coefficient (effect of kinetics, see Jouzel and Merlivat, 1984)
+                alpha_c(ig) = (alpha(ig)*satu(ig))/
+     &             ( (alpha(ig)*(Dv/Dv_hdo)*(satu(ig)-1.)) + 1.)
                else
                 alpha_c(ig) = 1.
Index: /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 2592)
+++ /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 2593)
@@ -1027,6 +1027,6 @@
 
             CALL hdo_surfex(ngrid,nlay,nq,ptimestep,
-     &                      zt,zq,pqsurf,
-     &                      old_h2o_vap,pdqsdif,dwatercap_dif,
+     &                      zt,pplay,zq,pqsurf,
+     &                      old_h2o_vap,qsat,pdqsdif,dwatercap_dif,
      &                      hdoflux)
             DO ig=1,ngrid