Index: trunk/LMDZ.COMMON/libf/evolution/ComputeCO2Adsorption.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/ComputeCO2Adsorption.F90	(revision 2850)
+++ 	(revision )
@@ -1,96 +1,0 @@
-   SUBROUTINE regolith_co2adsorption_old(ngrid,nslope,nsoil_PEM,ps,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,m_co2_completesoil,delta_mreg)
-#ifndef CPP_STD
-      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km
-      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad 
-      use comslope_mod, only : subslope_dist,def_slope_mean
-      IMPLICIT NONE
-! Input: 
- INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM
- REAL,INTENT(IN) :: ps(ngrid)                                      ! Average surface pressure [Pa]
- REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)               ! Mean Soil Temperature [K]
- REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope)                  ! Mean Thermal Inertia [USI]
- REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) !tendancies on the glaciers
- REAL,INTENT(IN) :: waterice(ngrid,nslope)                          ! water ice at the surface [KG/m^2]
- REAL,INTENT(IN) :: co2ice(ngrid,nslope)                            ! co2 ice at the surface [Kg/m^2] 
-
-! Outputs:
- REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)
- REAL,INTENT(INOUT) :: delta_mreg(ngrid)
-  
-! Constants:
-
- REAL ::  as = 1.9e4      ! Bulher & Piqueux 2020, assuming exchange only between 0 and 868m, last layer of the PEM
- REAL ::  delta = 5.749e1 ! Zent & Quinn, 1995
- REAL ::  beta =  -4.0711 ! Zent & Quinn, 1995
- REAL ::  gama =  0.2788  ! Zent & Quinn, 1995
- REAL ::  inertie_thresold = 2000. ! look for ice
- REAL ::  rho_regolith = 2000. ! density of the reoglith, buhler and piqueux 2021
-! Input: Local 
- INTEGER :: ig,islope,iloop,iref,k
- REAL :: dm_co2_regolith_slope(ngrid,nsoil_PEM,nslope)   ! elementary mass adsorded per mesh per slope
- INTEGER :: ispermanent_co2glaciers(ngrid,nslope)
- INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)
- REAL :: deltam_reg_complete(ngrid,n_1km,nslope)
- REAL :: deltam_reg_slope(ngrid,nslope)
-
-dm_co2_regolith_slope(:,:,:) = 0
-delta_mreg(:) = 0.
-!0. Initialization
-  do ig = 1,ngrid
-    do islope = 1,nslope
-     if((abs(tend_h2oglaciers(ig,islope)).gt.1e-5).and.(abs(waterice(ig,islope)).gt.0)) then
-        ispermanent_h2oglaciers(ig,islope) = 1
-     else
-        ispermanent_h2oglaciers(ig,islope) = 0
-     endif
-
-     if((abs(tend_co2glaciers(ig,islope)).gt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then
-        ispermanent_co2glaciers(ig,islope) = 1
-     else
-        ispermanent_co2glaciers(ig,islope) = 0
-     endif
-    enddo
-   enddo
-! 1. we compute the mass of co2 adsorded in each layer of the meshes  
-
- do ig = 1,ngrid
-  do islope = 1,nslope
-    do iloop = 1,n_1km
-     if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
-     dm_co2_regolith_slope(ig,iloop,islope) = rho_regolith*as*delta*(ps(ig)**gama)*(tsoil_PEM(ig,iloop,islope))**beta 
-     else
-        if(abs(m_co2_completesoil(ig,iloop,islope)).lt.1-10) then !!! we are at first call
-          dm_co2_regolith_slope(ig,iloop,islope) = rho_regolith*as*delta*(ps(ig)**gama)*(tsoil_PEM(ig,iloop,islope))**beta 
-        else ! no change: permanent ice stick the atoms of CO2
-          dm_co2_regolith_slope(ig,iloop,islope) = m_co2_completesoil(ig,iloop,islope)
-        endif
-     endif
-  enddo
- enddo
-enddo
-
-! 2. Check the exchange between the atmosphere and the regolith
-
-  do ig = 1,ngrid
-   delta_mreg(ig) = 0.
-   do islope = 1,nslope
-    deltam_reg_slope(ig,islope) = 0.
-    do iloop = 1,n_1km
-       if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
-             deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope)) &
-                                                   *(layer_PEM(iloop+1) - layer_PEM(iloop)) 
-       else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
-             deltam_reg_complete(ig,iloop,islope) = 0.
-       endif
-       deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
-    enddo
-   delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
-   enddo
-  enddo
-   m_co2_completesoil(:,:,:) = dm_co2_regolith_slope(:,:,:)
-
-
-!=======================================================================
-#endif
-      RETURN
-      END
Index: trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90	(revision 2855)
@@ -1,16 +1,24 @@
    SUBROUTINE computeice_table(ngrid,nslope,nsoil_PEM,rhowatersurf_ave,rhowatersoil_ave,ice_table)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Compute the ice table depth knowing the yearly average water
+!!! density at the surface and at depth.
+!!! Computations are made following the methods in Schorgofer et al., 2005
+!!! 
+!!! Author: LL
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
 #ifndef CPP_STD
-    USE comsoil_h, only:  inertiedat, volcapa
-    USE vertical_layers_mod, ONLY: ap,bp
-    USE comsoil_h_PEM, only: layer_PEM 
+    USE comsoil_h_PEM, only: layer_PEM                             ! Depth of the vertical grid 
     implicit none 
 
-    integer,intent(in) :: ngrid,nslope,nsoil_PEM
-    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)    ! surface temperature, timeseries [K]
-    real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope)    ! soil water density, yearly average [kg/m^3]
-    real,intent(inout) :: ice_table(ngrid,nslope)                 ! ice table [m]
-    real :: z1,z2
-    integer ig, islope,isoil
-    real :: diff_rho(nsoil_PEM)                    ! difference of vapor content
+    integer,intent(in) :: ngrid,nslope,nsoil_PEM                   ! Size of the physical grid, number of subslope, number of soil layer in the PEM 
+    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)              ! Water density at the surface, yearly averaged [kg/m^3]
+    real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope)    ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged  [kg/m^3]
+    real,intent(inout) :: ice_table(ngrid,nslope)                  ! ice table depth [m]
+    real :: z1,z2                                                  ! intermediate variables used when doing a linear interpolation between two depths to find the root 
+    integer ig, islope,isoil                                       ! loop variables
+    real :: diff_rho(nsoil_PEM)                                    ! difference of water vapor density between the surface and at depth [kg/m^3]
 
 
@@ -24,8 +32,8 @@
          enddo
 
-         if(diff_rho(1) > 0) then
+         if(diff_rho(1) > 0) then                                   ! ice is at the surface
            ice_table(ig,islope) = 0.
          else
-           do isoil = 1,nsoil_PEM -1
+           do isoil = 1,nsoil_PEM -1                                ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
              if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
                z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
@@ -38,10 +46,6 @@
         enddo
       enddo
-
-    
-
 !=======================================================================
       RETURN
-
 #endif
       END
Index: trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 2855)
@@ -2,24 +2,23 @@
 
 implicit none
-! nsoilmx : number of subterranean layers
-!EM: old soil routine:      integer, parameter :: nsoilmx = 10
-  integer, parameter :: nsoilmx_PEM = 27 
-  integer, parameter :: n_1km = 23
+  integer, parameter :: nsoilmx_PEM = 27               ! number of layers in the PEM
+  integer, parameter :: n_1km = 23                     ! index at which we have overcome 1km
 
-  real,save,allocatable,dimension(:) :: layer_PEM      ! soil layer depths
-  real,save,allocatable,dimension(:) :: mlayer_PEM     ! soil mid-layer depths
-  real,save,allocatable,dimension(:,:,:) :: TI_PEM ! soil thermal inertia
-    real,save,allocatable,dimension(:,:) :: inertiedat_PEM ! soil thermal inertia
+  real,save,allocatable,dimension(:) :: layer_PEM      ! soil layer depths   [m]
+  real,save,allocatable,dimension(:) :: mlayer_PEM     ! soil mid-layer depths [m]
+  real,save,allocatable,dimension(:,:,:) :: TI_PEM     ! soil thermal inertia [SI]
+    real,save,allocatable,dimension(:,:) :: inertiedat_PEM ! soil thermal inertia saved as reference for current climate [SI]
   ! variables (FC: built in firstcall in soil.F)
-  REAL,SAVE,ALLOCATABLE :: tsoil_PEM(:,:,:)       ! sub-surface temperatures (K)
-  real,save,allocatable :: mthermdiff_PEM(:,:)  ! (FC) mid-layer thermal diffusivity
-  real,save,allocatable :: thermdiff_PEM(:,:)   ! (FC) inter-layer thermal diffusivity
-  real,save,allocatable :: coefq_PEM(:)         ! (FC) q_{k+1/2} coefficients
-  real,save,allocatable :: coefd_PEM(:,:)       ! (FC) d_k coefficients
-  real,save,allocatable :: alph_PEM(:,:,:)        ! (FC) alpha_k coefficients
-  real,save,allocatable :: beta_PEM(:,:,:)        ! beta_k coefficients
-  real,save :: mu_PEM
-  real,parameter :: fluxgeo = 30e-3 !W/m^2
-  real, save, allocatable :: co2_adsorbded_phys(:,:,:)  ! co2 that is in the regolith (kg/m^2)
+  REAL,SAVE,ALLOCATABLE :: tsoil_PEM(:,:,:)       ! sub-surface temperatures [K]
+  real,save,allocatable :: mthermdiff_PEM(:,:)  ! (FC) mid-layer thermal diffusivity [SI]
+  real,save,allocatable :: thermdiff_PEM(:,:)   ! (FC) inter-layer thermal diffusivity [SI]
+  real,save,allocatable :: coefq_PEM(:)         ! (FC) q_{k+1/2} coefficients [SI]
+  real,save,allocatable :: coefd_PEM(:,:)       ! (FC) d_k coefficients [SI]
+  real,save,allocatable :: alph_PEM(:,:,:)      ! (FC) alpha_k coefficients [SI]
+  real,save,allocatable :: beta_PEM(:,:,:)      ! beta_k coefficients [SI]
+  real,save :: mu_PEM                           ! mu coefficient [SI]
+  real,parameter :: fluxgeo = 30e-3             ! Geothermal flux [W/m^2]
+  real, save, allocatable :: co2_adsorbded_phys(:,:,:)  ! co2 that is in the regolith [kg/m^2]
+  LOGICAL soil_pem  ! True by default, to run with the subsurface physic. Read in pem.def
 
 contains
Index: trunk/LMDZ.COMMON/libf/evolution/ini_soil_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/ini_soil_mod.F90	(revision 2850)
+++ 	(revision )
@@ -1,362 +1,0 @@
-      MODULE ini_soil_mod
-
-      IMPLICIT NONE
-      
-      CONTAINS    
-
-     subroutine ini_icetable(timelen,ngrid,nsoil_PEM, &
-               therm_i, timestep,tsurf_ave,tsoil_ave,tsurf_inst, tsoil_inst,q_co2,q_h2o,ps,ice_table)
-
-#ifndef CPP_STD
-
-      use vertical_layers_mod, only: ap,bp
-      use comsoil_h_PEM, only:  fluxgeo,layer_PEM,inertiedat_PEM
-      use comsoil_h,only: volcapa, nsoilmx
-#endif
-
-   implicit none
-
-!-----------------------------------------------------------------------
-!  Author: LL
-!  Purpose: Compute soil temperature using an implict 1st order scheme
-!  
-!  Note: depths of layers and mid-layers, soil thermal inertia and 
-!        heat capacity are commons in comsoil_PEM.h
-!        A convergence loop is added until equilibrium
-!-----------------------------------------------------------------------
-
-#include "dimensions.h"
-
-!-----------------------------------------------------------------------
-!  arguments
-!  ---------
-!  inputs:
-      integer,intent(in) :: timelen	! Time length in for time-series data
-      integer,intent(in) :: ngrid	! number of (horizontal) grid-points 
-      integer,intent(in) :: nsoil_PEM	! number of soil layers  in the PEM
-
-      real,intent(in) :: timestep	    ! time step
-      real,intent(in) :: tsurf_ave(ngrid)   ! surface temperature
-
-      real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
-      real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
-      real,intent(in) :: ps(ngrid,timelen)                        ! surface pressure [Pa]
-      real,intent(in) :: tsurf_inst(ngrid,timelen) ! soil (mid-layer) temperature
-
-! outputs:
-      real,intent(inout) :: tsoil_ave(ngrid,nsoil_PEM) ! soil (mid-layer) temperature. 
-      real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,timelen) ! soil (mid-layer) temperature
-      real,intent(out) :: ice_table(ngrid)                 ! ice table [m]
-      real,intent(inout) :: therm_i(ngrid,nsoil_PEM) ! thermal inertia
-
-! local variables:
-      integer ig,isoil,it,k,iref
-      REAL :: error_depth = 2.
-      REAL :: tsoil_saved(nsoil_PEM)
-      integer :: countmax = 20
-      integer :: countloop 
-      REAL :: tol_error = 0.1
-      REAL :: ice_inertia = 2120.
-      REAL :: alph_PEM(nsoil_PEM-1)
-      REAL :: beta_PEM(nsoil_PEM-1)
-      real :: rhoc(nsoil_PEM)
-      real :: volcapa_ice = 1.43e7
-
-      real :: k_soil(nsoil_PEM)
-      real :: d1(nsoil_PEM),d2(nsoil_PEM),d3(nsoil_PEM),re(nsoil_PEM)
-      real :: Tcol_saved(nsoil_PEM)
-      real :: tsoil_prev(nsoil_PEM)
-      real :: tsurf_prev
-      real :: icedepth_prev
-
-    real :: m_h2o = 18.01528E-3
-    real :: m_co2 = 44.01E-3  
-    real :: m_noco2 = 33.37E-3  
-    real :: A,B,z1,z2
-    real :: alpha = -6143.7 
-    real :: beta = 29.9074
-   
-    real,allocatable :: mass_mean(:)                            ! mean mass above the surface
-    real,allocatable :: zplev(:)                           ! pressure above the surface
-    real,allocatable :: pvapor(:)                               ! partial pressure above the surface
-   
-    real,allocatable :: rhovapor(:)
-    real :: rhovapor_avg                          ! mean vapor_density at the surface yearly averaged
-
-    real :: psv_surf
-    real,allocatable :: rho_soil(:)            ! water vapor in the soil
-    real,allocatable :: rho_soil_avg(:)                ! water vapor in the soil yearly averaged
-
-    real,allocatable :: diff_rho(:)                    ! difference of vapor content
-
-#ifndef CPP_STD
-
-    A =(1/m_co2 - 1/m_noco2)
-    B=1/m_noco2
-
-    ice_table(:) = 1.
-    do ig = 1,ngrid
-      error_depth = 1.
-      countloop = 0
-    
-      do while(( error_depth.gt.tol_error).and.(countloop.lt.countmax).and.(ice_table(ig).gt.-1e-20))
-    
-        countloop = countloop +1
-        Tcol_saved(:) = tsoil_ave(ig,:)
-
-!2.  Compute ice table
-
-! 2.1 Compute  water density at the surface, yearly averaged
-        allocate(mass_mean(timelen))
-!   1.1 Compute the partial pressure of vapor
-! a. the molecular mass into the column
-       mass_mean(:) = 1/(A*q_co2(ig,:) +B)
-! b. pressure level
-        allocate(zplev(timelen))
-        do it = 1,timelen
-         zplev(it) = ap(1) + bp(1)*ps(ig,it)
-        enddo
-
-! c. Vapor pressure
-        allocate(pvapor(timelen))
-        pvapor(:) = mass_mean(:)/m_h2o*q_h2o(ig,:)*zplev(:)
-        deallocate(zplev) 
-        deallocate(mass_mean)
-  
-!  d!  Check if there is frost at the surface and then compute the density
-!    at the surface
-        allocate(rhovapor(timelen))
-
-         do it = 1,timelen
-           psv_surf = exp(alpha/tsurf_inst(ig,it) +beta)
-           rhovapor(it) = min(psv_surf,pvapor(it))/tsurf_inst(ig,it)
-         enddo
-        deallocate(pvapor)
-        rhovapor_avg =SUM(rhovapor(:),1)/timelen                     
-        deallocate(rhovapor)
-
-! 2.2 Compute  water density at the soil layer, yearly averaged
-
-        allocate(rho_soil(timelen))
-        allocate(rho_soil_avg(nsoil_PEM))
-
-        do isoil = 1,nsoil_PEM
-          do it = 1,timelen 
-            rho_soil(it) = exp(alpha/tsoil_inst(ig,isoil,it) +beta)/tsoil_inst(ig,isoil,it)        
-          enddo
-         rho_soil_avg(isoil) = SUM(rho_soil(:),1)/timelen     
-        enddo
-        deallocate(rho_soil)
-         
-!2.3 Final: compute ice table
-        icedepth_prev = ice_table(ig)
-        ice_table(ig) = -1
-        allocate(diff_rho(nsoil_PEM))
-        do isoil = 1,nsoil_PEM
-          diff_rho(isoil) = rhovapor_avg -  rho_soil_avg(isoil)
-        enddo
-        deallocate(rho_soil_avg)
-
-        if(diff_rho(1) > 0) then
-          ice_table(ig) = 0.
-        else
-          do isoil = 1,nsoil_PEM -1
-            if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
-              z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
-              z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
-              ice_table(ig) = -z2/z1
-              exit
-            endif
-          enddo
-        endif
-    
-        deallocate(diff_rho)
-
-!3. Update Soil Thermal Inertia
-
-        if (ice_table(ig).gt. 0.) then
-          if (ice_table(ig).lt. 1e-10) then
-            do isoil = 1,nsoil_PEM
-              therm_i(ig,isoil)=ice_inertia
-            enddo
-          else 
-      ! 4.1 find the index of the mixed layer
-          iref=0 ! initialize iref
-          do k=1,nsoil_PEM ! loop on layers
-            if (ice_table(ig).ge.layer_PEM(k)) then
-              iref=k ! pure regolith layer up to here
-            else
-              ! correct iref was obtained in previous cycle
-            exit
-            endif
-          enddo
-
-      ! 4.2 Build the new ti
-         do isoil=1,iref
-           therm_i(ig,isoil) =inertiedat_PEM(ig,isoil) 
-         enddo
-
-         if (iref.lt.nsoil_PEM) then
-           if (iref.ne.0) then
-             ! mixed layer
-             therm_i(ig,iref+1)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
-              (((ice_table(ig)-layer_PEM(iref))/(inertiedat_PEM(ig,iref)**2))+ &
-                      ((layer_PEM(iref+1)-ice_table(ig))/(ice_inertia**2))))
-            
-           else ! first layer is already a mixed layer
-            ! (ie: take layer(iref=0)=0)
-            therm_i(ig,1)=sqrt((layer_PEM(1))/ &
-                          (((ice_table(ig))/(inertiedat_PEM(ig,1)**2))+ &
-                           ((layer_PEM(1)-ice_table(ig))/(ice_inertia**2))))
-           endif ! of if (iref.ne.0)        
-           ! lower layers of pure ice
-           do isoil=iref+2,nsoil_PEM
-             therm_i(ig,isoil)=ice_inertia
-           enddo
-         endif ! of if (iref.lt.(nsoilmx))
-       endif ! permanent glaciers
-
-       call soil_pem_1D(nsoil_PEM,.true.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM) 
-       call soil_pem_1D(nsoil_PEM,.false.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM) 
-
-       do it = 1,timelen 
-         tsoil_inst(ig,:,it) = tsoil_inst(ig,:,it) - (Tcol_saved(:) - tsoil_ave(ig,:))
-       enddo
-
-       error_depth = abs(icedepth_prev - ice_table(ig))      
-       endif
-
-     enddo
-
-     error_depth = 1.
-     countloop = 0
-   enddo
-
-#endif
- 
-      END SUBROUTINE ini_icetable
-      subroutine soil_pem_1D(nsoil,firstcall, &
-               therm_i,                          &
-               timestep,tsurf,tsoil,alph,beta)
-#ifndef CPP_STD
-      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
-                           mu_PEM,fluxgeo
-      use comsoil_h,only: volcapa
-#endif
-      implicit none
-
-!-----------------------------------------------------------------------
-!  Author: LL
-!  Purpose: Compute soil temperature using an implict 1st order scheme
-!  
-!  Note: depths of layers and mid-layers, soil thermal inertia and 
-!        heat capacity are commons in comsoil_PEM.h
-!        A convergence loop is added until equilibrium
-!-----------------------------------------------------------------------
-
-#include "dimensions.h"
-!#include "dimphys.h"
-
-!#include"comsoil.h"
-
-
-!-----------------------------------------------------------------------
-!  arguments
-!  ---------
-!  inputs:
-      integer,intent(in) :: nsoil	! number of soil layers  
-      logical,intent(in) :: firstcall ! identifier for initialization call 
-      real,intent(in) :: therm_i(nsoil) ! thermal inertia
-      real,intent(in) :: timestep	    ! time step
-      real,intent(in) :: tsurf   ! surface temperature
- 
-! outputs:
-      real,intent(inout) :: tsoil(nsoil) ! soil (mid-layer) temperature
-      real,intent(inout) :: alph(nsoil-1)
-      real,intent(inout) :: beta(nsoil-1)
-      
-! local variables:
-      integer ik
-      real :: thermdiff_PEM(nsoil-1)
-      real :: mthermdiff_PEM(0:nsoil-1)
-      real :: coefd_PEM(nsoil-1)
-      real :: coefq_PEM(0:nsoil-1)
-
-#ifndef CPP_STD
-
-! 0. Initialisations and preprocessing step
- if (firstcall) then
-    
-! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
-        do ik=0,nsoil-1
-	  mthermdiff_PEM(ik)=therm_i(ik+1)*therm_i(ik+1)/volcapa
-    
-	enddo
-
-! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
-        do ik=1,nsoil-1
-      thermdiff_PEM(ik)=((layer_PEM(ik)-mlayer_PEM(ik-1))*mthermdiff_PEM(ik) &
-                     +(mlayer_PEM(ik)-layer_PEM(ik))*mthermdiff_PEM(ik-1))  &
-                         /(mlayer_PEM(ik)-mlayer_PEM(ik-1))
-	enddo
-
-! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alpha_k and capcal
-      ! mu_PEM
-      mu_PEM=mlayer_PEM(0)/(mlayer_PEM(1)-mlayer_PEM(0))
-
-      ! q_{1/2}
-      coefq_PEM(0)=volcapa*layer_PEM(1)/timestep
-	! q_{k+1/2}
-        do ik=1,nsoil-1
-          coefq_PEM(ik)=volcapa*(layer_PEM(ik+1)-layer_PEM(ik)) 		 &
-                      /timestep
-	enddo
-
-	! d_k
-	do ik=1,nsoil-1
-	  coefd_PEM(ik)=thermdiff_PEM(ik)/(mlayer_PEM(ik)-mlayer_PEM(ik-1))
-	enddo
-	
-	! alph_PEM_{N-1}
-	alph(nsoil-1)=coefd_PEM(nsoil-1)/ 			 &
-                       (coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1))
-        ! alph_PEM_k
-        do ik=nsoil-2,1,-1
-	  alph(ik)=coefd_PEM(ik)/(coefq_PEM(ik)+coefd_PEM(ik+1)*    &
-                                   (1.-alph(ik+1))+coefd_PEM(ik))
-	enddo
-
-endif ! of if (firstcall)
-
-      IF (.not.firstcall) THEN
-!  2. Compute soil temperatures
-! First layer:
-        tsoil(1)=(tsurf+mu_PEM*beta(1)* &  
-                                  thermdiff_PEM(1)/mthermdiff_PEM(0))/ &
-                   (1.+mu_PEM*(1.0-alph(1))*&
-                    thermdiff_PEM(1)/mthermdiff_PEM(0))
-
-! Other layers:
-      do ik=1,nsoil-1
-          tsoil(ik+1)=alph(ik)*tsoil(ik)+beta(ik)
-       enddo
-      
-          ENDIF
-
-!  2. Compute beta_PEM coefficients (preprocessing for next time step)
-! Bottom layer, beta_PEM_{N-1}
-
-        beta(nsoil-1)=coefq_PEM(nsoil-1)*tsoil(nsoil)		 &
-                        /(coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1)) & 
-                 +  fluxgeo/(coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1))
-! Other layers
-      do ik=nsoil-2,1,-1
-	  beta(ik)=(coefq_PEM(ik)*tsoil(ik+1)+		 &
-                      coefd_PEM(ik+1)*beta(ik+1))/		 &
-                      (coefq_PEM(ik)+coefd_PEM(ik+1)*(1.0-alph(ik+1)) &
-                       +coefd_PEM(ik))
-      enddo
-#endif
-      end
-
-      END MODULE ini_soil_mod      
Index: trunk/LMDZ.COMMON/libf/evolution/interpolate_TIPEM_TIGCM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/interpolate_TIPEM_TIGCM.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/interpolate_TIPEM_TIGCM.F90	(revision 2855)
@@ -1,4 +1,12 @@
    subroutine interpolate_TIPEM_TIGCM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_PEM,TI_GCM)
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Transfer the thermal inertia from the PEM vertical  grid to the GCM vertical grid
+!!! 
+!!! 
+!!! Author: LL
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
@@ -10,16 +18,13 @@
 !  ---------
 !  inputs:
-      integer,intent(in) :: ngrid	! # of horizontal grid points
-      integer,intent(in) :: nslope	! # of subslope wihtin the mesh
-      integer,intent(in) :: nsoil_PEM	! # of soil layers in the PEM
-      integer,intent(in) :: nsoil_GCM	! # of soil layers in the GCM
-      real,intent(in) :: TI_PEM(ngrid,nsoil_PEM,nslope)	! # of soil layers in the PEM
-      real,intent(inout) :: TI_GCM(ngrid,nsoil_GCM,nslope)	! # of soil layers in the PEM
+      integer,intent(in) :: ngrid                                ! # of horizontal grid points
+      integer,intent(in) :: nslope                               ! # of subslope wihtin the mesh
+      integer,intent(in) :: nsoil_PEM                            ! # of soil layers in the PEM
+      integer,intent(in) :: nsoil_GCM                            ! # of soil layers in the GCM
+      real,intent(in) :: TI_PEM(ngrid,nsoil_PEM,nslope)          !  Thermal inertia in the PEM vertical grid [J/m^2/K/s^{1/2}]
+      real,intent(inout) :: TI_GCM(ngrid,nsoil_GCM,nslope)       ! Thermal inertia in the GCM vertical grid [J/m^2/K/s^{1/2}]
 
 !local variable
-      integer :: ig,islope,iloop
-
-
-
+      integer :: ig,islope,iloop                                 ! loop variables
      do ig = 1,ngrid
        do islope = 1,nslope
Index: trunk/LMDZ.COMMON/libf/evolution/iostart_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/iostart_PEM.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/iostart_PEM.F90	(revision 2855)
@@ -1,3 +1,11 @@
 MODULE iostart_PEM
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Read and write specific netcdf for the PEM 
+!!! 
+!!! 
+!!! Author: LL, inspired by iostart from the GCM
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
     IMPLICIT NONE
Index: trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90	(revision 2855)
@@ -1,11 +1,18 @@
 module lask_param_mod
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Define parameters from Laskar et al., 2004 evolution
+!!! 
+!!! Author: RV
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 implicit none
 
-  real,save,allocatable :: yearlask(:)
-  real,save,allocatable :: oblask(:)
-  real,save,allocatable :: exlask(:)
-  real,save,allocatable :: lsplask(:)
-  integer, save :: last_ilask
+  real,save,allocatable :: yearlask(:)          ! year before present from Laskar et al. Tab
+  real,save,allocatable :: oblask(:)            ! obliquity [deg]
+  real,save,allocatable :: exlask(:)            ! excentricity[deg]
+  real,save,allocatable :: lsplask(:)           ! ls perihelie [deg]
+  integer, save :: last_ilask                   ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
 
 contains
Index: trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90	(revision 2855)
@@ -12,6 +12,7 @@
 !=======================================================================
 !
-! Read initial confitions file
+! Purpose: Read in the data_GCM_Yr*.nc the number of time step
 !
+! Author: RV
 !=======================================================================
 
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 2855)
@@ -87,6 +87,5 @@
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
-      USE soil_evolution_mod, ONLY: soil_pem, soil_pem_CN
-      use comsoil_h_PEM, only: ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
+      use comsoil_h_PEM, only: soil_pem,ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
                               TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM,         & ! soil thermal inertia          
                               tsoil_PEM, mlayer_PEM,layer_PEM,                  & !number of subsurface layers, soil mid layer depths
@@ -175,12 +174,12 @@
       REAL , dimension(:,:), allocatable :: min_co2_ice_s_2        ! LON x LAT field : minimum of water ice at each point for the second year
 
-      REAL :: global_ave_press_GCM
-      REAL :: global_ave_press_old           ! physical point field : Global average pressure of initial/previous time step
-      REAL :: global_ave_press_new           ! physical point field : Global average pressure of current time step
-
-      REAL , dimension(:,:), allocatable ::  zplev_new
-      REAL , dimension(:,:), allocatable ::  zplev_gcm
-      REAL , dimension(:,:,:), allocatable ::  zplev_new_timeseries  ! same but with the time series
-      REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries   ! same but with the time series
+      REAL :: global_ave_press_GCM           ! constant: global average pressure retrieved in the GCM [Pa]
+      REAL :: global_ave_press_old           ! constant: Global average pressure of initial/previous time step
+      REAL :: global_ave_press_new           ! constant: Global average pressure of current time step
+
+      REAL , dimension(:,:), allocatable ::  zplev_new ! Physical x Atmospheric field : mass of the atmospheric  layers in the pem at current time step [kg/m^2]
+      REAL , dimension(:,:), allocatable ::  zplev_gcm ! same but retrieved from the gcm [kg/m^2]
+      REAL , dimension(:,:,:), allocatable ::  zplev_new_timeseries  ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
+      REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries   ! same but with the time series, for oldest time step
 
       LOGICAL :: STOPPING_water   ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
@@ -189,48 +188,43 @@
       LOGICAL :: STOPPING_1_co2   ! Logical : is there still water ice to sublimate?
 
-      REAL, dimension(:,:,:),allocatable  :: q_co2_GCM ! Initial amount of co2 in the first layer
-
-      real,save :: m_co2, m_noco2, A , B, mmean
-      real ,allocatable :: vmr_co2_gcm_phys(:,:) !(ngrid) ! co2 volume mixing ratio
-      real ,allocatable :: vmr_co2_pem_phys(:,:) !(ngrid) ! co2 volume mixing ratio
-      real ,allocatable :: q_h2o_GCM_phys(:,:)
-      real ,allocatable :: q_co2_GCM_phys(:,:)
-      real ,allocatable :: q_co2_PEM_phys(:,:)
-      REAL, ALLOCATABLE ::  ps_GCM(:,:,:)
-      REAL, ALLOCATABLE :: ps_GCM_yr1(:,:,:)
-      REAL, ALLOCATABLE ::  vmr_co2_gcm(:,:,:)
-      REAL, ALLOCATABLE ::  q_h2o_GCM(:,:,:)
-      REAL ,allocatable ::  q_h2o_PEM_phys(:,:)
-      integer :: timelen
-      REAL :: ave
-
-      REAL, ALLOCATABLE :: p(:,:)  !(ngrid,llmp1)
-      REAL :: extra_mass           ! Extra mass of a tracer if it is greater than 1
+      REAL, dimension(:,:,:),allocatable  :: q_co2_GCM ! Lon x Lat x Time : mass mixing ratio of co2 in the first layer [kg/kg]
+
+      real,save :: m_co2, m_noco2, A , B, mmean        ! Molar mass of co2, no co2 (Ar, ...), intermediate A, B for computations, mean molar mass of the layer [mol/kg]
+      real ,allocatable :: vmr_co2_gcm_phys(:,:) ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
+      real ,allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times  co2 volume mixing ratio used in the PEM
+      real ,allocatable :: q_h2o_GCM_phys(:,:)   ! Physics x Times h2o mass mixing ratio in the first layer from the GCM  [kg/kg]
+      real ,allocatable :: q_co2_GCM_phys(:,:)   ! Physics x Times co2 mass mixing ratio in the first layer from the GCM  [kg/kg]
+      real ,allocatable :: q_co2_PEM_phys(:,:)   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM  [kg/kg]
+      REAL, ALLOCATABLE ::  ps_GCM(:,:,:)        ! Lon x Lat x Times: surface pressure from the GCM [Pa]
+      REAL, ALLOCATABLE :: ps_GCM_yr1(:,:,:)     ! Lon x Lat x Times: surface pressure from the 1st year of the GCM [Pa]
+      REAL, ALLOCATABLE ::  vmr_co2_gcm(:,:,:)   ! Lon x Lat x Times: co2 volumemixing ratio retrieve from the gcm [m^3/m^3]
+      REAL, ALLOCATABLE ::  q_h2o_GCM(:,:,:)     ! Lon x Lat x Times: h2o volume mixing ratio retrieved from the GCM
+      REAL ,allocatable ::  q_h2o_PEM_phys(:,:)  ! Physics x Times: h2o mass mixing ratio computed in the PEM
+      integer :: timelen                         ! # time samples
+      REAL :: ave                                ! intermediate varibale to compute average
+
+      REAL, ALLOCATABLE :: p(:,:)  ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
+      REAL :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
+     REAL :: beta_clap_co2 = 3182.48 ! clapeyron's law for CO2
+     REAL :: alpha_clap_co2 = 23.3494 ! Clapeyron's law for CO2
+
 
 !!!!!!!!!!!!!!!!!!!!!!!! SLOPE
-      character*2 :: str2
-      REAL ,allocatable :: watercap_slope(:,:)                           !(ngrid,nslope)
-      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1        ! LON x LAT field : minimum of water ice at each point for the first year
-      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2        ! LON x LAT field : minimum of water ice at each point for the second year
-      REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_1        ! LON x LAT field : minimum of water ice at each point for the first year
-      REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_2        ! LON x LAT field : minimum of water ice at each point for the second year
-
-
-      REAL , dimension(:,:,:,:), allocatable :: co2_ice_GCM_slope        ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM
-      REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope     ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM
-
-
+      REAL ,allocatable :: watercap_slope(:,:)                           ! Physics x Nslope: watercap per slope 
+      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1        ! LON x LAT field : minimum of co2 ice at each point for the first year [kg/m^2]
+      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2        ! LON x LAT field : minimum of co2 ice at each point for the second year [kg/m^2]
+      REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_1        ! LON x LAT field : minimum of water ice at each point for the first year [kg/m^2]
+      REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_2        ! LON x LAT field : minimum of water ice at each point for the second year [kg/m^2]
+      REAL , dimension(:,:,:,:), allocatable :: co2_ice_GCM_slope        ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM [kg/m^2]
+      REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope     ! Physics x NSLOPE x Times field : co2 ice given by the GCM  [kg/m^2]
       REAL, dimension(:,:),allocatable  :: initial_co2_ice_sublim_slope    ! physical point field : Logical array indicating sublimating point of co2 ice
-      REAL, dimension(:,:),allocatable  :: initial_h2o_ice_slope           ! physical point field : Logical array indicating sublimating point of h2o ice
-      REAL, dimension(:,:),allocatable  :: initial_co2_ice_slope           ! physical point field : Logical array indicating sublimating point of co2 ice
-
-      REAL , dimension(:,:,:), allocatable :: tendencies_co2_ice_slope     ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year
-      REAL , dimension(:,:,:), allocatable :: tendencies_h2o_ice_slope     ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year
-      REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope   ! physical point field : Tendency of evolution of perenial co2 ice over a year
-      REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope_ini ! physical point field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
-      REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice_phys_slope   ! physical point field : Tendency of evolution of perenial co2 ice over a year
-
+      REAL, dimension(:,:),allocatable  :: initial_h2o_ice_slope           ! physical point field : Logical array indicating if there is water ice at initial state
+      REAL, dimension(:,:),allocatable  :: initial_co2_ice_slope           ! physical point field : Logical array indicating if there is co2 ice at initial state
+      REAL , dimension(:,:,:), allocatable :: tendencies_co2_ice_slope     ! LON x LAT x nslope field : Tendency of evolution of perenial co2 ice over a year
+      REAL , dimension(:,:,:), allocatable :: tendencies_h2o_ice_slope     ! LON x LAT x slope field : Tendency of evolution of perenial water ice over a year
+      REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope   ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
+      REAL, dimension(:,:),allocatable  :: tendencies_co2_ice_phys_slope_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
+      REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice_phys_slope   ! physical pointx slope  field : Tendency of evolution of perenial co2 ice over a year
       REAL,SAVE,ALLOCATABLE,DIMENSION(:) ::  co2_hmax                      ! Maximum height  for CO2 deposit on slopes (m)
-
       REAL, PARAMETER :: rho_co2 = 1600           ! CO2 ice density (kg/m^3)
       INTEGER :: iaval                            ! Index of the neighboord slope ()
@@ -241,53 +235,51 @@
 
      REAL, ALLOCATABLE :: tsurf_ave(:,:,:)          ! LON x LAT x SLOPE field : Averaged Surface Temperature [K]
-     REAL, ALLOCATABLE  :: tsurf_ave_phys(:,:)      ! IG x LAT x SLOPE field : Averaged Surface Temperature [K]
+     REAL, ALLOCATABLE  :: tsurf_ave_phys(:,:)      ! Physic x LAT x SLOPE field : Averaged Surface Temperature [K]
      REAL, ALLOCATABLE :: tsoil_ave(:,:,:,:)        ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
-     REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:)    ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
-     REAL, ALLOCATABLE :: tsoil_ave_phys_yr1(:,:,:) !IG x SLOPE field : Averaged Soil Temperature [K]
-
+     REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:)    ! LON x LAT x SLOPE field : Averaged Soil Temperature  during 1st year of the GCM [K]
+     REAL, ALLOCATABLE :: tsoil_ave_phys_yr1(:,:,:) ! Physics x SLOPE field : Averaged Soil Temperature  during 1st year [K]
      REAL, ALLOCATABLE :: TI_GCM(:,:,:,:)                  ! LON x LAT x SLOPE field : Averaged Thermal Inertia  [SI]
-     REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:)    !LONX LAT x SLOPE XTULES field : NOn averaged Surf Temperature [K]
-     REAL, ALLOCATABLE :: tsurf_phys_GCM_timeseries(:,:,:) !IG x SLOPE XTULES field : NOn averaged Surf Temperature [K]
+     REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:)    ! LON X LAT x SLOPE XTULES field : Surface Temperature in timeseries [K]
+     REAL, ALLOCATABLE :: tsurf_phys_GCM_timeseries(:,:,:) ! Physic x SLOPE XTULES field : NOn averaged Surf Temperature [K]
 
      REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
      REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:,:)    !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
-
      REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:,:)     ! LON x LAT x SLOPE field : Averaged Surface Temperature of the first year of the gcm [K]
-     REAL, ALLOCATABLE  :: tsurf_ave_phys_yr1(:,:) ! IG x LAT x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
-
-     REAL, PARAMETER :: year_step = 1
-     INTEGER :: year_iter        !Counter for the number of PEM iteration
-     INTEGER :: year_iter_max
-     REAL :: timestep
-
-     REAL, ALLOCATABLE :: inertiesoil(:,:)    ! Thermal inertia of the mesh for restart
-
-     REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:)  ! Averaged GCM Thermal Inertia  [SI]
-     REAL, ALLOCATABLE :: TI_GCM_start(:,:,:) ! Averaged GCM Thermal Inertia  [SI]
-
-     REAL,ALLOCATABLE  :: ice_depth(:,:)
-     REAL,ALLOCATABLE  :: TI_locslope(:,:)
-     REAL,ALLOCATABLE  :: Tsoil_locslope(:,:)
-     REAL,ALLOCATABLE  :: Tsurf_locslope(:)
-     REAL,ALLOCATABLE  :: alph_locslope(:,:)
-     REAL,ALLOCATABLE  :: beta_locslope(:,:)   
-     REAL,ALLOCATABLE  :: watersurf_density_timeseries(:,:,:,:)
-     REAL,ALLOCATABLE  :: watersoil_density_timeseries(:,:,:,:,:)
-     REAL,ALLOCATABLE  :: watersurf_density_phys_timeseries(:,:,:)
-     REAL,ALLOCATABLE  :: watersurf_density_phys_ave(:,:)
-     REAL,ALLOCATABLE  :: watersoil_density_phys_PEM_timeseries(:,:,:,:)
-     REAL,ALLOCATABLE  :: watersoil_density_phys_PEM_ave(:,:,:)
-     REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)
-     REAL, ALLOCATABLE :: delta_co2_adsorbded(:)
-     REAL :: totmass_adsorbded
-   real :: alpha_clap_h2o = -6143.7
-   real :: beta_clap_h2o = 28.9074
-
+     REAL, ALLOCATABLE  :: tsurf_ave_phys_yr1(:,:) ! Physic SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
+     REAL, ALLOCATABLE :: inertiesoil(:,:)    !Physic x Depth  Thermal inertia of the mesh for restart [SI]
+
+     REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:)  ! Physic x Depth x Slope Averaged GCM Thermal Inertia per slope  [SI]
+     REAL, ALLOCATABLE :: TI_GCM_start(:,:,:) ! Same but for the start
+
+     REAL,ALLOCATABLE  :: ice_depth(:,:)      ! Physic x SLope: Ice table depth [m]
+     REAL,ALLOCATABLE  :: TI_locslope(:,:)    ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
+     REAL,ALLOCATABLE  :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
+     REAL,ALLOCATABLE  :: Tsurf_locslope(:)   ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
+     REAL,ALLOCATABLE  :: alph_locslope(:,:)  ! Physic x Soil: Intermediate to compute Tsoil [1]
+     REAL,ALLOCATABLE  :: beta_locslope(:,:)  ! Physic x Soil : Intermediate tocompute Tsoil [K]
+     REAL,ALLOCATABLE  :: watersurf_density_timeseries(:,:,:,:) ! Lon x Lat x Slope x Times: water surface density, time series [kg/m^3]
+     REAL,ALLOCATABLE  :: watersoil_density_timeseries(:,:,:,:,:) ! Lon x Lat x Soil x Slope x Times water soil density, time series [kg /m^3]
+     REAL,ALLOCATABLE  :: watersurf_density_phys_timeseries(:,:,:) ! Physic  x Slope x Times, water surface density, time series [kg/m^3]
+     REAL,ALLOCATABLE  :: watersurf_density_phys_ave(:,:)          ! Physic  x Slope, water surface density, yearly averaged [kg/m^3]
+     REAL,ALLOCATABLE  :: watersoil_density_phys_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
+     REAL,ALLOCATABLE  :: watersoil_density_phys_PEM_ave(:,:,:)          ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
+     REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)                     ! Surface temperature saved from previous time step [K]
+     REAL, ALLOCATABLE :: delta_co2_adsorbded(:)                         ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
+     REAL :: totmass_adsorbded                                           ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
+     REAL :: alpha_clap_h2o = -6143.7                                    ! coeffcient to compute psat, from Murphie et Kood 2005 [K]
+     REAL :: beta_clap_h2o = 28.9074                                     ! coefficient to compute psat, from Murphie et Kood 2005 [1]
+     LOGICAL :: bool_sublim                                              ! logical to check if there is sublimation or not
+
+!! Some parameters for the PEM run
+    REAL, PARAMETER :: year_step = 1     ! timestep for the pem 
+    INTEGER ::  year_iter                ! number of iteration
+    INTEGER :: year_iter_max             ! maximum number of iterations before stopping
+    REAL :: timestep                     ! timestep [s] 
 #ifdef CPP_STD
 !     INTEGER :: nsplit_phys=1
 !     LOGICAL :: iflag_phys=.true.
-     REAL :: frost_albedo_threshold=0.05
-     REAL :: albedo_h2o_frost
-     REAL,ALLOCATABLE :: co2ice(:)
+     REAL :: frost_albedo_threshold=0.05             ! frost albedo threeshold to convert fresh frost to old ice
+     REAL :: albedo_h2o_frost                        ! albedo of h2o frost
+     REAL,ALLOCATABLE :: co2ice(:)                   ! Physics: co2 ice mesh averaged [kg/m^2]
 #endif
 
@@ -295,9 +287,5 @@
 
 ! Loop variable
-     LOGICAL :: bool_sublim
-     INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil
-     REAL :: beta = 3182.48
-     REAL :: alpha = 23.3494
-  
+     INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil 
 #ifndef CPP_STD
 ! Parallel variables
@@ -566,5 +554,5 @@
      print *, "Downloading data Y1..."
 
-     call read_data_GCM("data_GCM_Y1.nc", min_h2o_ice_s_1,min_co2_ice_s_1,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_yr1,timelen,min_co2_ice_slope_1,min_h2o_ice_slope_1,&   
+     call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm, min_h2o_ice_s_1,min_co2_ice_s_1,vmr_co2_gcm,ps_GCM_yr1,min_co2_ice_slope_1,min_h2o_ice_slope_1,&   
                        nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &      
                        watersurf_density_timeseries,watersoil_density_timeseries)
@@ -581,5 +569,5 @@
      print *, "Downloading data Y2"
 
-     call read_data_GCM("data_GCM_Y2.nc", min_h2o_ice_s_2,min_co2_ice_s_2,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM,timelen,min_co2_ice_slope_2,min_h2o_ice_slope_2, &
+     call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm ,min_h2o_ice_s_2,min_co2_ice_s_2,vmr_co2_gcm,ps_GCM,min_co2_ice_slope_2,min_h2o_ice_slope_2, &
                   nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &      
                   watersurf_density_timeseries,watersoil_density_timeseries)
@@ -815,6 +803,6 @@
 !---------------------------- Read the PEMstart --------------------- 
 
-      call pemetat0(.true.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
-      tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsurf_phys_GCM_timeseries,tsoil_phys_PEM_timeseries,&
+      call pemetat0(.false.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
+      tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsoil_phys_PEM_timeseries,&
       tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM, co2_adsorbded_phys,delta_co2_adsorbded,&
       watersurf_density_phys_ave,watersoil_density_phys_PEM_ave)
@@ -1071,5 +1059,5 @@
             do t=1,timelen
               if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-3) then
-                ave = ave + beta/(alpha-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
+                ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
               else
                 ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2855)
@@ -1,10 +1,9 @@
    SUBROUTINE pemetat0(startpem_file,filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
-                      tsurf_ave_yr1,tsurf_ave,q_co2,q_h2o,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
+                      tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
                       global_ave_pressure, m_co2_regolith_phys,deltam_co2_regolith_phys, watersurf_ave,watersoil_ave)
    
    use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
-   use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
+   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
    use comsoil_h, only:  volcapa,inertiedat 
-   use soil_evolution_mod, only: soil_pem,soil_pem_CN
    use adsorption_mod, only : regolith_co2adsorption
    USE temps_mod_evol, ONLY: year_PEM
@@ -12,55 +11,52 @@
   implicit none
  
-  character(len=*), intent(in) :: filename
-  LOGICAL,intent(in) :: startpem_file
-  character*8 :: fichnom
-  integer,intent(in) :: ngrid
-  integer,intent(in) :: nsoil_GCM
-  integer,intent(in) :: nsoil_PEM
-  integer,intent(in) :: nslope
-  integer,intent(in) :: timelen
-  real, intent(in) :: tsurf_ave_yr1(ngrid,nslope)                          ! surface temperature at the first year of GCM call
-  real,intent(in) :: tsurf_ave(ngrid,nslope)                  ! surface temperature at the current year
+  character(len=*), intent(in) :: filename                    ! name of the startfi_PEM.nc
+  LOGICAL,intent(in) :: startpem_file                         ! boolean to check if we read the startfile or not
+  integer,intent(in) :: ngrid                                 ! # of physical grid points
+  integer,intent(in) :: nsoil_GCM                             ! # of vertical grid points in the GCM
+  integer,intent(in) :: nsoil_PEM                             ! # of vertical grid points in the PEM
+  integer,intent(in) :: nslope                                ! # of sub-grid slopes
+  integer,intent(in) :: timelen                               ! # time samples
+  real, intent(in) :: tsurf_ave_yr1(ngrid,nslope)             ! surface temperature at the first year of GCM call [K]
+  real,intent(in) :: tsurf_ave_yr2(ngrid,nslope)              ! surface temperature at the second  year of GCM call [K]
   real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
   real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
-  real,intent(in) :: ps_inst(ngrid,timelen)                        ! surface pressure [Pa]
-  real,intent(in) :: tsurf_inst(ngrid,nslope,timelen) ! soil (mid-layer) temperature
-  real,intent(in) :: timestep	    ! time step
-  real,intent(in) :: tend_h2oglaciers(ngrid,nslope)
-  real,intent(in) :: tend_co2glaciers(ngrid,nslope)
-  real,intent(in) :: co2ice(ngrid,nslope)
-  real,intent(in) :: waterice(ngrid,nslope) 
-  real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope) 
-  real, intent(in) :: watersurf_ave(ngrid,nslope)     ! surface water ice density, yearly averaged  (kg/m^3)
-  real, intent(inout) :: watersoil_ave(ngrid,nsoil_PEM,nslope)    ! surface water ice density, yearly averaged (kg/m^3)
-  real, intent(in) :: global_ave_pressure                       ! mean average pressure on the planet
+  real,intent(in) :: ps_inst(ngrid,timelen)                   ! surface pressure [Pa]
+  real,intent(in) :: timestep                                 ! time step [s]
+  real,intent(in) :: tend_h2oglaciers(ngrid,nslope)           ! tendencies on h2o glaciers 
+  real,intent(in) :: tend_co2glaciers(ngrid,nslope)           ! tendencies on co2 glaciers
+  real,intent(in) :: co2ice(ngrid,nslope)                     ! co2 ice amount [kg/m^2]
+  real,intent(in) :: waterice(ngrid,nslope)                   ! water ice amount [kg/m^2]
+  real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope)   ! soil temperature during 1st year [K]
+  real, intent(in) :: watersurf_ave(ngrid,nslope)             ! surface water ice density, yearly averaged  (kg/m^3)
+  real, intent(inout) :: watersoil_ave(ngrid,nsoil_PEM,nslope)! surface water ice density, yearly averaged (kg/m^3)
+  real, intent(in) :: global_ave_pressure                     ! mean average pressure on the planet [Pa]
 ! outputs
 
-  real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) thermal inertia (SI)
-  real,intent(inout) :: tsoil_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) temperature (K)
-  real,intent(inout) :: ice_table(ngrid,nslope) ! (m)
-  real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,nslope,timelen) ! instantaneous soil (mid-layer) temperature (k)
-  real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of co2 adsorbed (kg/m^2)
-  real,intent(out) :: deltam_co2_regolith_phys(ngrid)                 ! mass of co2 that is exchanged due to adsorption desorption 
+  real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope)                ! soil (mid-layer) thermal inertia in the PEM grid [SI]
+  real,intent(inout) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)             ! soil (mid-layer) temperature [K]
+  real,intent(inout) :: ice_table(ngrid,nslope)                       ! Ice table depth [m]
+  real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,nslope,timelen)    ! instantaneous soil (mid-layer) temperature [K]
+  real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of co2 adsorbed [kg/m^2]
+  real,intent(out) :: deltam_co2_regolith_phys(ngrid)                 ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
 ! local
-   real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)
-   real :: TI_startPEM(ngrid,nsoil_PEM,nslope)
-   LOGICAL :: found
-   integer :: iloop,ig,islope,it,isoil
-   REAL :: TI_breccia = 750.
-   REAL :: TI_bedrock = 2300.
-   real :: kcond
-   real :: delta
-   CHARACTER*2 :: num
-   real :: tsoil_saved(ngrid,nsoil_PEM)
-   real :: tsoil_tmp_yr1(ngrid,nsoil_PEM,nslope)
-   real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope)
-   real :: tsoil_tmp(ngrid,nsoil_PEM,nslope)
-   real :: alph_tmp(ngrid,nsoil_PEM-1)
-   real :: beta_tmp(ngrid,nsoil_PEM-1)
-   real :: co2_ads_prev(ngrid)
-   real :: year_PEM_read
-   real :: alpha_clap_h2o = -6143.7
-   real :: beta_clap_h2o = 28.9074
+   real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)                     ! soil temperature saved in the start [K]
+   real :: TI_startPEM(ngrid,nsoil_PEM,nslope)                        ! soil thermal inertia saved in the start [SI]
+   LOGICAL :: found                                                   ! check if variables are found in the start
+   integer :: iloop,ig,islope,it,isoil                                ! index for loops
+   REAL :: TI_breccia = 750.                                          ! Thermal inertia of Breccia following Wood 2009 [SI]
+   REAL :: TI_bedrock = 2300.                                         ! Thermal inertia of Bedrock following Wood 2009 [SI]
+   real :: kcond                                                      ! Thermal conductivity, intermediate variable [SI]
+   real :: delta                                                      ! Depth of the interface regolith-breccia, breccia -bedrock [m]
+   CHARACTER*2 :: num                                                 ! intermediate string to read PEM start sloped variables
+   real :: tsoil_saved(ngrid,nsoil_PEM)                               ! saved soil temperature [K]
+   real :: tsoil_tmp_yr1(ngrid,nsoil_PEM,nslope)                      ! intermediate soil temperature during yr1[K]
+   real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope)                      ! intermediate soil temperature during yr 2 [K]
+   real :: alph_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computation []
+   real :: beta_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computatio []                               
+   real :: year_PEM_read                                              ! Year of the PEM previous run
+   real :: alpha_clap_h2o = -6143.7                                   ! Intermediate coefficient to compute psat using clapeyron law, Murphie et al. 2005 [K^-1]
+   real :: beta_clap_h2o = 28.9074                                    ! Intermediate coefficient to compute psat using clapeyron law, Murphie et al. 2005 [1]
+
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -68,8 +64,9 @@
 !!! Purpose: read start_pem. Need a specific iostart_PEM
 !!!
-!!! Order: 1. Thermal Inertia 
-!!!        2. Soil Temperature
-!!!        3. Ice table
-!!!        4. Mass of CO2 adsorbed
+!!! Order: 1. Previous year of the PEM run
+!!!        2. Thermal Inertia 
+!!!        3. Soil Temperature
+!!!        4. Ice table
+!!!        5. Mass of CO2 adsorbed
 !!!
 !!! /!\ This order must be respected !
@@ -212,5 +209,5 @@
     call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
     tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
-    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
+    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
 
 
@@ -419,13 +416,7 @@
      if(soil_pem) then
 
-    DO ig = 1,ngrid
-      DO islope = 1,nslope
-
-     write(*,*) 'ig,islope ,ice table=',ig,islope,ice_table(ig,islope)
-
-      ENDDO
-    ENDDO
-
-!! small test
+   
+
+!! Sanity check
     DO ig = 1,ngrid
       DO islope = 1,nslope
@@ -441,6 +432,5 @@
      endif!soil_pem
 
-      write(*,*) "construction ok, no nan" 
-     
+
 
    END SUBROUTINE
Index: trunk/LMDZ.COMMON/libf/evolution/pemredem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 2855)
@@ -1,3 +1,11 @@
 module pemredem
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Write specific netcdf restart for the PEM 
+!!! 
+!!! 
+!!! Author: LL, inspired by phyredem  from the GCM
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 implicit none
@@ -8,5 +16,5 @@
                          day_ini,time,nslope,def_slope,subslope_dist)
 ! create physics restart file and write time-independent variables
-  use comsoil_h_PEM, only: mlayer_PEM,fluxgeo,inertiedat_PEM
+  use comsoil_h_PEM, only: soil_pem,mlayer_PEM,fluxgeo,inertiedat_PEM
   use iostart_PEM, only : open_restartphy, close_restartphy, & 
                       put_var, put_field, length
@@ -69,7 +77,6 @@
                       put_var, put_field
 
-  use comsoil_h_PEM, only: mlayer_PEM,fluxgeo,inertiedat_PEM
+  use comsoil_h_PEM, only: mlayer_PEM,fluxgeo,inertiedat_PEM,soil_pem
   USE temps_mod_evol, ONLY: year_PEM
-  use soil_evolution_mod, only: soil_pem
                 
   implicit none
Index: trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2855)
@@ -2,5 +2,5 @@
 ! $Id $
 !
-SUBROUTINE read_data_GCM(fichnom,min_h2o_ice_s,min_co2_ice_s,iim_input,jjm_input,nlayer,vmr_co2_gcm,ps_GCM,timelen, & 
+SUBROUTINE read_data_GCM(fichnom,timelen, iim_input,jjm_input,min_h2o_ice_s,min_co2_ice_s,vmr_co2_gcm,ps_GCM, & 
              min_co2_ice_slope,min_h2o_ice_slope,nslope,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2_GCM,q_h2o_GCM,co2_ice_slope, &
              watersurf_density,watersoil_density)
@@ -10,5 +10,5 @@
                         nf90_inquire_dimension,nf90_close
       use comsoil_h, only: nsoilmx
-      USE soil_evolution_mod, ONLY: soil_pem
+      USE comsoil_h_PEM, ONLY: soil_pem
 
       IMPLICIT NONE
@@ -16,6 +16,7 @@
 !=======================================================================
 !
-! Read initial confitions file
+! Purpose: Read initial confitions file from the GCM
 !
+! Authors: RV & LL
 !=======================================================================
 
@@ -24,40 +25,34 @@
 !===============================================================================
 ! Arguments:
+
   CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
   INTEGER, INTENT(IN) :: timelen                   ! number of times stored in the file
-
-  INTEGER :: iim_input,jjm_input,nlayer,nslope
-  REAL, ALLOCATABLE ::  h2o_ice_s(:,:,:)                       ! h2o_ice_s of the concatenated file
-  REAL, ALLOCATABLE ::  co2_ice_s(:,:,:)                       ! co2_ice_s of the concatenated file
-
-  REAL, ALLOCATABLE ::  h2o_ice_s_slope(:,:,:,:)                       ! co2_ice_s of the concatenated file
-
-  REAL, INTENT(OUT) ::  min_h2o_ice_s(iim_input+1,jjm_input+1) ! Minimum of h2o_ice_s of the year
-  REAL, INTENT(OUT) ::  min_co2_ice_s(iim_input+1,jjm_input+1) ! Minimum of co2_ice_s of the year
-  REAL, INTENT(OUT) ::  min_co2_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
-  REAL, INTENT(OUT) ::  min_h2o_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
-  REAL, INTENT(OUT) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,timelen)      !!!!vmr_co2_phys_gcm(iim_input+1,jjm_input+1,timelen)
-  REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,timelen)
-  REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,timelen)
-  REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,timelen)
+  INTEGER :: iim_input,jjm_input,nslope            ! number of points in the lat x lon dynamical grid, number of subgrid slopes
+
+! Ouputs
+  REAL, INTENT(OUT) ::  min_h2o_ice_s(iim_input+1,jjm_input+1) ! Minimum of h2o ice, mesh averaged of the year  [kg/m^2]
+  REAL, INTENT(OUT) ::  min_co2_ice_s(iim_input+1,jjm_input+1) ! Minimum of co2 ice, mesh averaged  of the year [kg/m^2]
+  REAL, INTENT(OUT) ::  min_co2_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2 ice  per slope of the year [kg/m^2]
+  REAL, INTENT(OUT) ::  min_h2o_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of h2o ice per slope of the year [kg/m^2]
+  REAL, INTENT(OUT) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,timelen)      ! CO2 volume mixing ratio in the first layer  [mol/m^3]
+  REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,timelen)        ! H2O mass mixing ratio in the first layer [kg/m^3]
+  REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,timelen)        ! CO2 mass mixing ratio in the first layer [kg/m^3]
+  REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,timelen)          ! Surface Pressure [Pa]
+  REAL, INTENT(OUT) ::  co2_ice_slope(iim_input+1,jjm_input+1,nslope,timelen) ! co2 ice amount per  slope of the year [kg/m^2]
 
 !SOIL
-  REAL, INTENT(OUT) ::  tsurf_ave(iim_input+1,jjm_input+1,nslope) ! Average surface temperature of the concatenated file
-  REAL, INTENT(OUT) ::  tsoil_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average soil temperature of the concatenated file
-
-  REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen) ! Surface temperature of the concatenated file
-  REAL , INTENT(OUT) ::  tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file
-  REAL , INTENT(OUT) ::  watersurf_density(iim_input+1,jjm_input+1,nslope,timelen) ! Soil temperature of the concatenated file
-  REAL , INTENT(OUT) ::  watersoil_density(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file
-
-  REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Thermal Inertia  of the concatenated file
-  REAL, INTENT(OUT) ::  TI_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average Thermal Inertia  of the concatenated file
-  REAL, INTENT(OUT) ::  co2_ice_slope(iim_input+1,jjm_input+1,nslope,timelen) ! Minimum of co2_ice slope of the year
+  REAL, INTENT(OUT) ::  tsurf_ave(iim_input+1,jjm_input+1,nslope)         ! Average surface temperature of the concatenated file [K]
+  REAL, INTENT(OUT) ::  tsoil_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average soil temperature of the concatenated file [K]
+  REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen)                  ! Surface temperature of the concatenated file, time series [K]
+  REAL , INTENT(OUT) ::  tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen)         ! Soil temperature of the concatenated file, time series [K]
+  REAL , INTENT(OUT) ::  watersurf_density(iim_input+1,jjm_input+1,nslope,timelen)         ! Water density at the surface, time series [kg/m^3]
+  REAL , INTENT(OUT) ::  watersoil_density(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Water density in the soil layer, time series [kg/m^3]
+  REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen)                          ! Thermal Inertia  of the concatenated file, times series [SI]
+  REAL, INTENT(OUT) ::  TI_ave(iim_input+1,jjm_input+1,nsoilmx,nslope)                     ! Average Thermal Inertia  of the concatenated file [SI]
 !===============================================================================
 !   Local Variables 
-  CHARACTER(LEN=256) :: msg, var, modname
-  INTEGER,PARAMETER :: length=100
-  INTEGER :: iq, fID, vID, idecal
-  INTEGER :: ierr
+  CHARACTER(LEN=256) :: msg, var, modname               ! for reading
+  INTEGER :: iq, fID, vID, idecal                       ! for reading
+  INTEGER :: ierr                                       ! for reading
   CHARACTER(len=12) :: start_file_type="earth" ! default start file type
 
@@ -66,12 +61,15 @@
 
   INTEGER :: edges(4),corner(4)
-  INTEGER :: i,j,t
-  real,save :: m_co2, m_noco2, A , B, mmean
-
-  INTEGER :: islope
-  CHARACTER*2 :: num
+  INTEGER :: i,j,t                                                     ! loop variables
+  real,save :: m_co2, m_noco2, A , B, mmean                            ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
+
+  INTEGER :: islope                                                    ! loop for variables
+  CHARACTER*2 :: num                                                   ! for reading sloped variables
+  REAL, ALLOCATABLE ::  h2o_ice_s(:,:,:)                               ! h2o ice, mesh averaged, of the concatenated file [kg/m^2]
+  REAL, ALLOCATABLE ::  co2_ice_s(:,:,:)                               ! co2 ice, mesh averaged, of the concatenated file [kg/m^2]
+  REAL, ALLOCATABLE ::  h2o_ice_s_slope(:,:,:,:)                       ! h2o ice per slope of the concatenated file [kg/m^2]
 
 !-----------------------------------------------------------------------
-  modname="pemetat0"
+  modname="read_data_gcm"
 
       m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
Index: trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 2855)
@@ -1,4 +1,10 @@
       MODULE recomp_orb_param_mod
-
+!=======================================================================
+!
+! Purpose: Recompute orbit parameters based on values by Laskar et al.,
+! 2004
+!
+! Authors: RV 
+!=======================================================================
       IMPLICIT NONE
 
Index: trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90	(revision 2855)
@@ -1,3 +1,13 @@
 program reshape_XIOS_output
+
+!=======================================================================
+!
+! Purpose: Read XIOS files, and convert them into the correct GCM grid
+!          XIOS  longitudes start at -180 but stop before -180 (not duplicated)
+!          We basically add the last point, and complete the XIOS file. Looped
+!          over the two GCM runs
+!
+! Authors: RV & LL
+!=======================================================================
     use netcdf
     implicit none
Index: trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM.F90	(revision 2855)
@@ -33,7 +33,7 @@
       INTEGER, INTENT(IN) :: nsoil                  ! Number of soil layers 
       REAL :: icedepth                  ! Ice cover thickness (m)
-      REAL :: inert_h2o_ice = 800.
-      REAL :: rho_ice = 920.
-      REAL :: prev_thermi(ngrid,nsoil)
+      REAL :: inert_h2o_ice = 800.      ! surface water ice thermal inertia [SI]
+      REAL :: rho_ice = 920.            ! density of water ice [kg/m^3]
+      REAL :: prev_thermi(ngrid,nsoil)  ! previous thermal inertia [SI]
 !Inputs
 !------
Index: trunk/LMDZ.COMMON/libf/evolution/soil_evolution_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_evolution_mod.F90	(revision 2850)
+++ 	(revision )
@@ -1,139 +1,0 @@
-      MODULE soil_evolution_mod
-
-      IMPLICIT NONE
-
-      LOGICAL soil_pem  ! True by default, to run with the subsurface physic. Read in pem.def
-      
-      CONTAINS    
-
-     subroutine soil_pem_CN(ngrid,nsoil, &
-               therm_i, timestep,tsurf,tsurf_prev,tsoil)
-
-
-      use comsoil_h_PEM, only:   & fluxgeo,layer_PEM
-      use comsoil_h,only: volcapa
-      implicit none
-
-!-----------------------------------------------------------------------
-!  Author: LL
-!  Purpose: Compute soil temperature using an implict 1st order scheme
-!  
-!  Note: depths of layers and mid-layers, soil thermal inertia and 
-!        heat capacity are commons in comsoil_PEM.h
-!        A convergence loop is added until equilibrium
-!-----------------------------------------------------------------------
-
-#include "dimensions.h"
-
-!-----------------------------------------------------------------------
-!  arguments
-!  ---------
-!  inputs:
-      integer,intent(in) :: ngrid	! number of (horizontal) grid-points 
-      integer,intent(in) :: nsoil	! number of soil layers  
-      real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia
-      real,intent(in) :: timestep	    ! time step
-      real,intent(in) :: tsurf(ngrid)   ! surface temperature
-      real,intent(in) :: tsurf_prev(ngrid)   ! surface temperature at previous timestep
-
-! outputs:
-      real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
-
-! local variables:
-      integer ig,ik,isoil
-      REAL :: ice_inertia = 2120.
-      real :: alph_PEM(nsoil)
-      real :: beta_PEM(nsoil)
-      real :: rhoc(nsoil)
-      real :: volcapa_ice = 1.43e7
-
-      real :: k_soil(nsoil)
-      real :: d1(nsoil),d2(nsoil),d3(nsoil),re(nsoil)
-      real :: Tcol(nsoil)
-
-do ig = 1,ngrid
-
-   Tcol(:) = tsoil(ig,:)
-!0. Build thermal properties of the ground
-
-  do isoil = 1,nsoil
-    if(abs(ice_inertia-therm_i(ig,isoil)).lt.1e-6) then
-      rhoc(isoil) = volcapa_ice
-    else
-      rhoc(isoil) = volcapa
-    endif
-    k_soil(isoil) = therm_i(ig,isoil)*therm_i(ig,isoil)/rhoc(isoil)
-  enddo
-
-!1. Build Crank-Nicholson scheme for heat equation, following Schorgofer derivations
-!1.a all the points which are not impacted by boundary condition
-  do isoil = 2,nsoil-1
-    
-    alph_pem(isoil) = timestep/((rhoc(isoil) + rhoc(isoil+1))/2.)*k_soil(isoil+1)/((layer_PEM(isoil+1) - layer_PEM(isoil))*(layer_PEM(isoil+1) - layer_PEM(isoil-1)))
-    beta_pem(isoil) = timestep/((rhoc(isoil) + rhoc(isoil+1))/2.)*k_soil(isoil)/((layer_PEM(isoil) - layer_PEM(isoil-1))*(layer_PEM(isoil+1) - layer_PEM(isoil-1)))
-  enddo
-
-!1.b At the surface: T(z=0) = Tsurf.
-
-alph_pem(1) = k_soil(2)*timestep/(rhoc(1)*layer_PEM(2)*(layer_PEM(2)-layer_PEM(1)))
-beta_pem(1) = k_soil(1)*timestep/(rhoc(1)*layer_PEM(2)*layer_PEM(1))
-
-!1.c At the bottom: dT/dz = -Fgeo/k.
-beta_pem(nsoil) = timestep*k_soil(nsoil)/(2.*rhoc(nsoil)*(layer_PEM(isoil) - layer_PEM(isoil-1))**2)
-
-!1 .d set the tridiagonal with diagonals d1,d2,d3, and right element re: 
-
-
-  d1(:) = -beta_pem(:)   
-  d2(:) = 1. + alph_pem(:) + beta_pem(:)
-  d3(:) = -alph_pem(:)   
-  d2(nsoil) = 1. + beta_pem(nsoil)
- 
-  re(1)= alph_pem(1)*Tcol(2) + (1.-alph_pem(1)-beta_pem(1))*Tcol(1) + beta_pem(1)*(tsurf(ig)+tsurf_prev(ig))
-  do isoil=2,nsoil-1
-     re(isoil) = beta_pem(isoil)*Tcol(isoil-1) + (1.-alph_pem(isoil)-beta_pem(isoil))*Tcol(isoil) + alph_pem(isoil)*Tcol(isoil+1)
-  enddo
-  re(nsoil) = beta_pem(nsoil)*Tcol(nsoil-1) + (1.-beta_pem(nsoil))*Tcol(nsoil) + &
-  timestep/rhoc(nsoil)*fluxgeo/(layer_PEM(nsoil)-layer_PEM(nsoil-1)) 
-
-!2. Update by tridiagonal inversion 
-  call tridag(d1,d2,d3,re,Tcol,nsoil) 
-  tsoil(ig,:) = Tcol(:)
-  enddo
-
-
-      END SUBROUTINE soil_pem_CN
-
-!====================================
-
-      SUBROUTINE TRIDAG(A,B,C,R,U,N)
-! From Numerical Recipes in Fortran 77
-      INTEGER, INTENT(IN)::  N
-      INTEGER ,PARAMETER:: NMAX = 1000
-
-      REAL, INTENT(IN) :: A(N),B(N),C(N),R(N)
-      REAL,INTENT(INOUT) :: U(N)
-
-      INTEGER j
-      REAL bet,gam(NMAX)
-
-      IF(B(1).EQ.0.)PAUSE
-      BET=B(1)
-      U(1)=R(1)/BET
-      DO 11 J=2,N
-        GAM(J)=C(J-1)/BET
-        BET=B(J)-A(J)*GAM(J)
-        IF(BET.EQ.0.)PAUSE
-        U(J)=(R(J)-A(J)*U(J-1))/BET
-11    CONTINUE
-      DO 12 J=N-1,1,-1
-        U(J)=U(J)-GAM(J+1)*U(J+1)
-12    CONTINUE
-      RETURN
-      END  SUBROUTINE TRIDAG
-
-
-      END MODULE soil_evolution_mod
-
-
-      
Index: trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90	(revision 2855)
@@ -25,15 +25,15 @@
 !  ---------
 !  inputs:
-      integer,intent(in) :: ngrid	! number of (horizontal) grid-points 
-      integer,intent(in) :: nsoil	! number of soil layers  
-      logical,intent(in) :: firstcall ! identifier for initialization call 
-      real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia
-      real,intent(in) :: timestep	    ! time step
-      real,intent(in) :: tsurf(ngrid)   ! surface temperature
+      integer,intent(in) :: ngrid       ! number of (horizontal) grid-points 
+      integer,intent(in) :: nsoil       ! number of soil layers  
+      logical,intent(in) :: firstcall    ! identifier for initialization call 
+      real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia [SI]
+      real,intent(in) :: timestep       ! time step [s]
+      real,intent(in) :: tsurf(ngrid)   ! surface temperature [K]
  
 ! outputs:
-      real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
-      real,intent(inout) :: alph_PEM(ngrid,nsoil-1)
-      real,intent(inout) :: beta_PEM(ngrid,nsoil-1)
+      real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature [K]
+      real,intent(inout) :: alph_PEM(ngrid,nsoil-1) ! intermediate for computations [1]
+      real,intent(inout) :: beta_PEM(ngrid,nsoil-1) ! intermediate for computations [K]
       
 ! local variables:
Index: trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F	(revision 2850)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F	(revision 2855)
@@ -12,5 +12,5 @@
 
 !======================================================================
-!  Author: Ehouarn Millour (07/2006)
+!  Author: LL, based on work by  Ehouarn Millour (07/2006)
 !
 !  Purpose: Read and/or initialise soil depths and properties
@@ -22,17 +22,7 @@
 !
 !
-!  This subroutine reads from a NetCDF file (opened by the caller)
-!  of "startfi.nc" format.
 !  The various actions and variable read/initialized are:
-!  1. Check out the number of soil layers (in datafile); if it isn't equal
-!     to nsoil, then some interpolation will be required
-!     Also check if data in file "startfi.nc" is in older format (ie:
-!     thermal inertia was depth-independent; and there was no "depth"
-!     coordinate.
-!     Read/build layer (and midlayer) depths
-!  2. Read volumetric specific heat (or initialise it to default value)
-!  3. Read Thermal inertia
-!  4. Read soil temperatures
-!  5. Interpolate thermal inertia and temperature on the new grid, if
+!  1. Read/build layer (and midlayer) depth
+!  2. Interpolate thermal inertia and temperature on the new grid, if
 !     necessary
 !======================================================================
@@ -46,6 +36,6 @@
       integer,intent(in) :: nsoil_PEM	! # of soil layers in the PEM
       integer,intent(in) :: nsoil_GCM	! # of soil layers in the GCM
-      real,intent(in) :: TI_GCM(ngrid,nsoil_GCM,nslope)	! # of soil layers in the GCM
-      real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope)	! # of soil layers in the PEM
+      real,intent(in) :: TI_GCM(ngrid,nsoil_GCM,nslope)	!  Thermal inertia  in the GCM [SI]
+      real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope)	! Thermal inertia   in the PEM [SI]
 
 !======================================================================
