Index: /trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 3110)
+++ /trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 3111)
@@ -21,5 +21,5 @@
       use surfdat_h, only: albedo_h2o_cap,albedo_h2o_frost,
      &                     frost_albedo_threshold, inert_h2o_ice,
-     &                     frost_metam_threshold
+     &                     frost_metam_threshold,old_wsublimation_scheme
       use time_phylmdz_mod, only: ecritphy,day_step,iphysiq,ecritstart,
      &                            daysec,dtphys
@@ -890,4 +890,10 @@
      &            frost_albedo_threshold
 
+! TMP: old_wsublimation_scheme
+         write(*,*) "Old water sublimation scheme?"
+         old_wsublimation_scheme = .true.
+         call getin_p("old_wsublimation_scheme",old_wsublimation_scheme)
+         write(*,*) "old_wsublimation_scheme",old_wsublimation_scheme
+
 ! call Titus crocus line -- DEFAULT IS NONE
          write(*,*) "Titus crocus line ?"
Index: /trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90	(revision 3110)
+++ /trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90	(revision 3111)
@@ -20,5 +20,5 @@
     real, save, allocatable :: d_coef(:,:)  ! Diffusion coeficent
     LOGICAL,SAVE :: lag_layer ! does lag layer is present?
-!$OMP THREADPRIVATE(h2o_ice_depth,lag_co2_ice,albedo_perenialco2)
+!$OMP THREADPRIVATE(h2o_ice_depth,d_coef,lag_co2_ice,albedo_perenialco2)
 
     CONTAINS
Index: /trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90	(revision 3110)
+++ /trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90	(revision 3111)
@@ -31,8 +31,11 @@
   real,save :: z0_default ! default (constant over planet) surface roughness (m)
 
+  LOGICAL, SAVE :: old_wsublimation_scheme    ! TEMPORARY : TO USE THE OLD WATER SUBLIMATION SCHEME (i.e., using Cd instead of Ch), true by default
+
 !$OMP THREADPRIVATE(albedo_h2o_cap,albedo_h2o_frost,inert_h2o_ice,               &
 !$OMP                frost_albedo_threshold,frost_metam_threshold,TESice_Ncoef,  &
 !$OMP                TESice_Scoef,iceradius,dtemisice,                           &
-!$OMP                zmea,zstd,zsig,zgam,zthe,hmons,summit,base,z0,z0_default )
+!$OMP                zmea,zstd,zsig,zgam,zthe,hmons,summit,base,z0,z0_default, &
+!$OMP                old_wsublimation_scheme)
 
   !! mountain top dust flows
Index: /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3110)
+++ /trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3111)
@@ -20,5 +20,6 @@
      &                      igcm_hdo_vap, igcm_hdo_ice,
      &                      igcm_stormdust_mass, igcm_stormdust_number
-      use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness
+      use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness,
+     &                     old_wsublimation_scheme 
       USE comcstfi_h, ONLY: cpp, r, rcp, g, pi
       use watersat_mod, only: watersat
@@ -952,5 +953,4 @@
 c           the subtimestep.
            saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap)    
-
            DO ig=1,ngrid
             subtimestep = ptimestep/nsubtimestep(ig)
@@ -963,6 +963,11 @@
              zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
      &                     /float(nsubtimestep(ig))
-             zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
+             if(old_wsublimation_scheme) then
+               zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
      &                     /float(nsubtimestep(ig))
+             else
+               zb(1:ngrid,1)=zcdh(1:ngrid)*zb0(1:ngrid,1)
+     &                     /float(nsubtimestep(ig))
+             endif
              zb(1:ngrid,1)=dryness(1:ngrid)*zb(1:ngrid,1)
              
@@ -986,7 +991,11 @@
              zd(ig,1)=zb(ig,1)*z1(ig)
              zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
-
-             zdqsdif(ig)=rho(ig)*dryness(ig)*zcdv(ig)
-     &                      *(zq1temp(ig)-qsat(ig))
+             if(old_wsublimation_scheme) then
+                zdqsdif(ig)=rho(ig)*dryness(ig)*zcdv(ig)
+     &                           *(zq1temp(ig)-qsat(ig))
+             else
+                zdqsdif(ig)=rho(ig)*dryness(ig)*zcdh(ig)
+     &                           *(zq1temp(ig)-qsat(ig))
+             endif
 c             write(*,*)'subliming more than available frost:  qsurf!'
 
