Index: trunk/LMDZ.COMMON/libf/evolution/ComputeCO2Adsorption.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/ComputeCO2Adsorption.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/ComputeCO2Adsorption.F90	(revision 2794)
@@ -0,0 +1,97 @@
+   SUBROUTINE regolith_co2adsorption_old(ngrid,nslope,nsoil_PEM,ps,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,m_co2_completesoil,delta_mreg)
+      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(:,:,:)
+
+
+!=======================================================================
+      RETURN
+      END
Index: trunk/LMDZ.COMMON/libf/evolution/Example_ice_evol.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/Example_ice_evol.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/Example_ice_evol.F90	(revision 2794)
@@ -0,0 +1,27 @@
+         ! Update surface ice distribution to iterate to steady state if requested
+         if(ice_update)then
+
+            do ig=1,ngrid
+
+               delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
+               
+               ! add multiple years of evolution
+               qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep 
+
+               ! if ice has gone -ve, set to zero
+               if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
+                  qsurf_hist(ig,igcm_h2o_ice) = 0.0 
+               endif
+
+               ! if ice is seasonal, set to zero (NEW)
+               if(ice_min(ig).lt.0.01)then
+                  qsurf_hist(ig,igcm_h2o_ice) = 0.0 
+               endif
+
+            enddo
+
+            ! enforce ice conservation
+            ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*cell_area(:) )/SUM(cell_area(:))
+            qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
+
+         endif
Index: trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 2794)
@@ -0,0 +1,286 @@
+  module adsorption_mod  
+  implicit none
+  contains
+
+
+
+  subroutine regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen, ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, theta_h2o_adsorbded,m_h2o_adsorbed)
+
+ use vertical_layers_mod, ONLY: ap,bp
+ use comsoil_h_PEM, only: n_1km
+
+ implicit none
+
+
+ INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen ! size dimension
+ REAL,INTENT(IN) :: ps(ngrid,timelen)                   ! surface pressure (Pa)          
+ REAL,INTENT(IN) :: q_co2(ngrid,timelen)                ! Mass mixing ratio of co2 in the first layer (kg/kg)
+ REAL,INTENT(IN) :: q_h2o(ngrid,timelen)                ! Mass mixing ratio of H2o in the first layer (kg/kg)
+ REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope)      ! Soil Thermal inertia (J/K/^2/s^1/2)
+ REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)   ! Soil temperature (K)
+
+
+! output
+ REAL,INTENT(OUT) ::  m_h2o_adsorbed(ngrid,nsoil_PEM,nslope)     ! Density of h2o adsorbed (kg/m^3)
+ REAL,INTENT(OUT) :: theta_h2o_adsorbded(ngrid,nsoil_PEM,nslope) ! Fraction of the pores occupied by H2O molecules
+
+
+! constant
+ REAL :: Ko =  1.57e-8            ! Jackosky et al. 1997
+ REAL :: e = 2573.9               ! Jackosky et al. 1997
+ REAL :: mu = 0.48                ! Jackosky et al. 1997
+ REAL ::  inertie_thresold = 800. ! TI > 800 means cementation
+ real :: m_h2o = 18.01528E-3      ! Molecular weight of h2o (kg/mol)
+ real :: m_co2 = 44.01E-3         ! Molecular weight of co2 (kg/mol)
+ real :: m_noco2 = 33.37E-3       ! Molecular weight of non co2 (kg/mol)
+ REAL ::  rho_regolith = 2000.    ! density of the reoglith, Buhler & Piqueux 2021
+ real :: alpha_clapeyron = -6143.7! eq. (2) in Murphy & Koop 2005
+ real :: beta_clapeyron = 29.9074 ! eq. (2) in Murphy & Koop 2005
+ real :: mi = 2.84e-7             ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
+ real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
+
+
+
+
+
+! local variable
+ real :: A,B                                                   ! Used to compute the mean mass above the surface
+ real :: K                                                     ! Used to compute theta
+ real :: p_sat                                                 ! saturated vapor pressure of ice
+ integer ig,iloop, islope,isoil,it                             ! for loops
+ real,allocatable :: mass_mean(:,:)                            ! mean mass above the surface
+ real,allocatable :: zplev_mean(:,:)                           ! pressure above the surface
+ real,allocatable :: pvapor(:,:)                               ! partial pressure above the surface
+ real, allocatable :: pvapor_avg(:)                            ! yearly average vapor pressure above the surface
+
+
+
+! 0. Some initializations
+
+
+   
+     allocate(mass_mean(ngrid,timelen))
+     allocate(zplev_mean(ngrid,timelen))
+     allocate(pvapor(ngrid,timelen))
+     allocate(pvapor_avg(ngrid))
+
+      m_h2o_adsorbed(:,:,:) = 0.
+      theta_h2o_adsorbded(:,:,:) = 0.
+      A =(1/m_co2 - 1/m_noco2)
+      B=1/m_noco2
+! 1. Compute rho surface yearly averaged
+
+!   1.1 Compute the partial pressure of vapor
+!a. the molecular mass into the column
+     do ig = 1,ngrid
+       mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
+     enddo
+
+
+! b. pressure level
+     do it = 1,timelen
+       do ig = 1,ngrid
+         zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
+       enddo
+     enddo
+
+! c. Vapor pressure
+     pvapor(:,:) = mass_mean(:,:)/m_h2o*q_h2o(:,:)*zplev_mean(:,:)
+     pvapor_avg(:) = sum(pvapor(:,:),2)/timelen
+
+     deallocate(pvapor)
+     deallocate(zplev_mean)
+     deallocate(mass_mean)
+
+
+
+
+
+! 2. 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)  then
+       K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
+       theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1+K*pvapor_avg(ig)))**mu
+       m_h2o_adsorbed(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*mi*rho_regolith
+     else
+        p_sat =exp(alpha_clapeyron/tsoil_PEM(ig,iloop,islope) +beta_clapeyron) ! we assume fixed temperature in the ice ... not really:q good but ...
+        theta_h2o_adsorbded(ig,iloop,islope) = (K*p_sat/(1+K*p_sat))**mu
+        m_h2o_adsorbed(ig,iloop,islope) =as*theta_h2o_adsorbded(ig,iloop,islope)*mi*rho_regolith
+      endif
+ 
+   enddo
+  enddo
+ enddo
+
+ RETURN
+  end subroutine
+
+
+  SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_completesoil,delta_mreg)
+      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
+      use vertical_layers_mod, ONLY: ap,bp
+
+      IMPLICIT NONE
+! Input:  
+ INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen             ! size dimension
+ REAL,INTENT(IN) :: ps(ngrid,timelen)                               ! 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) :: q_co2(ngrid,timelen),q_h2o(ngrid,timelen)       ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg)
+ 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)   ! Density of co2 adsorbed (kg/m^3)
+ REAL,INTENT(INOUT) :: delta_mreg(ngrid)                            ! Difference density of co2 adsorbed (kg/m^3)
+
+
+  
+! Constants:
+
+ REAL :: alpha = 7.512e-6 ! Zent & Quinn 1995
+ REAL :: beta =  -1541.5  ! Zent & Quinn 1995
+ REAL ::  inertie_thresold = 800. ! TI > 800 means cementation
+ REAL ::  rho_regolith = 2000. ! density of the reoglith, buhler & piqueux 2021
+ real :: m_co2 = 44.01E-3      ! Molecular weight of co2 (kg/mol)
+ real :: m_noco2 = 33.37E-3    ! Molecular weight of h2o (kg/mol)
+ real :: m_theta = 4.27e-7     ! Mass of co2 per m^2 absorbed 
+ real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
+
+! Local          
+ real :: A,B                                             ! Used to compute the mean mass above the surface
+ INTEGER :: ig,islope,iloop,it                           ! for loops
+ REAL :: dm_co2_regolith_slope(ngrid,nsoil_PEM,nslope)   ! elementary mass adsorded per mesh per slope
+ INTEGER :: ispermanent_co2glaciers(ngrid,nslope)        ! Check if the glacier is permanent
+ INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)        ! Check if the glacier is permanent
+ REAL :: deltam_reg_complete(ngrid,n_1km,nslope)         ! Difference in the mass per slope and soil layer (kg/m^3)
+ REAL :: deltam_reg_slope(ngrid,nslope)                  ! Difference in the mass per slope  (kg/m^3)
+ REAL :: m_h2o_adsorbed(ngrid,nsoil_PEM,nslope)          ! Density of CO2 adsorbed (kg/m^3)
+ REAL :: theta_h2o_adsorbed(ngrid,nsoil_PEM,nslope)     ! Fraction of the pores occupied by H2O molecules
+!timelen array are allocated because heavy ...
+ real,allocatable :: mass_mean(:,:)                            ! mean mass above the surface
+ real,allocatable :: zplev_mean(:,:)                           ! pressure above the surface
+ real,allocatable :: pco2(:,:)                                  ! partial pressure above the surface
+ real, allocatable :: pco2_avg(:)                              ! yearly averaged
+
+
+
+! 0. Some initializations
+
+
+
+     allocate(mass_mean(ngrid,timelen))
+     allocate(zplev_mean(ngrid,timelen))
+     allocate(pco2(ngrid,timelen))
+     allocate(pco2_avg(ngrid))
+  
+      m_h2o_adsorbed(:,:,:) = 0.
+      A =(1/m_co2 - 1/m_noco2)
+      B=1/m_noco2
+
+     dm_co2_regolith_slope(:,:,:) = 0
+     delta_mreg(:) = 0.
+
+
+
+!0.1 Look at perenial ice
+  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
+
+
+!   0.2  Compute the partial pressure of CO2
+!a. the molecular mass into the column
+     do ig = 1,ngrid
+       mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
+     enddo
+
+
+! b. pressure level
+     do it = 1,timelen
+       do ig = 1,ngrid
+         zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
+       enddo
+     enddo
+
+! c. Vapor pressure
+     pco2(:,:) = mass_mean(:,:)/m_co2*q_co2(:,:)*zplev_mean(:,:)
+     pco2_avg(:) = sum(pco2(:,:),2)/timelen
+
+ 
+     deallocate(zplev_mean)
+     deallocate(mass_mean)
+     deallocate(pco2)
+
+
+! 1. Compute the fraction of the pores occupied by H2O
+ call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen, ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,theta_h2o_adsorbed, m_h2o_adsorbed)
+
+
+
+! 2.  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) = as*rho_regolith*m_theta*(1-theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
+                                             (alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
+     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) = as*rho_regolith*m_theta*(1-theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig) &
+                                                 /(alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope))) 
+        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(:,:,:)
+
+
+!=======================================================================
+      RETURN
+      END
+
+
+   end module
Index: trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod_slope.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod_slope.F90	(revision 2792)
+++ trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod_slope.F90	(revision 2794)
@@ -46,11 +46,17 @@
   ENDDO
 
+     print *, "jjm_input+1", jjm_input+1
+     print *, "iim_input+1", iim_input+1
+     print *, "nslope+1", nslope+1
+
 !  If the difference is too small; there is no evolution
   DO j=1,jjm_input+1
     DO i = 1, iim_input
        DO islope = 1, nslope
+!         print *, "tendencies_h2o_ice(i,j,islope)LAAA", tendencies_h2o_ice(i,j,islope)
          if(abs(tendencies_h2o_ice(i,j,islope)).LT.1.0E-10) then
             tendencies_h2o_ice(i,j,islope)=0.
          endif
+!         print *, "tendencies_h2o_ice(i,j,islope)HERE", tendencies_h2o_ice(i,j,islope)
        enddo
     ENDDO
@@ -71,4 +77,6 @@
   tendencies_h2o_ice_phys(ig0,:) = tendencies_h2o_ice(1,jjm_input+1,:)
 
+!  print *, "tendencies_h2o_ice_physze", tendencies_h2o_ice_phys(:,:)
+
 
 END SUBROUTINE compute_tendencies_slope
Index: trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90	(revision 2794)
@@ -0,0 +1,182 @@
+   SUBROUTINE computeice_table(timelen,ngrid,nslope,nsoil_GCM,nsoil_PEM,tsoil,tsurf,q_co2,q_h2o,ps,ice_table)
+    USE comsoil_h, only:  inertiedat, volcapa
+    USE vertical_layers_mod, ONLY: ap,bp
+    USE comsoil_h_PEM, only: layer_PEM 
+    implicit none 
+
+
+    integer,intent(in) :: timelen,ngrid,nslope,nsoil_PEM,nsoil_GCM
+    real,intent(in) :: tsoil(ngrid,nsoil_PEM,nslope,timelen)    ! soil temperature, timeseries [K]
+    real,intent(in) :: tsurf(ngrid,nslope,timelen)              ! surface temperature [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(ngrid,timelen)                        ! surface pressure [Pa]
+    real,intent(out) :: ice_table(ngrid,nslope)                 ! ice table [m]
+
+
+    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
+
+
+    integer ig, islope,isoil,it
+    real,allocatable :: mass_mean(:,:)                            ! mean mass above the surface
+    real,allocatable :: zplev_mean(:,:)                           ! pressure above the surface
+    real,allocatable :: pvapor(:,:)                               ! partial pressure above the surface
+
+    real,allocatable :: rhovapor(:,:,:)
+    real,allocatable :: 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
+
+
+
+     allocate(rhovapor(ngrid,nslope,timelen))
+     allocate(rhovapor_avg(ngrid,nslope))
+     allocate(pvapor(ngrid,timelen))
+
+     allocate(mass_mean(ngrid,timelen))
+     allocate(zplev_mean(ngrid,timelen))
+
+
+
+
+! 0. Some initializations
+
+      A =(1/m_co2 - 1/m_noco2)
+      B=1/m_noco2
+! 1. Compute rho surface yearly averaged
+
+!   1.1 Compute the partial pressure of vapor
+!a. the molecular mass into the column
+     do ig = 1,ngrid
+       mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
+     enddo
+
+
+! b. pressure level
+     do it = 1,timelen
+       do ig = 1,ngrid
+         zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
+       enddo
+     enddo
+
+! c. Vapor pressure
+     pvapor(:,:) = mass_mean(:,:)/m_h2o*q_h2o(:,:)*zplev_mean(:,:)
+    
+     deallocate(mass_mean)
+     deallocate(zplev_mean)
+  
+!    1.2   Check if there is frost at the surface and then compute the density
+!    at the surface
+     do ig = 1,ngrid
+       do islope = 1,nslope
+         do it = 1,timelen
+           psv_surf = exp(alpha/tsurf(ig,islope,it) +beta)
+          if ((isnan(psv_surf)).or.(isnan(pvapor(ig,it)))) then
+          write(*,*) 'pb vapor',ig,islope,it
+         stop
+         endif       
+         rhovapor(ig,islope,it) = min(psv_surf,pvapor(ig,it))/tsurf(ig,islope,it)
+         enddo
+       enddo
+     enddo
+     deallocate(pvapor)
+
+
+
+
+!    1.3 Density at the surface yearly averaged
+     rhovapor_avg(:,:) = SUM(rhovapor(:,:,:),3)/timelen
+
+
+     deallocate(rhovapor)
+
+
+
+
+! 2. Compute rho soil vapor
+   
+     allocate(rho_soil_avg(ngrid,nslope,nsoil_PEM))
+     allocate(rho_soil(ngrid,nslope,nsoil_PEM,timelen))
+
+     do ig = 1,ngrid
+       do islope = 1,nslope
+         do isoil = 1,nsoil_PEM
+            do it = 1,timelen 
+              rho_soil(ig,islope,isoil,it) = exp(alpha/tsoil(ig,isoil,islope,it) +beta)/tsoil(ig,isoil,islope,it)        
+             enddo
+           enddo
+       enddo
+     enddo
+
+
+
+    rho_soil_avg(:,:,:) = SUM( rho_soil(:,:,:,:),4)/timelen
+    deallocate(rho_soil)
+
+
+
+! 3. Computing ice table
+   
+    ice_table (:,:) = -1e4
+
+         allocate(diff_rho(ngrid,nslope,nsoil_PEM))
+
+         do isoil = 1,nsoil_PEM
+
+
+             diff_rho(:,:,isoil) = rhovapor_avg(:,:) - rho_soil_avg(:,:,isoil)
+!             write(*,*) 'diff =',ig,islope,isoil,diff_rho(ig,islope,isoil),rhovapor_avg(ig,islope) ,rho_soil_avg(ig,islope,isoil)
+         enddo
+    
+
+  deallocate(rhovapor_avg)
+      deallocate(rho_soil_avg)
+
+
+
+
+
+
+
+     do ig = 1,ngrid
+       do islope = 1,nslope
+         if(diff_rho(ig,islope,1) > 0) then
+           ice_table(ig,islope) = 0.
+         else
+           do isoil = 1,nsoil_PEM -1
+             if((diff_rho(ig,islope,isoil).lt.0).and.(diff_rho(ig,islope,isoil+1).gt.0.)) then
+               z1 = (diff_rho(ig,islope,isoil) - diff_rho(ig,islope,isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
+               z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(ig,islope,isoil+1)
+               ice_table(ig,islope) = -z2/z1
+               exit
+             endif
+           enddo
+          endif
+        enddo
+      enddo
+
+    
+    deallocate(diff_rho)
+
+
+
+
+
+
+!=======================================================================
+      RETURN
+
+
+
+
+
+      END
Index: trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 2794)
@@ -0,0 +1,66 @@
+module comsoil_h_PEM
+
+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
+
+  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
+  ! 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)
+
+contains
+
+  subroutine ini_comsoil_h_PEM(ngrid,nslope)
+  
+  implicit none
+  integer,intent(in) :: ngrid ! number of atmospheric columns
+  integer,intent(in) :: nslope ! number of slope within a mesh 
+
+    allocate(layer_PEM(nsoilmx_PEM)) !soil layer depths
+    allocate(mlayer_PEM(0:nsoilmx_PEM-1)) ! soil mid-layer depths
+    allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope)) ! soil thermal inertia
+    allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope)) ! soil temperatures
+    allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1))
+    allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1))
+    allocate(coefq_PEM(0:nsoilmx_PEM-1))
+    allocate(coefd_PEM(ngrid,nsoilmx_PEM-1))
+    allocate(alph_PEM(ngrid,nsoilmx_PEM-1,nslope))
+    allocate(beta_PEM(ngrid,nsoilmx_PEM-1,nslope))
+    allocate(inertiedat_PEM(ngrid,nsoilmx_PEM)) ! soil thermal inertia
+    allocate(co2_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
+  end subroutine ini_comsoil_h_PEM
+
+
+  subroutine end_comsoil_h_PEM
+
+  implicit none
+
+    if (allocated(layer_PEM)) deallocate(layer_PEM)
+    if (allocated(mlayer_PEM)) deallocate(mlayer_PEM)
+    if (allocated(TI_PEM)) deallocate(TI_PEM)
+    if (allocated(tsoil_PEM)) deallocate(tsoil_PEM)
+    if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM)
+    if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM)
+    if (allocated(coefq_PEM)) deallocate(coefq_PEM) 
+    if (allocated(coefd_PEM)) deallocate(coefd_PEM)
+    if (allocated(alph_PEM)) deallocate(alph_PEM)
+    if (allocated(beta_PEM)) deallocate(beta_PEM)
+    if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
+    if (allocated(co2_adsorbded_phys)) deallocate(co2_adsorbded_phys)
+  end subroutine end_comsoil_h_PEM
+
+end module comsoil_h_PEM
Index: trunk/LMDZ.COMMON/libf/evolution/criterion_co2_ice_stop_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/criterion_co2_ice_stop_mod.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/criterion_co2_ice_stop_mod.F90	(revision 2794)
@@ -0,0 +1,72 @@
+!
+! $Id $
+!
+SUBROUTINE criterion_co2_ice_stop(cell_area,initial_co2_ice,co2ice,STOPPING,ngrid,latitude,n_band_lat)
+
+  USE temps_mod_evol, ONLY: alpha_criterion
+
+      IMPLICIT NONE
+
+!=======================================================================
+!
+!  Routine that checks if the criterion to stop the PEM is reached
+!
+!=======================================================================
+
+!   arguments:
+!   ----------
+
+!   INPUT
+  INTEGER, intent(in) :: ngrid                  ! # of grid physical grid points 
+  REAL,    intent(in) :: cell_area(ngrid)       ! physical point field : Area of the cells
+  REAL,    intent(in) ::  co2ice(ngrid)          ! physical point field : Actual density of water ice
+  REAL,    intent(in) ::  latitude(ngrid)          ! physical point field : Latitude
+  REAL,    intent(in) ::  initial_co2_ice(n_band_lat)  ! Initial/Actual surface of water ice
+
+
+
+!   OUTPUT
+  LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
+
+!   local:
+!   -----
+  INTEGER :: i,j,n_band_lat                    ! Loop
+  REAL :: present_co2(n_band_lat)  ! Initial/Actual surface of water ice
+  REAL :: pi
+
+!=======================================================================
+
+      pi=4.D0*DATAN(1.D0)
+
+!   initialisation to false
+    STOPPING=.FALSE.
+
+     do j=1,n_band_lat
+        present_co2(j)=0.
+     enddo
+
+  do i=1,ngrid
+            j=floor((latitude(i)+(pi/2))/(pi)*n_band_lat)+1
+      if(j.GT.n_band_lat) then
+          j=n_band_lat
+      endif
+      present_co2(j)=present_co2(j)+co2ice(i)*cell_area(i)
+  enddo
+  
+!   check of the criterion
+  do j=1,n_band_lat
+    if(present_co2(j).LT.initial_co2_ice(j)*(1-alpha_criterion) .OR. &
+       present_co2(j).GT.initial_co2_ice(j)*(1+alpha_criterion)) then
+         STOPPING=.TRUE. 
+         print *, "j", j
+         print *, "present_co2(j)", present_co2(j)
+         print *, "initial_co2_ice(j)", initial_co2_ice(j)
+    endif
+  enddo
+
+END SUBROUTINE criterion_co2_ice_stop
+
+
+
+
+
Index: trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_slope.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_slope.F90	(revision 2792)
+++ trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_slope.F90	(revision 2794)
@@ -19,20 +19,20 @@
 
 !   INPUT
-  INTEGER, intent(in) :: ngrid,nslope                  ! # of grid physical grid points , # subslopes
-  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
+  INTEGER, intent(in) :: ngrid,nslope                  ! # of grid physical grid points 
+  REAL,    intent(in) :: cell_area(ngrid)       ! physical point field : Area of the cells
   REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
-  REAL,    intent(in) :: ini_surf                      ! Surface of sublimating ice in the GCM output
-  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope) ! Boolean to check if it is a sublimating area
-  REAL,    intent(in) :: global_ave_press_GCM          ! Average pressure in the GCM output
-  REAL,    intent(in) :: global_ave_press_new          ! Actual pressure
+  REAL,    intent(in) :: ini_surf
+  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
+  REAL,    intent(in) :: global_ave_press_GCM
+  REAL,    intent(in) :: global_ave_press_new
 
 
 !   OUTPUT
-  LOGICAL, intent(out) :: STOPPING                     ! Logical : is the criterion reached?
+  LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
 
 !   local:
 !   -----
-  INTEGER :: i,islope                                  ! Loop
-  REAL :: present_surf                                 ! Initial/Actual surface of water ice
+  INTEGER :: i,islope                    ! Loop
+  REAL :: present_surf  ! Initial/Actual surface of water ice
 
 !=======================================================================
@@ -41,13 +41,25 @@
     STOPPING=.FALSE.
 
-!   computation of the actual surface of sublimating ice
+!   computation of the actual surface
   present_surf=0.
   do i=1,ngrid
    do islope=1,nslope
       if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
+         print *, "i", i
+         print *, "initial_h2o_ice(i,islope)", initial_h2o_ice(i,islope)
+         print *, "qsurf(i,:)", qsurf(i,:)
+         print *, "cell_area(i)", cell_area(i)
+         print *, "present_surf",present_surf
          present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
       endif
    enddo
   enddo
+
+!  print *, "initial_h2o_ice", initial_h2o_ice
+!  print *, "qsurf", qsurf
+
+  print *, "present_surf", present_surf
+  print *, "ini_surf", ini_surf
+  print *, "ini_surf*0.8", ini_surf*(1-alpha_criterion)
   
 !   check of the criterion
Index: trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_water_slope.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_water_slope.F90	(revision 2792)
+++ trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_water_slope.F90	(revision 2794)
@@ -5,5 +5,5 @@
 
   USE temps_mod_evol, ONLY: alpha_criterion
-  use comslope_mod, ONLY: subslope_dist,nslope
+	  use comslope_mod, ONLY: subslope_dist,nslope
 
       IMPLICIT NONE
@@ -19,18 +19,18 @@
 
 !   INPUT
-  INTEGER, intent(in) :: ngrid                         ! # of grid physical grid points 
-  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
+  INTEGER, intent(in) :: ngrid                  ! # of grid physical grid points 
+  REAL,    intent(in) :: cell_area(ngrid)       ! physical point field : Area of the cells
   REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
-  REAL,    intent(in) :: ini_surf                      ! Surface of sublimating ice in the GCM output
-  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope) ! Boolean to check if it is a sublimating area
+  REAL,    intent(in) :: ini_surf
+  REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
 
 
 !   OUTPUT
-  LOGICAL, intent(out) :: STOPPING                     ! Logical : is the criterion reached?
+  LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
 
 !   local:
 !   -----
-  INTEGER :: i,islope                                  ! Loop
-  REAL :: present_surf                                 ! Initial/Actual surface of water ice
+  INTEGER :: i,islope                    ! Loop
+  REAL :: present_surf  ! Initial/Actual surface of water ice
 
 !=======================================================================
@@ -44,8 +44,20 @@
     do islope=1, nslope
       if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
+!         print *, "i", i
+!         print *, "initial_h2o_ice(i,islope)", initial_h2o_ice(i,islope)
+!         print *, "qsurf(i,islope)", qsurf(i,islope)
+!         print *, "cell_area(i)", cell_area(i)
+!         print *, "present_surf",present_surf
          present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
       endif
     enddo
   enddo
+
+!  print *, "initial_h2o_ice", initial_h2o_ice
+!  print *, "qsurf", qsurf
+
+!  print *, "present_surf", present_surf
+!  print *, "ini_surf", ini_surf
+!  print *, "ini_surf*0.8", ini_surf*(1-alpha_criterion)
   
 !   check of the criterion
Index: trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod_slope.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod_slope.F90	(revision 2792)
+++ trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod_slope.F90	(revision 2794)
@@ -20,10 +20,11 @@
 !   INPUT
 
-  INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope      ! # of grid points along longitude/latitude/ total
-  REAL, intent(in) ::  cell_area(ngrid)                         ! Area of each cell
+  INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope   ! # of grid points along longitude/latitude/ total
+!  REAL, intent(in) ::  tendencies_h2o_ice_phys(ngrid) ! physical point field : Evolution of perenial ice over one year
+  REAL, intent(in) ::  cell_area(ngrid)
 
 !   OUTPUT
-  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                   ! physical point field : Previous and actual density of water ice
-  LOGICAL :: STOPPING                                           ! Not used
+  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                ! physical point field : Previous and actual density of water ice
+  LOGICAL :: STOPPING
   REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
 
@@ -33,9 +34,12 @@
 
   INTEGER :: i,j,ig0,islope                                  ! loop variable
+!  REAL :: not_budget, budget
+  REAL :: pos_tend, neg_tend, real_coefficient,negative_part
+  REAL ::  new_tendencies(ngrid)
 
   STOPPING=.false.
 
 
-! Evolution of the co2 ice for each physical point
+! Evolution of the water ice for each physical point
   do i=1,ngrid
     do islope=1,nslope
Index: trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90	(revision 2792)
+++ trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90	(revision 2794)
@@ -6,17 +6,12 @@
 
   USE temps_mod_evol, ONLY: dt_pem
-  use comslope_mod,   ONLY: subslope_dist
+  	  use comslope_mod, ONLY: subslope_dist
 
       IMPLICIT NONE
 
 !=======================================================================
-!  
-!  Routine that compute the evolution of the water ice.
-!  
-!  We suppose that the water exchange uniquely form one point to another 
-!  (water ice disappear from a sublimating point and increase at another 
-!  accumulation point). Therfor the total positive and negative tendencies 
-!  must be equal to conserve the total amount of water
-!  
+!
+!  Routine that compute the evolution of the water ice
+!
 !=======================================================================
 
@@ -26,10 +21,11 @@
 !   INPUT
 
-  INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope      ! # of grid points along longitude/latitude/ total
-  REAL, intent(in) ::  cell_area(ngrid)                         ! Area of each cell
+  INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope   ! # of grid points along longitude/latitude/ total
+!  REAL, intent(in) ::  tendencies_h2o_ice_phys(ngrid) ! physical point field : Evolution of perenial ice over one year
+  REAL, intent(in) ::  cell_area(ngrid)
 
 !   OUTPUT
-  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                   ! physical point field : Previous and actual density of water ice
-  LOGICAL :: STOPPING						! Boolean to alert that there is no sublimating point anymore
+  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                ! physical point field : Previous and actual density of water ice
+  LOGICAL :: STOPPING
   REAL, intent(inout) ::  tendencies_h2o_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
 
@@ -39,12 +35,13 @@
 
   INTEGER :: i,j,ig0,islope                                  ! loop variable
-  REAL :: pos_tend, neg_tend, real_coefficient               ! The integral of positive/negative tendencies overall the planet. Ratio
-  REAL :: negative_part                                      ! After applying the tendencies, some point can be negative, we sum all of them in this variable
-  REAL :: new_tendencies(ngrid,nslope)                       ! Tendencies after taking into account the non physical negative part above
+!  REAL :: not_budget, budget
+  REAL :: pos_tend, neg_tend, real_coefficient,negative_part
+  REAL ::  new_tendencies(ngrid,nslope)
 
 
 !=======================================================================
 
-! Integral of the positive and negative tendencies
+!  budget=sum(qsurf(:))
+
   pos_tend=0.
   neg_tend=0.
@@ -62,5 +59,6 @@
   enddo
 
-! These two quantities must be equal to conserve water, therfor we apply a ratio to them
+!  print *, "pos_tend", pos_tend
+!  print *, "neg_tend", neg_tend
 
   if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then
@@ -68,4 +66,5 @@
        do islope=1,nslope
        if(tendencies_h2o_ice_phys(i,islope).LT.0) then
+!          print *, "pos_tend/neg_tend", pos_tend/neg_tend
           new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)*(pos_tend/neg_tend)
        else
@@ -75,4 +74,5 @@
      enddo
   elseif(neg_tend.LT.pos_tend .and. neg_tend.GT.0) then
+!          print *, "neg_tend/pos_tend", neg_tend/pos_tend
      do i=1,ngrid
        do islope=1,nslope
@@ -85,5 +85,7 @@
      enddo
   elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then 
-      call criterion_ice_stop_water_slope(cell_area,1,qsurf(:,:)*0,STOPPING,ngrid,cell_area)
+
+!      call criterion_ice_stop(cell_area,1,qsurf*0.,STOPPING,ngrid,cell_area)
+      call criterion_ice_stop_water_slope(cell_area,1,qsurf(:,:)*0.,STOPPING,ngrid,cell_area)
       do i=1,ngrid
          do islope=1,nslope
@@ -93,19 +95,38 @@
   endif
 
+ 
+  
+  negative_part = 0.
+
 
 ! Evolution of the water ice for each physical point
   do i=1,ngrid
     do islope=1, nslope
+!      qsurf(i)=qsurf(i)+tendencies_h2o_ice_phys(i)*dt_pem
       qsurf(i,islope)=qsurf(i,islope)+new_tendencies(i,islope)*dt_pem
-      if (qsurf(i,islope).lt.0) then         ! If we have negative number, we put them to zero add sum them
+!      budget=budget+tendencies_h2o_ice_phys(i)*dt_pem
+      if (qsurf(i,islope).lt.0) then
+!        not_budget=not_budget+qsurf(i)
+!       print *, "NNqsurf(i,islope)", qsurf(i,islope)
+!       print *, "NNnew_tendencies(i,islope)", new_tendencies(i,islope)
+!       print *, "NNtendencies_h2o_ice_phys(i,islope)", tendencies_h2o_ice_phys(i,islope)
         negative_part=negative_part-qsurf(i,islope)*cell_area(i)*subslope_dist(i,islope)
         qsurf(i,islope)=0.
         tendencies_h2o_ice_phys(i,islope)=0.
+!        print *, "NNineg", i
+      endif
+      if(qsurf(i,islope).NE.qsurf(i,islope)) then
+!          print *, "qsurf(i,islope)",qsurf(i,islope)
+!          print *, "new_tendencies",new_tendencies(i,islope)
+!          print *, "tendencies_h2o_ice_phys",tendencies_h2o_ice_phys(i,islope)
+!          print *, "i", i
+!          print *,"islope",islope
       endif
     enddo
   enddo
 
-! This negative part must be put somewhere else to conserve water, therfor we compute new tendencies
+!  print *, "negative_part", negative_part
   real_coefficient=negative_part/pos_tend
+!  print *, "real_coefficient", real_coefficient
 
   do i=1,ngrid
@@ -118,3 +139,8 @@
 
 
+
+! Conservation of water ice
+!  qsurf(:)=qsurf(:)*budget/sum(qsurf(:))
+
+
 END SUBROUTINE evol_h2o_ice_s_slope
Index: trunk/LMDZ.COMMON/libf/evolution/ini_soil_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/ini_soil_mod.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/ini_soil_mod.F90	(revision 2794)
@@ -0,0 +1,417 @@
+      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)
+
+
+
+      use vertical_layers_mod, only: ap,bp
+      use comsoil_h_PEM, only:  fluxgeo,layer_PEM,inertiedat_PEM
+      use comsoil_h,only: volcapa, nsoilmx
+
+   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) :: 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
+
+
+      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
+ 
+
+      END SUBROUTINE ini_icetable
+      subroutine soil_pem_1D(nsoil,firstcall, &
+               therm_i,                          &
+               timestep,tsurf,tsoil,alph,beta)
+
+
+      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
+                           mu_PEM,fluxgeo
+      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"
+!#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)
+
+! 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
+
+
+
+      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 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/interpolate_TIPEM_TIGCM.F90	(revision 2794)
@@ -0,0 +1,33 @@
+   subroutine interpolate_TIPEM_TIGCM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_PEM,TI_GCM)
+
+
+
+      implicit none
+
+
+!======================================================================
+!  arguments
+!  ---------
+!  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
+
+!local variable
+      integer :: ig,islope,iloop
+
+
+
+     do ig = 1,ngrid
+       do islope = 1,nslope
+         do iloop = 1,nsoil_GCM
+           TI_GCM(ig,iloop,islope) = TI_PEM(ig,iloop,islope)
+         enddo
+       enddo
+     enddo
+
+
+  end subroutine interpolate_TIPEM_TIGCM
Index: trunk/LMDZ.COMMON/libf/evolution/iostart_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/iostart_PEM.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/iostart_PEM.F90	(revision 2794)
@@ -0,0 +1,1040 @@
+MODULE iostart_PEM
+
+    IMPLICIT NONE
+PRIVATE
+    INTEGER,SAVE :: nid_start ! NetCDF file identifier for startfi.nc file
+    INTEGER,SAVE :: nid_restart ! NetCDF file identifier for restartfi.nc file
+    
+    ! restartfi.nc file dimension identifiers: (see open_restartphy())
+    INTEGER,SAVE :: idim1 ! "index" dimension
+    INTEGER,SAVE :: idim2 ! "physical_points" dimension
+    INTEGER,SAVE :: idim3 ! "subsurface_layers" dimension
+    INTEGER,SAVE :: idim4 ! "nlayer_plus_1" dimension
+    INTEGER,SAVE :: idim5 ! "number_of_advected_fields" dimension
+    INTEGER,SAVE :: idim6 ! "nlayer" dimension
+    INTEGER,SAVE :: idim7 ! "Time" dimension
+    INTEGER,SAVE :: idim8 ! "nslope" dimension
+    INTEGER,SAVE :: idim9 ! "inter slope" dimension
+    INTEGER,SAVE :: timeindex ! current time index (for time-dependent fields)
+    INTEGER,PARAMETER :: length=100 ! size of tab_cntrl array
+    
+    INTERFACE get_field
+      MODULE PROCEDURE Get_field_r1,Get_field_r2,Get_field_r3
+    END INTERFACE get_field
+    
+    INTERFACE get_var
+      MODULE PROCEDURE get_var_r0,Get_var_r1,Get_var_r2,Get_var_r3
+    END INTERFACE get_var
+
+    INTERFACE put_field
+      MODULE PROCEDURE put_field_r1,put_field_r2,put_field_r3
+    END INTERFACE put_field
+
+    INTERFACE put_var
+      MODULE PROCEDURE put_var_r0,put_var_r1,put_var_r2,put_var_r3
+    END INTERFACE put_var
+
+    PUBLIC nid_start, length
+    PUBLIC get_field,get_var,put_field,put_var
+    PUBLIC inquire_dimension, inquire_dimension_length
+    PUBLIC inquire_field, inquire_field_ndims
+    PUBLIC open_startphy,close_startphy,open_restartphy,close_restartphy
+    
+CONTAINS
+
+  SUBROUTINE open_startphy(filename)
+  USE netcdf, only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE, nf90_strerror
+  USE mod_phys_lmdz_para, only: is_master, bcast
+  IMPLICIT NONE
+    CHARACTER(LEN=*) :: filename
+    INTEGER          :: ierr
+
+    IF (is_master) THEN
+      ierr = NF90_OPEN (filename, NF90_NOWRITE, nid_start)
+      IF (ierr.NE.NF90_NOERR) THEN
+        write(*,*)'open_startphy: problem opening file '//trim(filename)
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("open_startphy","Cannot open file",1)
+      ENDIF
+    ENDIF
+    
+    CALL bcast(nid_start) ! tell all procs about nid_start
+  
+  END SUBROUTINE open_startphy
+
+  SUBROUTINE close_startphy
+  USE netcdf, only: NF90_CLOSE
+  USE mod_phys_lmdz_para, only: is_master
+  IMPLICIT NONE
+    INTEGER          :: ierr
+
+    IF (is_master) THEN
+        ierr = NF90_CLOSE (nid_start)
+    ENDIF
+
+  END SUBROUTINE close_startphy
+
+
+  FUNCTION inquire_field(Field_name)
+  ! check if a given field is present in the input file
+  USE netcdf, only: NF90_INQ_VARID, NF90_NOERR
+  USE mod_phys_lmdz_para, only: is_master, bcast
+  IMPLICIT NONE
+    CHARACTER(LEN=*),INTENT(IN) :: Field_name
+    LOGICAL :: inquire_field
+    INTEGER :: varid
+    INTEGER :: ierr
+    
+    IF (is_master) THEN
+      ierr=NF90_INQ_VARID(nid_start,Field_name,varid)
+      IF (ierr==NF90_NOERR) THEN
+        Inquire_field=.TRUE.
+      ELSE
+        Inquire_field=.FALSE.
+      ENDIF
+    ENDIF
+
+    CALL bcast(inquire_field)
+
+  END FUNCTION inquire_field
+
+
+  FUNCTION inquire_field_ndims(Field_name)
+  ! give the number of dimensions of "Field_name" stored in the input file 
+  USE netcdf, only: nf90_inq_varid, nf90_inquire_variable, &
+                    NF90_NOERR, nf90_strerror
+  USE mod_phys_lmdz_para, only: is_master, bcast
+  IMPLICIT NONE
+    CHARACTER(LEN=*),INTENT(IN) :: Field_name
+    INTEGER :: inquire_field_ndims
+    INTEGER :: varid
+    INTEGER :: ierr
+    
+    IF (is_master) THEN
+      ierr=nf90_inq_varid(nid_start,Field_name,varid)
+      ierr=nf90_inquire_variable(nid_start,varid,&
+                                  ndims=inquire_field_ndims)
+      IF (ierr.NE.NF90_NOERR) THEN
+        write(*,*)'inquire_field_ndims: problem obtaining ndims of '&
+                  //trim(field_name)
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("inquire_field_ndims","Failed to get ndims",1)
+      ENDIF
+    ENDIF
+
+    CALL bcast(inquire_field_ndims)
+
+  END FUNCTION inquire_field_ndims
+
+
+  FUNCTION inquire_dimension(Field_name)
+  ! check if a given dimension is present in the input file
+  USE netcdf, only: nf90_inq_dimid, NF90_NOERR
+  USE mod_phys_lmdz_para, only: is_master, bcast
+  IMPLICIT NONE
+    CHARACTER(LEN=*),INTENT(IN) :: Field_name
+    LOGICAL :: inquire_dimension
+    INTEGER :: varid
+    INTEGER :: ierr
+    
+    IF (is_master) THEN
+      ierr=NF90_INQ_DIMID(nid_start,Field_name,varid)
+      IF (ierr==NF90_NOERR) THEN
+        Inquire_dimension=.TRUE.
+      ELSE
+        Inquire_dimension=.FALSE.
+      ENDIF
+    ENDIF
+
+    CALL bcast(inquire_dimension)
+
+  END FUNCTION inquire_dimension
+
+  FUNCTION inquire_dimension_length(Field_name)
+  ! give the length of the "Field_name" dimension stored in the input file
+  USE netcdf, only: nf90_inquire_dimension, nf90_inq_dimid, &
+                    NF90_NOERR, nf90_strerror
+  USE mod_phys_lmdz_para, only: is_master, bcast
+  IMPLICIT NONE
+    CHARACTER(LEN=*),INTENT(IN) :: Field_name
+    INTEGER :: inquire_dimension_length
+    INTEGER :: varid
+    INTEGER :: ierr
+    
+    IF (is_master) THEN
+      ierr=nf90_inq_dimid(nid_start,Field_name,varid)
+      ierr=nf90_inquire_dimension(nid_start,varid,&
+                                  len=inquire_dimension_length)
+      IF (ierr.NE.NF90_NOERR) THEN
+        write(*,*)'inquire_field_length: problem obtaining length of '&
+                  //trim(field_name)
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("inquire_field_ndims","Failed to get length",1)
+      ENDIF
+    ENDIF
+
+    CALL bcast(inquire_dimension_length)
+
+  END FUNCTION inquire_dimension_length
+
+
+
+  SUBROUTINE Get_Field_r1(field_name,field,found,timeindex)
+  ! For a surface field
+  use mod_grid_phy_lmdz, only: klon_glo ! number of atmospheric columns (full grid)
+  IMPLICIT NONE
+    CHARACTER(LEN=*),INTENT(IN)    :: Field_name
+    REAL,INTENT(INOUT)               :: Field(:)
+    LOGICAL,INTENT(OUT),OPTIONAL   :: found 
+    INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
+
+    integer :: edges(4), corners(4)
+
+    edges(1)=klon_glo
+    edges(2:4)=1
+    corners(1:4)=1
+    if (PRESENT(timeindex)) then
+      corners(2)=timeindex
+    endif
+
+    IF (PRESENT(found)) THEN
+      CALL Get_field_rgen(field_name,field,1,corners,edges,found)
+    ELSE
+      CALL Get_field_rgen(field_name,field,1,corners,edges)
+    ENDIF
+      
+  END SUBROUTINE Get_Field_r1
+  
+  SUBROUTINE Get_Field_r2(field_name,field,found,timeindex)
+  ! For a "3D" horizontal-vertical field
+  use mod_grid_phy_lmdz, only: klon_glo ! number of atmospheric columns (full grid)
+  IMPLICIT NONE
+    CHARACTER(LEN=*),INTENT(IN)    :: Field_name
+    REAL,INTENT(INOUT)               :: Field(:,:)
+    LOGICAL,INTENT(OUT),OPTIONAL   :: found 
+    INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
+
+    integer :: edges(4), corners(4)
+
+    edges(1)=klon_glo
+    edges(2)=size(field,2)
+    edges(3:4)=1
+    corners(1:4)=1
+    if (PRESENT(timeindex)) then
+      corners(3)=timeindex
+    endif
+    
+    IF (PRESENT(found)) THEN
+      CALL Get_field_rgen(field_name,field,size(field,2),&
+                          corners,edges,found)
+    ELSE
+      CALL Get_field_rgen(field_name,field,size(field,2),&
+                          corners,edges)
+    ENDIF
+
+      
+  END SUBROUTINE Get_Field_r2
+  
+  SUBROUTINE Get_Field_r3(field_name,field,found,timeindex)
+  ! for a "4D" field surf/alt/??
+  use mod_grid_phy_lmdz, only: klon_glo ! number of atmospheric columns (full grid)
+  IMPLICIT NONE
+    CHARACTER(LEN=*),INTENT(IN)    :: Field_name
+    REAL,INTENT(INOUT)               :: Field(:,:,:)
+    LOGICAL,INTENT(OUT),OPTIONAL   :: found 
+    INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
+
+    integer :: edges(4), corners(4)
+
+    edges(1)=klon_glo
+    edges(2)=size(field,2)
+    edges(3)=size(field,3)
+    edges(4)=1
+    corners(1:4)=1
+    if (PRESENT(timeindex)) then
+      corners(4)=timeindex
+    endif
+    
+    IF (PRESENT(found)) THEN
+      CALL Get_field_rgen(field_name,field,size(field,2)*size(field,3),&
+                          corners,edges,found)
+    ELSE
+      CALL Get_field_rgen(field_name,field,size(field,2)*size(field,3),&
+                          corners,edges)
+    ENDIF
+      
+  END SUBROUTINE Get_Field_r3
+  
+  SUBROUTINE Get_field_rgen(field_name,field,field_size, &
+                            corners,edges,found)
+  USE netcdf
+  USE dimphy
+  USE mod_grid_phy_lmdz
+  USE mod_phys_lmdz_para
+  IMPLICIT NONE
+    CHARACTER(LEN=*) :: Field_name
+    INTEGER          :: field_size
+    REAL             :: field(klon,field_size)
+    INTEGER,INTENT(IN) :: corners(4)
+    INTEGER,INTENT(IN) :: edges(4)
+    LOGICAL,OPTIONAL :: found
+    
+    REAL    :: field_glo(klon_glo,field_size)
+    LOGICAL :: tmp_found
+    INTEGER :: varid
+    INTEGER :: ierr
+    
+    IF (is_master) THEN
+  
+      ierr=NF90_INQ_VARID(nid_start,Field_name,varid)
+      
+      IF (ierr==NF90_NOERR) THEN
+        CALL body(field_glo)
+        tmp_found=.TRUE.
+      ELSE
+        tmp_found=.FALSE.
+      ENDIF
+    
+    ENDIF
+    
+    CALL bcast(tmp_found)
+
+    IF (tmp_found) THEN
+      CALL scatter(field_glo,field)
+    ENDIF
+    
+    IF (PRESENT(found)) THEN
+      found=tmp_found
+    ELSE
+      IF (.NOT. tmp_found) THEN
+        PRINT*, 'get_field_rgen: Field <'//field_name//'> not found'
+        CALL abort_physic("get_field_rgen","Failed to get field",1)
+      ENDIF
+    ENDIF
+ 
+    
+    CONTAINS
+     
+     SUBROUTINE body(field_glo)
+       REAL :: field_glo(klon_glo*field_size)
+         ierr=NF90_GET_VAR(nid_start,varid,field_glo,corners,edges)
+         IF (ierr/=NF90_NOERR) THEN
+           ! La variable exist dans le fichier mais la lecture a echouee. 
+           PRINT*, 'get_field_rgen: Failed reading <'//field_name//'>'
+
+!           IF (field_name=='CLWCON' .OR. field_name=='RNEBCON' .OR. field_name=='RATQS') THEN
+!              ! Essaye de lire le variable sur surface uniqument, comme fait avant
+!              field_glo(:)=0.
+!              ierr=NF90_GET_VAR(nid_start,varid,field_glo(1:klon_glo))
+!              IF (ierr/=NF90_NOERR) THEN
+!                 PRINT*, 'phyetat0: Lecture echouee aussi en 2D pour <'//field_name//'>'
+!                 CALL abort
+!              ELSE
+!                 PRINT*, 'phyetat0: La variable <'//field_name//'> lu sur surface seulement'!, selon ancien format, le reste mis a zero'
+!              END IF
+!           ELSE
+              CALL abort_physic("get_field_rgen","Failed to read field",1)
+!           ENDIF
+         ENDIF
+
+     END SUBROUTINE body
+
+  END SUBROUTINE Get_field_rgen
+
+
+  SUBROUTINE get_var_r0(var_name,var,found)
+  ! Get a scalar from input file
+  IMPLICIT NONE  
+    CHARACTER(LEN=*),INTENT(IN)  :: var_name
+    REAL,INTENT(INOUT)             :: var
+    LOGICAL,OPTIONAL,INTENT(OUT) :: found
+
+    REAL                         :: varout(1)
+    
+    IF (PRESENT(found)) THEN
+      CALL Get_var_rgen(var_name,varout,size(varout),found)
+    ELSE
+      CALL Get_var_rgen(var_name,varout,size(varout))
+    ENDIF
+    var=varout(1)
+ 
+  END SUBROUTINE get_var_r0
+
+  SUBROUTINE get_var_r1(var_name,var,found)
+  ! Get a vector from input file
+  IMPLICIT NONE  
+    CHARACTER(LEN=*),INTENT(IN)  :: var_name
+    REAL,INTENT(INOUT)             :: var(:)
+    LOGICAL,OPTIONAL,INTENT(OUT) :: found
+    
+    IF (PRESENT(found)) THEN
+      CALL Get_var_rgen(var_name,var,size(var),found)
+    ELSE
+      CALL Get_var_rgen(var_name,var,size(var))
+    ENDIF
+  
+  END SUBROUTINE get_var_r1
+
+  SUBROUTINE get_var_r2(var_name,var,found)
+  ! Get a 2D field from input file
+  IMPLICIT NONE  
+    CHARACTER(LEN=*),INTENT(IN)  :: var_name
+    REAL,INTENT(OUT)             :: var(:,:)
+    LOGICAL,OPTIONAL,INTENT(OUT) :: found
+    
+    IF (PRESENT(found)) THEN
+      CALL Get_var_rgen(var_name,var,size(var),found)
+    ELSE
+      CALL Get_var_rgen(var_name,var,size(var))
+    ENDIF
+  
+  END SUBROUTINE get_var_r2
+
+  SUBROUTINE get_var_r3(var_name,var,found)
+  ! Get a 3D field frominput file
+  IMPLICIT NONE  
+    CHARACTER(LEN=*),INTENT(IN)  :: var_name
+    REAL,INTENT(INOUT)             :: var(:,:,:)
+    LOGICAL,OPTIONAL,INTENT(OUT) :: found
+    
+    IF (PRESENT(found)) THEN
+      CALL Get_var_rgen(var_name,var,size(var),found)
+    ELSE
+      CALL Get_var_rgen(var_name,var,size(var))
+    ENDIF
+  
+  END SUBROUTINE get_var_r3
+
+  SUBROUTINE Get_var_rgen(var_name,var,var_size,found)
+  USE netcdf
+  USE dimphy
+  USE mod_grid_phy_lmdz
+  USE mod_phys_lmdz_para
+  IMPLICIT NONE
+    CHARACTER(LEN=*) :: var_name
+    INTEGER          :: var_size
+    REAL             :: var(var_size)
+    LOGICAL,OPTIONAL :: found
+    
+    LOGICAL :: tmp_found
+    INTEGER :: varid
+    INTEGER :: ierr
+    
+    IF (is_mpi_root .AND. is_omp_root) THEN
+  
+      ierr=NF90_INQ_VARID(nid_start,var_name,varid)
+      
+      IF (ierr==NF90_NOERR) THEN
+        ierr=NF90_GET_VAR(nid_start,varid,var)
+        IF (ierr/=NF90_NOERR) THEN
+          PRINT*, 'phyetat0: Failed loading <'//trim(var_name)//'>'
+          CALL abort_physic("get_var_rgen","Failed to read variable",1)
+        ENDIF
+        tmp_found=.TRUE.
+      ELSE
+        tmp_found=.FALSE.
+      ENDIF
+    
+    ENDIF
+    
+    CALL bcast(tmp_found)
+
+    IF (tmp_found) THEN
+      CALL bcast(var)
+    ENDIF
+    
+    IF (PRESENT(found)) THEN
+      found=tmp_found
+    ELSE
+      IF (.NOT. tmp_found) THEN
+        PRINT*, 'phyetat0: Variable <'//trim(var_name)//'> not found'
+        CALL abort_physic("get_var_rgen","Failed to read variable",1)
+      ENDIF
+    ENDIF
+
+  END SUBROUTINE Get_var_rgen
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  SUBROUTINE open_restartphy(filename)
+  USE netcdf, only: NF90_CREATE, NF90_CLOBBER, NF90_64BIT_OFFSET, &
+                    NF90_NOERR, nf90_strerror, &
+                    NF90_PUT_ATT, NF90_GLOBAL, NF90_DEF_DIM, &
+                    NF90_UNLIMITED, NF90_ENDDEF, &
+                    NF90_WRITE, NF90_OPEN
+  USE mod_phys_lmdz_para, only: is_master
+  USE mod_grid_phy_lmdz, only: klon_glo
+  USE dimphy, only: klev, klevp1
+  USE tracer_mod, only: nqmx
+  USE comsoil_h_PEM, only:  nsoilmx_PEM
+  USE comslope_mod, only: nslope
+  IMPLICIT NONE
+    CHARACTER(LEN=*),INTENT(IN) :: filename
+    INTEGER                     :: ierr
+    LOGICAL,SAVE :: already_created=.false.
+    
+    IF (is_master) THEN
+      if (.not.already_created) then
+        ! At the very first call, create the file
+        ierr=NF90_CREATE(filename,IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), &
+                          nid_restart)
+        IF (ierr/=NF90_NOERR) THEN
+          write(*,*)'open_restartphy: problem creating file '//trim(filename)
+          write(*,*)trim(nf90_strerror(ierr))
+          CALL abort_physic("open_restartphy","Failed creating file",1)
+        ENDIF
+        already_created=.true.
+      else
+        ! Just open the file
+        ierr=NF90_OPEN(filename,NF90_WRITE,nid_restart)
+        IF (ierr/=NF90_NOERR) THEN
+          write(*,*)'open_restartphy: problem opening file '//trim(filename)
+          write(*,*)trim(nf90_strerror(ierr))
+          CALL abort_physic("open_restartphy","Failed opening file",1)
+        ENDIF
+        return
+      endif ! of if (.not.already_created)
+
+      ierr=NF90_PUT_ATT(nid_restart,NF90_GLOBAL,"title",&
+                        "Physics start file")
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'open_restartphy: problem writing title '
+        write(*,*)trim(nf90_strerror(ierr))
+      ENDIF
+
+      ierr=NF90_DEF_DIM(nid_restart,"index",length,idim1)
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'open_restartphy: problem defining index dimension '
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("open_restartphy","Failed defining index",1)
+      ENDIF
+      
+      ierr=NF90_DEF_DIM(nid_restart,"physical_points",klon_glo,idim2)
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'open_restartphy: problem defining physical_points dimension '
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("open_restartphy","Failed defining physical_points",1)
+      ENDIF
+      
+      ierr=NF90_DEF_DIM(nid_restart,"subsurface_layers",nsoilmx_PEM,idim3)
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'open_restartphy: problem defining subsurface_layers dimension '
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("open_restartphy","Failed defining subsurface_layers",1)
+      ENDIF
+      
+      ierr=NF90_DEF_DIM(nid_restart,"nlayer_plus_1",klevp1,idim4)
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'open_restartphy: problem defining nlayer_plus_1 dimension '
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("open_restartphy","Failed defining nlayer_plus_1",1)
+      ENDIF
+      
+      if (nqmx.ne.0) then
+        ierr=NF90_DEF_DIM(nid_restart,"number_of_advected_fields",nqmx,idim5)
+      else
+        ! pretend nqmx=1 because 0 size implies unlimited dimension for NetCDF
+        ierr=NF90_DEF_DIM(nid_restart,"number_of_advected_fields",1,idim5)
+      endif
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'open_restartphy: problem defining number_of_advected_fields dimension '
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("open_restartphy","Failed defining number_of_advected_fields",1)
+      ENDIF
+
+      ierr=NF90_DEF_DIM(nid_restart,"nlayer",klev,idim6)
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'open_restartphy: problem defining nlayer dimension '
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("open_restartphy","Failed defining nlayer",1)
+      ENDIF
+      
+      ierr=NF90_DEF_DIM(nid_restart,"Time",NF90_UNLIMITED,idim7)
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'open_restartphy: problem defining Time dimension '
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("open_restartphy","Failed defining Time",1)
+      ENDIF
+
+      ierr=NF90_DEF_DIM(nid_restart,"nslope",nslope,idim8)
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'phyredem: problem defining nslope dimension'
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL ABORT
+      ENDIF
+
+      ierr=NF90_DEF_DIM(nid_restart,"inter slope",nslope+1,idim9)
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'phyredem: problem defining inter slope dimension'
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL ABORT
+      ENDIF
+
+
+      ierr=NF90_ENDDEF(nid_restart)
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'open_restartphy: problem ending definition mode '
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("open_restartphy","Failed ending definition mode",1)
+      ENDIF
+    ENDIF
+
+  END SUBROUTINE open_restartphy
+
+  SUBROUTINE close_restartphy
+  USE netcdf, only: NF90_CLOSE
+  USE mod_phys_lmdz_para, only: is_master
+  IMPLICIT NONE
+    INTEGER          :: ierr
+
+    IF (is_master) THEN
+      ierr = NF90_CLOSE (nid_restart)
+    ENDIF
+ 
+  END SUBROUTINE close_restartphy
+
+  SUBROUTINE put_field_r1(field_name,title,field,time)
+  ! For a surface field
+  IMPLICIT NONE
+  CHARACTER(LEN=*),INTENT(IN)    :: field_name
+  CHARACTER(LEN=*),INTENT(IN)    :: title
+  REAL,INTENT(IN)                :: field(:)
+  REAL,OPTIONAL,INTENT(IN)       :: time
+  
+  IF (present(time)) THEN
+    ! if timeindex is present, it is a time-dependent variable
+    CALL put_field_rgen(field_name,title,field,1,time)
+  ELSE
+    CALL put_field_rgen(field_name,title,field,1)
+  ENDIF
+  
+  END SUBROUTINE put_field_r1
+
+  SUBROUTINE put_field_r2(field_name,title,field,time)
+  ! For a "3D" horizontal-vertical field
+  IMPLICIT NONE
+  CHARACTER(LEN=*),INTENT(IN)    :: field_name
+  CHARACTER(LEN=*),INTENT(IN)    :: title
+  REAL,INTENT(IN)                :: field(:,:)
+  REAL,OPTIONAL,INTENT(IN)       :: time
+  
+  IF (present(time)) THEN
+    ! if timeindex is present, it is a time-dependent variable
+    CALL put_field_rgen(field_name,title,field,size(field,2),time)
+  ELSE
+    CALL put_field_rgen(field_name,title,field,size(field,2))
+  ENDIF
+  
+  END SUBROUTINE put_field_r2
+
+  SUBROUTINE put_field_r3(field_name,title,field,time)
+  ! For a "4D" field surf/alt/??
+  IMPLICIT NONE
+  CHARACTER(LEN=*),INTENT(IN)    :: field_name
+  CHARACTER(LEN=*),INTENT(IN)    :: title
+  REAL,INTENT(IN)                :: field(:,:,:)
+  REAL,OPTIONAL,INTENT(IN)       :: time
+  
+  IF (present(time)) THEN
+    ! if timeindex is present, it is a time-dependent variable
+    CALL put_field_rgen(field_name,title,field,size(field,2)*size(field,3),&
+                        time)
+  ELSE  
+    CALL put_field_rgen(field_name,title,field,size(field,2)*size(field,3))
+  ENDIF
+  
+  END SUBROUTINE put_field_r3
+  
+  SUBROUTINE put_field_rgen(field_name,title,field,field_size,time)
+  USE netcdf
+  USE dimphy
+  USE comsoil_h_PEM, only:  nsoilmx_PEM
+  USE comslope_mod, ONLY: nslope
+  USE mod_grid_phy_lmdz
+  USE mod_phys_lmdz_para
+  IMPLICIT NONE
+  CHARACTER(LEN=*),INTENT(IN)    :: field_name
+  CHARACTER(LEN=*),INTENT(IN)    :: title
+  INTEGER,INTENT(IN)             :: field_size
+  REAL,INTENT(IN)                :: field(klon,field_size)
+  REAL,OPTIONAL,INTENT(IN)       :: time
+  
+  REAL                           :: field_glo(klon_glo,field_size)
+  INTEGER                        :: ierr
+  INTEGER                        :: nvarid
+  INTEGER                        :: idim
+   
+    CALL gather(field,field_glo)
+    
+    IF (is_master) THEN
+
+      IF (field_size==1) THEN
+        ! input is a 1D "surface field" array
+        if (.not.present(time)) then ! for a time-independent field
+          ierr=NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                            (/idim2/),nvarid)
+#else
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                            (/idim2/),nvarid)
+#endif
+          if (ierr.ne.NF90_NOERR) then
+            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
+            write(*,*)trim(nf90_strerror(ierr))
+          endif
+          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+          ierr=NF90_ENDDEF(nid_restart)
+          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo)
+        else
+          ! check if the variable has already been defined:
+          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
+          if (ierr/=NF90_NOERR) then ! variable not found, define it
+            ierr=NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                              (/idim2,idim7/),nvarid)
+#else
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                              (/idim2,idim7/),nvarid)
+#endif
+            if (ierr.ne.NF90_NOERR) then
+              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
+              write(*,*)trim(nf90_strerror(ierr))
+            endif
+            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+            ierr=NF90_ENDDEF(nid_restart)
+          endif
+          ! Write the variable
+          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
+                            start=(/1,timeindex/))
+        endif ! of if (.not.present(timeindex))
+
+      ELSE IF (field_size==klev) THEN
+        ! input is a 2D "atmospheric field" array
+        if (.not.present(time)) then ! for a time-independent field
+          ierr=NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                            (/idim2,idim6/),nvarid)
+#else
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                            (/idim2,idim6/),nvarid)
+#endif
+          if (ierr.ne.NF90_NOERR) then
+            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
+            write(*,*)trim(nf90_strerror(ierr))
+          endif
+          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+          ierr=NF90_ENDDEF(nid_restart)
+          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo)
+        else
+          ! check if the variable has already been defined:
+          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
+          if (ierr/=NF90_NOERR) then ! variable not found, define it
+            ierr=NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                              (/idim2,idim6,idim7/),nvarid)
+#else
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                              (/idim2,idim6,idim7/),nvarid)
+#endif
+            if (ierr.ne.NF90_NOERR) then
+              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
+              write(*,*)trim(nf90_strerror(ierr))
+            endif
+            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+            ierr=NF90_ENDDEF(nid_restart)
+          endif
+          ! Write the variable
+          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
+                            start=(/1,1,timeindex/))
+        endif ! of if (.not.present(time))
+
+      ELSE IF (field_size==klevp1) THEN
+        ! input is a 2D "interlayer atmospheric field" array
+        if (.not.present(time)) then ! for a time-independent field
+          ierr=NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                            (/idim2,idim4/),nvarid)
+#else
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                            (/idim2,idim4/),nvarid)
+#endif
+          if (ierr.ne.NF90_NOERR) then
+            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
+            write(*,*)trim(nf90_strerror(ierr))
+          endif
+          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+          ierr = NF90_ENDDEF(nid_restart)
+          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
+        else
+          ! check if the variable has already been defined:
+          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
+          if (ierr/=NF90_NOERR) then ! variable not found, define it
+            ierr=NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                              (/idim2,idim4,idim7/),nvarid)
+#else
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                              (/idim2,idim4,idim7/),nvarid)
+#endif
+            if (ierr.ne.NF90_NOERR) then
+              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
+              write(*,*)trim(nf90_strerror(ierr))
+            endif
+            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+            ierr=NF90_ENDDEF(nid_restart)
+          endif
+          ! Write the variable
+          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
+                            start=(/1,1,timeindex/))
+        endif ! of if (.not.present(timeindex))
+
+      ELSE IF (field_size==nsoilmx_PEM) THEN
+        ! input is a 2D "subsurface field" array
+        if (.not.present(time)) then ! for a time-independent field
+          ierr = NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                            (/idim2,idim3/),nvarid)
+#else
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                            (/idim2,idim3/),nvarid)
+#endif
+          if (ierr.ne.NF90_NOERR) then
+            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
+            write(*,*)trim(nf90_strerror(ierr))
+          endif
+          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+          ierr = NF90_ENDDEF(nid_restart)
+          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
+        else
+          ! check if the variable has already been defined:
+          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
+          if (ierr/=NF90_NOERR) then ! variable not found, define it
+            ierr=NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                              (/idim2,idim3,idim7/),nvarid)
+#else
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                              (/idim2,idim3,idim7/),nvarid)
+#endif
+           if (ierr.ne.NF90_NOERR) then
+              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
+              write(*,*)trim(nf90_strerror(ierr))
+            endif
+            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+            ierr=NF90_ENDDEF(nid_restart)
+          endif
+          ! Write the variable
+          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
+                            start=(/1,1,timeindex/))
+
+        endif ! of if (.not.present(time))
+
+
+      ELSE IF (field_size==nslope) THEN
+        ! input is a 2D "subsurface field" array
+        if (.not.present(time)) then ! for a time-independent field
+          ierr = NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                            (/idim2,idim8/),nvarid)
+#else
+          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                            (/idim2,idim8/),nvarid)
+#endif
+          if (ierr.ne.NF90_NOERR) then
+            write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
+            write(*,*)trim(nf90_strerror(ierr))
+          endif
+          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+          ierr = NF90_ENDDEF(nid_restart)
+          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
+        else
+          ! check if the variable has already been defined:
+          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
+          if (ierr/=NF90_NOERR) then ! variable not found, define it
+            ierr=NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                              (/idim2,idim8,idim7/),nvarid)
+#else
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                              (/idim2,idim8,idim7/),nvarid)
+#endif
+           if (ierr.ne.NF90_NOERR) then
+              write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
+              write(*,*)trim(nf90_strerror(ierr))
+            endif
+            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+            ierr=NF90_ENDDEF(nid_restart)
+          endif
+          ! Write the variable
+          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
+                            start=(/1,1,timeindex/))
+
+        endif ! of if (.not.present(time))
+
+      ELSE
+        PRINT *, "Error phyredem(put_field_rgen) : wrong dimension for ",trim(field_name)
+        write(*,*) "  field_size =",field_size
+        CALL abort_physic("put_field_rgen","wrong field dimension",1)
+      ENDIF
+
+      ! Check the writting of field to file went OK
+      if (ierr.ne.NF90_NOERR) then
+        write(*,*) " Error phyredem(put_field_rgen) : failed writing ",trim(field_name)
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("put_field_rgen","Failed writing field",1)
+      endif
+
+    ENDIF ! of IF (is_master)
+    
+  END SUBROUTINE put_field_rgen  
+  
+  SUBROUTINE put_var_r0(var_name,title,var)
+  ! Put a scalar in file
+   IMPLICIT NONE
+     CHARACTER(LEN=*),INTENT(IN) :: var_name
+     CHARACTER(LEN=*),INTENT(IN) :: title
+     REAL,INTENT(IN)             :: var
+     REAL                        :: varin(1)
+     
+     varin(1)=var
+     
+     CALL put_var_rgen(var_name,title,varin,size(varin))
+
+  END SUBROUTINE put_var_r0
+
+
+  SUBROUTINE put_var_r1(var_name,title,var)
+  ! Put a vector in file
+   IMPLICIT NONE
+     CHARACTER(LEN=*),INTENT(IN) :: var_name
+     CHARACTER(LEN=*),INTENT(IN) :: title
+     REAL,INTENT(IN)             :: var(:)
+     
+     CALL put_var_rgen(var_name,title,var,size(var))
+
+  END SUBROUTINE put_var_r1
+ 
+  SUBROUTINE put_var_r2(var_name,title,var)
+  ! Put a 2D field in file
+   IMPLICIT NONE
+     CHARACTER(LEN=*),INTENT(IN) :: var_name
+     CHARACTER(LEN=*),INTENT(IN) :: title
+     REAL,INTENT(IN)             :: var(:,:)
+     
+     CALL put_var_rgen(var_name,title,var,size(var))
+
+  END SUBROUTINE put_var_r2     
+  
+  SUBROUTINE put_var_r3(var_name,title,var)
+  ! Put a 3D field in file
+   IMPLICIT NONE
+     CHARACTER(LEN=*),INTENT(IN) :: var_name
+     CHARACTER(LEN=*),INTENT(IN) :: title
+     REAL,INTENT(IN)             :: var(:,:,:)
+     
+     CALL put_var_rgen(var_name,title,var,size(var))
+
+  END SUBROUTINE put_var_r3
+
+  SUBROUTINE put_var_rgen(var_name,title,var,var_size)
+  USE netcdf, only: NF90_REDEF, NF90_DEF_VAR, NF90_ENDDEF, NF90_PUT_VAR, &
+                    NF90_FLOAT, NF90_DOUBLE, &
+                    NF90_PUT_ATT, NF90_NOERR, nf90_strerror, &
+                    nf90_inq_dimid, nf90_inquire_dimension, NF90_INQ_VARID
+  USE comsoil_h_PEM, only:  nsoilmx_PEM
+  USE comslope_mod, only: nslope
+  USE mod_phys_lmdz_para, only: is_master
+  IMPLICIT NONE
+     CHARACTER(LEN=*),INTENT(IN) :: var_name
+     CHARACTER(LEN=*),INTENT(IN) :: title
+     INTEGER,INTENT(IN)          :: var_size
+     REAL,INTENT(IN)             :: var(var_size)
+     
+     INTEGER :: ierr
+     INTEGER :: nvarid
+     INTEGER :: idim1d
+     logical,save :: firsttime=.true.
+         
+    IF (is_master) THEN
+
+      IF (var_name=="Time") THEN
+        ! Very specific case of "Time" variable
+        if (firsttime) then
+          ! Create the "Time variable"
+          ierr=NF90_REDEF(nid_restart)
+#ifdef NC_DOUBLE
+          ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_DOUBLE,&
+                            (/idim7/),nvarid)
+#else
+          ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_FLOAT,&
+                            (/idim7/),nvarid)
+#endif
+          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+          ierr=NF90_ENDDEF(nid_restart)
+          
+          firsttime=.false.
+        endif
+        ! Append "time" value
+        ! get current length of "Time" dimension
+        ierr=nf90_inq_dimid(nid_restart,var_name,idim1d)
+        ierr=nf90_inquire_dimension(nid_restart,idim1d,len=timeindex)
+        timeindex=timeindex+1
+        ierr=NF90_INQ_VARID(nid_restart,var_name,nvarid)
+        ierr=NF90_PUT_VAR(nid_restart,nvarid,var,&
+                           start=(/timeindex/))
+        IF (ierr/=NF90_NOERR) THEN
+          write(*,*)'put_var_rgen: problem writing Time'
+          write(*,*)trim(nf90_strerror(ierr))
+          CALL abort_physic("get_var_rgen","Failed to write Time",1)
+        ENDIF
+        return ! nothing left to do
+      ELSEIF (var_size==length) THEN
+        ! We know it is a "controle" kind of 1D array
+        idim1d=idim1
+      ELSEIF (var_size==nsoilmx_PEM) THEN
+        ! We know it is an  "mlayer" kind of 1D array
+        idim1d=idim3
+      ELSEIF (var_size==nslope+1) THEN
+        ! We know it is an  "inter slope" kind of 1D array
+        idim1d=idim9
+      ELSE 
+        PRINT *, "put_var_rgen error : wrong dimension"
+        write(*,*) "  var_size =",var_size
+        CALL abort_physic("get_var_rgen","Wrong variable dimension",1)
+
+      ENDIF ! of IF (var_size==length) THEN
+
+      ! Swich to NetCDF define mode
+      ierr=NF90_REDEF (nid_restart)
+      ! Define the variable
+#ifdef NC_DOUBLE
+      ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_DOUBLE,(/idim1d/),nvarid)
+#else
+      ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_FLOAT,(/idim1d/),nvarid)
+#endif
+      ! Add a "title" attribute
+      IF (LEN_TRIM(title)>0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
+      ! Swich out of define mode
+      ierr=NF90_ENDDEF(nid_restart)
+      ! Write variable to file
+      ierr=NF90_PUT_VAR(nid_restart,nvarid,var)
+      IF (ierr/=NF90_NOERR) THEN
+        write(*,*)'put_var_rgen: problem writing '//trim(var_name)
+        write(*,*)trim(nf90_strerror(ierr))
+        CALL abort_physic("get_var_rgen","Failed writing variable",1)
+      ENDIF
+    ENDIF ! of IF (is_master)
+    
+  END SUBROUTINE put_var_rgen     
+
+END MODULE iostart_PEM
Index: trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90	(revision 2794)
@@ -0,0 +1,75 @@
+!
+! $Id $
+!
+SUBROUTINE nb_time_step_GCM(fichnom,timelen)
+
+      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
+                        nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
+                        nf90_inquire_dimension,nf90_close
+
+      IMPLICIT NONE
+
+!=======================================================================
+!
+! Read initial confitions file
+!
+!=======================================================================
+
+  include "dimensions.h"
+
+!===============================================================================
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
+!===============================================================================
+!   Local Variables 
+  CHARACTER(LEN=256) :: msg, var, modname
+  INTEGER :: iq, fID, vID, idecal
+  INTEGER :: ierr
+
+  INTEGER :: timelen ! number of times stored in the file
+!-----------------------------------------------------------------------
+  modname="nb_time_step_GCM"
+
+!  Open initial state NetCDF file
+  var=fichnom
+  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
+
+      ierr = nf90_inq_varid (fID, "temps", vID)
+      IF (ierr .NE. nf90_noerr) THEN
+        write(*,*)"read_data_GCM: Le champ <temps> est absent"
+        write(*,*)"read_data_GCM: J essaie <Time>"
+        ierr = nf90_inq_varid (fID, "Time", vID)
+        IF (ierr .NE. nf90_noerr) THEN
+           write(*,*)"read_data_GCM: Le champ <Time> est absent"
+           write(*,*)trim(nf90_strerror(ierr))
+           CALL ABORT_gcm("read_data_GCM", "", 1)
+        ENDIF
+        ! Get the length of the "Time" dimension
+        ierr = nf90_inq_dimid(fID,"Time",vID)
+        ierr = nf90_inquire_dimension(fID,vID,len=timelen)
+      ELSE   
+        ! Get the length of the "temps" dimension
+        ierr = nf90_inq_dimid(fID,"temps",vID)
+        ierr = nf90_inquire_dimension(fID,vID,len=timelen)
+      ENDIF
+
+  CALL err(NF90_CLOSE(fID),"close",fichnom)
+
+  CONTAINS
+
+
+SUBROUTINE err(ierr,typ,nam)
+  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
+  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
+  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
+  IF(ierr==NF90_NoERR) RETURN
+  SELECT CASE(typ)
+    CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
+    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
+    CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
+    CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
+  END SELECT
+  CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
+END SUBROUTINE err
+
+END SUBROUTINE nb_time_step_GCM
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 2792)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 2794)
@@ -4,13 +4,18 @@
 ! I   Initialisation
 !    I_a READ run.def , run_pem.def
-!    I_b READ starfi_0.nc and start_0.nc
+!    I_b READ starfi_0.nc
+!    I_c READ GCM data and convert to the physical grid 
+!    I_d Compute tendencies & Save initial situation
 
 ! II  Run
-!    II_a Compute tendencies
-!    II_b Save initial situation
-!    II_c Time loop
+!    II_a update pressure, ice and tracers 
+!    II_b CO2 glaciers flows 
+!    II_c Update surface and soil temperatures
+!    II_d Update the tendencies and stopping criterion
 
 ! III Output
-!    III_a Write startfi.nc
+!    III_a Recomp tendencies for the start and startfi
+!    III_b Write start and starfi.nc
+!    III_c Write start_pem
 
 !------------------------
@@ -21,18 +26,21 @@
 !module needed for INITIALISATION
       use phyetat0_mod, only: phyetat0
-      use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer
+      use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa
       use surfdat_h, only: tsurf, co2ice, emis,&
       &                    qsurf,watercap, ini_surfdat_h, &
                            albedodat, zmea, zstd, zsig, zgam, zthe, &
-      &                    hmons, summit, base
+                           hmons, summit, base,albedo_h2o_frost, &
+                          frost_albedo_threshold,emissiv
       use dimradmars_mod, only: totcloudfrac, albedo, & 
                                 ini_dimradmars_mod
       use turb_mod, only: q2, wstar, ini_turb_mod
       use dust_param_mod, only: tauscaling, ini_dust_param_mod
+!      use co2cloud_mod, only: mem_Mccn_co2, mem_Mh2o_co2,&
+!      &                        mem_Nccn_co2,ini_co2cloud
       use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
                         nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
                         nf90_inquire_dimension,nf90_close
       use phyredem, only: physdem0, physdem1
-      use tracer_mod, only: noms,nqmx ! tracer names
+      use tracer_mod, only: noms,nqmx,igcm_h2o_ice ! tracer names
 
 ! For phyredem :
@@ -49,7 +57,8 @@
       USE infotrac
       USE temps_mod_evol, ONLY: nyear, dt_pem
-      USE vertical_layers_mod, ONLY: ap,bp
-
-      use comslope_mod, ONLY: nslope,def_slope,def_slope_mean, &
+!     USE vertical_layers_mod, ONLY: ap,bp
+      USE comcstfi_h, only: pi 
+
+      USE comslope_mod, ONLY: nslope,def_slope,def_slope_mean, &
                            subslope_dist,co2ice_slope, &
                            tsurf_slope,tsoil_slope,fluxgrd_slope,&
@@ -59,26 +68,41 @@
                            iflat
 
-
+      USE geometry_mod, only: longitude_deg,latitude_deg 
+
+      use pemredem, only:  pemdem1
+      USE soil_evolution_mod, ONLY: soil_pem_CN
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
+
+      use comsoil_h_PEM, only: ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
+                              TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM, &        ! soil thermal inertia          
+                              tsoil_PEM ,       &        !number of subsurface layers
+                              mlayer_PEM,layer_PEM, &       ! soil mid layer depths
+                              fluxgeo, &  ! geothermal flux
+                              co2_adsorbded_phys ! mass of co2 in the regolith
+      use adsorption_mod, only : regolith_co2adsorption
   IMPLICIT NONE
 
   include "dimensions.h"
   include "paramet.h"
+
+  INTEGER ngridmx
+  PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
+
   include "comdissnew.h"
   include "comgeom.h"
   include "iniprint.h"
-
 ! Same variable's name as in the GCM
 
-      INTEGER :: ngrid      !Number of physical grid points
-      INTEGER :: nlayer     !Number of vertical layer
-      INTEGER :: nq         !Number of tracer
-      INTEGER :: day_ini    !First day of the simulation
-      REAL :: pday          !Physical day
-      REAL :: time_phys     !Same as GCM
-      REAL :: ptimestep     !Same as GCM
-      REAL :: ztime_fin     !Same as GCM
+  INTEGER :: ngrid      !Number of physical grid points
+  INTEGER :: nlayer     !Number of vertical layer
+  INTEGER :: nq         !Number of tracer
+  INTEGER :: day_ini    !First day of the simulation
+  REAL :: pday          !Physical day
+  REAL :: time_phys     !Same as GCM
+  REAL :: ptimestep     !Same as GCM
+  REAL :: ztime_fin     !Same as GCM
 
 ! Variable for reading start.nc
-      character (len = *), parameter :: FILE_NAME_start = "start_evol.nc" !Name of the file used for initialsing the PEM
+      character (len = *), parameter :: FILE_NAME_start = "start_0.nc" !Name of the file used for initialsing the PEM
   !   variables dynamiques
   REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
@@ -86,5 +110,8 @@
   REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
   REAL ps(ip1jmp1)                       ! pression  au sol
-  REAL, dimension(:),allocatable :: ps_phys !(ngrid) pression  au sol
+  REAL, dimension(:),allocatable :: ps_phys !(ngrid)                       ! pression  au sol
+  REAL, dimension(:,:),allocatable :: ps_phys_timeseries !(ngrid x timelen) ! pression  au sol instantannées
+  REAL, dimension(:,:),allocatable :: ps_phys_timeseries_yr1 !(ngrid x timelen) ! pression  au sol instantannées for the first year of the gcm
+
   REAL masse(ip1jmp1,llm)                ! masse d'air
   REAL phis(ip1jmp1)                     ! geopotentiel au sol
@@ -93,8 +120,9 @@
 ! Variable for reading starfi.nc
 
-      character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM
+      character (len = *), parameter :: FILE_NAME = "startfi_0.nc" !Name of the file used for initialsing the PEM
       integer :: ncid, varid,status                                !Variable for handling opening of files
       integer :: phydimid, subdimid, nlayerdimid, nqdimid          !Variable ID for Netcdf files
       integer :: lonvarid, latvarid, areavarid,sdvarid             !Variable ID for Netcdf files
+      integer :: apvarid,bpvarid				   !Variable ID for Netcdf files
 
 ! Variable for reading starfi.nc and writting restartfi.nc
@@ -102,4 +130,6 @@
       REAL, dimension(:),allocatable :: longitude                  !Longitude read in FILE_NAME and written in restartfi 
       REAL, dimension(:),allocatable  :: latitude                  !Latitude read in FILE_NAME and written in restartfi 
+      REAL, dimension(:),allocatable :: ap                         !Coefficient ap read in FILE_NAME_start and written in restart 
+      REAL, dimension(:),allocatable :: bp                         !Coefficient bp read in FILE_NAME_start and written in restart
       REAL, dimension(:),allocatable  :: cell_area                 !Cell_area read in FILE_NAME and written in restartfi 
       REAL :: Total_surface                                        !Total surface of the planet
@@ -113,16 +143,10 @@
       REAL, dimension(:),allocatable  :: tendencies_co2_ice_phys   ! physical point field : Tendency of evolution of perenial co2 ice over a year
 
-! Variable for initial situation of water ice
-
       REAL :: ini_surf                                             ! Initial surface of sublimating water ice
-      REAL :: ini_surf_h2o                                         ! Initial surface of sublimating water ice
+      REAL :: ini_surf_h2o                                             ! Initial surface of sublimating water ice
       REAL, dimension(:),allocatable  :: initial_h2o_ice           ! physical point field : Logical array indicating sublimating point
-
-! Variable for initial situation of co2 ice
 
       REAL :: ini_surf_co2                                         ! Initial surface of sublimating co2 ice
       REAL, dimension(:),allocatable  :: initial_co2_ice           ! physical point field : Logical array indicating sublimating point of co2 ice
-
-! Variable for needed to compute tendencies for water ice
 
       REAL , dimension(:,:), allocatable :: min_h2o_ice_s_1        ! LON x LAT field : minimum of water ice at each point for the first year
@@ -130,69 +154,149 @@
       REAL , dimension(:,:,:), allocatable :: h2o_ice_first_last_day     ! LON x LAT x 2 field : First and Last value of a GCM year simulation of water ice
 
-! Variable for needed to compute tendencies for co2 ice
-
       REAL , dimension(:,:), allocatable :: min_co2_ice_s_1        ! LON x LAT field : minimum of water ice at each point for the first year
       REAL , dimension(:,:), allocatable :: min_co2_ice_s_2        ! LON x LAT field : minimum of water ice at each point for the second year
       REAL , dimension(:,:,:), allocatable :: co2_ice_first_last_day     ! LON x LAT x 2 field : First and Last value of a GCM year simulation of water ice
 
-      REAL :: global_ave_press_GCM           ! physical point field : Global average pressure in the GCM output
+      REAL, dimension(:),allocatable  :: local_old_press           ! physical point field : Local pressure of initial/previous time step
+      REAL, dimension(:),allocatable  :: local_new_press           ! physical point field : Local pressure of current time step
+
+      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 ! ngrid x nlayer+1 : Pressure level at the end of the iteration
-      REAL , dimension(:,:), allocatable :: ZPLEV_OLD ! ngrid x nlayer+1 : Pressure level at the beginning of the iteration
-
-      INTEGER :: year_iter        !Counter for the number of PEM iteration
-      INTEGER :: year             ! Year read in the startfi files and written back in it counting the year iterations
-
+      REAL , dimension(:,:), allocatable ::  zplev_new
+      REAL , dimension(:,:), allocatable :: zplev_old
+      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 :: tot_co2_atm,tot_var_co2_atm
+
+
+      INTEGER :: year_iter  !Counter for the number of PEM iteration
       LOGICAL :: STOPPING_water   ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
       LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate?
-
-      LOGICAL :: STOPPING_co2     ! Logical : is the criterion (% of change in the surface of sublimating co2 ice) reached?
-      LOGICAL :: STOPPING_1_co2   ! Logical : is there still co2 ice to sublimate?
-
-      REAL, dimension(:),allocatable  :: q_co2_GCM ! Initial amount of co2 in the first layer
-      REAL ps_GCM(ip1jmp1)                         ! pression  au sol donné par le GCM
-
-      real ,allocatable :: vmr_co2_gcm_phys(:,:)   !(ngrid x Timelen) : co2 volume mixing ratio given by the GCM 
-      REAL, ALLOCATABLE ::  vmr_co2_gcm(:,:,:)     ! iim+1 x jjm+1 x Timelen : co2 volume mixing ratio given by the GCM
-
-      REAL, ALLOCATABLE ::  ps_GCM_1(:,:,:)        ! iim+1 x jjm+1 x Timelen : Surface pressure in the first GCM year
-      REAL, ALLOCATABLE ::  ps_GCM_2(:,:,:)        ! iim+1 x jjm+1 x Timelen : Surface pressure in the second GCM year
-
-      integer :: timelen                           ! Number of time point in the GCM diagfi output
-
-      REAL, ALLOCATABLE :: p(:,:)                  !(ngrid,llmp1) : Pressure
+      LOGICAL :: STOPPING_co2   ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
+      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(:,:)
+      real ,allocatable :: q_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 :: deficit_mass ! Extra mass of a tracer if it is lower than 0
+
+
+
 
 !!!!!!!!!!!!!!!!!!!!!!!! SLOPE
-      REAL ,allocatable :: watercap_slope(:,:)                           ! (ngrid,nslope)
-      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1        ! LON x LAT field : minimum of co2 ice at each point for the first year
-      REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2        ! LON x LAT field : minimum of co2 ice at each point for the second year
+      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, 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  :: initial_h2o_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, PARAMETER :: co2_hmax = 10 			  ! 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 ()
+      REAL, PARAMETER :: rho_co2 = 1600           ! CO2 ice density (kg/m^3)
+      INTEGER :: iaval                            ! Index of the neighboord slope ()
       REAL , dimension(:,:), allocatable :: flag_co2flow(:,:)   !(ngrid,nslope)          ! To flag where there is a glacier flow
       REAL , dimension(:), allocatable :: flag_co2flow_mesh(:)  !(ngrid)          ! To flag where there is a glacier flow
 
 
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL 
+
+
+     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 :: 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 :: 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 :: 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 :: tol_soil = 1.e-6
+     INTEGER :: countmax = 15
+     INTEGER :: countloop
+     REAL, ALLOCATABLE :: tsoil_phys_PEM_saved(:,:,:)
+     REAL :: error_tsoil
+
+
+     LOGICAL :: firstcall
+!    REAL, PARAMETER :: daysec=88775.              ! duree du sol (s)  ~88775 s
+     REAL, PARAMETER :: year_day = 669
+     REAL, PARAMETER :: year_step = 1
+     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  :: q_h2o_PEM_phys_ave(:) ! averaged water vapor content
+     REAL,ALLOCATABLE  :: interp_coef(:)
+
+     REAL,ALLOCATABLE  :: ice_depth(:,:)
+     REAL,ALLOCATABLE  :: TI_locslope(:,:)
+     REAL,ALLOCATABLE  :: Tsoil_locslope(:,:)
+     REAL,ALLOCATABLE  :: Tsoil_locslope_CN(:,:)
+     REAL,ALLOCATABLE  :: Tsurf_locslope(:)
+     REAL,ALLOCATABLE  :: Tsurf_locslope_prev(:)
+     REAL,ALLOCATABLE  :: alph_locslope(:,:)
+     REAL,ALLOCATABLE  :: beta_locslope(:,:)   
+     REAL :: kcond
+     REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)
+     REAL, ALLOCATABLE :: delta_co2_adsorbded(:)
+     REAL :: totmass_adsorbded
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 ! Loop variable
-     INTEGER :: i,j,ig0,l,ig,nnq,t,islope
-! Constant
-     REAL :: pi,beta,alpha
-
+     LOGICAL :: bool_sublim
+     INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop
+     REAL :: beta = 3182.48
+     REAL :: alpha = 23.3494
+  
 ! Parallel variables
       is_sequential=.true.
@@ -205,7 +309,14 @@
       time_phys=0. !test
 
-!------------------------
-
-!=======================================================
+!      n_band_lat=18
+
+      m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
+      m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
+      A =(1/m_co2 - 1/m_noco2)
+      B=1/m_noco2
+
+
+
+!------------------------
 
 ! I   Initialisation
@@ -219,4 +330,5 @@
       CALL conf_evol( 99, .TRUE. )
 
+
 !------------------------
 
@@ -262,14 +374,28 @@
 
 !----------------------------READ start.nc --------------------- 
-
-    call infotrac_init
+    
+     call infotrac_init
 
      allocate(q(ip1jmp1,llm,nqtot))
-
      CALL dynetat0(FILE_NAME_start,vcov,ucov, &
                     teta,q,masse,ps,phis, time_0)
 
-     CALL iniconst 
+     CALL iniconst !new
      CALL inigeom
+     allocate(ap(nlayer+1))
+     allocate(bp(nlayer+1))
+     status =nf90_open(FILE_NAME_start, NF90_NOWRITE, ncid)
+     status = nf90_inq_varid(ncid, "ap", apvarid)
+     status = nf90_get_var(ncid, apvarid, ap)
+     status = nf90_inq_varid(ncid, "bp", bpvarid)
+     status = nf90_get_var(ncid, bpvarid, bp)
+     status =nf90_close(ncid)
+     
+        
+   
+    daysec=88775.
+    dtphys=0
+
+ 
 
     CALL iniphysiq(iim,jjm,llm, &
@@ -278,8 +404,9 @@
           rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
           iflag_phys)
-
-     year=day_ini
-
-
+   DO nnq=1,nqtot
+    if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
+   ENDDO 
+
+ 
 !----------------------------reading ---------------------
 
@@ -287,5 +414,10 @@
 ! First we read the initial state (starfi.nc)
 
+
        allocate(watercap_slope(ngrid,nslope))
+       allocate(TI_GCM_start(ngrid,nsoilmx,nslope))
+       allocate(q_h2o_PEM_phys_ave(ngrid))
+      allocate(inertiesoil(ngrid,nsoilmx))
+
 
          CALL phyetat0 (FILE_NAME,0,0, &
@@ -294,10 +426,27 @@
               tsurf,tsoil,albedo,emis,   &
               q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar,     &
-              watercap,nslope,tsurf_slope,                       &
+              watercap,inertiesoil,nslope,tsurf_slope,           &
               tsoil_slope,co2ice_slope,def_slope,def_slope_mean, &
-              subslope_dist,albedo_slope,emiss_slope,            &
+              subslope_dist,albedo_slope,emiss_slope, TI_GCM_start,     &
               qsurf_slope,watercap_slope)
 
-! Slope : Definition of the flat slope
+
+
+     deallocate(TI_GCM_start) !not used then
+
+     DO i=1,ngrid
+         DO nnq=1,nqtot
+           DO islope=1,nslope
+             if(qsurf_slope(i,nnq,islope).LT.0) then
+                qsurf_slope(i,nnq,islope)=0.
+             endif
+           enddo
+         enddo
+       enddo
+
+!     Define some slope statistics
+
+
+      
 
        iflat=1
@@ -306,66 +455,242 @@
            abs(def_slope_mean(iflat)))THEN
            iflat = islope
-         ENDIF
+         ENDIF     
+
        ENDDO
        PRINT*,'Flat slope for islope = ',iflat
        PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
 
-     allocate(flag_co2flow(ngrid,nslope))
-     allocate(flag_co2flow_mesh(ngrid))
-
        flag_co2flow(:,:) = 0.     
        flag_co2flow_mesh(:) = 0.
 
-
-! Then we read the evolution of water and co2 ice over the first year of the GCM run, saving only the minimum value
-
+!------------------------
+
+! I   Initialisation
+!    I_a READ run.def , run_pem.def
+!    I_b READ starfi_0.nc
+!    I_c READ GCM data and convert to the physical grid 
+
+!------------------------
+
+! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the GCM run, saving only the minimum value
+
+
+     call nb_time_step_GCM("concat_year_one.nc",timelen)
      allocate(min_h2o_ice_s_1(iim+1,jjm+1))
      allocate(min_co2_ice_s_1(iim+1,jjm+1))
-
-     allocate(ps_GCM_1(iim+1,jjm+1,2676))
-
+     allocate(vmr_co2_gcm(iim+1,jjm+1,timelen))
+     allocate(q_h2o_GCM(iim+1,jjm+1,timelen))
+     allocate(q_co2_GCM(iim+1,jjm+1,timelen))
+     allocate(ps_GCM(iim+1,jjm+1,timelen))
+     allocate(ps_GCM_yr1(iim+1,jjm+1,timelen))
      allocate(min_co2_ice_slope_1(iim+1,jjm+1,nslope))
      allocate(min_h2o_ice_slope_1(iim+1,jjm+1,nslope))
-
-     allocate(vmr_co2_gcm(iim+1,jjm+1,2676))
-
-     call pemetat0 ("concat_year_one.nc", min_h2o_ice_s_1,min_co2_ice_s_1,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_1,timelen,min_co2_ice_slope_1,min_h2o_ice_slope_1,nslope)
-
-! Then we read the evolution of water and co2 ice over the second year of the GCM run, saving only the minimum value
+     allocate(tsurf_ave(iim+1,jjm+1,nslope))
+     allocate(tsurf_ave_yr1(iim+1,jjm+1,nslope))
+     allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope))
+     allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope))
+     allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope))
+     allocate(tsurf_ave_phys(ngrid,nslope))
+     allocate(tsurf_ave_phys_yr1(ngrid,nslope))
+     allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope))
+     allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
+     allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
+
+     allocate(tsurf_GCM_timeseries(iim+1,jjm+1,nslope,timelen))
+     allocate(tsurf_phys_GCM_timeseries(ngrid,nslope,timelen))
+     allocate(co2_ice_GCM_phys_slope(ngrid,nslope,timelen))
+     allocate(co2_ice_GCM_slope(iim+1,jjm+1,nslope,timelen))
+     allocate(Tsurfave_before_saved(ngrid,nslope))
+     allocate(tsoil_phys_PEM_saved(ngrid,nsoilmx_PEM,nslope))
+     allocate(delta_co2_adsorbded(ngrid))
+     allocate(q_phys(ngrid,llm,nqtot))
+
+
+
+     call read_data_GCM("concat_year_one.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,&   
+                       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)
+
+
+! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the GCM run, saving only the minimum value
 
      allocate(min_h2o_ice_s_2(iim+1,jjm+1))
      allocate(min_co2_ice_s_2(iim+1,jjm+1))
-
      allocate(min_co2_ice_slope_2(iim+1,jjm+1,nslope))
      allocate(min_h2o_ice_slope_2(iim+1,jjm+1,nslope))
 
-     allocate(ps_GCM_2(iim+1,jjm+1,2676))
-
-     call pemetat0 ("concat_year_two.nc", min_h2o_ice_s_2,min_co2_ice_s_2,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_2,timelen,min_co2_ice_slope_2,min_h2o_ice_slope_2,nslope)
-
-! We read the evolutaion of water ice over the first year of the GCM run, saving the first and last state
+     call read_data_GCM("concat_year_two.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, &
+                  nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope)
+
+
+
+
+! The variables in the dynamic grid are transfered to the physical grid
+
 
      allocate(vmr_co2_gcm_phys(ngrid,timelen))
-
-  vmr_co2_gcm_phys(1,:)=vmr_co2_gcm(1,1,:)
-
-  ig0 = 2
-  DO j=2,jjm
-    DO i = 1, iim
-       vmr_co2_gcm_phys(ig0,:)  =vmr_co2_gcm(i,j,:)
-       ig0= ig0 + 1
-    ENDDO
-  ENDDO
-
-  vmr_co2_gcm_phys(ig0,:) = vmr_co2_gcm(1,jjm+1,:)
-
-!=======================================================
-
-! II  Run
-!    II_a Compute tendencies
-
-!------------------------
-
-!---------------------------- RUN ---------------------
+     allocate(vmr_co2_pem_phys(ngrid,timelen))
+     allocate(q_h2o_GCM_phys(ngrid,timelen))
+     allocate(q_h2o_PEM_phys(ngrid,timelen))
+     allocate(q_co2_GCM_phys(ngrid,timelen))
+     allocate(q_co2_PEM_phys(ngrid,timelen))
+
+
+     CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,vmr_co2_gcm,vmr_co2_gcm_phys)
+     CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_h2o_GCM,q_h2o_GCM_phys)
+     CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_co2_GCM,q_co2_GCM_phys)
+
+     do nnq=1,nqtot
+
+        call gr_dyn_fi(llm,iip1,jjp1,ngridmx,q(:,:,nnq),q_phys(:,:,nnq))
+     enddo
+
+
+     deallocate(vmr_co2_gcm)
+     deallocate(q_h2o_GCM)
+     deallocate(q_co2_GCM)
+
+
+    q_co2_PEM_phys(:,:)=  q_co2_GCM_phys(:,:)
+    q_h2o_PEM_phys(:,:)=  q_h2o_GCM_phys(:,:)
+
+
+     allocate(ps_phys(ngrid))
+     allocate(ps_phys_timeseries(ngrid,timelen))
+     allocate(ps_phys_timeseries_yr1(ngrid,timelen))
+
+
+
+
+
+     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys)
+
+
+! for the time series:
+
+
+
+    call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries)
+    call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1)
+    deallocate(ps_GCM)
+    deallocate(ps_GCM_yr1)
+
+
+!!!!!!!!!!!!!!!!!!!!! Initialisation for soil_PEM
+
+      call end_comsoil_h_PEM
+      call ini_comsoil_h_PEM(ngrid,nslope)
+      timestep=year_day*daysec/year_step
+
+
+      allocate(ice_depth(ngrid,nslope))
+      allocate(TI_GCM_phys(ngrid,nsoilmx,nslope))
+
+      DO islope = 1,nslope
+          CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,TI_GCM(:,:,:,islope),TI_GCM_phys(:,:,islope))
+      ENDDO
+
+
+      call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM)
+
+      deallocate(TI_GCM)
+
+
+      DO islope = 1,nslope
+
+        DO l=1,nsoilmx
+           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope))
+        ENDDO
+        DO l=nsoilmx+1,nsoilmx_PEM
+           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope))
+        ENDDO
+      ENDDO
+
+
+
+       deallocate(tsoil_ave_yr1)
+
+
+
+
+
+
+
+
+
+
+      DO islope = 1,nslope
+
+        DO l=1,nsoilmx
+           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope))
+        ENDDO
+        DO l=nsoilmx+1,nsoilmx_PEM
+           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope))
+        ENDDO
+      ENDDO
+     
+
+
+       deallocate(tsoil_ave)
+
+
+      DO islope = 1,nslope
+       
+        DO l=1,nsoilmx
+          DO t=1,timelen
+           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
+          ENDDO
+        ENDDO
+      ENDDO
+
+      deallocate(tsoil_GCM_timeseries)
+     
+ 
+
+
+       CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys)
+       
+       CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1)
+
+       deallocate(tsurf_ave)
+       deallocate(tsurf_ave_yr1)
+
+   
+
+
+      DO islope = 1,nslope
+
+        DO t=1,timelen
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t))
+        enddo
+      enddo
+      deallocate(tsurf_GCM_timeseries)
+
+
+  DO islope = 1,nslope
+
+        DO t=1,timelen
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t))
+        enddo
+      enddo
+
+
+
+      deallocate(co2_ice_GCM_slope)
+
+
+
+!------------------------
+
+! I   Initialisation
+!    I_a READ run.def , run_pem.def
+!    I_b READ starfi_0.nc
+!    I_c READ GCM data and convert to the physical grid 
+!    I_d Compute tendencies & Save initial situation
+
+!------------------------
+
+
+
+
 
 !----- Compute tendencies
@@ -380,4 +705,6 @@
      allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope))
      allocate(tendencies_co2_ice_phys_slope(ngrid,nslope))
+     allocate(tendencies_co2_ice_phys_slope_ini(ngrid,nslope))
+
 
      allocate(tendencies_h2o_ice_slope(iim+1,jjm+1,nslope))
@@ -389,38 +716,27 @@
              min_h2o_ice_s_2,iim,jjm,ngrid,tendencies_h2o_ice_phys)
 
-!  Compute the tendencies of the evolution of co2 ice over the years
-
       call compute_tendencies(tendencies_co2_ice,min_co2_ice_s_1,&
              min_co2_ice_s_2,iim,jjm,ngrid,tendencies_co2_ice_phys)
 
-!  Compute the tendencies of the evolution of water ice over the years with slopes
-
       call compute_tendencies_slope(tendencies_co2_ice_slope,min_co2_ice_slope_1,&
              min_co2_ice_slope_2,iim,jjm,ngrid,tendencies_co2_ice_phys_slope,nslope)
 
-!  Compute the tendencies of the evolution of co2 ice over the years with slopes
+
+      tendencies_co2_ice_phys_slope_ini(:,:)=tendencies_co2_ice_phys_slope(:,:)
 
       call compute_tendencies_slope(tendencies_h2o_ice_slope,min_h2o_ice_slope_1,&
              min_h2o_ice_slope_2,iim,jjm,ngrid,tendencies_h2o_ice_phys_slope,nslope)
 
-!------------------------
-
-! II  Run
-!    II_a Compute tendencies
-!    II_b Save initial situation
-
-!------------------------
-
+     
 !----- Save initial situation
 
      allocate(initial_h2o_ice(ngrid))
-
      allocate(initial_co2_ice(ngrid))
-     allocate(q_co2_GCM(ngrid))
-
+     allocate(initial_co2_ice_sublim_slope(ngrid,nslope))
      allocate(initial_co2_ice_slope(ngrid,nslope))
      allocate(initial_h2o_ice_slope(ngrid,nslope))
-
-! We save the places where water/co2 ice is sublimating
+     year_iter=0
+
+! We save the places where water ice is sublimating
   do i=1,ngrid
       if (tendencies_h2o_ice_phys(i).LT.0) then
@@ -432,4 +748,9 @@
 
       if (tendencies_co2_ice_phys_slope(i,islope).LT.0) then
+         initial_co2_ice_sublim_slope(i,islope)=1.
+      else
+         initial_co2_ice_sublim_slope(i,islope)=0.         
+      endif
+      if (co2ice_slope(i,islope).GT.0) then
          initial_co2_ice_slope(i,islope)=1.
       else
@@ -449,5 +770,4 @@
   ini_surf_h2o=0.
   Total_surface=0.
-
   do i=1,ngrid
     if (initial_h2o_ice(i).GT.0.5) then
@@ -455,5 +775,5 @@
     endif
     do islope=1,nslope
-      if (initial_co2_ice_slope(i,islope).GT.0.5) then
+      if (initial_co2_ice_sublim_slope(i,islope).GT.0.5) then
          ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
       endif
@@ -465,118 +785,232 @@
   enddo
 
-! We put the pressure on the physical grid
-
-  allocate(ps_phys(ngrid))
-
-  ig0 = iim+1
-  ps_phys(1)=ps(ig0)
-  ig0=ig0+1
-
-  DO i = 2, ngrid-1
-     if(modulo(ig0,iim+1).eq.0 .and. i.gt.3) then
-       ig0=ig0+1
-     endif
-     ps_phys(i) = ps(ig0)
-     ig0= ig0 + 1
-  ENDDO
-
-  ps_phys(ngrid)=ps(ig0)
-
-! We compute the global average pressure of the initial state
-
-  global_ave_press_old=0.
-  do i=1,ngrid
-    global_ave_press_old=global_ave_press_old+cell_area(i)*ps_phys(i)/Total_surface
-  enddo
-
-  global_ave_press_GCM=global_ave_press_old
-
-! We compute the initial pressure level
-
+     print *, "ini_surf_co2=", ini_surf_co2
+     print *, "ini_surf=", ini_surf
+     print *, "ini_surf_h2o=", ini_surf_h2o
+
+
+     allocate(local_old_press(ngrid))
+     allocate(local_new_press(ngrid))
+    
      allocate(zplev_new(ngrid,nlayer+1))
      allocate(zplev_old(ngrid,nlayer+1))
 
+
+     global_ave_press_old=0.
+     do i=1,ngrid
+       global_ave_press_old=global_ave_press_old+cell_area(i)*ps_phys(i)/Total_surface
+     enddo
+
+     global_ave_press_GCM=global_ave_press_old
+     print *, "global_ave_press_old", global_ave_press_old
+     
+
+
+!----- Time loop
+     print *, "Before timeloop"
+     allocate(flag_co2flow(ngrid,nslope))
+     allocate(flag_co2flow_mesh(ngrid))
+
+     firstcall=.TRUE.
+
+
+
+
+
+!------------------------
+
+! I   Initialisation
+!    I_a READ run.def , run_pem.def
+!    I_b READ starfi_0.nc
+!    I_c READ GCM data and convert to the physical grid 
+!    I_d Compute tendencies & Save initial situation
+!    I_e Read startfi_PEM
+
+!------------------------
+
+
+! now we complete with the PEMstart
+
+      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_yr1,ps_phys_timeseries,tsurf_phys_GCM_timeseries,tsoil_phys_PEM_timeseries,&
+      tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), co2_adsorbded_phys,delta_co2_adsorbded)
+totmass_adsorbded = 0.
+     do ig = 1,ngrid
+      do islope =1, nslope
+       do l = 1,nsoilmx_PEM - 1
+
+          totmass_adsorbded = totmass_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
+          subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
+          cell_area(ig)
+       enddo
+
+      enddo
+
+     enddo
+
+     write(*,*) "tot mass in the regolith", totmass_adsorbded
+     stop
+     deallocate(tsurf_ave_phys_yr1)
+     deallocate(ps_phys_timeseries_yr1)
+     deallocate(tsoil_ave_phys_yr1)    
+
+
+!--------------------------- END INITIALISATION ---------------------
+
+
+
+
+!---------------------------- RUN ---------------------
+
+
+
+
+!------------------------
+
+! II  Run
+!    II_a update pressure,ice and tracers 
+
+!------------------------
+
+    
+     do while (.true.)
+! II.a.1. global pressure
+     global_ave_press_new=global_ave_press_old
+
+     do i=1,ngrid
+       do islope=1,nslope
+           global_ave_press_new=global_ave_press_new-g*cell_area(i)*tendencies_co2_ice_phys_slope(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
+       enddo
+     
+      global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
+     enddo
+! II.a.2. tracers
+     allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
         DO l=1,nlayer+1
-         DO ig=1,ip1jmp1
-          zplev_old(ig,l) = ap(l)  + bp(l)*ps(ig)
+         DO ig=1,ngrid
+          zplev_old(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
+          zplev_old_timeseries(ig,l,:) = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
          ENDDO
         ENDDO
 
-
-!=======================================================
+     do i=1,ip1jmp1
+
+       ps(i)=ps(i)*global_ave_press_new/global_ave_press_old
+
+     enddo
+     do i = 1,ngrid
+       ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_old
+       ps_phys_timeseries(i,:) = ps_phys_timeseries(i,:)*global_ave_press_new/global_ave_press_old
+     enddo
+
+     allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
+
+        do l=1,nlayer+1
+         do ig=1,ngrid
+          zplev_new(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
+          zplev_new_timeseries(ig,l,:)  = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
+         enddo
+        enddo
+
+
+      print *, "1 is okay"
+
+        DO nnq=1,nqtot
+        if (noms(nnq).NE."co2") then
+          DO l=1,llm-1
+            DO ig=1,ngrid
+               q(ig,l,nnq)=q(ig,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))
+            ENDDO
+            q(:,llm,nnq)=q(:,llm-1,nnq)
+          ENDDO
+        else
+          DO l=1,llm-1
+            DO ig=1,ngrid
+                q(ig,l,nnq)=q(ig,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))  &
+                                + (   (zplev_new(ig,l)-zplev_new(ig,l+1))  -       &
+                                      (zplev_old(ig,l)-zplev_old(ig,l+1))     )  / &
+                                      (zplev_new(ig,l)-zplev_new(ig,l+1))
+            ENDDO
+           q(:,llm,nnq)=q(:,llm-1,nnq)
+          ENDDO
+        endif
+        ENDDO
+
+
+ DO nnq=1,nqtot
+       DO ig=1,ngrid
+         DO l=1,llm-1
+          if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number" .or. noms(nnq).NE."ccn_number") ) then
+              extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
+              q(ig,l,nnq)=1.
+              q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
+          endif
+          if(q(ig,l,nnq).LT.0) then
+              q(ig,l,nnq)=1E-30
+          endif
+         ENDDO
+       enddo
+   enddo
+
+
+
+
+  l=1
+  DO ig=1,ngrid
+      DO t = 1, timelen
+
+         q_h2o_PEM_phys(ig,t)=q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
+         if(q_h2o_PEM_phys(ig,t).LT.0) then
+           q_h2o_PEM_phys(ig,t)=1E-30
+         endif
+         if(q_h2o_PEM_phys(ig,t).GT.1) then
+           q_h2o_PEM_phys(ig,t)=1.
+         endif
+
+   enddo
+  enddo
+
+
+
+  DO ig=1,ngrid
+      DO t = 1, timelen
+         q_co2_PEM_phys(ig,t)=q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))  &
+                                + (   (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))  -       &
+                                      (zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))     )  / &
+                                      (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
+         if (q_co2_PEM_phys(ig,t).LT.0) then
+              q_co2_PEM_phys(ig,t)=1E-30
+         elseif (q_co2_PEM_phys(ig,t).GT.1) then
+              q_co2_PEM_phys(ig,t)=1.
+         endif
+         mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
+         vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
+      ENDDO
+  ENDDO
+   
+     deallocate(zplev_new_timeseries)
+     deallocate(zplev_old_timeseries)
+
+
+! II.a.3. Evolution of the ice
+
+
+     call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
+      print *, "4 is okay"
+
+     call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
+      print *, "5 is okay"
+
 
 !------------------------
 
 ! II  Run
-!    II_a Compute tendencies
-!    II_b Save initial situation
-!    II_c Time loop
-
-!------------------------
-
-
-!----- Time loop
-
-     year_iter=0
-
-     do while (.true.)
-
-! Compute the new global average pressure given the tendencies of co2 condensation
-     global_ave_press_new=global_ave_press_old
-     do i=1,ngrid
-       do islope=1,nslope
-           global_ave_press_new=global_ave_press_new-dt_pem*g*cell_area(i)*tendencies_co2_ice_phys_slope(i,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
-       enddo
-     enddo
-
-! Compute the new surface pressure at each physical point
-     do i=1,ngrid
-       ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_old
-     enddo
-
-! Compute the new surface pressure at each dynamical point
-     do i=1,ip1jmp1
-       ps(i)=ps(i)*global_ave_press_new/global_ave_press_old
-     enddo
-
-! Compute the new pressure level at each physical point
-
-     DO l=1,nlayer+1
-       DO ig=1,ip1jmp1
-        zplev_new(ig,l) = ap(l)  + bp(l)*ps(ig)
-       ENDDO
-     ENDDO
-
-! Compute the new values of tracer given the change of pressure and pressure level
-
-     DO nnq=1,nqtot      ! For each tracer
-      if (nnq.NE.1) then ! If it is not co2
-        DO l=1,llm-1     ! For each level
-          DO i=1,ip1jmp1 ! Each point
-             q(i,l,nnq)=q(i,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))
-          ENDDO          ! i loop
-          q(:,llm,nnq)=q(:,llm-1,nnq)
-        ENDDO            ! l loop
-      else               ! If it is co2
-        DO l=1,llm-1     ! For each level
-          DO i=1,ip1jmp1 ! Each point
-              q(i,l,nnq)=q(i,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))  &
-                              + (   (zplev_new(ig,l)-zplev_new(ig,l+1))  -       &
-                                    (zplev_old(ig,l)-zplev_old(ig,l+1))     )  / &
-                                    (zplev_new(ig,l)-zplev_new(ig,l+1))
-          ENDDO          ! i loop
-         q(:,llm,nnq)=q(:,llm-1,nnq)
-        ENDDO            ! l loop
-      endif
-     ENDDO               ! nnq loop
-
-
-!We apply the tendency on the ice
-
-     call evol_h2o_ice_s_slope(qsurf_slope(:,6,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
-
-     call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
-
-! FLOW : Make the co2 ice flow
+!    II_a update pressure, ice and tracers 
+!    II_b CO2 glaciers flows 
+
+!------------------------
+
+    
+
 
        DO ig = 1,ngrid
@@ -619,44 +1053,231 @@
              flag_co2flow_mesh(ig) = 1.
             ENDIF ! co2ice > hmax
-          ENDIF ! iflat
+          ENDIF ! iflattsoil_lope
         ENDDO !islope 
        ENDDO !ig
 
-
-! End of FLOW
-
-! We recompute the tendencies of co2 given the new surface pressure
-
-    call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,vmr_co2_gcm_phys,ps_GCM_2,global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
-
+      print *, "6 is okay"
+
+!------------------------
+
+! II  Run
+!    II_a update pressure, ice and tracers 
+!    II_b CO2 glaciers flows 
+!    II_c Update surface and soil temperatures
+
+!------------------------
+
+
+! II_c.1 Update Tsurf
+      bool_sublim=0
+      Tsurfave_before_saved(:,:) = tsurf_ave_phys(:,:)
+       DO ig = 1,ngrid
+        DO islope = 1,nslope
+           if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-5) THEN !co2ice disappeared, look for closest point without co2ice
+
+
+              if(latitude_deg(ig).gt.0) then
+                do ig_loop=ig,ngrid
+                  DO islope_loop=islope,iflat,-1
+                     if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
+                            tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
+                            bool_sublim=1
+                            exit
+                     endif
+                  enddo
+                  if (bool_sublim.eq.1) then
+                   exit
+                  endif
+                enddo
+              else
+               do ig_loop=ig,1,-1
+                  DO islope_loop=islope,iflat
+                     if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
+                            tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
+                            bool_sublim=1
+                            exit
+                     endif
+                  enddo
+                  if (bool_sublim.eq.1) then
+                   exit
+                  endif
+                enddo
+
+              endif
+              initial_co2_ice_slope(ig,islope)=0
+
+              if ((co2ice_slope(ig,islope).lt.1e-5).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then
+                 
+               albedo_slope(ig,1,islope) = albedo_h2o_frost
+               albedo_slope(ig,2,islope) = albedo_h2o_frost
+
+              else
+                            
+               albedo_slope(ig,1,islope) = albedodat(ig)
+               albedo_slope(ig,2,islope) = albedodat(ig) 
+               emiss_slope(ig,islope) = emissiv         
+
+              endif
+
+           elseif( co2ice_slope(ig,islope).GT. 1E-5) THEN !Put tsurf as tcond co2
+  
+             ave=0.
+             do t=1,timelen
+              if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-5) then
+                 ave = ave + beta/(alpha-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
+              else
+                 ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
+              endif
+             enddo
+
+             tsurf_ave_phys(ig,islope)=ave/timelen
+           endif
+        enddo
+       enddo
+
+      print *, "7 is okay"
+   
+     do t = 1,timelen
+     tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:)) 
+     enddo
+! for the start
+     do ig = 1,ngrid
+      do islope = 1,nslope
+     tsurf_slope(ig,islope) =  tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope))
+     enddo
+     enddo
+
+! II_c.2 Update soil temperature
+
+allocate(TI_locslope(ngrid,nsoilmx_PEM))
+allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
+allocate(Tsurf_locslope(ngrid))
+allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
+allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
+
+
+
+
+print *,"8 is ok"
+
+! Soil averaged
+  do islope = 1,nslope
+     TI_locslope(:,:) = TI_PEM(:,:,islope)
+     Tsoil_locslope(:,:) = tsoil_PEM(:,:,islope)
+     Tsurf_locslope(:) = tsurf_ave_phys(:,islope)
+     alph_locslope(:,:) = alph_PEM(:,:,islope)
+     beta_locslope(:,:) = beta_PEM(:,:,islope)
+     call soil_pem(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
+     call soil_pem(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
+     do t = 1,timelen
+       tsoil_phys_PEM_timeseries(:,:,islope,t) =  tsoil_phys_PEM_timeseries(:,:,islope,t) + (Tsoil_locslope(:,:)- tsoil_PEM(:,:,islope))
+     enddo
+     TI_PEM(:,:,islope) = TI_locslope(:,:)
+     tsoil_PEM(:,:,islope) = Tsoil_locslope(:,:)
+     tsurf_ave_phys(:,islope)  =  Tsurf_locslope(:)
+     alph_PEM(:,:,islope) = alph_locslope(:,:)
+     beta_PEM(:,:,islope) = beta_locslope(:,:)
+  enddo
+
+
+      print *, "9 is okay"
+
+
+!!!!!! Check for the time series
+
+   
+      do ig = 1,ngrid
+       do islope = 1,nslope
+         do iloop = 1,nsoilmx_PEM
+ 
+         if(isnan(tsoil_PEM(ig,iloop,islope))) then
+          write(*,*) "is nan tsoil", ig, iloop, islope, tsoil_PEM(ig,iloop,islope)
+           stop 
+         endif
+         enddo
+        enddo
+      enddo
+
+
+write(*,*) "before ice table, no stop"
+
+! II_c.3 Update the ice table
+     call computeice_table(timelen,ngrid,nslope,nsoilmx,nsoilmx_PEM, tsoil_phys_PEM_timeseries,tsurf_phys_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & 
+                           ps_phys_timeseries, ice_depth)
+       
+      print *, "10 is okay"
+! II_c.3 Update the soil thermal properties
+      call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),ps_phys, &
+        cell_area,ice_depth,TI_PEM)
+
+! II_c.4 Update the mass of the regolith adsorbded
+   
+  
+   call regolith_co2adsorption(ngrid,nslope,nsoilmx_PEM,timelen,ps_phys_timeseries,tsoil_PEM,TI_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope, & 
+                               co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), q_co2_PEM_phys,q_h2o_PEM_phys,co2_adsorbded_phys,delta_co2_adsorbded)
+
+
+!------------------------
+
+! II  Run
+!    II_a update pressure, ice and tracers 
+!    II_b CO2 glaciers flows 
+!    II_c Update surface and soil temperatures
+!    II_d Update the tendencies and stopping criterion
+
+!------------------------
+
+
+      call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,tendencies_co2_ice_phys_slope_ini,vmr_co2_gcm_phys,vmr_co2_pem_phys,ps_phys_timeseries,&     
+                               global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
+
+          print *, "12 is okay"
       year_iter=year_iter+dt_pem
 
-! We check if we should stop : the criterion is based on the area of ice that is sublimating
-
-      call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,6,:),STOPPING_water,ngrid,initial_h2o_ice)
-
-      call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_slope,global_ave_press_GCM,global_ave_press_new,nslope)
+!We check with we should stop
+
+
+      call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
+
+
+          print *, "13 is okay"
+      call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
+
+ print *, "criterion ok"
+
 
       print *, "STOPPING_water", STOPPING_water, "STOPPING1_water", STOPPING_1_water
       print *, "STOPPING_co2", STOPPING_co2, "STOPPING1_co2", STOPPING_1_co2
       print *, "year_iter=", year_iter
-
-      year=year+dt_pem
-
-! We check if the orbits parameters have changed too much
-      call orbit_param(year)
-
       if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2)  then
         exit
       endif
 
-! Preparing for next iteration
       global_ave_press_old=global_ave_press_new
       zplev_old=zplev_new
 
 
-      enddo
-
-! We put the correct thickness of co2ice given the slope
+      enddo  ! big loop of the pem
+
+
+!---------------------------- END RUN PEM ---------------------
+
+
+
+
+!---------------------------- OUTPUT ---------------------
+
+
+
+
+!------------------------
+
+! III Output
+!    III_a Recomp tendencies for the start and startfi
+
+!------------------------
+
+
+! III_a.1 Ice update
       DO ig = 1,ngrid
           co2ice(ig) = 0.
@@ -668,9 +1289,8 @@
       ENDDO ! of DO ig=1,ngrid
 
-! We put the correct thickness of water ice given the slope
-      DO ig = 1,ngrid 
-          qsurf(ig,6) = 0.
+      DO ig = 1,ngrid !!!!! TO DOC GHANGE IF QSURF ?ü!?!?!?
+          qsurf(ig,igcm_h2o_ice) = 0.
              DO islope = 1,nslope
-                 qsurf(ig,6) = qsurf(ig,6) + qsurf_slope(ig,6,islope) &
+                 qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
                             * subslope_dist(ig,islope) /          &
                            cos(pi*def_slope_mean(islope)/180.)
@@ -680,12 +1300,67 @@
 
 
+
+! III_a.2 Tsurf and Tsoil update
+
+
+   call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
+   tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
+
+
+      DO ig = 1,ngrid
+          tsurf(ig) = 0.
+          
+             DO islope = 1,nslope
+                 tsurf(ig) = tsurf(ig) + tsurf_slope(ig,islope) &
+                            * subslope_dist(ig,islope)      
+             ENDDO
+      ENDDO ! of DO ig=1,ngrid
+
+      DO ig = 1,ngrid
+         DO iloop = 1,nsoilmx
+          tsoil(ig,iloop) = 0.
+          inertiesoil(ig,iloop) = 0.
+             DO islope = 1,nslope
+                  tsoil(ig,iloop) =  tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
+                            * subslope_dist(ig,islope)   
+                  inertiesoil(ig,iloop) =  inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) &
+                            * subslope_dist(ig,islope)         
+             ENDDO
+          ENDDO
+      ENDDO ! of DO ig=1,ngrid
+
+! III_a.3 Surface optical properties
+
+    DO ig = 1,ngrid
+       DO l = 1,2
+        albedo(ig,l) =0. 
+        DO islope = 1,nslope
+          albedo(ig,l)= albedo(ig,l)+albedo_slope(ig,l,islope) &
+                         *subslope_dist(ig,islope) 
+        ENDDO
+       ENDDO
+      ENDDO
+
+
+    DO ig = 1,ngrid
+
+        emis(ig) =0. 
+        DO islope = 1,nslope
+          emis(ig)= emis(ig)+emiss_slope(ig,islope) &
+                         *subslope_dist(ig,islope) 
+        ENDDO
+       ENDDO
+      
+
+
 !------------------------
 
 ! III Output
-!    III_a Write startfi.nc
-
-!------------------------
-
-!----------------------------WRITE restart.nc ---------------------
+!    III_a Recomp tendencies for the start and startfi
+!    III_b Write start and starfi.nc
+
+!------------------------
+
+! III_b.1 WRITE restart.nc 
 
       ptimestep=iphysiq*daysec/REAL(day_step)/nsplit_phys
@@ -693,44 +1368,66 @@
       ztime_fin=0.
 
+         print *, "RVip1jmp1", ip1jmp1
      allocate(p(ip1jmp1,nlayer+1))
-     CALL pression (ip1jmp1,ap,bp,ps,p)
-     CALL massdair(p,masse)
-
-     day_ini=year
-
-     CALL dynredem0("restart_evol.nc", day_ini, phis)
-
-     CALL dynredem1("restart_evol.nc",   &
+         print *, "ngrid", ngrid
+         print *, "nlayer", nlayer
+          CALL pression (ip1jmp1,ap,bp,ps,p)
+     print *, "masse 1", masse(1,:)
+          CALL massdair(p,masse)
+     print *, "masse 2", masse(1,:)
+
+     do nnq=1,nqtot
+        call gr_fi_dyn(llm,ngridmx,iip1,jjp1,q_phys(:,:,nnq),q(:,:,nnq))
+     enddo
+
+
+      CALL dynredem0("restart_evol.nc", day_ini, phis)
+
+      CALL dynredem1("restart_evol.nc",   &
       		time_0,vcov,ucov,teta,q,masse,ps)
 
 
-!----------------------------WRITE restartfi.nc ---------------------
-
-     pday=year
-
-     call physdem0("restartfi_evol.nc",longitude,latitude, &
-                     nsoilmx,ngrid,nlayer,nq,              &
-                     ptimestep,pday,0.,cell_area,          &
-                     albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
-                     hmons,summit,base,nslope,def_slope,   &
-                     subslope_dist)
-
-    call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,   &
+! III_b.2 WRITE restartfi.nc 
+
+      
+             call physdem0("restartfi_evol.nc",longitude,latitude, &
+                        nsoilmx,ngrid,nlayer,nq,              &
+                        ptimestep,pday,0.,cell_area,          &
+                        albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
+                        hmons,summit,base,nslope,def_slope,   &
+                        subslope_dist)
+
+
+          call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,  &
                      ptimestep,ztime_fin,                        &
                      tsurf,tsoil,co2ice,albedo,emis,             &
                      q2,qsurf,tauscaling,totcloudfrac,wstar,     &
-                     watercap,nslope,co2ice_slope,               &
+                     watercap,inertiesoil,nslope,co2ice_slope,   &
                      tsurf_slope,tsoil_slope, albedo_slope,      &
-                     emiss_slope,qsurf_slope,watercap_slope)
-
-
-      print *, "RV : So far so good"
-
-
-
-
-      deallocate(longitude)
-      deallocate(latitude)
-      deallocate(cell_area)
+                     emiss_slope,qsurf_slope,watercap_slope, TI_GCM_phys)
+
+
+!------------------------
+
+! III Output
+!    III_a Recomp tendencies for the start and startfi
+!    III_b Write start and starfi.nc
+!    III_c Write start_pem
+
+!------------------------
+
+  
+
+         call pemdem1("restartfi_PEM.nc",ztime_fin,nsoilmx_PEM,ngrid,nslope , &
+                      tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys)
+
+
+ 
+
+      print *, "LL & RV : So far so good"
+
+
+
+
 
 
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2792)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2794)
@@ -1,225 +1,479 @@
-!
-! $Id $
-!
-SUBROUTINE pemetat0(fichnom,min_h2o_ice_s,min_co2_ice_s,iim_input,jjm_input,nlayer,vmr_co2_gcm,ps_GCM,timelen,min_co2_ice_slope,min_h2o_ice_slope,nslope)
-
-      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
-                        nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
-                        nf90_inquire_dimension,nf90_close
-
-      IMPLICIT NONE
-
-!=======================================================================
-!
-! Read initial confitions file
-!
-!=======================================================================
-
-  include "dimensions.h"
-
-!===============================================================================
-! Arguments:
-
-! INPUT:
-  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
-
-  INTEGER :: iim_input,jjm_input,nlayer,nslope
-
-! LOCAL:
-
-  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(:,:,:,:)               ! h2o_ice_s slope of the concatenated file
-  REAL, ALLOCATABLE ::  co2_ice_s_slope(:,:,:,:)               ! co2_ice_s slope of the concatenated file
-
-  REAL, ALLOCATABLE ::  q_co2_GCM(:,:,:)                       ! vmr of co2 of the first layer in the GCM year
-
-! OUTPUT:
-
-  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_h2o_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope 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) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,2676)         ! vmr of co2 of the first layer in the GCM year
-
-  REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,2676)  ! Surface pressure in the GCM year
-
-!===============================================================================
-!   Local Variables 
-  CHARACTER(LEN=256) :: msg, var, modname
-  INTEGER,PARAMETER :: length=100
-  INTEGER :: iq, fID, vID, idecal
-  INTEGER :: ierr
-  CHARACTER(len=12) :: start_file_type="earth" ! default start file type
-
-  REAL,ALLOCATABLE :: time(:) ! times stored in start
-  INTEGER :: timelen ! number of times stored in the file
-  INTEGER :: indextime ! index of selected time
-
-  INTEGER :: edges(4),corner(4)
-  INTEGER :: i,j,t
-  real,save :: m_co2, m_noco2, A , B, mmean
-
-  INTEGER :: islope
-  CHARACTER*2 :: num
-
-
-!-----------------------------------------------------------------------
-  modname="pemetat0"
-
-      m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
-      m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
-      A =(1/m_co2 - 1/m_noco2)
-      B=1/m_noco2
-
-!  Open initial state NetCDF file
-  var=fichnom
-  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
-
-      ierr = nf90_inq_varid (fID, "temps", vID)
-      IF (ierr .NE. nf90_noerr) THEN
-        write(*,*)"pemetat0: Le champ <temps> est absent"
-        write(*,*)"pemetat0: J essaie <Time>"
-        ierr = nf90_inq_varid (fID, "Time", vID)
-        IF (ierr .NE. nf90_noerr) THEN
-           write(*,*)"pemetat0: Le champ <Time> est absent"
-           write(*,*)trim(nf90_strerror(ierr))
-           CALL ABORT_gcm("pemetat0", "", 1)
-        ENDIF
-        ! Get the length of the "Time" dimension
-        ierr = nf90_inq_dimid(fID,"Time",vID)
-        ierr = nf90_inquire_dimension(fID,vID,len=timelen)
-      ELSE   
-        ! Get the length of the "temps" dimension
-        ierr = nf90_inq_dimid(fID,"temps",vID)
-        ierr = nf90_inquire_dimension(fID,vID,len=timelen)
-      ENDIF
-
-      allocate(co2_ice_s(iim+1,jjm+1,timelen))
-      allocate(q_co2_GCM(iim+1,jjm+1,timelen))
-
-      allocate(co2_ice_s_slope(iim+1,jjm+1,nslope,timelen))
-      allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
-
-
-
-! Get all the needed variable of the concatenated file
-  CALL get_var3("h2o_ice_s"   ,h2o_ice_s)
-  CALL get_var3("co2ice"   ,co2_ice_s)
-  CALL get_var3("co2_cropped"   ,q_co2_GCM)
-  CALL get_var3("ps"   ,ps_GCM)
-
-! Slope
+   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_yr1_inst,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, m_co2_regolith_phys,deltam_co2_regolith_phys)
+   
+   use iostart_PEM, only:  open_startphy, close_startphy, get_field 
+   use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
+   use comsoil_h, only:  volcapa,inertiedat 
+   use ini_soil_mod, only:  ini_icetable
+   use soil_evolution_mod, only: soil_pem_CN
+   use adsorption_mod, only : regolith_co2adsorption
+  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
+  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_yr1_inst(ngrid,timelen)
+  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) 
+
+! 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 
+! 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(nsoil_PEM-1)
+   real :: beta_tmp(nsoil_PEM-1)
+   real :: co2_ads_prev(ngrid) 
+   real :: ps_ave_yr1(ngrid)
+   real :: ps_ave_yr2(ngrid)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: read start_pem. Need a specific iostart_PEM
+!!!
+!!! Order: 1. Thermal Inertia 
+!!!        2. Soil Temperature
+!!!        3. Ice table
+!!!        4. Mass of CO2 adsorbed
+!!!
+!!! /!\ This order must be respected !
+!!! Author: LL
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+! 0. Some initializations
+ ps_ave_yr1(:) = sum(ps_yr1_inst(:,:),2)/timelen
+ ps_ave_yr2(:) = sum(ps_inst(:,:),2)/timelen
+
+!1. Run
+
+
+
+
+if (startpem_file) then
+   ! open pem initial state file:
+   call open_startphy(filename)
+ 
+!1. Thermal Inertia
+! a. General case
 
 DO islope=1,nslope
   write(num,fmt='(i2.2)') islope
-  call get_var3("co2ice_slope"//num,co2_ice_s_slope(:,:,islope,:))
-ENDDO
+  call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
+   if(.not.found) then
+      write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
+      write(*,*)'will reconstruct the values of TI_PEM'
+
+      do ig = 1,ngrid
+
+          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then 
+!!! transition
+             delta = 50.
+             TI_PEM(ig,nsoil_GCM+1,islope) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
+            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
+                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
+
+          do iloop=nsoil_GCM+2,n_1km
+            TI_PEM(ig,iloop,islope) = TI_breccia
+         enddo      
+         else ! we keep the high ti values
+           do iloop=nsoil_GCM+1,n_1km
+                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
+           enddo 
+
+         endif ! TI PEM and breccia comparison
+
+
+!! transition
+             delta = 1000.
+             TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
+            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
+                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
+          do iloop=n_1km+2,nsoil_PEM
+            TI_PEM(ig,iloop,islope) = TI_bedrock
+         enddo   
+      enddo
+
+   else
+     do iloop = nsoil_GCM+1,nsoil_PEM
+       TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
+     enddo
+   endif ! not found
+ENDDO ! islope
+
+print *,'PEMETAT0: THERMAL INERTIA DONE'
+
+
+
+! b. Special case for inertiedat, inertiedat_PEM
+call get_field("inertiedat_PEM",inertiedat_PEM,found)
+if(.not.found) then
+ do iloop = 1,nsoil_GCM
+   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
+ enddo
+
+  
+!!! zone de transition
+delta = 50.
+do ig = 1,ngrid
+if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then 
+inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
+            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
+                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
+
+
+ do iloop = nsoil_GCM+2,n_1km 
+       inertiedat_PEM(ig,iloop) = TI_breccia
+  enddo
+
+else
+   do iloop=nsoil_GCM+1,n_1km
+      inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
+   enddo 
+endif ! comparison ti breccia
+enddo!ig
+
+!!! zone de transition
+delta = 1000.
+do ig = 1,ngrid
+inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
+            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
+                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
+enddo
+
+
+
+ do iloop = n_1km+2, nsoil_PEM
+   do ig = 1,ngrid
+      inertiedat_PEM(ig,iloop) = TI_bedrock
+   enddo
+ enddo
+endif ! not found
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!2. Soil Temperature
 
 DO islope=1,nslope
   write(num,fmt='(i2.2)') islope
-  call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_slope(:,:,islope,:))
+   call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
+  if(.not.found) then
+      write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
+      write(*,*)'will reconstruct the values of Tsoil'
+      do ig = 1,ngrid
+        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
+        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))   
+
+       do iloop=nsoil_GCM+2,n_1km
+            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
+      enddo
+        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
+        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1))  
+
+       do iloop=n_1km+2,nsoil_PEM
+            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
+      enddo
+    enddo
+
+   else
+! predictor corrector: restart from year 1 of the GCM and build the evolution of
+! tsoil at depth
+
+    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
+    call soil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
+    call soil_pem(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(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
+
+
+     do iloop = nsoil_GCM+1,nsoil_PEM
+       tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
+     enddo
+   endif
+
+    do it = 1,timelen
+        do isoil = nsoil_GCM+1,nsoil_PEM
+        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
+        enddo
+     enddo
+
+
+
+    
 ENDDO
-
-! Compute the minimum over the year for each point
-  min_h2o_ice_s(:,:)=minval(h2o_ice_s,3)
-  min_co2_ice_s(:,:)=minval(co2_ice_s,3)
-
-  min_co2_ice_slope(:,:,:)=minval(co2_ice_s_slope,4)
-  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
-
-! By definition, a density is positive, we get rid of the negative values
-  DO i=1,iim+1
-    DO j = 1, jjm+1
-       if (min_co2_ice_s(i,j).LT.0) then
-          min_h2o_ice_s(i,j)  = 0.
-          min_co2_ice_s(i,j)  = 0.
-       endif
-       DO islope=1,nslope
-          if (min_co2_ice_slope(i,j,islope).LT.0) then
-            min_co2_ice_slope(i,j,islope)  = 0.
-          endif
-          if (min_h2o_ice_slope(i,j,islope).LT.0) then
-            min_h2o_ice_slope(i,j,islope)  = 0.
-          endif
-       ENDDO
-    ENDDO
-  ENDDO
-
-! Compute the volume mixing ratio of co2 in the first layer
-
-  DO i=1,iim+1
-    DO j = 1, jjm+1
-      DO t = 1, timelen
-         if (q_co2_GCM(i,j,t).LT.0) then
-              q_co2_GCM(i,j,t)=1E-10
-         elseif (q_co2_GCM(i,j,t).GT.1) then
-              q_co2_GCM(i,j,t)=1.
+print *,'PEMETAT0: SOIL TEMP  DONE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!3. Ice Table
+
+ 
+   call get_field("ice_table",ice_table,found)
+   if(.not.found) then
+      write(*,*)'PEM settings: failed loading <Ice Table>'
+      write(*,*)'will reconstruct the values of ice table'
+
+
+      call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:),  tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
+
+
+   else 
+   ! update ice table
+     call computeice_table(timelen,ngrid,nslope,nsoil_PEM,tsoil_inst,tsurf_inst,q_co2,q_h2o, ps_inst, ice_table)
+
+
+   endif
+
+print *,'PEMETAT0: ICE TABLE  DONE'
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+!4. CO2 Adsorption
+DO islope=1,nslope
+  write(num,fmt='(i2.2)') islope
+   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
+   if(.not.found) then
+      write(*,*)'PEM settings: failed loading <m_co2_regolith_phys>'
+      write(*,*)'will reconstruct the values of co2 adsorbded'
+   m_co2_regolith_phys(:,:,:) = 0.
+   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
+
+   deltam_co2_regolith_phys(:) = 0.
+   exit
+   endif
+ENDDO
+
+if (found) then
+   DO iloop = 1,nsoil_GCM
+      tsoil_tmp_yr1(:,iloop,:) = tsoil_PEM_yr1(:,iloop,:)
+
+   ENDDO  
+
+   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
+
+
+endif
+print *,'PEMETAT0: CO2 adsorption done  '
+
+
+
+call close_startphy
+
+
+else !No startfi, let's build all by hand
+
+
+
+!a) Thermal inertia
+   do islope = 1,nslope
+      do ig = 1,ngrid
+
+          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then 
+!!! transition
+             delta = 50.
+             TI_PEM(ig,nsoil_GCM+1,islope) =sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
+            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
+                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
+
+          do iloop=nsoil_GCM+2,n_1km
+            TI_PEM(ig,iloop,islope) = TI_breccia
+         enddo      
+         else ! we keep the high ti values
+           do iloop=nsoil_GCM+1,n_1km
+                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
+           enddo 
+
          endif
-         mmean=1/(A*q_co2_GCM(i,j,t) +B)
-         vmr_co2_gcm(i,j,t) = q_co2_GCM(i,j,t)*mmean/m_co2
+
+
+!! transition
+             delta = 1000.
+             TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
+            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
+                      ((layer_PEM(n_1km+1)-delta)/(TI_breccia**2))))
+          do iloop=n_1km+2,nsoil_PEM
+            TI_PEM(ig,iloop,islope) = TI_bedrock
+         enddo   
+      enddo
+
+enddo
+
+
+ do iloop = 1,nsoil_GCM
+   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
+ enddo
+!!! zone de transition
+delta = 50.
+do ig = 1,ngrid
+      if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then 
+inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
+            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
+                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
+
+
+ do iloop = nsoil_GCM+2,n_1km 
+
+       inertiedat_PEM(ig,iloop) = TI_breccia
+   
+ enddo
+else 
+   do iloop = nsoil_GCM+1,n_1km 
+       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
+    enddo
+
+endif
+enddo
+
+
+!!! zone de transition
+delta = 1000.
+do ig = 1,ngrid
+inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
+            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
+                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
+enddo
+
+
+
+ do iloop = n_1km+2, nsoil_PEM
+   do ig = 1,ngrid
+      inertiedat_PEM(ig,iloop) = TI_bedrock
+   enddo
+ enddo
+
+print *,'PEMETAT0: TI  DONE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+!b) Soil temperature 
+    do islope = 1,nslope
+
+     write(*,*) "islope=",islope
+     do ig = 1,ngrid
+        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
+        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))
+
+       do iloop=nsoil_GCM+2,n_1km
+            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
+      enddo
+        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
+        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1))
+
+       do iloop=n_1km+2,nsoil_PEM
+            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
+            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
+      enddo
+!     write(*,*) "ig, islope, T=", ig,islope,tsoil_PEM(ig,:,islope)
+     enddo
+    
+    do it = 1,timelen
+        do isoil = nsoil_GCM+1,nsoil_PEM
+        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
+        enddo
+     enddo
+enddo
+print *,'PEMETAT0: TSOIL DONE  '
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!c) Ice table    
+      do islope = 1,nslope
+      call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:), &  
+       tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
+      enddo
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+print *,'PEMETAT0: Ice table DONE  '
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!d) Regolith adsorbed
+
+   m_co2_regolith_phys(:,:,:) = 0.
+   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
+   deltam_co2_regolith_phys(:) = 0.
+
+print *,'PEMETAT0: CO2 adsorption done  '
+
+
+
+
+
+
+
+
+
+
+endif ! of if (startphy_file)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    DO ig = 1,ngrid
+      DO islope = 1,nslope
+
+     write(*,*) 'ig,islope ,ice table=',ig,islope,ice_table(ig,islope)
+
       ENDDO
     ENDDO
-  ENDDO
-
-
-
-
-
-
-  CONTAINS
-
-SUBROUTINE check_dim(n1,n2,str1,str2)
-  INTEGER,          INTENT(IN) :: n1, n2
-  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
-  CHARACTER(LEN=256) :: s1, s2
-  IF(n1/=n2) THEN
-    s1='value of '//TRIM(str1)//' ='
-    s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
-    WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(s1),n1,TRIM(s2),n2
-    CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
-  END IF
-END SUBROUTINE check_dim
-
-
-SUBROUTINE get_var1(var,v)
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-END SUBROUTINE get_var1
-
-
-SUBROUTINE get_var3(var,v) ! on U grid
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:,:,:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-
-END SUBROUTINE get_var3
-
-SUBROUTINE get_var4(var,v) 
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:,:,:,:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-END SUBROUTINE get_var4
-
-SUBROUTINE err(ierr,typ,nam)
-  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
-  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
-  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
-  IF(ierr==NF90_NoERR) RETURN
-  SELECT CASE(typ)
-    CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
-    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
-    CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
-    CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
-  END SELECT
-  CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
-END SUBROUTINE err
-
-END SUBROUTINE pemetat0
+
+
+
+!! small test
+    DO ig = 1,ngrid
+      DO islope = 1,nslope
+        DO iloop = 1,nsoil_PEM
+           if(isnan(tsoil_PEM(ig,iloop,islope))) then 
+         write(*,*) "failed nan construction", ig, iloop, islope
+           stop
+           endif
+        ENDDO
+      ENDDO
+     ENDDO
+
+      write(*,*) "construction ok, no nan" 
+     
+
+   END SUBROUTINE
Index: trunk/LMDZ.COMMON/libf/evolution/pemredem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 2794)
@@ -0,0 +1,126 @@
+module pemredem
+
+implicit none
+
+contains
+
+subroutine pemdem0(filename,lonfi,latfi,cell_area,nsoil_PEM,ngrid, &
+                         day_ini,time,nslope,def_slope,subslope_dist)
+! create physics restart file and write time-independent variables
+!  use comsoil_h, only: inertiedat, volcapa, mlayer
+  use comsoil_h_PEM, only: mlayer_PEM,fluxgeo,inertiedat_PEM
+  use iostart_PEM, only : open_restartphy, close_restartphy, & 
+                      put_var, put_field, length
+  use mod_grid_phy_lmdz, only : klon_glo
+  use planete_h, only: aphelie, emin_turb, lmixmin, obliquit, &
+                       peri_day, periheli, year_day
+  use comcstfi_h, only: g, mugaz, omeg, rad, rcp
+  use time_phylmdz_mod, only: daysec
+  implicit none
+ 
+  character(len=*), intent(in) :: filename
+  real,intent(in) :: lonfi(ngrid)
+  real,intent(in) :: latfi(ngrid)
+  integer,intent(in) :: nsoil_PEM
+  integer,intent(in) :: ngrid
+  real,intent(in) :: day_ini
+  real,intent(in) :: time
+  real, intent(in) :: cell_area(ngrid) !boundaries for bining of the slopes
+  integer,intent(in) :: nslope
+  real, intent(in) :: def_slope(nslope+1) !boundaries for bining of the slopes
+  real, intent(in) :: subslope_dist(ngrid,nslope) !undermesh statistics
+
+
+  
+  ! Create physics start file
+  call open_restartphy(filename)
+  
+  ! Write the mid-layer depths
+  call put_var("soildepth","Soil mid-layer depth",mlayer_PEM)
+  
+  ! Write longitudes
+  call put_field("longitude","Longitudes of physics grid",lonfi)
+  
+  ! Write latitudes
+  call put_field("latitude","Latitudes of physics grid",latfi)
+  
+  ! Write mesh areas
+  call put_field("area","Mesh area",cell_area)
+     
+  ! Multidimensionnal variables (nogcm undermesh slope statistics)
+  call put_var("def_slope","slope criterium stages",def_slope)
+  call put_field("subslope_dist","under mesh slope distribution",subslope_dist)
+
+  call put_var("FluxGeo","Geothermal Flux (W/m^2)",fluxgeo)
+
+  ! Close file
+  call close_restartphy
+  
+end subroutine pemdem0
+
+subroutine pemdem1(filename,time,nsoil_PEM,ngrid,nslope, &
+                    tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table,m_co2_regolith)
+  ! write time-dependent variable to restart file
+  use iostart_PEM, only : open_restartphy, close_restartphy, & 
+                      put_var, put_field
+
+  use comsoil_h_PEM, only: mlayer_PEM,fluxgeo,inertiedat_PEM
+                
+  implicit none
+  
+  include "callkeys.h"
+  
+  character(len=*),intent(in) :: filename
+  integer,intent(in) :: nsoil_PEM
+  integer,intent(in) :: ngrid
+  real,intent(in) :: time
+  real,intent(IN) :: tsoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
+  real,intent(IN) :: inertiesoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
+  real,intent(IN) :: ice_table(ngrid,nslope) !under mesh bining according to slope
+  integer,intent(IN) :: nslope
+  real,intent(in) :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
+  integer :: islope
+  CHARACTER*2 :: num  
+
+ 
+  integer :: iq
+  character(len=30) :: txt ! to store some text
+  ! indexes of water vapour & water ice tracers (if any):
+
+  
+
+  ! Open file
+  call open_restartphy(filename)
+  
+  ! First variable to write must be "Time", in order to correctly
+  ! set time counter in file
+  call put_var("Time","Temps de simulation",time)
+  
+ 
+  
+  ! Multidimensionnal variables (undermesh slope statistics)
+DO islope=1,nslope
+  write(num,fmt='(i2.2)') islope
+  call put_field("tsoil_PEM_slope"//num,"Soil temperature by slope type", &
+                 tsoil_slope_PEM(:,:,islope),time)
+  call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type", &
+                 inertiesoil_slope_PEM(:,:,islope),time)
+
+  call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith", &
+                    m_co2_regolith(:,:,islope), time)
+
+ENDDO
+
+  call put_field("ice_table","Depth of ice table", &
+                 ice_table,time)
+
+  call put_field("inertiedat_PEM","Thermal inertie of PEM ", &
+                 inertiedat_PEM,time)
+
+
+  ! Close file
+  call close_restartphy
+  
+end subroutine pemdem1
+
+end module pemredem
Index: trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2794)
@@ -0,0 +1,307 @@
+!
+! $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, & 
+             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)
+
+      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
+                        nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
+                        nf90_inquire_dimension,nf90_close
+      use comsoil_h, only: nsoilmx
+
+      IMPLICIT NONE
+
+!=======================================================================
+!
+! Read initial confitions file
+!
+!=======================================================================
+
+  include "dimensions.h"
+
+!===============================================================================
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
+
+  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, ALLOCATABLE ::  vmr_co2_gcm(:,:,:)                     !!!!vmr_co2_phys_gcm(iim_input+1,jjm_input+1,timelen)
+  REAL, INTENT(OUT) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,2676)                     !!!!vmr_co2_phys_gcm(iim_input+1,jjm_input+1,timelen)
+!  REAL, ALLOCATABLE ::  q_h2o_GCM(:,:,:)
+  REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,2676)
+  REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,2676)
+!  REAL, ALLOCATABLE ::  q_co2_GCM(:,:,:)
+  REAL, ALLOCATABLE ::  q1_co2_GCM(:,:,:)
+!  real, INTENT(OUT) ::  vmr_co2_phys_gcm(:,:)                  !!!!vmr_co2_gcm(ngrid,timelen)
+!  REAL, ALLOCATABLE ::  ps_GCM(:,:,:)
+  REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,2676)
+
+
+!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,2676) ! Surface temperature of the concatenated file
+  REAL , INTENT(OUT) ::  tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,2676) ! Soil temperature of the concatenated file
+
+  REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,2676) ! 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,2676) ! Minimum of co2_ice slope of the year
+!===============================================================================
+!   Local Variables 
+  CHARACTER(LEN=256) :: msg, var, modname
+  INTEGER,PARAMETER :: length=100
+  INTEGER :: iq, fID, vID, idecal
+  INTEGER :: ierr
+  CHARACTER(len=12) :: start_file_type="earth" ! default start file type
+
+  REAL,ALLOCATABLE :: time(:) ! times stored in start
+  INTEGER :: timelen ! number of times stored in the file
+  INTEGER :: indextime ! index of selected time
+
+  INTEGER :: edges(4),corner(4)
+  INTEGER :: i,j,t
+  real,save :: m_co2, m_noco2, A , B, mmean
+
+  INTEGER :: islope
+  CHARACTER*2 :: num
+
+
+!-----------------------------------------------------------------------
+  modname="pemetat0"
+
+      m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
+      m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
+      A =(1/m_co2 - 1/m_noco2)
+      B=1/m_noco2
+
+!  Open initial state NetCDF file
+  var=fichnom
+  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
+
+      ierr = nf90_inq_varid (fID, "temps", vID)
+      IF (ierr .NE. nf90_noerr) THEN
+        write(*,*)"pemetat0: Le champ <temps> est absent"
+        write(*,*)"pemetat0: J essaie <Time>"
+        ierr = nf90_inq_varid (fID, "Time", vID)
+        IF (ierr .NE. nf90_noerr) THEN
+           write(*,*)"pemetat0: Le champ <Time> est absent"
+           write(*,*)trim(nf90_strerror(ierr))
+     print *, "ICIIII0"
+           CALL ABORT_gcm("pemetat0", "", 1)
+        ENDIF
+        ! Get the length of the "Time" dimension
+     print *, "ICIIIITIME"
+        ierr = nf90_inq_dimid(fID,"Time",vID)
+        ierr = nf90_inquire_dimension(fID,vID,len=timelen)
+      ELSE   
+        ! Get the length of the "temps" dimension
+     print *, "ICIIIITEMPS"
+        ierr = nf90_inq_dimid(fID,"temps",vID)
+        ierr = nf90_inquire_dimension(fID,vID,len=timelen)
+      ENDIF
+
+     print *, "ICIIII"
+
+      allocate(co2_ice_s(iim+1,jjm+1,timelen))
+
+     print *, "ICIIIIAAAA"
+
+!      allocate(q_co2_GCM(iim+1,jjm+1,timelen))
+
+     print *, "ICIIIIBBBB"
+
+!      allocate(q_h2o_GCM(iim+1,jjm+1,timelen))
+
+     print *, "ICIIIICCCC"
+
+      allocate(q1_co2_GCM(iim+1,jjm+1,timelen))
+
+     print *, "ICIIII2"
+
+
+
+      allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
+      allocate(h2o_ice_s(iim+1,jjm+1,timelen))
+          print *, "ICIIII3"
+
+! Get h2o_ice_s of the concatenated file
+  CALL get_var3("h2o_ice_s"   ,h2o_ice_s)
+
+  print *, "A"
+
+  CALL get_var3("co2ice"   ,co2_ice_s)
+  CALL get_var3("co2_cropped"   ,q_co2_GCM)
+  CALL get_var3("h2o_cropped"   ,q_h2o_GCM)
+
+  print *, "B"
+
+  CALL get_var3("ps"   ,ps_GCM)
+
+  print *, "C"
+
+  print *, "nslope", nslope
+
+DO islope=1,nslope
+  write(num,fmt='(i2.2)') islope
+  call get_var3("co2ice_slope"//num,co2_ice_slope(:,:,islope,:))
+ENDDO
+
+  print *, "co2ice_slope"
+
+DO islope=1,nslope
+  write(num,fmt='(i2.2)') islope
+  call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_slope(:,:,islope,:))
+ENDDO
+
+  print *, "h2o_ice_s_slope"
+
+DO islope=1,nslope
+  write(num,fmt='(i2.2)') islope
+  call get_var3("tsurf_slope"//num,tsurf_gcm(:,:,islope,:))
+ENDDO
+
+  print *, "tsurf_slope"
+
+DO islope=1,nslope
+  write(num,fmt='(i2.2)') islope
+  call get_var4("tsoil_slope"//num,tsoil_gcm(:,:,:,islope,:))
+ENDDO
+
+  print *, "tsoil_slope"
+
+DO islope=1,nslope
+  write(num,fmt='(i2.2)') islope
+  call get_var4("inertiesoil_slope"//num,TI_gcm(:,:,:,islope,:))
+ENDDO
+
+  print *, "inertiesoil_slope"
+
+
+
+
+
+! Compute the minimum over the year for each point
+  min_h2o_ice_s(:,:)=minval(h2o_ice_s,3)
+  min_co2_ice_s(:,:)=minval(co2_ice_s,3)
+
+  min_co2_ice_slope(:,:,:)=minval(co2_ice_slope,4)
+  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
+
+!Compute averages
+
+!  DO i=1,timelen
+    tsurf_ave(:,:,:)=SUM(tsurf_gcm(:,:,:,:),4)/timelen
+    tsoil_ave(:,:,:,:)=SUM(tsoil_gcm(:,:,:,:,:),5)/timelen
+    TI_ave(:,:,:,:)=SUM(TI_gcm(:,:,:,:,:),5)/timelen
+!  ENDDO
+
+
+
+
+! By definition, a density is positive, we get rid of the negative values
+  DO i=1,iim+1
+    DO j = 1, jjm+1
+       if (min_co2_ice_s(i,j).LT.0) then
+          min_h2o_ice_s(i,j)  = 0.
+          min_co2_ice_s(i,j)  = 0.
+       endif
+       DO islope=1,nslope
+          if (min_co2_ice_slope(i,j,islope).LT.0) then
+            min_co2_ice_slope(i,j,islope)  = 0.
+          endif
+          if (min_h2o_ice_slope(i,j,islope).LT.0) then
+            min_h2o_ice_slope(i,j,islope)  = 0.
+          endif
+       ENDDO
+    ENDDO
+  ENDDO
+
+  DO i=1,iim+1
+    DO j = 1, jjm+1
+      DO t = 1, timelen
+!         q1_co2_GCM(i,j,t)=q_co2_GCM(i,j,t)
+!         ps_GCM(i,j,t)=ps(i)
+!c             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
+         if (q_co2_GCM(i,j,t).LT.0) then
+              q_co2_GCM(i,j,t)=1E-10
+         elseif (q_co2_GCM(i,j,t).GT.1) then
+              q_co2_GCM(i,j,t)=1.
+         endif
+         if (q_h2o_GCM(i,j,t).LT.0) then
+              q_h2o_GCM(i,j,t)=1E-30
+         elseif (q_h2o_GCM(i,j,t).GT.1) then
+              q_h2o_GCM(i,j,t)=1.
+         endif
+         mmean=1/(A*q_co2_GCM(i,j,t) +B)
+         vmr_co2_gcm(i,j,t) = q_co2_GCM(i,j,t)*mmean/m_co2
+      ENDDO
+    ENDDO
+  ENDDO
+
+
+
+
+
+
+  CONTAINS
+
+SUBROUTINE check_dim(n1,n2,str1,str2)
+  INTEGER,          INTENT(IN) :: n1, n2
+  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
+  CHARACTER(LEN=256) :: s1, s2
+  IF(n1/=n2) THEN
+    s1='value of '//TRIM(str1)//' ='
+    s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
+    WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(s1),n1,TRIM(s2),n2
+    CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
+  END IF
+END SUBROUTINE check_dim
+
+
+SUBROUTINE get_var1(var,v)
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+END SUBROUTINE get_var1
+
+
+SUBROUTINE get_var3(var,v) ! on U grid
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:,:,:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+
+END SUBROUTINE get_var3
+
+SUBROUTINE get_var4(var,v) 
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:,:,:,:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+END SUBROUTINE get_var4
+
+SUBROUTINE err(ierr,typ,nam)
+  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
+  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
+  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
+  IF(ierr==NF90_NoERR) RETURN
+  SELECT CASE(typ)
+    CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
+    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
+    CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
+    CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
+  END SELECT
+  CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
+END SUBROUTINE err
+
+END SUBROUTINE read_data_gcm
Index: trunk/LMDZ.COMMON/libf/evolution/read_data_humidity.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/read_data_humidity.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/read_data_humidity.F90	(revision 2794)
@@ -0,0 +1,139 @@
+!
+! $Id $
+!
+SUBROUTINE read_data_humidity(fichnom,iim_input,jjm_input,nlayer,timelen, nslope,q_h2o_GCM,q_co2_GCM,ps_GCM,tsurf_gcm)
+
+      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
+                        nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
+                        nf90_inquire_dimension,nf90_close
+      use comsoil_h, only: nsoilmx
+
+      IMPLICIT NONE
+
+!=======================================================================
+!
+! Read initial confitions file to compute the ice table
+!
+!=======================================================================
+
+  include "dimensions.h"
+
+!===============================================================================
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
+
+  INTEGER, INTENT(IN) :: iim_input,jjm_input,nlayer,nslope,timelen
+
+! atmosphere
+
+  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)
+
+
+! Surface
+  REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen) ! Surface temperature of the concatenated file
+
+
+
+!===============================================================================
+!   Local Variables 
+  CHARACTER(LEN=256) :: msg, var, modname
+  INTEGER,PARAMETER :: length=100
+  INTEGER :: iq, fID, vID, idecal
+  INTEGER :: ierr
+  CHARACTER(len=12) :: start_file_type="earth" ! default start file type
+
+  REAL,ALLOCATABLE :: time(:) ! times stored in start
+  INTEGER :: indextime ! index of selected time
+
+  INTEGER :: edges(4),corner(4)
+  INTEGER :: i,j,t
+  real,save :: m_co2, m_noco2, A , B, mmean
+
+  INTEGER :: islope
+  CHARACTER*2 :: num
+
+
+!-----------------------------------------------------------------------
+  modname="pemetat0"
+
+      m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
+      m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
+      A =(1/m_co2 - 1/m_noco2)
+      B=1/m_noco2
+
+!  Open initial state NetCDF file
+  var=fichnom
+  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
+
+  CALL get_var3("co2_cropped" ,q_co2_GCM)
+  CALL get_var3("h2o_cropped" ,q_h2o_GCM)
+  CALL get_var3("ps"   ,ps_GCM)
+
+
+  DO islope=1,nslope
+    write(num,fmt='(i2.2)') islope
+    call get_var3("tsurf_slope"//num,tsurf_gcm(:,:,islope,:))
+  ENDDO
+
+  print *, "tsurf_slope"
+
+
+
+
+
+  CONTAINS
+
+SUBROUTINE check_dim(n1,n2,str1,str2)
+  INTEGER,          INTENT(IN) :: n1, n2
+  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
+  CHARACTER(LEN=256) :: s1, s2
+  IF(n1/=n2) THEN
+    s1='value of '//TRIM(str1)//' ='
+    s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
+    WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(s1),n1,TRIM(s2),n2
+    CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
+  END IF
+END SUBROUTINE check_dim
+
+
+SUBROUTINE get_var1(var,v)
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+END SUBROUTINE get_var1
+
+
+SUBROUTINE get_var3(var,v) ! on U grid
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:,:,:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+
+END SUBROUTINE get_var3
+
+SUBROUTINE get_var4(var,v) 
+  CHARACTER(LEN=*), INTENT(IN)  :: var
+  REAL,             INTENT(OUT) :: v(:,:,:,:)
+  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
+  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
+END SUBROUTINE get_var4
+
+SUBROUTINE err(ierr,typ,nam)
+  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
+  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
+  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
+  IF(ierr==NF90_NoERR) RETURN
+  SELECT CASE(typ)
+    CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
+    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
+    CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
+    CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
+  END SELECT
+  CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
+END SUBROUTINE err
+
+END SUBROUTINE read_data_humidity
Index: trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90	(revision 2792)
+++ trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90	(revision 2794)
@@ -2,5 +2,5 @@
 ! $Id $
 !
-SUBROUTINE recomp_tend_co2_slope(tendencies_co2_ice_phys,vmr_co2_gcm,ps_GCM_2,global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
+SUBROUTINE recomp_tend_co2_slope(tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,vmr_co2_gcm,vmr_co2_pem,ps_GCM_2,global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
 
       IMPLICIT NONE
@@ -16,13 +16,17 @@
 
 !   INPUT
-  INTEGER, intent(in) :: timelen,ngrid,nslope                    ! # of timepoint if diagfi, # grid point, # subslope
-  REAL, INTENT(in) ::  vmr_co2_gcm(ngrid,timelen)                ! physical point x timelen field : Volume mixing ratio of co2 in the first layer
-  REAL, intent(in) :: ps_GCM_2(ngrid,timelen)                    ! physical point x timelen field : Surface pressure in the GCM
-  REAL, intent(in) :: global_ave_press_GCM                       ! Average pressure in the GCM output
-  REAL, intent(in) :: global_ave_press_new                       ! Actual average pressure
+  INTEGER, intent(in) :: timelen,ngrid,nslope
+  REAL, INTENT(in) ::  vmr_co2_gcm(ngrid,timelen)                ! physical point field : Volume mixing ratio of co2 in the first layer
+  REAL, INTENT(in) ::  vmr_co2_pem(ngrid,timelen)                ! physical point field : Volume mixing ratio of co2 in the first layer
+  REAL, intent(in) :: ps_GCM_2(ngrid,timelen)                 ! physical point field : Surface pressure in the GCM
+!  REAL, INTENT(in) ::  q_co2_GCM(ngrid)                ! physical point field : Density of co2 in the first layer
+!  REAL, intent(in) :: ps_GCM(ngrid)                 ! physical point field : Density of co2 in the first layer
+  REAL, intent(in) :: global_ave_press_GCM
+  REAL, intent(in) :: global_ave_press_new
+  REAL, intent(in) ::  tendencies_co2_ice_phys_ini(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
 
 
 !   OUTPUT
-  REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope)  ! physical point field : Evolution of perenial ice over one year
+  REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
 
 
@@ -30,6 +34,6 @@
 !   ----
 
-  INTEGER :: i,t,islope                                          ! Loop variables 
-  REAL :: eps, sigma, L, beta, alpha, coef, ave                  ! Constant
+  INTEGER :: i,t,islope
+  REAL :: eps, sigma, L, beta, alpha, coef, ave
 
   eps=0.95
@@ -41,17 +45,29 @@
   coef=669*24*3600*eps*sigma/L
 
+  print *, "coef", coef
+  print *, "global_ave_press_GCM", global_ave_press_GCM
+  print *, "global_ave_press_new", global_ave_press_new
+
 ! Evolution of the water ice for each physical point
   do i=1,ngrid
     do islope=1,nslope
        ave=0.
-       if (tendencies_co2_ice_phys(i,islope).LT. 1E-5 .and. tendencies_co2_ice_phys(i,islope).GT.-1E-5) then
-       else
-         do t=1,timelen
-           ave=ave+(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100)))**4  &
-              -(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)*global_ave_press_GCM/global_ave_press_new/100)))**4
-         enddo
-         tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys(i,islope)+coef*ave/timelen
-       endif
+    do t=1,timelen
+
+!       write(*,*)'i,t=',i,t,islope, alpha,beta,ave,vmr_co2_gcm(i,t),vmr_co2_pem(i,t),ps_GCM_2(i,t),global_ave_press_GCM,global_ave_press_new
+       ave=ave+(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100)))**4  &
+              -(beta/(alpha-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*global_ave_press_GCM/global_ave_press_new/100)))**4
     enddo
+!    print *, "i", i
+!    print *, "tendencies_co2_ice_phys_ini bef", tendencies_co2_ice_phys_ini(i,islope)
+!    tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys_ini(i,islope)+coef*ave/timelen
+    tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys_ini(i,islope)-coef*ave/timelen
+
+!    print *, "tendencies after", tendencies_co2_ice_phys(i,islope)
+!    print *, "ave", ave
+!    print *, "timelen", timelen
+!    print *, "vmr_co2_pem(i,t)*ps_GCM_2(i,t)", vmr_co2_pem(i,t)*ps_GCM_2(i,t)
+!    print *, "-------------------"
+  enddo
   enddo
 
Index: trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM.F90	(revision 2794)
@@ -0,0 +1,99 @@
+      SUBROUTINE soil_TIfeedback_PEM(ngrid,nsoil,icecover,   newtherm_i)
+
+!      use tracer_mod, only:  rho_ice
+      use comsoil_h_PEM, only: layer_PEM, inertiedat_PEM
+
+      IMPLICIT NONE
+
+!=======================================================================
+!   Description :
+!       Surface water ice / Thermal inertia feedback.
+!
+!   When surface water-ice is thick enough, this routine creates a new
+!   soil thermal inertia with three different layers :
+!   - One layer of surface water ice (the thickness is given
+!     by the variable icecover (in kg of ice per m2) and the thermal
+!     inertia is prescribed by inert_h2o_ice (see surfdat_h)); 
+!   - A transitional layer of mixed thermal inertia;
+!   - A last layer of regolith below the ice cover whose thermal inertia
+!     is equal to inertiedat.
+!
+!  To use the model :
+!      SET THE tifeedback LOGICAL TO ".true." in callphys.def.
+!
+!  Author: Adapted from J.-B. Madeleine Mars 2008 ( Updated November 2012) by LL, 2022
+!=======================================================================
+ 
+
+!Local variables
+!---------------
+
+      INTEGER :: ig                     ! Grid point (ngrid)
+      INTEGER :: ik                     ! Grid point (nsoil)
+      INTEGER :: iref                   ! Ice/Regolith boundary index
+      INTEGER, INTENT(IN) :: ngrid                  ! Number of horizontal grid points 
+      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)
+!Inputs
+!------
+
+      REAL ,INTENT(IN):: icecover(ngrid)         ! tracer on the surface (kg.m-2)
+                                       
+                                        ! water ice
+!Outputs
+!-------
+
+      REAL,INTENT(INOUT) :: newtherm_i(ngrid,nsoil)    ! New soil thermal inertia
+
+
+      prev_thermi(:,:) = newtherm_i(:,:)
+
+!Creating the new soil thermal inertia table
+!-------------------------------------------
+      DO ig=1,ngrid
+!      Calculating the ice cover thickness
+       
+        icedepth=icecover(ig)/rho_ice
+       
+!      If the ice cover is too thick or watercaptag=true,
+!        the entire column is changed :
+        IF (icedepth.ge.layer_PEM(nsoil)) THEN
+          DO ik=1,nsoil
+               newtherm_i(ig,ik)=max(inert_h2o_ice,prev_thermi(ig,ik))
+          ENDDO 
+!      We neglect the effect of a very thin ice cover :
+        ELSE IF (icedepth.lt.layer_PEM(1)) THEN
+          DO ik=1,nsoil
+               newtherm_i(ig,ik)=inertiedat_PEM(ig,ik)
+          ENDDO 
+        ELSE
+!        Ice/regolith boundary index :
+          iref=1
+!        Otherwise, we find the ice/regolith boundary:
+          DO ik=1,nsoil-1
+              IF ((icedepth.ge.layer_PEM(ik)).and. (icedepth.lt.layer_PEM(ik+1))) THEN
+                  iref=ik+1
+                  EXIT
+              ENDIF
+          ENDDO
+!        And we change the thermal inertia:
+          DO ik=1,iref-1
+            newtherm_i(ig,ik)=max(inert_h2o_ice,prev_thermi(ig,ik))
+          ENDDO
+!        Transition (based on the equations of thermal conduction):
+          newtherm_i(ig,iref)=sqrt( (layer_PEM(iref)-layer_PEM(iref-1)) / &
+           ( ((icedepth-layer_PEM(iref-1))/newtherm_i(ig,iref-1)**2) + &
+             ((layer_PEM(iref)-icedepth)/inertiedat_PEM(ig,ik)**2) ) )
+!        Underlying regolith:
+          DO ik=iref+1,nsoil
+              newtherm_i(ig,ik)=inertiedat_PEM(ig,ik)
+          ENDDO
+        ENDIF ! icedepth
+      ENDDO ! ig
+
+!=======================================================================
+      RETURN
+      END
Index: trunk/LMDZ.COMMON/libf/evolution/soil_evolution_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_evolution_mod.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_evolution_mod.F90	(revision 2794)
@@ -0,0 +1,154 @@
+      MODULE soil_evolution_mod
+
+
+      IMPLICIT NONE
+      
+      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"
+!#include "dimphys.h"
+
+!#include"comsoil.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 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90	(revision 2794)
@@ -0,0 +1,162 @@
+      subroutine soil_pem(ngrid,nsoil,firstcall, &
+               therm_i,                          &
+               timestep,tsurf,tsoil,alph_PEM,beta_PEM)
+
+
+      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
+                          mthermdiff_PEM, thermdiff_PEM, coefq_PEM, &
+                          coefd_PEM, mu_PEM,fluxgeo
+      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"
+!#include "dimphys.h"
+
+!#include"comsoil.h"
+
+
+!-----------------------------------------------------------------------
+!  arguments
+!  ---------
+!  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
+ 
+! 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)
+
+
+
+
+      
+! local variables:
+      integer ig,ik
+
+
+     
+
+! 0. Initialisations and preprocessing step
+ if (firstcall) then
+    
+! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
+      do ig=1,ngrid
+        do ik=0,nsoil-1
+	  mthermdiff_PEM(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
+    
+	enddo
+      enddo
+
+
+! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
+      do ig=1,ngrid
+        do ik=1,nsoil-1
+      thermdiff_PEM(ig,ik)=((layer_PEM(ik)-mlayer_PEM(ik-1))*mthermdiff_PEM(ig,ik) &
+                     +(mlayer_PEM(ik)-layer_PEM(ik))*mthermdiff_PEM(ig,ik-1))  &
+                         /(mlayer_PEM(ik)-mlayer_PEM(ik-1))
+
+	enddo
+      enddo
+
+! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_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
+
+      do ig=1,ngrid
+	! d_k
+	do ik=1,nsoil-1
+	  coefd_PEM(ig,ik)=thermdiff_PEM(ig,ik)/(mlayer_PEM(ik)-mlayer_PEM(ik-1))
+	enddo
+	
+	! alph_PEM_{N-1}
+	alph_PEM(ig,nsoil-1)=coefd_PEM(ig,nsoil-1)/ 			 &
+                       (coefq_PEM(nsoil-1)+coefd_PEM(ig,nsoil-1))
+        ! alph_PEM_k
+        do ik=nsoil-2,1,-1
+	  alph_PEM(ig,ik)=coefd_PEM(ig,ik)/(coefq_PEM(ik)+coefd_PEM(ig,ik+1)*    &
+                                   (1.-alph_PEM(ig,ik+1))+coefd_PEM(ig,ik))
+	enddo
+
+
+      enddo ! of do ig=1,ngrid
+
+      
+  
+endif ! of if (firstcall)
+
+
+
+      IF (.not.firstcall) THEN
+!  2. Compute soil temperatures
+! First layer:
+
+
+
+
+      do ig=1,ngrid
+        tsoil(ig,1)=(tsurf(ig)+mu_PEM*beta_PEM(ig,1)* &  
+                                  thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
+                   (1.+mu_PEM*(1.0-alph_PEM(ig,1))*&
+                    thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
+     
+
+  
+
+
+! Other layers:
+      do ik=1,nsoil-1
+
+          tsoil(ig,ik+1)=alph_PEM(ig,ik)*tsoil(ig,ik)+beta_PEM(ig,ik)
+
+ 
+       enddo
+      
+     enddo
+          ENDIF
+
+
+
+
+!  2. Compute beta_PEM coefficients (preprocessing for next time step)
+! Bottom layer, beta_PEM_{N-1}
+      do ig=1,ngrid
+        beta_PEM(ig,nsoil-1)=coefq_PEM(nsoil-1)*tsoil(ig,nsoil)		 &
+                        /(coefq_PEM(nsoil-1)+coefd_PEM(ig,nsoil-1)) & 
+                 +  fluxgeo/(coefq_PEM(nsoil-1)+coefd_PEM(ig,nsoil-1))
+      enddo
+! Other layers
+      do ik=nsoil-2,1,-1
+        do ig=1,ngrid
+	  beta_PEM(ig,ik)=(coefq_PEM(ik)*tsoil(ig,ik+1)+		 &
+                      coefd_PEM(ig,ik+1)*beta_PEM(ig,ik+1))/		 &
+                      (coefq_PEM(ik)+coefd_PEM(ig,ik+1)*(1.0-alph_PEM(ig,ik+1)) &
+                       +coefd_PEM(ig,ik))
+	enddo
+      enddo
+
+
+
+      end
+
Index: trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F	(revision 2794)
@@ -0,0 +1,105 @@
+      subroutine soil_settings_PEM(ngrid,nslope,nsoil_PEM,nsoil_GCM,
+     &   TI_GCM,TI_PEM)
+
+
+!      use netcdf
+      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM
+      use comsoil_h, only: inertiedat,layer,mlayer, volcapa
+      use iostart, only: inquire_field_ndims, get_var, get_field,
+     &                   inquire_field, inquire_dimension_length
+
+      implicit none
+
+!======================================================================
+!  Author: Ehouarn Millour (07/2006)
+!
+!  Purpose: Read and/or initialise soil depths and properties
+!
+! Modifications: Aug.2010 EM : use NetCDF90 to load variables (enables using
+!                      r4 or r8 restarts independently of having compiled
+!                      the GCM in r4 or r8)
+!                June 2013 TN : Possibility to read files with a time axis
+!
+!
+!  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
+!     necessary
+!======================================================================
+
+!======================================================================
+!  arguments
+!  ---------
+!  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_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
+
+!======================================================================
+! local variables:
+      integer ig,iloop,islope	! loop counters
+      logical found
+
+      real alpha,lay1 ! coefficients for building layers
+      real xmin,xmax ! to display min and max of a field
+
+
+      
+!======================================================================
+
+! 1. Depth coordinate
+! -------------------
+
+! 1.4 Build mlayer(), if necessary
+!      if (interpol) then
+      ! default mlayer distribution, following a power law:
+      !  mlayer(k)=lay1*alpha**(k-1/2)
+        lay1=2.e-4
+	alpha=2
+        do iloop=0,nsoil_PEM-1
+	  mlayer_PEM(iloop)=lay1*(alpha**(iloop-0.5))
+	enddo
+!      endif
+
+! 1.5 Build layer(); following the same law as mlayer()
+      ! Assuming layer distribution follows mid-layer law:
+      ! layer(k)=lay1*alpha**(k-1)
+      lay1=sqrt(mlayer_PEM(0)*mlayer_PEM(1))
+      alpha=mlayer_PEM(1)/mlayer_PEM(0)
+      do iloop=1,nsoil_PEM
+        layer_PEM(iloop)=lay1*(alpha**(iloop-1))
+      enddo
+
+
+
+
+! 2. Thermal inertia (note: it is declared in comsoil_h)
+! ------------------
+
+
+      do ig = 1,ngrid
+        do islope = 1,nslope
+          do iloop = 1,nsoil_GCM
+            TI_PEM(ig,iloop,islope) = TI_GCM(ig,iloop,islope)
+          enddo
+          if(nsoil_PEM.gt.nsoil_GCM) then
+           do iloop = nsoil_GCM+1,nsoil_PEM
+             TI_PEM(ig,iloop,islope) = TI_GCM(ig,nsoil_GCM,islope)
+           enddo
+          endif
+        enddo
+      enddo
+      end subroutine soil_settings_PEM
Index: trunk/LMDZ.COMMON/libf/evolution/update_soil.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/update_soil.F90	(revision 2794)
+++ trunk/LMDZ.COMMON/libf/evolution/update_soil.F90	(revision 2794)
@@ -0,0 +1,183 @@
+   SUBROUTINE update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,ps_new,&
+                          cellarea,ice_depth,TI_PEM)
+
+
+ USE comsoil_h, only:  inertiedat, volcapa
+ USE comsoil_h_PEM, only: layer_PEM,n_1km,inertiedat_PEM
+ USE vertical_layers_mod, ONLY: ap,bp
+ implicit none
+! Input: 
+ INTEGER,INTENT(IN) ::  ngrid, nslope, nsoil_PEM
+ REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope)
+ REAL,INTENT(IN) :: ps_new(ngrid)
+ REAL,INTENT(IN) :: cellarea(ngrid)
+ REAL,INTENT(IN) :: co2ice(ngrid,nslope)
+ REAL,INTENT(IN) :: waterice(ngrid,nslope)
+ REAL,INTENT(in) :: ice_depth(ngrid,nslope)
+ 
+
+! Outputs:
+
+ REAL,INTENT(INOUT) :: TI_PEM(ngrid,nsoil_PEM,nslope)
+  
+! Constants:
+
+ REAL ::  alpha = 0.2
+ REAL ::  beta = 1.08e7
+ REAL ::  To = 273.15
+ REAL ::  R =  8.314
+ REAL ::  L =  51058.
+ REAL ::  inertie_thresold = 800. ! look for ice
+ REAL ::  inertie_co2glaciers = 2120 ! Mellon et al. 2000
+ REAL ::  inertie_averaged = 250 ! Mellon et al. 2000
+ REAL ::  ice_inertia = 2120  ! Inertia of ice
+ REAL ::  surfaceice_inertia = 800  ! Inertia of ice
+ REAL ::  P610 = 610.0
+ REAL ::  m_h2o = 18.01528E-3
+ REAL ::  m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
+ REAL ::  m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
+ REAL ::  A,B,mmean             
+
+! Local variables:
+
+ INTEGER :: ig,islope,iloop,iref,k
+ REAL :: regolith_inertia(ngrid,nslope) ! TI of the regolith 
+ REAL :: d(ngrid,nsoil_PEM,nslope)
+ REAL :: p_avg_new
+ REAL :: Total_surface
+ INTEGER :: ispermanent_co2glaciers(ngrid,nslope)
+ INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)
+
+
+! 0. Initialisation
+
+
+
+
+  Total_surface = sum(cellarea)
+  p_avg_new = 0.
+  do ig = 1,ngrid
+    p_avg_new = p_avg_new + ps_new(ig)*cellarea(ig)/Total_surface
+  enddo
+
+
+
+  do ig = 1,ngrid
+    do islope = 1,nslope
+     if((abs(tend_h2oglaciers(ig,islope)).lt.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)).lt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then
+        ispermanent_co2glaciers(ig,islope) = 1
+     else
+        ispermanent_co2glaciers(ig,islope) = 0
+     endif
+    enddo
+
+  enddo
+
+
+ ispermanent_co2glaciers(:,:) = 0
+ ispermanent_h2oglaciers(:,:) = 0
+
+
+!  1.Ice TI feedback
+
+!  do ig=1,ngrid
+!    do islope=1,nslope
+!      if ((ispermanent_co2glaciers(ig,islope).eq.1).or.(ispermanent_h2oglaciers(ig,islope).eq.1)) then
+!       do iloop = 1,n_1km
+!         TI_PEM(ig,iloop,islope) = surfaceice_inertia
+!       enddo
+!      endif
+!    enddo
+!   enddo
+    do islope = 1,nslope
+     call soil_TIfeedback_PEM(ngrid,nsoil_PEM,waterice(:,islope),   TI_PEM(:,:,islope))
+    enddo
+
+! 2. Modification of the regolith thermal inertia.
+
+
+
+  do ig=1,ngrid
+    do islope=1,nslope
+     do iloop =1,n_1km
+       d(ig,iloop,islope) = ((inertiedat_PEM(ig,iloop)*inertiedat_PEM(ig,iloop))/(volcapa*alpha*P610**0.6))**(-1/(0.11*log10(P610/beta)))
+   if(TI_PEM(ig,iloop,islope).lt.inertie_thresold) then !           we are modifying the regolith properties, not ice
+          TI_PEM(ig,iloop,islope) = sqrt(volcapa*alpha*(p_avg_new**0.6)* d(ig,iloop,islope)**(-0.11*log10(p_avg_new/beta)))
+
+   endif 
+     enddo
+   enddo
+enddo
+
+
+
+
+
+
+!  3. Build new TI for the PEM
+! a) For the regolith
+  do ig=1,ngrid
+    do islope=1,nslope
+
+  if (ice_depth(ig,islope).gt. -1.e-10) then
+! 3.0 FIrst if permanent ice
+   if (ice_depth(ig,islope).lt. 1e-10) then
+      do iloop = 1,nsoil_PEM
+       TI_PEM(ig,iloop,islope)=max(ice_inertia,inertiedat_PEM(ig,iloop))
+      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_depth(ig,islope).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 iloop=1,iref
+         TI_PEM(ig,iloop,islope) =TI_PEM(ig,iloop,islope)
+      enddo
+      if (iref.lt.nsoil_PEM) then
+        if (iref.ne.0) then
+          ! mixed layer
+           TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
+            (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ &
+                      ((layer_PEM(iref+1)-ice_depth(ig,islope))/(ice_inertia**2))))
+        else ! first layer is already a mixed layer
+          ! (ie: take layer(iref=0)=0)
+          TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
+                          (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
+                           ((layer_PEM(1)-ice_depth(ig,islope))/(ice_inertia**2))))
+        endif ! of if (iref.ne.0)        
+        ! lower layers of pure ice
+        do iloop=iref+2,nsoil_PEM
+          TI_PEM(ig,iloop,islope)=ice_inertia
+        enddo
+      endif ! of if (iref.lt.(nsoilmx))
+      endif ! permanent glaciers
+      endif ! depth > 0
+     enddo !islope
+    enddo !ig
+
+
+
+
+
+
+
+!=======================================================================
+      RETURN
+      END
