Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 2377)
+++ trunk/LMDZ.MARS/README	(revision 2378)
@@ -3104,2 +3104,8 @@
 Minor fix to r2362 (addition of CO2 microphysics): make CO2 conservation
 tests only if CO2 microphysics is on.
+
+== 18/06/2020 == LR
+HDO
+Correction of an error in newstart for inihdo.
+Other minor corrections for HDO cycle.
+Transition from fractionation coefficients from Merlivat et al. 1967 to Lamb et al. 2017
Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 2377)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 2378)
@@ -1010,20 +1010,16 @@
         !D/H (SMOW) = 155.76e-6 so HDO/H2O is twice that
 
-          DO l=1,llm
-            DO j=1,jjp1
-              DO i=1,iip1-1
-              q(i,j,l,igcm_hdo_vap)= q(i,j,l,igcm_h2o_vap) * DoverH
-              q(i,j,l,igcm_hdo_ice)= q(i,j,l,igcm_h2o_ice) * DoverH
-              ENDDO
-              ! We want to have the very same value at lon -180 and lon 180
-              q(iip1,j,l,igcm_hdo_vap) = q(1,j,l,igcm_hdo_vap)
-              q(iip1,j,l,igcm_hdo_ice) = q(1,j,l,igcm_hdo_ice)
-            ENDDO
-          ENDDO
-
            do ig=1,ngridmx
            qsurf(ig,igcm_h2o_ice)=max(0.,qsurf(ig,igcm_h2o_ice))
-           qsurf(ig,igcm_hdo_ice)= DoverH * qsurf(ig,igcm_h2o_ice)
            end do
+
+        ! Update the hdo tracers
+         q(1:iip1,1:jjp1,1:llm,igcm_hdo_vap)
+     &      =q(1:iip1,1:jjp1,1:llm,igcm_h2o_vap)* DoverH
+         q(1:iip1,1:jjp1,1:llm,igcm_hdo_ice)
+     &      =q(1:iip1,1:jjp1,1:llm,igcm_h2o_ice)* DoverH
+
+         qsurf(1:ngridmx,igcm_hdo_ice)
+     &      =qsurf(1:ngridmx,igcm_h2o_ice)*DoverH
 
 
Index: trunk/LMDZ.MARS/libf/phymars/hdo_surfex_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/hdo_surfex_mod.F	(revision 2377)
+++ trunk/LMDZ.MARS/libf/phymars/hdo_surfex_mod.F	(revision 2378)
@@ -105,6 +105,6 @@
 
                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
+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
                else
                 alpha_c(ig) = 1.
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2377)
+++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2378)
@@ -609,4 +609,13 @@
          end do
        endif
+
+       if (txt.eq."hdo_ice") then
+         do ig=1,ngrid
+          if (qsurf(ig,iq).lt.0.0) then
+             qsurf(ig,iq) = 0.0
+          end if
+         end do
+       endif
+
       enddo
      endif ! of if (nq.ge.1)
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2377)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2378)
@@ -3032,5 +3032,5 @@
 
                 do ig=1,ngrid
-                if (qsurf(ig,igcm_h2o_ice).gt.1e-30) then
+                if (qsurf(ig,igcm_h2o_ice).gt.qperemin) then
                     DoH_surf(ig) = 0.5*( qsurf(ig,igcm_hdo_ice)/
      &                  qsurf(ig,igcm_h2o_ice) )/155.76e-6
Index: trunk/LMDZ.MARS/libf/phymars/simpleclouds.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/simpleclouds.F	(revision 2377)
+++ trunk/LMDZ.MARS/libf/phymars/simpleclouds.F	(revision 2378)
@@ -172,6 +172,6 @@
 
                 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
+c               alpha_c(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
+                alpha_c(ig,l) = exp(13525./zt(ig,l)**2.-5.59d-2)  !Lamb
                 else
                 alpha_c(ig,l) = 1.d0
