Index: trunk/LMDZ.MARS/changelog.txt
===================================================================
--- trunk/LMDZ.MARS/changelog.txt	(revision 3252)
+++ trunk/LMDZ.MARS/changelog.txt	(revision 3253)
@@ -4541,2 +4541,6 @@
 == 04/03/2024 == LL
 Introduce the possibility to prescribe (in 1D) the soil water ice content when running with adsorption/desoprtion
+
+== 05/03/2024 == LL
+Following -r 3098; cleaning of vdfic for the management of subsurface water ice.
+Fixing some errors (wrong interpolation to compute the water ice temperature, wrong boundary conditions to compute qvap(1))
Index: trunk/LMDZ.MARS/deftank/field_def_physics_mars.xml
===================================================================
--- trunk/LMDZ.MARS/deftank/field_def_physics_mars.xml	(revision 3252)
+++ trunk/LMDZ.MARS/deftank/field_def_physics_mars.xml	(revision 3253)
@@ -533,6 +533,14 @@
                    long_name="co2 in the first layer"
                    unit="kg/kg" />
-
-
+            <!-- Subsurface ice interactions -->
+            <field id="zdqsdif_ssi_frost_tot"
+                   long_name="Flux between frost and subsurface ice"
+                   unit="kg.m-2.s-1" />
+            <field id="zdqsdif_ssi_atm_tot"
+                   long_name="Flux between atmosphere and subsurface ice"
+                   unit="kg.m-2.s-1" />
+             <field id="zdqsdif_ssi_tot"
+                   long_name="Total flux echanged with subsurface ice"
+                   unit="kg.m-2.s-1" />        
 
             <!-- CO2 cycle -->
@@ -1255,8 +1263,5 @@
             <field id="waterdensity_soil_slope07"
                    long_name="waterdensity_soil of slope 07"
-                   unit="kg.m-3" />
-            <field id="zdqsdif_ssi_frost"
-                   long_name="Flux between frost and subsurface"
-                   unit="kg.m-2.s-1" />
+                   unit="kg.m-3" />          
             <field id="subtimestep"
                    long_name="vdifc substimestep length"
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3252)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3253)
@@ -42,5 +42,5 @@
 use nonoro_gwd_ran_mod,       only: du_nonoro_gwd, dv_nonoro_gwd
 use conf_phys_mod,            only: conf_phys
-use paleoclimate_mod,         only: paleoclimate, ini_paleoclimate_h, end_paleoclimate_h
+use paleoclimate_mod,         only: paleoclimate, ini_paleoclimate_h, end_paleoclimate_h, h2o_ice_depth
 use comslope_mod,             only: nslope, subslope_dist, ini_comslope_h, end_comslope_h
 use co2condens_mod,           only: CO2cond_ps
@@ -541,5 +541,5 @@
     call getin("subsurface_ice_depth",ice_depth)
     write(*,*) " subsurface_ice_depth = ",ice_depth
-
+    h2o_ice_depth(:,:) = ice_depth
     z0(1) = z0_default ! default value for roughness
     write(*,*) 'Surface roughness length z0 (m)?'
@@ -668,5 +668,4 @@
 
 if (.not. therestartfi) qsoil = 0.
-
 
 if (.not. therestartfi) then
@@ -716,5 +715,5 @@
         inertiedat(1,:) = inertiedat(1,1) ! soil thermal inertia
     endif ! ice_depth > 0
-
+    
     do isoil = 1,nsoil
         inertiesoil(1,isoil,:) = inertiedat(1,isoil)
Index: trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3252)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3253)
@@ -128,11 +128,5 @@
       REAL zcst1
       REAL zu2(ngrid)
-      REAL Tice(ngrid,nslope) ! subsurface temperature where ice is located.
-      REAL qeq(ngrid,nslope) ! saturation water vapor in the subsurface
-      REAL dist_up(ngrid,nslope) !distance from ice to layer above
-      REAL dist_down(ngrid,nslope) !distance from ice to layer down
-      REAL dist_sum(ngrid,nslope) ! sum of distance
-      REAL zdqsdif_ssi(ngrid,nslope) !SSI flux
-      REAL zdqsdif_ssi_frost(ngrid,nslope) !SSI-frost interaction
+  
 
       EXTERNAL SSUM,SCOPY
@@ -172,6 +166,5 @@
       REAL rho(ngrid) ! near surface air density
       REAL qsat(ngrid)
-      REAL qsat2(ngrid,nslope)
-      REAL resist(ngrid,nslope) !subsurface ice flux reduction coef
+ 
 
       REAL hdoflux(ngrid,nslope)       ! value of vapour flux of HDO
@@ -244,6 +237,15 @@
       LOGICAL :: virtual
 
-
-
+!! Subsurface ice interactions
+      REAL Tice(ngrid,nslope)                  ! subsurface temperature where ice is located (K)
+      REAL qsat_ssi(ngrid,nslope)              ! qsat(Tice) (kg/kg)
+      REAL resist(ngrid,nslope)                ! subsurface ice flux reduction coef (1)
+      REAL zdqsdif_ssi_atm(ngrid,nslope)       ! SSI - atmosphere flux (kg/m^2/s^-1) at a given subtimestep
+      REAL zdqsdif_ssi_frost(ngrid,nslope)     ! SSI - frost flux (kg/m^2/s^-1) at a given subtimestep
+      REAL zdqsdif_ssi_atm_tot(ngrid,nslope)   ! SSI - atmosphere flux (kg/m^2/s^-1) summed over all subtimestep
+      REAL zdqsdif_ssi_frost_tot(ngrid,nslope) ! SSI - frost flux (kg/m^2/s^-1) summed over all subtimestep
+      REAL zdqsdif_ssi_tot(ngrid,nslope)       ! Total flux of the interactions with SSI kg/m^2/s^-1)
+      
+      REAL :: tol_frost = 1e-4 ! tolerence for frost thicnkess (kg/m^2) to avoid numerical noise effect
 c    ** un petit test de coherence
 c       --------------------------
@@ -306,7 +308,13 @@
       zq1temp_regolith(1:ngrid)=0 
       zdqsdif_tot(1:ngrid)=0
-      h2o_ice_depth(1:ngrid,1:nslope)=1
       virtual = .false.
       pore_icefraction(:,:,:) = 0.
+      zdqsdif_ssi_atm(:,:) = 0.
+      zdqsdif_ssi_frost(:,:) = 0.
+      zdqsdif_ssi_tot(:,:) = 0.
+      zdqsdif_ssi_atm_tot(:,:) = 0.
+      zdqsdif_ssi_frost_tot(:,:) = 0.
+      
+           
 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
@@ -983,8 +991,6 @@
 c --------- h2o_vap --------------------------------
 
-
-c      Traitement de la vapeur d'eau h2o_vap
-c      Utilisation d'un sous pas de temps afin
-c      de decrire le flux de chaleur latente 
+c Treatement of the water frost
+c We use a subtimestep to take into account the release of latent heat
          
         do iq=1,nq  
@@ -993,6 +999,4 @@
           DO islope = 1,nslope
            DO ig=1,ngrid
-
-             
 
              zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/
@@ -1001,25 +1005,23 @@
      &                       cos(pi*def_slope_mean(islope)/180.)
            ENDDO ! ig=1,ngrid
-
-c          make_tsub : sous pas de temps adaptatif
-c          la subroutine est a la fin du fichier
+           
+c             Computation of the subtimestep
               call make_tsub(ngrid,pdtsrf(:,islope),zqsurf,
      &                    ptimestep,dtmax,watercaptag,
-     &                    nsubtimestep)
-c           Calculation for turbulent exchange with the surface (for ice)
+     &                    nsubtimestep) 
+c           Calculation for turbulent exchange (rho Cd,h U (qatm - qsat(Tsurf)) with the surface (for ice)
 c           initialization of ztsrf, which is surface temperature in
 c           the subtimestep.
            saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap)    
            DO ig=1,ngrid
-!           nsubtimestep(ig)=1 !for debug
-            subtimestep = ptimestep/nsubtimestep(ig)
+c           nsubtimestep(ig)=1 !for debug
+           subtimestep = ptimestep/nsubtimestep(ig)
                           call write_output('subtimestep',
      &                'vdifc substimestep length','s',subtimestep)
             ztsrf(ig)=ptsrf(ig,islope)  !  +pdtsrf(ig)*subtimestep
             zq_tmp_vap(ig,:,:) =zq(ig,:,:)
-c           Debut du sous pas de temps
+c           Beginning of the subtimestep
             DO tsub=1,nsubtimestep(ig)
-              if(tsub.eq.nsubtimestep(ig)) writeoutput = .true.
-c           C'est parti ! 
+             if(tsub.eq.nsubtimestep(ig)) writeoutput = .true.
              zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
      &                     /float(nsubtimestep(ig))
@@ -1061,6 +1063,8 @@
 
              zdqsdif_tot(ig) = zdqsdif_surf(ig)
-!!! Subsurface exchange
-! Check for subsurface exchanges
+! --------------------------------------------------------------------------------------------------------------------------------              
+! We consider here the possible interactions between the subsurface and the atmosphere in the case of the surface is free of frost 
+! when computing the complete adsorption/desorption model
+! --------------------------------------------------------------------------------------------------------------------------------     
              if(.not.watercaptag(ig)) then
                 if (((-(zdqsdif_surf(ig))*
@@ -1074,6 +1078,8 @@
                   exchange = .false.
              endif
-                  zdqsdif_tot(ig) = zdqsdif_surf(ig)
-
+
+! --------------------------------------------------------------------------------------------------------------------------------              
+! If one consider adsorption, all the fluxes to/from the surface/subsurface/atmosphere are computed here
+! --------------------------------------------------------------------------------------------------------------------------------   
 
              if (adsorption_soil) then  
@@ -1087,5 +1093,4 @@
      &                     zq1temp_regolith(ig),
      &                     pore_icefraction(ig,:,islope))
-
 
 
@@ -1101,192 +1106,120 @@
                  endif  ! of "if not.watercaptag"
               endif ! adsorption          
-
-             if(.not.watercaptag(ig).and.(.not.adsorption_soil)) then
-               if ((-zdqsdif_tot(ig)*subtimestep)
-     &            .gt.(zqsurf(ig))) then
-                
+! --------------------------------------------------------------------------------------------------------------------------------------             
+! Here we do the same, but without computing the complete adsorpption/desorption. Note that it work only if one does not use adsorption
+! If no subsurface ice, then the models computes the surface flux/water vapor flux as usual
+! --------------------------------------------------------------------------------------------------------------------------------------     
+                     
+! ------------------------------------------------------------------------------------------------------------------------------------------             
+! First, We consider here the possible interactions between the subsurface and the atmosphere in the case of the surface is free of frost
+! ------------------------------------------------------------------------------------------------------------------------------------------ 
+
+              if(.not.watercaptag(ig).and.(.not.adsorption_soil)) then
+                 if ((-zdqsdif_tot(ig)*subtimestep) 
+     &                .ge.(zqsurf(ig))) then
+     
+                    zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep
+!                    zqsurf(ig) = 0
+                    if((h2o_ice_depth(ig,islope).gt.0)
+     &                  .and.lag_layer) then
+                  ! Atm <-> subsurface exchange, we need to update the exchange coefficient zb by a factor 1/1+R; R = zice*Cd,h*/D and add the flux from the subsurface
+                  
+                         if(old_wsublimation_scheme) then
+                            resist(ig,islope)=h2o_ice_depth(ig,islope)
+     &                          *zcdv(ig,islope)/d_coef(ig,islope)
+                         else
+                            resist(ig,islope)=h2o_ice_depth(ig,islope)
+     &                          *zcdh(ig,islope)/d_coef(ig,islope)                         
+                         endif
+                         
+                         zb(ig,1)=zb(ig,1)*1/(1+resist(ig,islope)) ! change zb to account subsurface ice
+                  ! Now we add the flux from the subsurface ice : rho Cd,h U*(1/1+R) * (q1-qsat_ssi(Tice))
+                    
+                        call compute_Tice(nsoil, ptsoil(ig,:,islope),
+     &                                    ztsrf(ig),
+     &                                    h2o_ice_depth(ig,islope),
+     &                                    Tice(ig,islope))        ! compute ice temperature
     
-                 !EV subsurface ice
-                IF(h2o_ice_depth(ig,islope) .gt. 0 .and. lag_layer)
-     &           then 
-                zdqsdif_tot(ig)=
-     &                        -zqsurf(ig)/subtimestep
-                zqsurf(ig)=0
-
-               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)
-              resist(ig,1)=(1+(h2o_ice_depth(ig,islope)*zcdh(ig,islope)
-     &                           /d_coef(ig,islope)))
-              !write(*,*)'R=',resist(ig,islope)
-              !write(*,*)'zice=',h2o_ice_depth(ig,islope)
-              !write(*,*)'D=',d_coef(ig,islope)
-              !write(*,*)'zcdh=',zcdh(ig)
-              zb(ig,1)=zb(ig,1)/resist(ig,1) ! change zb to account subsurface ice
-              !vdifc algorithem !!!!!needs to change  reseist io to the mean
-              !beacuse of the slopes!!!!!
-              z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
-             zc(ig,nlay)=za(ig,nlay)*zq_tmp_vap(ig,nlay,iq)*z1(ig)
-             zd(ig,nlay)=zb(ig,nlay)*z1(ig)
-             DO ilay=nlay-1,2,-1
-                 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
-     $            zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
-                 zc(ig,ilay)=(za(ig,ilay)*zq_tmp_vap(ig,ilay,iq)+
-     $            zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
-                 zd(ig,ilay)=zb(ig,ilay)*z1(ig)
-             ENDDO
-             z1(ig)=1./(za(ig,1)+zb(ig,1)+
-     $          zb(ig,2)*(1.-zd(ig,2)))
-             zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
-     $        zb(ig,2)*zc(ig,2)) * z1(ig)
-             zd(ig,1)=zb(ig,1)*z1(ig)
-             zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
-             zdqsdif_ssi(ig,1)=rho(ig)*dryness(ig)*zcdv(ig,islope)    
-     &                      *(zq1temp(ig)-qeq(ig,islope))
-                      call write_output('zdq_ssi',
-     &                '','',zdqsdif_ssi(ig,1))
-                       call write_output('qeq',
-     &                '','',qeq(ig,1))
-                       call write_output('q1',
-     &                '','',zq1temp(ig))
-             !write(*,*)'qeq=',qeq(ig,1)
-             !write(*,*)'q1=',zq1temp(ig) 
-             !write(*,*)'zdqsdif_ssi=',zdqsdif_ssi(ig,1)
-             !I should some all the interactions with the SSI accoridng
-             !to the statistics and the evaluate
+                        call watersat(1,Tice(ig,islope),pplev(ig,1),
+     &                                qsat_ssi(ig,islope))
+                 
+                         qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
+     &                                       *qsat_ssi(ig,islope)
+                  ! Flux from the subsurface      
+                        if(old_wsublimation_scheme) then
+                           zdqsdif_ssi_atm(ig,islope)=rho(ig)
+     &                     *dryness(ig)*zcdv(ig,islope)
+     &                     *1/(1+resist(ig,islope))    
+     &                     *(zq1temp(ig)-qsat_ssi(ig,islope))
+                        else
+                           zdqsdif_ssi_atm(ig,islope)=rho(ig)*
+     &                     *dryness(ig) *zcdh(ig,islope) 
+     &                     *1/(1+resist(ig,islope))   
+     &                     *(zq1temp(ig)-qsat_ssi(ig,islope))                         
+                        endif
+                  ! And now we solve correctly the system      
+                        z1(ig)=1./(za(ig,1)+zb(ig,1)+
+     &                         zb(ig,2)*(1.-zd(ig,2)))
+                        zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
+     &                           zb(ig,2)*zc(ig,2) + 
+     &                           (-zdqsdif_tot(ig)) *subtimestep)
+     &                           * z1(ig)
+                        zd(ig,1)=zb(ig,1)*z1(ig)
+                        zq1temp(ig)=zc(ig,1)+ zd(ig,1)
+     &                              *qsat_ssi(ig,islope)
+            
+
+
+                    else 
+                  ! No atm <-> subsurface exchange, we do it the usual way
+                        zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep
+                        z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
+                        zc(ig,1)=(za(ig,1)*
+     &                      zq_tmp_vap(ig,1,igcm_h2o_vap)+
+     &                      zb(ig,2)*zc(ig,2) +
+     &                      (-zdqsdif_tot(ig)) *subtimestep) *z1(ig)
+                            zq1temp(ig)=zc(ig,1)
+                    endif ! Of subsurface <-> atmosphere exchange
+                 endif ! sublimating more than available frost & surface - frost exchange 
+             endif !if .not.watercaptag(ig) 
+
+! --------------------------------------------------------------------------------------------------------------------------------              
+! We check possible frost subsurface ice interaction: since there is no subsurface water ice mass reservoir represented (yet), 
+! we do not include their effect on the mass of surface frost.
+! -------------------------------------------------------------------------------------------------------------------------------- 
+
+             if((h2o_ice_depth(ig,islope).gt.0).and.lag_layer
+     &                      .and.(.not.adsorption_soil)) then
+             ! First case: still frost at the surface but no watercaptag
+                 if(((watercaptag(ig))).or.
+     &                (((-zdqsdif_tot(ig)*subtimestep) 
+     &                .lt.(zqsurf(ig)))
+     &                .and. (zqsurf(ig).gt.tol_frost))) then
+                  ! Still frost at the surface:  we consider the possibility to have subsurface <-> frost exchange
+                  ! The flux between frost and ssi is D/zice *(qsat(Tsurf)-qsat_ssi(Tice))    
+                    call compute_Tice(nsoil, ptsoil(ig,:,islope),
+     &                                ztsrf(ig),
+     &                                h2o_ice_depth(ig,islope),
+     &                                Tice(ig,islope))          ! compute ice temperature
+    
+                    call watersat(1,Tice(ig,islope),pplev(ig,1),
+     &                            qsat_ssi(ig,islope))
+                 
+                    qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
+     &                                       *qsat_ssi(ig,islope)
+     
+                    zdqsdif_ssi_frost(ig,islope)= d_coef(ig,islope)
+     &                                 /h2o_ice_depth(ig,islope)
+     &                                 *rho(ig)*dryness(ig)
+     &                                 *(qsat(ig)-qsat_ssi(ig,islope))
+
+! Line to comment for now since we don't have a mass of subsurface frost in our computation (otherwise, we would not conserve the H2O mass in the system)               
+                    zdqsdif_tot(ig) = zdqsdif_tot(ig) - 
+     &                        zdqsdif_ssi_frost(ig,islope) 
+                 endif ! watercaptag or frost at the surface
+             endif ! interaction frost <-> subsurface ice
              
-             !write(*,*) "zdq_ssi*t", zdqsdif_ssi(ig)*subtimestep
-             if (zdqsdif_ssi(ig,1)<0) then
-               !zdqsdif_surf(ig)=zdqsdif_ssi(ig,1)-(zqsurf(ig)/subtimestep)
-               zdqsdif_tot(ig)=zdqsdif_tot(ig)+zdqsdif_ssi(ig,1)
-                call write_output('zdq_zdqssi',
-     &                '','',zdqsdif_tot(ig)+zdqsdif_ssi(ig,1))
-             endif
-              ELSE
-c             pdqsdif > 0 : ice condensing
-c             pdqsdif < 0 : ice subliming
-c             write(*,*) "subliming more than available frost:  qsurf!"
-                  zdqsdif_tot(ig)=
-     &                        -zqsurf(ig)/subtimestep
-c                write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
-                 z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
-                 zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,igcm_h2o_vap)+
-     $           zb(ig,2)*zc(ig,2) +
-     $           (-zdqsdif_tot(ig)) *subtimestep) *z1(ig)
-                 zq1temp(ig)=zc(ig,1)
-                 ENDIF !if h2o_ice_depth>0 and lag_layer
-               endif  !if .not.watercaptag(ig) 
-             endif ! if sublim more than surface
-
-              ! subsurface ice <--> frost interaction !EV
-             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))
-     &          *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))
-               !needs to change to the mean of eq
-               zdqsdif_tot(ig)=zdqsdif_tot(ig)-zdqsdif_ssi_frost(ig,1) 
-                !!!!! zdsqdif_ssi_frosst need to be changed to an
-                !average
-              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))
-     &          *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))
-               zdqsdif_tot(ig)=zdqsdif_tot(ig)-zdqsdif_ssi_frost(ig,1)
-               !needs to change to the mean of eq
-              ENDIF
-!              call write_output('subtimestep',
-!     &                'vdifc substimestep length','s',subtimestep)
-             ! ENDDO !subsurface ice subslope
-
-
 
 c             Starting upward calculations for water :
-c             Actualisation de h2o_vap dans le premier niveau
              zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
 c             Take into account the H2O latent heat impact on the surface temperature
@@ -1307,9 +1240,13 @@
               zqsurf(ig)= zqsurf(ig)+(
      &                       zdqsdif_tot(ig))*subtimestep
-              if (zqsurf(ig)<0 .and.
-     &        (.not.watercaptag(ig))) then
+              if (zqsurf(ig)<0 .and.(.not.watercaptag(ig))) then
                       zqsurf(ig)=0
               endif
-
+              zdqsdif_ssi_atm_tot(ig,islope) = 
+     &                              zdqsdif_ssi_atm_tot(ig,islope)
+     &                            + zdqsdif_ssi_atm(ig,islope)
+              zdqsdif_ssi_frost_tot(ig,islope) = 
+     &                                zdqsdif_ssi_frost_tot(ig,islope) 
+     &                              + zdqsdif_ssi_frost(ig,islope)
 c             Monitoring instantaneous latent heat flux in W.m-2 :
               zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+
@@ -1326,5 +1263,5 @@
                endif
 
-c             Fin du sous pas de temps
+c             End of the subtimestep
             ENDDO ! tsub=1,nsubtimestep
 
@@ -1339,4 +1276,8 @@
      &                       cos(pi*def_slope_mean(islope)/180.))
      &                       /ptimestep
+     
+            zdqsdif_ssi_tot(ig,islope) = 
+     &                        zdqsdif_ssi_atm_tot(ig,islope) 
+     &                      + zdqsdif_ssi_frost_tot(ig,islope)
 c if subliming more than qsurf(ice) and on watercaptag, water
 c sublimates from watercap reservoir (dwatercap_dif is <0)
@@ -1459,14 +1400,18 @@
      &                     "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('zdq_subtimestep',
-!     &          'Actual flux zdqsdif_surf*subtimestep',
-!     &          'kg.m-2',zdqsdif_tot(:)*subtimestep)
-!         call write_output('zdq_end',
-!     &          'Flux after all contributions',
-!     &          'kg.m-2.s-1',zdqsdif_tot(:))
+
+         call write_output('zdqsdif_ssi_frost_tot',
+     &     'Flux between frost and subsurface ice','kg.m-2.s-1',
+     &                zdqsdif_ssi_frost_tot(:,iflat))
+
+         call write_output('zdqsdif_ssi_atm_tot',
+     &     'Flux between atmosphere and subsurface ice','kg.m-2.s-1',
+     &                zdqsdif_ssi_atm_tot(:,iflat))
+
+         call write_output('zdqsdif_ssi_tot',
+     &     'Total flux echange with subsurface ice','kg.m-2.s-1',
+     &                zdqsdif_ssi_tot(:,iflat))
+
+
 C       Diagnostic output for HDO
 !        if (hdo) then
@@ -1553,3 +1498,53 @@
 
       END SUBROUTINE make_tsub
+      
+      
+c====================================
+
+      SUBROUTINE compute_Tice(nsoil, ptsoil, ptsurf, ice_depth, Tice)
+
+c Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells.
+      use comsoil_h, only: layer, mlayer
+
+      implicit none 
+      integer,intent(in) :: nsoil            ! Number of soil layers
+      real,intent(in) :: ptsoil(nsoil)       ! Soil temperature (K)
+      real,intent(in) :: ptsurf              ! Soil temperature (K)
+      real,intent(in) :: ice_depth           ! Ice depth (m)
+      real,intent(out) :: Tice               ! Ice temperatures (K)
+
+c Local
+      integer :: ik                          ! Loop variables 
+      integer :: indexice                    ! Index of the ice
+
+c Code:
+      indexice = -1
+      if(ice_depth.lt.mlayer(0)) then
+          indexice = 0.
+      else 
+          do ik = 0,nsoil-2 ! go through all the layers to find the ice locations
+             if((mlayer(ik).le.ice_depth).and.
+     &          (mlayer(ik+1).gt.ice_depth)) then
+                    indexice = ik+1
+                    exit
+             endif
+          enddo
+      endif
+      
+      if(indexice.lt.0) then
+          call abort_physic("vdifc - compute Tice",
+     &         "subsurface ice is below the last soil layer",1)
+      else
+          if(indexice .ge. 1) then ! Linear inteprolation between soil temperature
+              Tice = (ptsoil(indexice)-ptsoil(indexice+1))
+     &               /(mlayer(indexice-1)-mlayer(indexice))
+     &               *(ice_depth-mlayer(indexice)) + ptsoil(indexice+1)
+          else ! Linear inteprolation between the 1st soil temperature and the surface temperature
+              Tice = (ptsoil(1) - ptsurf)/mlayer(0)
+     &               *(ice_depth-mlayer(0)) + ptsoil(1)
+          endif ! index ice >=0
+      endif !indexice <0
+      
+
+      END SUBROUTINE compute_Tice
       END MODULE vdifc_mod