Index: trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90	(revision 3612)
+++ trunk/LMDZ.MARS/libf/phymars/paleoclimate_mod.F90	(revision 3619)
@@ -16,4 +16,5 @@
     real,    allocatable, dimension(:,:) :: h2o_ice_depth         ! Thickness of the lag before H2O ice [m]
     real,    allocatable, dimension(:,:) :: lag_co2_ice           ! Thickness of the lag before CO2 ice [m]
+    real,    allocatable, dimension(:,:) :: zdqsdif_ssi_tot       ! Total flux of the interactions with SSI (kg/m^2/s^-1)
     real,    allocatable, dimension(:,:) :: d_coef                ! Diffusion coefficent
     logical                              :: lag_layer             ! Does lag layer is present?
@@ -34,4 +35,5 @@
 allocate(h2o_ice_depth(ngrid,nslope))
 allocate(lag_co2_ice(ngrid,nslope))
+allocate(zdqsdif_ssi_tot(ngrid,nslope))
 allocate(d_coef(ngrid,nslope))
 
@@ -45,4 +47,5 @@
 if (allocated(d_coef)) deallocate(d_coef)
 if (allocated(h2o_ice_depth)) deallocate(h2o_ice_depth)
+if (allocated(zdqsdif_ssi_tot)) deallocate(zdqsdif_ssi_tot)
 if (allocated(lag_co2_ice)) deallocate(lag_co2_ice)
 
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 3612)
+++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 3619)
@@ -31,5 +31,5 @@
 use comsoil_h,           only: flux_geo
 use comslope_mod,        only: nslope, major_slope
-use paleoclimate_mod,    only: paleoclimate, h2o_ice_depth, lag_co2_ice, d_coef
+use paleoclimate_mod,    only: paleoclimate, h2o_ice_depth, lag_co2_ice, d_coef, zdqsdif_ssi_tot
 use comcstfi_h,          only: pi
 use geometry_mod,        only: latitude
@@ -845,4 +845,12 @@
      h2o_ice_depth(:,:) = -1.
    endif
+   
+! Total flux with SSI
+   call get_field("zdqsdif_ssi_tot",zdqsdif_ssi_tot,found,indextime)
+   if (.not.found) then
+     write(*,*) "phyetat0: Failed loading <zdqsdif_ssi_tot> : ", &
+                          "<zdqsdif_ssi_tot> is set as -1 (no subsurface ice)"
+     zdqsdif_ssi_tot(:,:) = 0.
+   endif
 
 ! Diffusion coeficent
@@ -873,4 +881,5 @@
     h2o_ice_depth = -1.
     lag_co2_ice = -1.
+    zdqsdif_ssi_tot = 0.
     d_coef = 4.e-4
     perennial_co2ice = 0.
@@ -880,4 +889,5 @@
    h2o_ice_depth = -1.
    lag_co2_ice = -1.
+   zdqsdif_ssi_tot = 0.
    d_coef = 4.e-4
    perennial_co2ice = 0.
Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 3612)
+++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 3619)
@@ -28,5 +28,5 @@
   use time_phylmdz_mod, only: daysec
   use comslope_mod, ONLY: nslope
-  USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice, d_coef
+  USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice, d_coef, zdqsdif_ssi_tot
 
   implicit none
@@ -193,4 +193,5 @@
     call put_field("h2o_ice_depth","Depth to the fisrt H2O ice",h2o_ice_depth)
     call put_field("lag_co2_ice","Depth of the CO2 lags",lag_co2_ice)
+    call put_field("zdqsdif_ssi_tot","Total flux with SSI",zdqsdif_ssi_tot)
     call put_field("d_coef","Diffusion coefficent",d_coef)
   endif
Index: trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3612)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 3619)
@@ -34,5 +34,6 @@
      &                      subslope_dist,major_slope,iflat
       use microphys_h, only: To
-      use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer
+      use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer,
+     &                            zdqsdif_ssi_tot
       use comsoil_h, only: layer, mlayer,adsorption_soil
       use vdif_cd_mod, only: vdif_cd
@@ -243,5 +244,4 @@
       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 = 2e-3 ! tolerence for frost thicnkess (kg/m^2) to avoid numerical noise effect