Index: trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 2406)
+++ trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 2407)
@@ -516,9 +516,4 @@
      &        ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rhocloud,
      &        zqi(1,1,iq),wq,beta)
-c             Surface Tendencies
-c             ~~~~~~~~~~
-              do ig=1,ngrid 
-                pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
-              end do
             else
 c             water ice sedimentation
@@ -527,53 +522,52 @@
      &        ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rho_q(iq),
      &        zqi(1,1,iq),wq,beta)
-c             Surface tendencies
-c             ~~~~~~~~~~
-              do ig=1,ngrid 
-                pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
+            endif ! of if (microphys)
+c           Surface tendencies
+c           ~~~~~~~~~~
+            do ig=1,ngrid 
+              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
+            end do
+c           Special case of isotopes
+c           ~~~~~~~~~~
+            !MVals: for now only water can have an isotope, a son ("fils"), and it has to be hdo
+            if (nqfils(iq).gt.0) then
+              if (iq.eq.igcm_h2o_ice) then
+               iq2=igcm_hdo_ice
+              else
+               call abort_physic("callsedim_mod","invalid isotope",1)
+              endif 
+              !MVals: input parameters in vlz_fi for hdo
+              do l=1,nlay
+               do ig=1,ngrid
+                if (zq0(ig,l,iq).gt.qperemin) then
+                 Ratio0(ig,l)=zq0(ig,l,iq2)/zq0(ig,l,iq)
+                else
+                 Ratio0(ig,l)=0.
+                endif
+                Ratio(ig,l)=Ratio0(ig,l)
+                masseq(ig,l)=max(masse(ig,l)*zq0(ig,l,iq),masseqmin)
+                w(ig,l)=wq(ig,l) !MVals: very important: hdo is transported by h2o (see vlsplt_p.F: correction bugg 15 mai 2015)
+               enddo !ig=1,ngrid
+              enddo !l=1,nlay
+              !MVals: no need to enter newsedim as the transporting mass w has been already calculated
+              call vlz_fi(ngrid,nlay,Ratio,2.,masseq,w,wq)
+              zqi(:,nlay,iq2)=zqi(:,nlay,iq)*Ratio0(:,nlay)
+              do l=1,nlay-1
+               do ig=1,ngrid
+                newmasse=max(masseq(ig,l)+w(ig,l+1)-w(ig,l),masseqmin)
+                Ratio(ig,l)=(Ratio0(ig,l)*masseq(ig,l)
+     &                       +wq(ig,l+1)-wq(ig,l))/newmasse                
+                zqi(ig,l,iq2)=zqi(ig,l,iq)*Ratio(ig,l)    
+               enddo
+              enddo !l=1,nlay-1
+              !MVals: hdo surface tendency
+              do ig=1,ngrid
+               if (w(ig,1).gt.masseqmin) then
+                 pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*(wq(ig,1)/w(ig,1))
+               else
+                 pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*Ratio0(ig,1)
+               endif
               end do
-c             Special case of isotopes
-c             ~~~~~~~~~~
-              !MVals: Loop over the sons ("fils")
-              if (nqfils(iq).gt.0) then
-                if (iq.eq.igcm_h2o_ice) then
-                 iq2=igcm_hdo_ice
-                else
-                 call abort_physic("callsedim_mod","invalid isotope",1)
-                endif 
-                !MVals: input paramters in vlz_fi for hdo
-                do l=1,nlay
-                 do ig=1,ngrid
-                  if (zq0(ig,l,iq).gt.qperemin) then
-                   Ratio0(ig,l)=zq0(ig,l,iq2)/zq0(ig,l,iq)
-                  else
-                   Ratio0(ig,l)=0.
-                  endif
-                  Ratio(ig,l)=Ratio0(ig,l)
-                  masseq(ig,l)=max(masse(ig,l)*zq0(ig,l,iq),masseqmin)
-                  w(ig,l)=wq(ig,l) !MVals: very important: hdo is transported by h2o (see vlsplt_p.F: correction bugg 15 mai 2015)
-                 enddo !ig=1,ngrid
-                enddo !l=1,nlay
-                !MVals: no need to enter newsedim as the transporting mass w has been already calculated
-                call vlz_fi(ngrid,nlay,Ratio,2.,masseq,w,wq)
-                zqi(:,nlay,iq2)=zqi(:,nlay,iq)*Ratio0(:,nlay)
-                do l=1,nlay-1
-                 do ig=1,ngrid
-                  newmasse=max(masseq(ig,l)+w(ig,l+1)-w(ig,l),masseqmin)
-                  Ratio(ig,l)=(Ratio0(ig,l)*masseq(ig,l)
-     &                         +wq(ig,l+1)-wq(ig,l))/newmasse                
-                  zqi(ig,l,iq2)=zqi(ig,l,iq)*Ratio(ig,l)    
-                 enddo
-                enddo !l=1,nlay-1
-                !MVals: hdo surface tendency
-                do ig=1,ngrid
-                 if (w(ig,1).gt.masseqmin) then
-                   pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*(wq(ig,1)/w(ig,1))
-                 else
-                   pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*Ratio0(ig,1)
-                 endif
-                end do
-              endif !(nqfils(iq).gt.0)
-            endif ! of if (microphys)
-
+            endif !(nqfils(iq).gt.0)
 c -----------------------------------------------------------------
 c         GENERAL CASE
Index: trunk/LMDZ.MARS/libf/phymars/growthrate.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/growthrate.F	(revision 2406)
+++ trunk/LMDZ.MARS/libf/phymars/growthrate.F	(revision 2407)
@@ -1,3 +1,3 @@
-      subroutine growthrate(temp,pmid,psat,rcrystal,res)
+      subroutine growthrate(temp,pmid,psat,rcrystal,res,Dv)
 
       use tracer_mod, only: rho_ice
@@ -34,5 +34,5 @@
 c     Output
       REAL res      ! growth resistance (res=Rk+Rd)
-
+      REAL Dv       ! water vapor diffusion coefficient
 
 c   local:
@@ -41,5 +41,5 @@
       REAL k,Lv                 
       REAL knudsen           ! Knudsen number (gas mean free path/particle radius)
-      REAL afactor,Dv,lambda       ! Intermediate computations for growth rate
+      REAL afactor,lambda       ! Intermediate computations for growth rate
       REAL Rk,Rd
       
Index: trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F	(revision 2406)
+++ trunk/LMDZ.MARS/libf/phymars/improvedclouds_mod.F	(revision 2407)
@@ -14,5 +14,7 @@
      &                      igcm_h2o_ice, igcm_dust_mass,
      &                      igcm_dust_number, igcm_ccn_mass,
-     &                      igcm_ccn_number
+     &                      igcm_ccn_number,
+     &                      igcm_hdo_vap,igcm_hdo_ice,
+     &                      qperemin
       use conc_mod, only: mmean
       use comcstfi_h, only: pi, cpp
@@ -82,4 +84,7 @@
       REAL cste
       REAL dMice           ! mass of condensed ice
+      REAL dMice_hdo       ! mass of condensed HDO ice
+      REAL alpha(ngrid,nlay) ! HDO equilibrium fractionation coefficient (Saturation=1)
+      REAL alpha_c(ngrid,nlay) ! HDO real fractionation coefficient
 !      REAL sumcheck
       REAL*8 ph2o          ! Water vapor partial pressure (Pa)
@@ -103,4 +108,5 @@
 
       REAL res      ! Resistance growth
+      REAL Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
       
 
@@ -364,5 +370,5 @@
 c       get resistance growth
         call growthrate(zt(ig,l),pplay(ig,l),
-     &             real(ph2o/satu),rice(ig,l),res)
+     &             real(ph2o/satu),rice(ig,l),res,Dv)
 
         res_out(ig,l) = res
@@ -391,5 +397,52 @@
        subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
           
-          
+c Special case of the isotope of water HDO    
+       if (hdo) then
+         !! condensation
+         if (dMice.gt.0.0) then
+         !! do we use fractionation?
+           if (hdofrac) then
+             !! Calculation of the HDO vapor coefficient
+             Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,l)/(pi*mhdo/nav) )
+     &          * kbz * zt(ig,l) / 
+     &          ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo) 
+     &          * sqrt(1.+mhdo/mco2) )
+             !! Calculation of the fractionnation coefficient at equilibrium
+             alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
+c             alpha = exp(13525./zt(ig,l)**2.-5.59d-2)  !Lamb
+             !! Calculation of the 'real' fractionnation coefficient
+             alpha_c(ig,l) = (alpha(ig,l)*satu)/
+     &          ( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.)
+c             alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics
+           else
+             alpha_c(ig,l) = 1.d0
+           endif
+           if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
+              dMice_hdo= 
+     &          dMice*alpha_c(ig,l)*
+     &     ( zq0(ig,l,igcm_hdo_vap)
+     &             /zq0(ig,l,igcm_h2o_vap) )
+           else
+             dMice_hdo=0.
+           endif
+         !! sublimation
+         else
+           if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
+             dMice_hdo= 
+     &            dMice*
+     &     ( zq0(ig,l,igcm_hdo_ice)
+     &             /zq0(ig,l,igcm_h2o_ice) )
+           else
+             dMice_hdo=0.
+           endif
+         endif !if (dMice.gt.0.0)
+
+       dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
+       dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
+
+       zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
+       zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo
+
+       endif ! if (hdo)        
      
 !=============================================================
@@ -405,4 +458,9 @@
      &                            + zq(ig,l,igcm_h2o_ice)
             zq(ig,l,igcm_h2o_ice) = 0.
+            if (hdo) then
+              zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) 
+     &                            + zq(ig,l,igcm_hdo_ice)
+              zq(ig,l,igcm_hdo_ice) = 0.
+            endif
 c           Dust particles
             zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
@@ -429,4 +487,12 @@
      &   (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
      &    zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
+        if (hdo) then
+          subpdqcloud(1:ngrid,1:nlay,igcm_hdo_vap) =
+     &     (zq(1:ngrid,1:nlay,igcm_hdo_vap) - 
+     &      zq0(1:ngrid,1:nlay,igcm_hdo_vap))/microtimestep
+          subpdqcloud(1:ngrid,1:nlay,igcm_hdo_ice) =
+     &     (zq(1:ngrid,1:nlay,igcm_hdo_ice) -
+     &      zq0(1:ngrid,1:nlay,igcm_hdo_ice))/microtimestep
+        endif
         subpdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
      &   (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
Index: trunk/LMDZ.MARS/libf/phymars/microphys.h
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/microphys.h	(revision 2406)
+++ trunk/LMDZ.MARS/libf/phymars/microphys.h	(revision 2407)
@@ -18,4 +18,6 @@
 !     Molecular weight of H2O (kg.mol-1)
       DOUBLE PRECISION, PARAMETER :: mh2o = 18.d-3
+!     Molecular weight of HDO (kg.mol-1)
+      DOUBLE PRECISION, PARAMETER :: mhdo = 19.d-3
 !     Molecular weight of CO2 (kg.mol-1)
       DOUBLE PRECISION, PARAMETER :: mco2 = 44.d-3
@@ -27,4 +29,6 @@
 !     Effective H2O gas molecular radius (m)
       DOUBLE PRECISION, PARAMETER :: molh2o = 1.2d-10
+!     Effective HDO gas molecular radius (m)
+      DOUBLE PRECISION, PARAMETER :: molhdo = 1.2d-10
 !     Surface tension of ice/vapor (N.m)
       DOUBLE PRECISION, PARAMETER :: sigh2o = 0.12
