Index: trunk/LMDZ.MARS/libf/phymars/callkeys.h
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 2312)
@@ -16,6 +16,7 @@
      &   ,co2useh2o,meteo_flux,CLFvaryingCO2,spantCO2,CLFvarying        &
      &   ,satindexco2,rdstorm,slpwind,calllott_nonoro                   &
-     &   ,latentheat_surfwater,gwd_convective_source,startphy_file
-
+     &   ,latentheat_surfwater,gwd_convective_source,startphy_file      &
+     &   ,hdo,hdofrac
+     
       COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
      &   ,dustbin,nltemodel,nircorr,solvarmod,solvaryear,dustinjection
@@ -78,9 +79,10 @@
       logical CLFvarying
       logical water
+      logical hdo
+      logical hdofrac
       logical microphys
       logical photochem
       integer nltemodel
       integer nircorr
-
 
       character(len=100) dustiropacity
Index: trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 2312)
@@ -17,5 +17,6 @@
      &                      rho_dust, rho_q, radius, varian,
      &                      igcm_ccn_mass, igcm_ccn_number,
-     &                      igcm_h2o_ice, nuice_sed, nuice_ref,
+     &                      igcm_h2o_ice, igcm_hdo_ice,
+     &                      nuice_sed, nuice_ref,
      &                      igcm_ccnco2_mass,igcm_ccnco2_number,
      &                      igcm_co2_ice, igcm_stormdust_mass, 
@@ -499,5 +500,6 @@
 c -----------------------------------------------------------------
            else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
-     &       .or. (iq .eq. igcm_h2o_ice)) then
+     &       .or. (iq .eq. igcm_h2o_ice)
+     &       .or. (iq .eq. igcm_hdo_ice)) then
             if (microphys) then
               ! water ice sedimentation
Index: trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2312)
@@ -541,4 +541,15 @@
          call getin_p("water",water)
          write(*,*) " water = ",water
+! hdo
+         write(*,*) "Compute hdo cycle ?"
+         hdo=.false. ! default value
+         call getin_p("hdo",hdo)
+         write(*,*) " hdo = ",hdo
+
+         write(*,*) "Use fractionation for hdo?"
+         hdofrac=.true. ! default value
+         call getin_p("hdofrac",hdofrac)
+         write(*,*) " hdofrac = ",hdofrac
+
 
 ! sub-grid cloud fraction: fixed
@@ -618,4 +629,10 @@
      &          "water requires tracer",1)
          endif
+
+         if (hdo.and..not.water) then
+           print*,'if hdo is used, water should be used too'
+           stop
+         endif
+
          
 ! water ice clouds effective variance distribution for sedimentaion       
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 2312)
@@ -117,4 +117,6 @@
       integer :: nq=1 ! number of tracers
       real :: latitude(1), longitude(1), cell_area(1)
+      integer iqh2ovap
+      integer iqh2oice
 
       character*2 str2
@@ -294,4 +296,5 @@
           ! WATER VAPOUR
           if (txt.eq."h2o_vap") then
+             iqh2ovap=iq !remember index for water vap
             !look for a "profile_h2o_vap" input file
             open(91,file='profile_h2o_vap',status='old',
@@ -309,4 +312,5 @@
           ! WATER ICE
           if (txt.eq."h2o_ice") then
+                iqh2oice = iq !remember index for water ice
             !look for a "profile_h2o_vap" input file
             open(91,file='profile_h2o_ice',status='old',
@@ -322,4 +326,44 @@
             close(91)
           endif ! of if (txt.eq."h2o_ice")
+
+          ! HDO VAPOUR
+          if (txt.eq."hdo_vap") then
+            !look for a "profile_hdo_vap" input file
+            open(91,file='profile_hdo_vap',status='old',
+     &       form='formatted',iostat=ierr)
+            if (ierr.eq.0) then
+              read(91,*) qsurf(iq)
+              do ilayer=1,nlayer
+                read(91,*) q(ilayer,iq)
+              enddo
+            else
+              write(*,*) "No profile_hdo_vap file!"
+              do ilayer=1,nlayer
+                q(ilayer,iq) =  q(ilayer,iqh2ovap)*2*155.76e-6*5 
+              enddo
+            endif
+            close(91)
+          endif ! of if (txt.eq."hdo_ice")
+          ! HDO ICE
+          if (txt.eq."hdo_ice") then
+            !look for a "profile_hdo_vap" input file
+            open(91,file='profile_hdo_ice',status='old',
+     &       form='formatted',iostat=ierr)
+            if (ierr.eq.0) then
+              read(91,*) qsurf(iq)
+              do ilayer=1,nlayer
+                read(91,*) q(ilayer,iq)
+              enddo
+            else
+              write(*,*) "No profile_hdo_ice file!"
+                qsurf(iq) = qsurf(iqh2oice) * 2*155.76e-6*5
+               do ilayer=1,nlayer
+                q(ilayer,iq) = q(ilayer,iqh2oice) * 2*155.76e-6*5 
+              enddo
+            endif
+            close(91)
+          endif ! of if (txt.eq."hdo_ice")
+
+
           ! DUST
           !if (txt(1:4).eq."dust") then
Index: trunk/LMDZ.MARS/libf/phymars/hdo_surfex_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/hdo_surfex_mod.F	(revision 2312)
+++ trunk/LMDZ.MARS/libf/phymars/hdo_surfex_mod.F	(revision 2312)
@@ -0,0 +1,148 @@
+      MODULE hdo_surfex_mod
+
+      IMPLICIT NONE
+
+      CONTAINS
+
+      subroutine hdo_surfex(ngrid,nlay,nq,ptimestep,
+     &                      zt,zq,pqsurf,
+     &                      old_h2o_vap,pdqsdif,h2oflux,
+     &                      hdoflux)
+
+      use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice,
+     &                      igcm_hdo_vap, igcm_hdo_ice
+
+      use surfdat_h, only: watercaptag
+      use geometry_mod, only: longitude_deg,latitude_deg
+
+      implicit none
+c------------------------------------------------------------------
+c               Routine to compute the fluxes between air and surface
+c               for HDO, based of the fluxes for H2O
+c           L. Rossi.; M. Vals 2019
+c------------------------------------------------------------------
+      include "callkeys.h"
+
+c------------------------------------------------------------------
+c     Arguments:
+c     ---------
+c     Inputs:
+      INTEGER, INTENT(IN) :: ngrid,nlay
+      INTEGER, INTENT(IN) :: nq                 ! nombre de traceurs
+
+      REAL, INTENT(IN) :: ptimestep             ! pas de temps physique (s)
+      REAL, INTENT(IN) :: zt(ngrid,nlay)       ! local value of temperature
+      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(INOUT) :: pdqsdif(ngrid,nq)    ! tendance towards surface 
+                                 !   (kg/kg.s-1)
+      REAL, INTENT(IN) :: h2oflux(ngrid)       ! value of vapour flux of H2O
+                                      ! same sign as pdqsdif
+
+c     Output:
+      REAL, INTENT(OUT) :: hdoflux(ngrid)       ! value of vapour flux of HDO
+
+c------------------------------------------------------------------
+c     Local variables:
+
+      REAL alpha_c(ngrid)  ! fractionation factor
+      REAL extrasublim ! sublimation in excess of surface ice
+      REAL tmpratio(ngrid)   ! D/H ratio in flux to surf
+
+      INTEGER ig,l
+
+      REAL DoH_vap(ngrid)
+
+c-----------------------------------------------------------------------
+c    Calculation of the fluxes for HDO
+        
+        alpha_c(1:ngrid)=0.
+
+        DO ig=1,ngrid
+
+            !! IF Sublimation
+            if (h2oflux(ig).le.0.) then
+
+               if (pqsurf(ig,igcm_h2o_ice).gt.1e-16) then
+                pdqsdif(ig,igcm_hdo_ice) =
+     &            pdqsdif(ig,igcm_h2o_ice)*
+     &             (pqsurf(ig,igcm_hdo_ice)/
+     &             pqsurf(ig,igcm_h2o_ice) )
+               else
+                pdqsdif(ig,igcm_hdo_ice) = 0.
+               endif
+
+                pdqsdif(ig,igcm_hdo_ice)=
+     &             max(pdqsdif(ig,igcm_hdo_ice),
+     &            -pqsurf(ig,igcm_hdo_ice)/ptimestep)
+
+                hdoflux(ig) = pdqsdif(ig,igcm_hdo_ice)
+
+             if(watercaptag(ig)) then
+
+              !if we sublimate more than qsurf
+              if ((-h2oflux(ig)*ptimestep)
+     &           .gt.pqsurf(ig,igcm_h2o_ice)) then
+
+C               This is how much we sublimate in excess of
+C               pqsurf for H2O                        
+                  extrasublim=
+     &              ( pqsurf(ig,igcm_h2o_ice) /ptimestep 
+     &                   +h2oflux(ig))
+c               extrasublim has same sign conv. as pdqsdif
+ 
+C               hdoflux(ig) is the flux of HDO from atm. to surf.
+c               The D/H of the old ice is supposed to be 5 SMOW
+c               We need D/H of the flux to be 5, so we need
+c               extrasublim* 5 * 2 * 155.76e-6 (=1 SMOW)
+                    hdoflux(ig)= hdoflux(ig)
+     &                       +(extrasublim*(2.*155.76e-6)*5.)
+!    &                       +extrasublim
+!                 hdoflux(ig) = h2oflux(ig) !test
+              endif
+                endif ! watercap
+
+            else ! condensation
+
+               if (hdofrac) then !do we use fractionation?
+                alpha_c(ig) = exp(16288./zt(ig,1)**2.-9.34e-2)
+c               alpha_c = exp(13525./zt(ig,l)**2.-5.59e-2) !Lamb
+               else
+                alpha_c(ig) = 1.
+               endif
+
+               if (old_h2o_vap(ig).gt.1.e-16) then
+                         pdqsdif(ig,igcm_hdo_ice)=
+     &                      alpha_c(ig)*pdqsdif(ig,igcm_h2o_ice)*
+     &                      (zq(ig,1,igcm_hdo_vap)/
+     &                          old_h2o_vap(ig))
+                else
+                 pdqsdif(ig,igcm_hdo_ice)= 0.
+                endif
+
+               if (hdofrac) then !do we use fractionation?
+                pdqsdif(ig,igcm_hdo_ice)= 
+     &      min(  pdqsdif(ig,igcm_hdo_ice), 
+     &           (zq(ig,1,igcm_hdo_vap)/ptimestep) )
+               endif
+
+                hdoflux(ig)=pdqsdif(ig,igcm_hdo_ice)
+
+            endif !sublimation
+
+        ENDDO ! of DO ig=1,ngrid
+
+c           CALL WRITEDIAGFI(ngrid,'extrasublim',
+c    &                       'extrasublimation',
+c    &                       ' ',2,tmpratio)
+c           CALL WRITEDIAGFI(ngrid,'alpha_c_s',
+c    &                       'alpha_c_s',
+c    &                       ' ',2,alpha_c)
+
+       return
+      end subroutine hdo_surfex
+c------------------------------------------------------------------
+
+      end module hdo_surfex_mod
Index: trunk/LMDZ.MARS/libf/phymars/initracer.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 2312)
@@ -70,4 +70,6 @@
       igcm_h2o_vap=0
       igcm_h2o_ice=0
+      igcm_hdo_vap=0
+      igcm_hdo_ice=0
       igcm_stormdust_mass=0
       igcm_stormdust_number=0
@@ -345,4 +347,9 @@
           count=count+1
         endif
+        if (noms(iq).eq."hdo_vap") then
+          igcm_hdo_vap=iq
+          mmol(igcm_hdo_vap)=18. !19
+          count=count+1
+        endif
         if (noms(iq).eq."co2_ice") then
           igcm_co2_ice=iq
@@ -353,4 +360,9 @@
           igcm_h2o_ice=iq
           mmol(igcm_h2o_ice)=18.
+          count=count+1
+        endif
+        if (noms(iq).eq."hdo_ice") then
+          igcm_hdo_ice=iq
+          mmol(igcm_hdo_ice)=18. !19
           count=count+1
         endif
@@ -403,4 +415,13 @@
       endif
 
+      ! Additional test required for HDO
+      ! We need to compute some things for H2O before HDO
+      if (hdo) then
+        if (igcm_h2o_vap.gt.igcm_hdo_vap) then
+           write(*,*) "Tracer H2O must be initialized before HDO"
+           STOP
+        endif
+      endif
+
 c------------------------------------------------------------
 c     Initialize tracers .... (in tracer_mod)
@@ -429,5 +450,6 @@
 c       iq=1: Q mass mixing ratio, iq=2: N number mixing ratio
 
-        if( (nq.lt.2).or.(water.and.(nq.lt.4)) ) then
+        if( (nq.lt.2).or.(water.and.(nq.lt.4))
+     *       .or.(hdo.and.(nq.lt.6) )) then
           write(*,*)'initracer: nq is too low : nq=', nq
           write(*,*)'water= ',water,' doubleq= ',doubleq   
@@ -594,4 +616,23 @@
 
       end if  ! (water)
+
+c     Initialization for hdo vapor
+c     ------------------------------
+      if (hdo) then
+         radius(igcm_hdo_vap)=0.
+         alpha_lift(igcm_hdo_vap) =0.
+         alpha_devil(igcm_hdo_vap)=0.
+         if(water.and.(nq.ge.2)) then
+           radius(igcm_hdo_ice)=3.e-6
+           rho_q(igcm_hdo_ice)=rho_ice
+           alpha_lift(igcm_hdo_ice) =0.
+           alpha_devil(igcm_hdo_ice)=0.
+         elseif(hdo.and.(nq.lt.6)) then
+            write(*,*) 'nq is too low : nq=', nq
+            write(*,*) 'hdo= ',hdo
+         endif
+
+      end if  ! (hdo)
+
       
 ! Initialisation for CO2 clouds
@@ -683,4 +724,21 @@
          endif
        endif
+
+       if (hdo) then
+       ! verify that we indeed have hdo_vap and hdo_ice tracers
+         if (igcm_hdo_vap.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use hdo option without ",
+     &                "an hdo_vap tracer !"
+           stop
+         endif
+         if (igcm_hdo_ice.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use hdo option without ",
+     &                "an hdo_ice tracer !"
+           stop
+         endif
+       endif
+
 
        if (co2clouds) then
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2312)
@@ -293,5 +293,4 @@
             minval(base), maxval(base)
 
- 
 ! Time axis
 ! obtain timestart from run.def
@@ -567,4 +566,13 @@
       write(*,*) 'iq = ', iq
     endif
+
+    if (hdo) then
+    if (txt.eq."hdo_vap") then
+      txt="hdo_ice"
+      write(*,*) 'phyetat0: loading surface tracer', &
+                           ' hdo_ice instead of hdo_vap'
+    endif
+    endif !hdo
+
     if (startphy_file) then
        call get_field(txt,qsurf(:,iq),found,indextime)
Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2312)
@@ -189,4 +189,6 @@
   integer :: i_h2o_vap=0
   integer :: i_h2o_ice=0
+  integer :: i_hdo_vap=0
+  integer :: i_hdo_ice=0
 
   
@@ -245,5 +247,14 @@
       i_h2o_ice=iq
     endif
+
+  ! look for HDO vapour & HDO ice tracers (if any)
+    if (noms(iq).eq."hdo_vap") then
+      i_hdo_vap=iq
+    endif
+    if (noms(iq).eq."hdo_ice") then
+      i_hdo_ice=iq
+    endif
   enddo
+
   
   if (nq.gt.0) then
@@ -264,4 +275,19 @@
         endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap)
       endif ! of if (txt.eq."h2o_vap")
+
+      if (txt.eq."hdo_vap") then
+        write(*,*)"physdem1: skipping HDO vapour tracer"
+        if (i_hdo_ice.eq.i_hdo_vap) then
+          ! then there is no "water ice" tracer; but still
+          ! there is some water ice on the surface
+          write(*,*)"          writing HDO ice instead"
+          txt="hdo_ice"
+        else
+          ! there is a "water ice" tracer which has been / will be
+          ! delt with in due time
+          cycle
+        endif ! of if (igcm_hdo_ice.eq.igcm_hdo_vap)
+      endif ! of if (txt.eq."hdo_vap")
+
       call put_field(trim(txt),"tracer on surface",qsurf(:,iq),time)
     enddo
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2312)
@@ -27,4 +27,5 @@
       use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice,
      &                      igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
+     &                      igcm_hdo_vap, igcm_hdo_ice,
      &                      igcm_ccn_mass, igcm_ccn_number,
      &                      igcm_ccnco2_mass, igcm_ccnco2_number,
@@ -87,4 +88,5 @@
       use wxios, only: wxios_context_init, xios_context_finalize
 #endif
+
 
       IMPLICIT NONE
@@ -395,4 +397,9 @@
                                !   Qabs instead of Qext
                                !   (direct comparison with TES)
+      REAL mtotD(ngrid)          ! Total mass of HDO vapor (kg/m2)
+      REAL icetotD(ngrid)        ! Total mass of HDO ice (kg/m2)
+      REAL DoH_vap(ngrid,nlayer) !D/H ratio
+      REAL DoH_ice(ngrid,nlayer) !D/H ratio
+      REAL DoH_surf(ngrid) !D/H ratio surface
                                
       REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s)
@@ -418,4 +425,5 @@
 c Test 1d/3d scavenging
       real h2otot(ngrid)
+      real hdotot(ngrid)
       REAL satu(ngrid,nlayer)  ! satu ratio for output
       REAL zqsat(ngrid,nlayer) ! saturation
@@ -701,4 +709,6 @@
       call update_xios_timestep
 #endif      
+
+c     Initialize various variables
 c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       pdv(:,:)=0
@@ -714,5 +724,5 @@
       dsotop(:,:)=0.
       dwatercap(:)=0
-
+      
 #ifdef DUSTSTORM
       pq_tmp(:,:,:)=0
@@ -773,4 +783,5 @@
          ENDDO
       ENDDO
+
 #ifndef MESOSCALE
 c-----------------------------------------------------------------------
@@ -886,4 +897,5 @@
            !! callradite for background dust in the case of slope wind entrainment
            nohmons=.true.
+
 c          Radiative transfer
 c          ------------------
@@ -1200,4 +1212,5 @@
            endif ! end of if (dustinjection.gt.0)
 
+
 c-----------------------------------------------------------------------
 c    4. Gravity wave and subgrid scale topography drag :
@@ -1263,4 +1276,5 @@
              pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1)
          ENDIF
+
 c        Calling vdif (Martian version WITH CO2 condensation)
          CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
@@ -1289,4 +1303,5 @@
             ENDDO
           ENDDO
+
           if (tracer) then
            DO iq=1, nq
@@ -1369,5 +1384,5 @@
           q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
         ENDDO      
-    
+
         if (tracer) then
         DO iq=1,nq
@@ -1388,4 +1403,5 @@
         lmax_th_out(:)=0.
       end if
+
 c-----------------------------------------------------------------------
 c   7. Dry convective adjustment:
@@ -1483,5 +1499,5 @@
      &          zdtcloud(1:ngrid,1:nlayer)
            endif
-           
+
 ! increment water vapour and ice atmospheric tracers tendencies
            pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = 
@@ -1491,4 +1507,14 @@
      &       pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + 
      &       zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice)
+
+           if (hdo) then
+! increment HDO vapour and ice atmospheric tracers tendencies
+           pdq(1:ngrid,1:nlayer,igcm_hdo_vap) = 
+     &       pdq(1:ngrid,1:nlayer,igcm_hdo_vap) + 
+     &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_vap)
+           pdq(1:ngrid,1:nlayer,igcm_hdo_ice) = 
+     &       pdq(1:ngrid,1:nlayer,igcm_hdo_ice) + 
+     &       zdqcloud(1:ngrid,1:nlayer,igcm_hdo_ice)
+           endif !hdo
 
 ! increment dust and ccn masses and numbers
@@ -1896,11 +1922,11 @@
             ENDDO
            ENDDO
-	   	   
-	   DO iq=1, nq
+
+          DO iq=1, nq
             DO ig=1,ngrid
               dqsurf(ig,iq)=dqsurf(ig,iq)+zdqsc(ig,iq)
             ENDDO  ! (ig)
            ENDDO    ! (iq)
-	   
+
          ENDIF ! of IF (tracer)
 
@@ -2248,4 +2274,5 @@
      &                dsotop(ig,:)*tauscaling(ig)
               enddo
+         
            if(doubleq) then
               do ig=1,ngrid 
@@ -2325,4 +2352,10 @@
              rave(:)=0
              tauTES(:)=0
+
+             IF (hdo) then
+                 mtotD(:)=0
+                 icetotD(:)=0
+             ENDIF !hdo
+
              do ig=1,ngrid 
                do l=1,nlayer
@@ -2333,4 +2366,13 @@
      &                        zq(ig,l,igcm_h2o_ice) * 
      &                        (zplev(ig,l) - zplev(ig,l+1)) / g
+                 IF (hdo) then
+                 mtotD(ig) = mtotD(ig) + 
+     &                      zq(ig,l,igcm_hdo_vap) * 
+     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
+                 icetotD(ig) = icetotD(ig) + 
+     &                        zq(ig,l,igcm_hdo_ice) * 
+     &                        (zplev(ig,l) - zplev(ig,l+1)) / g
+                 ENDIF !hdo
+
 c                Computing abs optical depth at 825 cm-1 in each
 c                  layer to simulate NEW TES retrieval
@@ -2887,4 +2929,56 @@
      &             3,zq(:,:,igcm_h2o_vap))
 
+            if (hdo) then
+            vmr=zq(1:ngrid,1:nlayer,igcm_hdo_ice)
+     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_ice)
+            CALL WRITEDIAGFI(ngrid,'vmr_hdoice','hdo ice vmr',
+     &                       'mol/mol',3,vmr)
+            vmr=zq(1:ngrid,1:nlayer,igcm_hdo_vap)
+     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_vap)
+            CALL WRITEDIAGFI(ngrid,'vmr_hdovap','hdo vap vmr',
+     &                       'mol/mol',3,vmr)
+            call WRITEDIAGFI(ngrid,'hdo_ice','hdo_ice','kg/kg',
+     &             3,zq(:,:,igcm_hdo_ice))
+            call WRITEDIAGFI(ngrid,'hdo_vap','hdo_vap','kg/kg',
+     &             3,zq(:,:,igcm_hdo_vap))
+
+            CALL WRITEDIAGFI(ngrid,'mtotD',
+     &                       'total mass of HDO vapor',
+     &                       'kg/m2',2,mtotD)
+            CALL WRITEDIAGFI(ngrid,'icetotD',
+     &                       'total mass of HDO ice',
+     &                       'kg/m2',2,icetotD)
+
+C           Calculation of the D/H ratio
+            do l=1,nlayer
+                do ig=1,ngrid
+                if (zq(ig,l,igcm_h2o_vap).gt.1e-16) then 
+                    DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/
+     &              zq(ig,l,igcm_h2o_vap) )*1./(2.*155.76e-6)
+                else
+                    DoH_vap(ig,l) = 0.
+                endif
+                enddo
+            enddo
+
+            do l=1,nlayer
+                do ig=1,ngrid
+                if (zq(ig,l,igcm_h2o_ice).gt.1e-16) then
+                    DoH_ice(ig,l) = ( zq(ig,l,igcm_hdo_ice)/
+     &                  zq(ig,l,igcm_h2o_ice) )/(2.*155.76e-6)
+                else
+                    DoH_ice(ig,l) = 0.
+                endif
+                enddo
+            enddo
+
+            CALL WRITEDIAGFI(ngrid,'DoH_vap',
+     &                       'D/H ratio in vapor',
+     &                       ' ',3,DoH_vap)  
+            CALL WRITEDIAGFI(ngrid,'DoH_ice',
+     &                       'D/H ratio in ice',
+     &                       '',3,DoH_ice)
+
+            endif !hdo
 
 !A. Pottier
@@ -2907,4 +3001,23 @@
      &                       'surface h2o_ice',
      &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
+            if (hdo) then
+            call WRITEDIAGFI(ngrid,'hdo_ice_s',
+     &                       'surface hdo_ice',
+     &                       'kg.m-2',2,qsurf(1,igcm_hdo_ice))
+
+                do ig=1,ngrid
+                if (qsurf(ig,igcm_h2o_ice).gt.1e-30) then
+                    DoH_surf(ig) = 0.5*( qsurf(ig,igcm_hdo_ice)/
+     &                  qsurf(ig,igcm_h2o_ice) )/155.76e-6
+                else
+                    DoH_surf(ig) = 0.
+                endif
+                enddo
+
+            call WRITEDIAGFI(ngrid,'DoH_surf',
+     &                       'surface D/H',
+     &                       '',2,DoH_surf)
+            endif ! hdo
+
             CALL WRITEDIAGFI(ngrid,'albedo',
      &                         'albedo',
@@ -3430,4 +3543,9 @@
            icetot = 0
            h2otot = qsurf(1,igcm_h2o_ice)
+           if (hdo) THEN
+           mtotD = 0
+           icetotD = 0
+           hdotot = qsurf(1,igcm_hdo_ice)
+           ENDIF !hdo
 
            do l=1,nlayer
@@ -3436,6 +3554,16 @@
              icetot = icetot +  zq(1,l,igcm_h2o_ice) 
      &                 * (zplev(1,l) - zplev(1,l+1)) / g
+               if (hdo) THEN
+                 mtotD = mtotD +  zq(1,l,igcm_hdo_vap) 
+     &                 * (zplev(1,l) - zplev(1,l+1)) / g
+                 icetotD = icetotD +  zq(1,l,igcm_hdo_ice) 
+     &                 * (zplev(1,l) - zplev(1,l+1)) / g
+               ENDIF !hdo
            end do
            h2otot = h2otot+mtot+icetot
+               IF (hdo) then
+                   hdotot = hdotot+mtotD+icetotD
+               ENDIF ! hdo
+                
 
              CALL WRITEDIAGFI(ngrid,'h2otot',
@@ -3448,4 +3576,44 @@
      &                         'icetot',
      &                         'kg/m2',0,icetot)
+
+             IF (hdo) THEN
+             CALL WRITEDIAGFI(ngrid,'mtotD',
+     &                         'mtotD',
+     &                         'kg/m2',0,mtotD)
+             CALL WRITEDIAGFI(ngrid,'icetotD',
+     &                         'icetotD',
+     &                         'kg/m2',0,icetotD)
+             CALL WRITEDIAGFI(ngrid,'hdotot',
+     &                         'hdotot',
+     &                         'kg/m2',0,hdotot)
+
+C           Calculation of the D/H ratio
+            do l=1,nlayer
+                if (zq(1,l,igcm_h2o_vap).gt.1e-16) then
+                    DoH_vap(1,l) = 0.5*( zq(1,l,igcm_hdo_vap)/
+     &                zq(1,l,igcm_h2o_vap) )/155.76e-6
+                else
+                    DoH_vap(1,l) = 0.
+                endif
+            enddo
+
+            do l=1,nlayer
+                if (zq(1,l,igcm_h2o_ice).gt.1e-16) then
+                    DoH_ice(1,l) = 0.5*( zq(1,l,igcm_hdo_ice)/
+     &                  zq(1,l,igcm_h2o_ice) )/155.76e-6
+                else
+                    DoH_ice(1,l) = 0.
+                endif
+            enddo
+
+            CALL WRITEDIAGFI(ngrid,'DoH_vap',
+     &                       'D/H ratio in vapor',
+     &                       ' ',1,DoH_vap)  
+            CALL WRITEDIAGFI(ngrid,'DoH_ice',
+     &                       'D/H ratio in ice',
+     &                       '',1,DoH_ice)
+
+             ENDIF !Hdo
+
 
            if (scavenging) then
@@ -3527,4 +3695,11 @@
      &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
      &                          +zdqcloud(1,:,igcm_h2o_vap))
+        IF (hdo) THEN
+        call WRITEDIAGFI(ngrid,'zdqcloud_iceD','cloud ice hdo',
+     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_hdo_ice))
+        call WRITEDIAGFI(ngrid,'zdqcloud_vapD','cloud vap hdo',
+     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_hdo_vap))
+
+        ENDIF ! hdo
      
         call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
Index: trunk/LMDZ.MARS/libf/phymars/simpleclouds.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/simpleclouds.F	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/simpleclouds.F	(revision 2312)
@@ -5,7 +5,9 @@
       USE updaterad
       USE watersat_mod, ONLY: watersat
-      use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice
+      use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice,
+     &                      igcm_hdo_vap, igcm_hdo_ice
       USE comcstfi_h
       use dimradmars_mod, only: naerkind
+
       implicit none
 c------------------------------------------------------------------
@@ -77,8 +79,10 @@
       REAL ccntyp(ngrid,nlay)
                                       ! Typical dust number density (#/kg)
+      REAL alpha_c(ngrid,nlay) !!MARGAUX: alpha_c as a spatial variable
+
 c     CCN reduction factor
 c      REAL, PARAMETER :: ccn_factor = 4.5  !! comme TESTS_JB // 1. avant
-      
-
+      REAL DoH_vap(ngrid,nlay)
+      REAL DoH_ice(ngrid,nlay)
 c-----------------------------------------------------------------------
 c    1. initialisation
@@ -100,13 +104,26 @@
           zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004 
           zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
+
+          if (hdo) then
+          zq(ig,l,igcm_hdo_vap)=
+     &      pq(ig,l,igcm_hdo_vap)+pdq(ig,l,igcm_hdo_vap)*ptimestep
+          zq(ig,l,igcm_hdo_vap)=max(zq(ig,l,igcm_hdo_vap),1e-30) ! FF 12/2004 
+          zq0(ig,l,igcm_hdo_vap)=zq(ig,l,igcm_hdo_vap)
+
+          zq(ig,l,igcm_hdo_ice)=
+     &      pq(ig,l,igcm_hdo_ice)+pdq(ig,l,igcm_hdo_ice)*ptimestep
+          zq(ig,l,igcm_hdo_ice)=max(zq(ig,l,igcm_hdo_ice),1e-30) ! FF 12/2004 
+          zq0(ig,l,igcm_hdo_ice)=zq(ig,l,igcm_hdo_ice)
+
+          endif !hdo
         enddo
       enddo
 
-
       pdqcloud(1:ngrid,1:nlay,1:nq)=0
       pdtcloud(1:ngrid,1:nlay)=0
 
+      alpha_c(1:ngrid,1:nlay)=0.
+
 c     ----------------------------------------------
-c
 c
 c     Rapport de melange a saturation dans la couche l : -------
@@ -122,18 +139,19 @@
 
           if (zq(ig,l,igcm_h2o_vap).ge.zqsat(ig,l))then  !  Condensation
-            dzq=zq(ig,l,igcm_h2o_vap)-zqsat(ig,l)                
+            dzq=zq(ig,l,igcm_h2o_vap)-zqsat(ig,l)
+
           elseif(zq(ig,l,igcm_h2o_vap).lt.zqsat(ig,l))then  ! Sublimation
             dzq=-min(zqsat(ig,l)-zq(ig,l,igcm_h2o_vap),
      &               zq(ig,l,igcm_h2o_ice))
           endif
-
+            
 c         Water Mass change
 c         ~~~~~~~~~~~~~~~~~
           zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)+dzq
           zq(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)-dzq
-          
-
+        
         enddo ! of do ig=1,ngrid
       enddo ! of do l=1,nlay
+
 
 c     Tendance finale
@@ -145,8 +163,58 @@
           pdqcloud(ig,l,igcm_h2o_ice) =
      &      (zq(ig,l,igcm_h2o_ice) - zq0(ig,l,igcm_h2o_ice))/ptimestep
+
           lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
           pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
-        end do
-      end do
+
+          if (hdo) then
+            if (pdqcloud(ig,l,igcm_h2o_ice).gt.0.0) then !condens
+
+                if (hdofrac) then ! do we use fractionation?
+                alpha_c(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
+c               alpha_c(ig,l) = exp(13525./zt(ig,l)**2.-5.59d-2)  !Lamb
+                else
+                alpha_c(ig,l) = 1.d0
+                endif
+               
+                if (zq0(ig,l,igcm_h2o_vap).gt.1.e-16) then
+                  pdqcloud(ig,l,igcm_hdo_ice)= 
+     &              pdqcloud(ig,l,igcm_h2o_ice)*alpha_c(ig,l)*
+     &         ( zq0(ig,l,igcm_hdo_vap)
+     &                 /zq0(ig,l,igcm_h2o_vap) )
+                else
+                   pdqcloud(ig,l,igcm_hdo_ice)= 0.0
+                endif
+
+                pdqcloud(ig,l,igcm_hdo_ice) =
+     &            min(pdqcloud(ig,l,igcm_hdo_ice),
+     &              zq0(ig,l,igcm_hdo_vap)/ptimestep)
+
+                pdqcloud(ig,l,igcm_hdo_vap)= 
+     &               -pdqcloud(ig,l,igcm_hdo_ice)       
+
+          else  ! sublimation
+
+             if (zq0(ig,l,igcm_h2o_ice).gt.1.e-16) then 
+                pdqcloud(ig,l,igcm_hdo_ice)= 
+     &               pdqcloud(ig,l,igcm_h2o_ice)*
+     &      ( zq0(ig,l,igcm_hdo_ice)
+     &              /zq0(ig,l,igcm_h2o_ice) )
+             else
+                pdqcloud(ig,l,igcm_hdo_ice)= 0.
+             endif
+
+              pdqcloud(ig,l,igcm_hdo_ice) =
+     &          max(pdqcloud(ig,l,igcm_hdo_ice),
+     &            -zq0(ig,l,igcm_hdo_ice)/ptimestep)
+
+              pdqcloud(ig,l,igcm_hdo_vap)= 
+     &             -pdqcloud(ig,l,igcm_hdo_ice)        
+
+            endif ! condensation/sublimation
+
+          endif ! hdo
+
+        enddo ! of do ig=1,ngrid
+      enddo ! of do l=1,nlay
 
 c     ice crystal radius
@@ -158,4 +226,9 @@
       end do
 
+c     if (hdo) then
+c           CALL WRITEDIAGFI(ngrid,'alpha_c',
+c    &                       'alpha_c',
+c    &                       ' ',3,alpha_c) 
+c     endif !hdo
 c------------------------------------------------------------------
       return
Index: trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90	(revision 2312)
@@ -53,4 +53,6 @@
       integer,save :: igcm_h2o_vap ! water vapour
       integer,save :: igcm_h2o_ice ! water ice
+      integer,save :: igcm_hdo_vap ! hdo vapour
+      integer,save :: igcm_hdo_ice ! hdo ice
       integer,save :: igcm_co2_ice ! co2 ice
 
Index: trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 2312)
@@ -18,4 +18,5 @@
      &                      igcm_dust_submicron, igcm_h2o_vap,
      &                      igcm_h2o_ice, alpha_lift,
+     &                      igcm_hdo_vap, igcm_hdo_ice,
      &                      igcm_stormdust_mass, igcm_stormdust_number
       use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness
@@ -24,4 +25,6 @@
       use turb_mod, only: turb_resolved, ustar, tstar
       use compute_dtau_mod, only: ti_injection_sol,tf_injection_sol
+      use hdo_surfex_mod, only: hdo_surfex
+
       IMPLICIT NONE
 
@@ -140,4 +143,8 @@
       REAL qsat(ngrid)
 
+      REAL hdoflux(ngrid)       ! value of vapour flux of HDO
+      REAL h2oflux(ngrid)       ! value of vapour flux of H2O
+      REAL old_h2o_vap(ngrid)   ! traceur d'eau avant traitement
+
       REAL kmixmin
 
@@ -167,4 +174,7 @@
 
       REAL,INTENT(OUT) :: sensibFlux(ngrid)
+
+!!MARGAUX
+      REAL DoH_vap(ngrid,nlay)
 
 c    ** un petit test de coherence
@@ -443,5 +453,5 @@
 c       donc les entrees sont /zcdv/ pour la condition a la limite sol
 c       et /zkv/ = Ku
- 
+
       zb(1:ngrid,2:nlay)=zkv(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
       zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
@@ -786,5 +796,5 @@
           zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
 
-          if ((water).and.(iq.eq.igcm_h2o_vap)) then 
+          if ((water).and.(iq.eq.igcm_h2o_vap)) then
 c            This line is required to account for turbulent transport 
 c            from surface (e.g. ice) to mid-layer of atmosphere:
@@ -811,5 +821,6 @@
           ENDDO
 
-          if (water.and.(iq.eq.igcm_h2o_ice)) then
+          if ( (water.and.(iq.eq.igcm_h2o_ice)) 
+     $      .or. (hdo.and.(iq.eq.igcm_hdo_ice)) ) then
             ! special case for water ice tracer: do not include
             ! h2o ice tracer from surface (which is set when handling
@@ -821,4 +832,17 @@
      $         zb(ig,2)*zc(ig,2)) *z1(ig)
             ENDDO
+         else if (hdo.and.(iq.eq.igcm_hdo_vap)) then
+            CALL hdo_surfex(ngrid,nlay,nq,ptimestep,
+     &                      zt,zq,pqsurf,
+     &                      old_h2o_vap,pdqsdif,h2oflux,
+     &                      hdoflux)
+            DO ig=1,ngrid
+                z1(ig)=1./(za(ig,1)+zb(ig,1)+
+     $           zb(ig,2)*(1.-zd(ig,2)))
+                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
+     $         zb(ig,2)*zc(ig,2) +
+     $        (-hdoflux(ig)) *ptimestep) *z1(ig)  !tracer flux from surface
+            ENDDO
+
           else ! general case
             DO ig=1,ngrid
@@ -829,4 +853,5 @@
      $        (-pdqsdif(ig,iq)) *ptimestep) *z1(ig)  !tracer flux from surface
             ENDDO
+
           endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
   
@@ -834,4 +859,5 @@
 c           Calculation for turbulent exchange with the surface (for ice)
             DO ig=1,ngrid
+              old_h2o_vap(ig)=zq(ig,1,igcm_h2o_vap)
               zd(ig,1)=zb(ig,1)*z1(ig)
               zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
@@ -843,4 +869,9 @@
 
             DO ig=1,ngrid
+            IF (hdo) then !if hdo, need to treat polar caps differently
+              h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice) ! h2oflux is
+                                                     ! uncorrected waterflux
+            ENDIF !hdo
+
               if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
      &             .gt.(pqsurf(ig,igcm_h2o_ice))) then
@@ -857,4 +888,7 @@
      &                         -pqsurf(ig,igcm_h2o_ice)/ptimestep
                    dwatercap_dif(ig) = 0.0
+                   if (hdo) then
+                     h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice)
+                   endif
 
                 endif   
@@ -910,5 +944,14 @@
         enddo ! of do iq=1,nq
       end if ! of if(tracer)
-      
+
+C       Diagnostic output for HDO
+        if (hdo) then
+            CALL WRITEDIAGFI(ngrid,'hdoflux',
+     &                       'hdoflux',
+     &                       ' ',2,hdoflux)  
+            CALL WRITEDIAGFI(ngrid,'h2oflux',
+     &                       'h2oflux',
+     &                       ' ',2,h2oflux)
+        endif
 
 c-----------------------------------------------------------------------
@@ -958,3 +1001,6 @@
       END SUBROUTINE vdifc
 
+c====================================
+
+
       END MODULE vdifc_mod
Index: trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F	(revision 2311)
+++ trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F	(revision 2312)
@@ -19,4 +19,5 @@
       USE watersat_mod, ONLY: watersat
       use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
+     &                      igcm_hdo_vap, igcm_hdo_ice,
      &                      igcm_dust_mass, igcm_dust_number,
      &                      igcm_ccn_mass, igcm_ccn_number,
@@ -124,4 +125,5 @@
       REAL :: mincloud ! min cloud frac
       LOGICAL, save :: flagcloud=.true. 
+
 c    ** un petit test de coherence
 c       --------------------------
@@ -312,4 +314,12 @@
      &          sum_subpdq(ig,l,igcm_h2o_vap)
      &        + pdq(ig,l,igcm_h2o_vap)
+            IF (hdo) THEN
+            sum_subpdq(ig,l,igcm_hdo_ice) = 
+     &          sum_subpdq(ig,l,igcm_hdo_ice)
+     &        + pdq(ig,l,igcm_hdo_ice)
+            sum_subpdq(ig,l,igcm_hdo_vap) = 
+     &          sum_subpdq(ig,l,igcm_hdo_vap)
+     &        + pdq(ig,l,igcm_hdo_vap)
+            ENDIF !hdo
           ENDDO
         ENDDO      
@@ -361,4 +371,14 @@
      &          sum_subpdq(ig,l,igcm_h2o_vap)
      &        + subpdqcloud(ig,l,igcm_h2o_vap)
+
+            IF (hdo) THEN
+            sum_subpdq(ig,l,igcm_hdo_ice) =
+     &          sum_subpdq(ig,l,igcm_hdo_ice)
+     &        + subpdqcloud(ig,l,igcm_hdo_ice)
+            sum_subpdq(ig,l,igcm_hdo_vap) =
+     &          sum_subpdq(ig,l,igcm_hdo_vap)
+     &        + subpdqcloud(ig,l,igcm_hdo_vap)
+            ENDIF ! hdo
+
           ENDDO
         ENDDO
@@ -397,4 +417,12 @@
      &        sum_subpdq(ig,l,igcm_h2o_vap)/real(imicro) 
      &       - pdq(ig,l,igcm_h2o_vap)
+            IF (hdo) THEN
+            pdqcloud(ig,l,igcm_hdo_ice) = 
+     &        sum_subpdq(ig,l,igcm_hdo_ice)/real(imicro) 
+     &       - pdq(ig,l,igcm_hdo_ice)
+            pdqcloud(ig,l,igcm_hdo_vap) = 
+     &        sum_subpdq(ig,l,igcm_hdo_vap)/real(imicro) 
+     &       - pdq(ig,l,igcm_hdo_vap)
+            ENDIF !hdo
          ENDDO
        ENDDO
@@ -478,4 +506,5 @@
         DO l=1,nlay
          DO ig=1,ngrid
+
           IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
      &       (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice)) 
@@ -484,4 +513,9 @@
      &     - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
            pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
+             if (hdo) then
+           pdqcloud(ig,l,igcm_hdo_ice) = 
+     &     - pq(ig,l,igcm_hdo_ice)/ptimestep - pdq(ig,l,igcm_hdo_ice)
+           pdqcloud(ig,l,igcm_hdo_vap) = -pdqcloud(ig,l,igcm_hdo_ice)
+             endif
           ENDIF 
           IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
@@ -491,8 +525,13 @@
      &     - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
            pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
+             if (hdo) then
+           pdqcloud(ig,l,igcm_hdo_vap) = 
+     &     - pq(ig,l,igcm_hdo_vap)/ptimestep - pdq(ig,l,igcm_hdo_vap)
+           pdqcloud(ig,l,igcm_hdo_ice) = -pdqcloud(ig,l,igcm_hdo_vap)
+             endif
           ENDIF
-         ENDDO
-        ENDDO
-
+
+         ENDDO
+        ENDDO
 
 c------Update the ice and dust particle size "rice" for output or photochemistry
@@ -627,4 +666,11 @@
      &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
      &      igcm_h2o_vap))
+      if (hdo) then
+      call WRITEDIAGFI(ngrid,"pdqiceD","pdqiceD apres microphysique"
+     &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,igcm_hdo_ice))
+      call WRITEDIAGFI(ngrid,"pdqvapD","pdqvapD apres microphysique"
+     &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
+     &      igcm_hdo_vap))
+      endif
       call WRITEDIAGFI(ngrid,"pdqccn2","pdqcloudccn apres microphysique"
      &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
