Index: /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 2886)
+++ /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 2887)
@@ -25,6 +25,7 @@
   real,save,allocatable :: beta(:,:)        ! beta_k coefficients
   real,save :: mu
+  real,save,allocatable :: flux_geo(:)       ! Geothermal Flux (W/m^2)
 
-!$OMP THREADPRIVATE(tsoil,mthermdiff,thermdiff,coefq,coefd,alph,beta,mu )
+!$OMP THREADPRIVATE(tsoil,mthermdiff,thermdiff,coefq,coefd,alph,beta,mu,flux_geo)
 
 contains
@@ -47,5 +48,5 @@
     allocate(alph(ngrid,nsoilmx-1))
     allocate(beta(ngrid,nsoilmx-1))
-
+    allocate(flux_geo(ngrid))
  
   end subroutine ini_comsoil_h
@@ -66,5 +67,5 @@
     if (allocated(alph)) deallocate(alph)
     if (allocated(beta)) deallocate(beta)
-
+    if (allocated(flux_geo)) deallocate(flux_geo)
   end subroutine end_comsoil_h
 
Index: /trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2886)
+++ /trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2887)
@@ -24,5 +24,5 @@
   use dust_param_mod, only: dustscaling_mode
   USE ioipsl_getin_p_mod, ONLY : getin_p
-
+  use comsoil_h, only: flux_geo
   implicit none
   
@@ -644,4 +644,6 @@
   ! as well as thermal inertia and volumetric heat capacity
   call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
+else
+    flux_geo(:) = 0.
 endif ! of if (startphy_file)
 
Index: /trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2886)
+++ /trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2887)
@@ -162,5 +162,5 @@
   use dust_rad_adjust_mod, only: dust_rad_adjust_prev,dust_rad_adjust_next
   use dust_param_mod, only: dustscaling_mode
-
+  use comsoil_h,only: flux_geo
   implicit none
   
@@ -314,4 +314,7 @@
   endif
 
+
+  ! Geothermal Flux
+     call put_field('flux_geo','Geothermal flux',flux_geo,time)
   ! Close file
   call close_restartphy
Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2886)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2887)
@@ -3598,6 +3598,6 @@
 
          ! Write soil temperature
-!        call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
-!    &                     3,tsoil)
+        call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
+     &                     3,tsoil)
          ! Write surface temperature
 !        call writediagsoil(ngrid,"tsurf","Surface temperature","K",
@@ -3628,4 +3628,6 @@
 c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
 c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
+        call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
+     &                     3,tsoil)
 
 ! THERMALS STUFF 1D
@@ -3977,4 +3979,6 @@
      &                   log(zplay(1,nlayer)/zplay(1,nlayer-1))*
      &                   rnew(1,nlayer)*tmean/g
+
+
 #endif
 
Index: /trunk/LMDZ.MARS/libf/phymars/soil.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/soil.F	(revision 2886)
+++ /trunk/LMDZ.MARS/libf/phymars/soil.F	(revision 2887)
@@ -6,5 +6,5 @@
       use comsoil_h, only: layer, mlayer, volcapa,
      &                     mthermdiff, thermdiff, coefq,
-     &                     coefd, alph, beta, mu
+     &                     coefd, alph, beta, mu,flux_geo
       use surfdat_h, only: watercaptag, inert_h2o_ice
 
@@ -141,6 +141,8 @@
       do ig=1,ngrid
         beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil)
-     &                   /(coefq(nsoil-1)+coefd(ig,nsoil-1))
+     &            /(coefq(nsoil-1)+coefd(ig,nsoil-1))
+     &            +flux_geo(ig)/(coefq(nsoil-1) + coefd(ig,nsoil-1))
       enddo
+   
 ! Other layers
       do ik=nsoil-2,1,-1
Index: /trunk/LMDZ.MARS/libf/phymars/soil_settings.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/soil_settings.F	(revision 2886)
+++ /trunk/LMDZ.MARS/libf/phymars/soil_settings.F	(revision 2887)
@@ -2,5 +2,5 @@
 
 !      use netcdf
-      use comsoil_h, only: layer, mlayer, inertiedat, volcapa
+      use comsoil_h, only: layer, mlayer, inertiedat, volcapa,flux_geo
       use iostart, only: inquire_field_ndims, get_var, get_field,
      &                   inquire_field, inquire_dimension_length
@@ -435,6 +435,21 @@
 	deallocate(oldtsoil)
       endif ! of if (interpol)
-      
-! 6. Report min and max values of soil temperatures and thermal inertias
+
+
+      
+! 6. Load Geothermal Flux
+! ----------------------------------------------------------------------
+
+
+      call get_field("flux_geo",flux_geo,found,timeindex=indextime)
+      if (.not.found) then
+           write(*,*) "soil_settings: Failed loading <flux_geo>"
+           write(*,*) "soil_settings: Will put  <flux_geo> to 0."
+           flux_geo(:) = 0.
+      endif
+
+
+      
+! 7. Report min and max values of soil temperatures and thermal inertias
 ! ----------------------------------------------------------------------