Index: trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3122)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3124)
@@ -296,5 +296,5 @@
       zq1temp_regolith(1:ngrid)=0 
       zdqsdif_tot(1:ngrid)=0
-!      h2o_ice_depth(1:ngrid,1:nslope)=5
+      !h2o_ice_depth(1:ngrid,1:nslope)=1
 c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
 c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
@@ -963,4 +963,5 @@
            saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap)    
            DO ig=1,ngrid
+            !nsubtimestep(ig)=1 for debug
             subtimestep = ptimestep/nsubtimestep(ig)
             ztsrf(ig)=ptsrf(ig,islope)  !  +pdtsrf(ig)*subtimestep
@@ -1152,4 +1153,37 @@
              IF(h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer
      &        .and. zqsurf(ig) .gt. 0) then
+                            DO ik=0,nsoil-2 ! go through all the layers to find the ice locations
+                IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and.
+     &           (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN
+                 lice = ik+1
+                EXIT
+                ENDIF
+               ENDDO !of subsurface loop
+                IF (lice .gt. 1) then !calculate the distance from the layers
+                 dist_up(ig,islope)=(h2o_ice_depth(ig,islope)
+     &                              -mlayer(lice-1))
+                 dist_down(ig,islope)=(mlayer(lice)
+     &                                -h2o_ice_depth(ig,islope))
+                 dist_sum(ig,islope)=dist_up(ig,islope)
+     &                               +dist_down(ig,islope)
+                 Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp
+     &                          *ptsoil(ig,lice-1,islope)
+     &           /dist_sum(ig,islope))+
+     &           (dist_down(ig,islope)*ptsoil(ig,lice,islope)
+     &           /dist_sum(ig,islope))
+                        ELSE
+                                Tice(ig,islope)=ptsoil(ig,1,islope)
+                        ENDIF
+
+             call watersat(1,Tice(ig,1),pplev(ig,1)
+     &                     ,qsat2(ig,1))
+              qeq(ig,1)=(ztsrf(ig)/Tice(ig,1))
+     &                        *qsat2(ig,1)
+             ! write(*,*)'icedep=',h2o_ice_depth(ig,1)
+             ! write(*,*)'qeq=',qeq(ig,1)
+             ! write(*,*)'d=',d_coef(ig,1)
+             ! write(*,*)'qsat=',qsat(ig)
+             ! write(*,*)'dry=',dryness(ig)
+             ! write(*,*)'rho=',rho(ig)
                zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)
      &                         /h2o_ice_depth(ig,1))
@@ -1161,4 +1195,32 @@
               ELSEIF (h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer
      &        .and. watercaptag(ig)) then
+                  DO ik=0,nsoil-2 ! go through all the layers to find the ice locations
+                    IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and.
+     &           (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN
+                 lice = ik+1
+                EXIT
+                ENDIF
+               ENDDO !of subsurface loop
+                IF (lice .gt. 1) then !calculate the distance from the layers
+                 dist_up(ig,islope)=(h2o_ice_depth(ig,islope)
+     &                              -mlayer(lice-1))
+                 dist_down(ig,islope)=(mlayer(lice)
+     &                                -h2o_ice_depth(ig,islope))
+                 dist_sum(ig,islope)=dist_up(ig,islope)
+     &                               +dist_down(ig,islope)
+                 Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp
+     &                          *ptsoil(ig,lice-1,islope)
+     &           /dist_sum(ig,islope))+
+     &           (dist_down(ig,islope)*ptsoil(ig,lice,islope)
+     &           /dist_sum(ig,islope))
+                        ELSE
+                                Tice(ig,islope)=ptsoil(ig,1,islope)
+                        ENDIF
+
+             call watersat(1,Tice(ig,1),pplev(ig,1)
+     &                     ,qsat2(ig,1))
+              qeq(ig,1)=(ztsrf(ig)/Tice(ig,1))
+     &                        *qsat2(ig,1)
+
               zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)
      &                         /h2o_ice_depth(ig,1))
@@ -1167,7 +1229,4 @@
                !needs to change to the mean of eq
               ENDIF
-              call write_output('zdqsdif_ssi_frost',
-     &                'Flux between frost and subsurface','kg.m-2.s-1',
-     &                zdqsdif_ssi_frost(ig,1))
               call write_output('subtimestep',
      &                'vdifc substimestep length','s',subtimestep)
@@ -1211,5 +1270,6 @@
 
 c             Fin du sous pas de temps
-            ENDDO ! tsub=1,nsubtimestep 
+            ENDDO ! tsub=1,nsubtimestep
+
 c             Integration of subtimestep temp and water budget :
 c             (btw could also compute the post timestep temp and ice
@@ -1342,10 +1402,16 @@
      &                     "Ground ice latent heat flux",
      &                     "W.m-2",surf_h2o_lh(:,iflat))
+         call write_output('zdqsdif_ssi_frost',
+     &          'Flux between frost and subsurface','kg.m-2.s-1',
+     &                zdqsdif_ssi_frost(:,1))
+         call write_output('subtimestep',
+     &          'vdifc substimestep length','s',subtimestep)
+
          call write_output('zdq_subtimestep',
-     &                     'Actual flux zdqsdif_surf*subtimestep',
-     &                     'kg.m-2',zdqsdif_surf(:)*subtimestep)
+     &          'Actual flux zdqsdif_surf*subtimestep',
+     &          'kg.m-2',zdqsdif_surf(:)*subtimestep)
          call write_output('zdq_end',
-     &                     'Flux after all contributions',
-     &                     'kg.m-2.s-1',zdqsdif_surf(:))
+     &          'Flux after all contributions',
+     &          'kg.m-2.s-1',zdqsdif_surf(:))
 C       Diagnostic output for HDO
 !        if (hdo) then