Index: unk/LMDZ.COMMON/libf/evolution/abort_pem.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/abort_pem.F	(revision 3075)
+++ 	(revision )
@@ -1,55 +1,0 @@
-!
-! $Id: abort_pem.F 
-!
-c
-c
-      SUBROUTINE abort_pem(modname, message, ierr)
-     
-#ifdef CPP_IOIPSL
-      USE IOIPSL
-#else
-! if not using IOIPSL, we still need to use (a local version of) getin_dump
-      USE ioipsl_getincom
-#endif
-
-#ifdef CPP_XIOS
-    ! ug Pour les sorties XIOS
-      USE wxios
-#endif
-
-#include "iniprint.h"
- 
-C
-C Stops the simulation cleanly, closing files and printing various
-C comments
-C
-C  Input: modname = name of calling program
-C         message = stuff to print
-C         ierr    = severity of situation ( = 0 normal )
-
-      character(len=*), intent(in):: modname
-      integer, intent(in):: ierr
-      character(len=*), intent(in):: message
-
-      write(lunout,*) 'in abort_pem'
-
-#ifdef CPP_XIOS
-    !Fermeture propre de XIOS
-      CALL wxios_close()
-#endif
-
-#ifdef CPP_IOIPSL
-      call histclo
-      call restclo
-#endif
-      call getin_dump
-      write(lunout,*) 'Stopping in ', modname
-      write(lunout,*) 'Reason = ',message
-      if (ierr .eq. 0) then
-        write(lunout,*) 'Everything is cool'
-        stop
-      else
-        write(lunout,*) 'Houston, we have a problem ', ierr
-        stop 1
-      endif
-      END
Index: /trunk/LMDZ.COMMON/libf/evolution/abort_pem_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/abort_pem_mod.F90	(revision 3076)
+++ /trunk/LMDZ.COMMON/libf/evolution/abort_pem_mod.F90	(revision 3076)
@@ -0,0 +1,59 @@
+MODULE abort_pem_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE abort_pem(modname,message,ierr)
+
+#ifdef CPP_IOIPSL
+    use IOIPSL
+#else
+    ! if not using IOIPSL, we still need to use (a local version of) getin_dump
+    use ioipsl_getincom
+#endif
+
+#ifdef CPP_XIOS
+    use wxios ! ug Pour les sorties XIOS
+#endif
+
+implicit none
+
+#include "iniprint.h"
+
+! Stop the simulation cleanly, closing files and printing various comments
+! Input: modname = name of calling program
+!        message = stuff to print
+!        ierr    = severity of situation (= 0 normal)
+character(len = *), intent(in) :: modname
+integer,            intent(in) :: ierr
+character(len = *), intent(in) :: message
+
+!----- Code
+write(lunout,*) 'in abort_pem'
+
+#ifdef CPP_XIOS
+    CALL wxios_close() ! Fermeture propre de XIOS
+#endif
+
+#ifdef CPP_IOIPSL
+    call histclo
+    call restclo
+#endif
+
+call getin_dump
+write(lunout,*) 'Stopping in ', modname
+write(lunout,*) 'Reason = ', message
+if (ierr == 0) then
+    write(lunout,*) 'Everything is cool'
+    stop
+else
+    write(lunout,*) 'Houston, we have a problem ', ierr
+    stop 1
+endif
+
+END SUBROUTINE abort_pem
+
+END MODULE abort_pem_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3075)
+++ /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3076)
@@ -89,2 +89,10 @@
 == 05/10/2023 == JBC
 Correction of a bug: the variable 'g' was not correctly initialized in 1D. A little of code cleaning.
+
+== 06/10/2023 == JBC
+Big cleaning/improvements of the PEM:
+    - Conversion of "abort_pem.F" and "soil_settings_PEM.F" into Fortran 90;
+    - Transformation of every PEM subroutines into module;
+    - Rewriting of many subroutines with modern Fortran syntax;
+    - Correction of a bug in "pem.F90" when calling 'recomp_tend_co2_slope'. The arguments were given in disorder and emissivity was missing;
+    - Update of "launch_pem.sh" in deftank.
Index: unk/LMDZ.COMMON/libf/evolution/compute_tendencies_slope.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_slope.F90	(revision 3075)
+++ 	(revision )
@@ -1,30 +1,0 @@
-SUBROUTINE compute_tendencies_slope(ngrid,nslope,min_ice_Y1,min_ice_Y2,tendencies_ice)
-
-implicit none
-
-!=======================================================================
-!
-!  Compute the tendencies of the evolution of water ice over the years
-!
-!=======================================================================
-
-!   arguments:
-!   ----------
-!   INPUT
-integer,                       intent(in) :: ngrid, nslope ! # of grid points along longitude/latitude/ total
-real, dimension(ngrid,nslope), intent(in) :: min_ice_Y1    ! LON x LAT field : minimum of water ice at each point for the first year
-real, dimension(ngrid,nslope), intent(in) :: min_ice_Y2    ! LON x LAT field : minimum of water ice at each point for the second year
-
-!   OUTPUT
-real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! physical point field : difference between the minima = evolution of perenial ice
-
-!=======================================================================
-
-!  We compute the difference
-tendencies_ice = min_ice_Y2 - min_ice_Y1
-
-!  If the difference is too small; there is no evolution
-where (abs(tendencies_ice) < 1.e-10) tendencies_ice = 0.
-
-END SUBROUTINE compute_tendencies_slope
-
Index: /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_slope_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_slope_mod.F90	(revision 3076)
+++ /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_slope_mod.F90	(revision 3076)
@@ -0,0 +1,39 @@
+MODULE compute_tendencies_slope_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE compute_tendencies_slope(ngrid,nslope,min_ice_Y1,min_ice_Y2,tendencies_ice)
+
+implicit none
+
+!=======================================================================
+!
+!  Compute the tendencies of the evolution of water ice over the years
+!
+!=======================================================================
+
+!   arguments:
+!   ----------
+!   INPUT
+integer,                       intent(in) :: ngrid, nslope ! # of grid points along longitude/latitude/ total
+real, dimension(ngrid,nslope), intent(in) :: min_ice_Y1    ! LON x LAT field : minimum of water ice at each point for the first year
+real, dimension(ngrid,nslope), intent(in) :: min_ice_Y2    ! LON x LAT field : minimum of water ice at each point for the second year
+
+!   OUTPUT
+real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! physical point field : difference between the minima = evolution of perenial ice
+!=======================================================================
+
+!  We compute the difference
+tendencies_ice = min_ice_Y2 - min_ice_Y1
+
+!  If the difference is too small; there is no evolution
+where (abs(tendencies_ice) < 1.e-10) tendencies_ice = 0.
+
+END SUBROUTINE compute_tendencies_slope
+
+END MODULE compute_tendencies_slope_mod
+
Index: /trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3075)
+++ /trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3076)
@@ -8,29 +8,33 @@
 !=======================================================================
 
-IMPLICIT NONE
+implicit none
 
-CONTAINS
+!=======================================================================
+contains
+!=======================================================================
 
-  SUBROUTINE conf_pem(i_myear,n_myears)
+SUBROUTINE conf_pem(i_myear,n_myears)
 
 #ifdef CPP_IOIPSL
-  use IOIPSL, only: getin
+    use IOIPSL,          only: getin
 #else
-  ! if not using IOIPSL, we still need to use (a local version of) getin
-  use ioipsl_getincom, only: getin
+    ! if not using IOIPSL, we still need to use (a local version of) getin
+    use ioipsl_getincom, only: getin
 #endif
-  
-  use time_evol_mod,  only: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, Max_iter_pem, &
-                            evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
-  use comsoil_h_pem,  only: soil_pem, fluxgeo, water_reservoir_nom, depth_breccia, depth_bedrock, reg_thprop_dependp 
-  use adsorption_mod, only: adsorption_pem
-  use glaciers_mod,   only: co2glaciersflow, h2oglaciersflow
-  use ice_table_mod,  only: icetable_equilibrium, icetable_dynamic 
 
-  integer, intent(out) :: i_myear, n_myears ! Current simulated Martian year and maximum number of Martian years to be simulated
+use time_evol_mod,  only: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, Max_iter_pem, &
+                          evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
+use comsoil_h_pem,  only: soil_pem, fluxgeo, water_reservoir_nom, depth_breccia, depth_bedrock, reg_thprop_dependp
+use adsorption_mod, only: adsorption_pem
+use glaciers_mod,   only: co2glaciersflow, h2oglaciersflow
+use ice_table_mod,  only: icetable_equilibrium, icetable_dynamic
 
-  character(len=20), parameter :: modname ='conf_pem'
-  integer                      :: ierr
-  integer                      :: year_earth_bp_ini ! Initial year (in Earth years) of the simulation of the PEM defined in run.def
+implicit none
+
+integer, intent(out) :: i_myear, n_myears ! Current simulated Martian year and maximum number of Martian years to be simulated
+
+character(20), parameter :: modname ='conf_pem'
+integer                  :: ierr
+integer                  :: year_earth_bp_ini ! Initial year (in Earth years) of the simulation of the PEM defined in run.def
 
 !PEM parameters
@@ -48,106 +52,106 @@
 
 ! #---------- Orbital parameters
-  evol_orbit_pem=.false.
-  call getin('evol_orbit_pem',evol_orbit_pem)
+evol_orbit_pem = .false.
+call getin('evol_orbit_pem',evol_orbit_pem)
 
-  year_bp_ini = 0.
-  call getin('year_earth_bp_ini',year_earth_bp_ini)
-  year_bp_ini = year_earth_bp_ini/convert_years
+year_bp_ini = 0.
+call getin('year_earth_bp_ini',year_earth_bp_ini)
+year_bp_ini = year_earth_bp_ini/convert_years
 
-  var_obl = .true.
-  call getin('var_obl',var_obl)
-  write(*,*) 'Does obliquity vary ?',var_obl
+var_obl = .true.
+call getin('var_obl',var_obl)
+write(*,*) 'Does obliquity vary ?',var_obl
 
-  var_ecc = .true.
-  call getin('var_ecc',var_ecc)
-  write(*,*) 'Does excentricity vary ?',var_ecc
+var_ecc = .true.
+call getin('var_ecc',var_ecc)
+write(*,*) 'Does excentricity vary ?',var_ecc
 
-  var_lsp = .true.
-  call getin('var_lsp',var_lsp)
-  write(*,*) 'Does Ls peri vary ?',var_lsp
+var_lsp = .true.
+call getin('var_lsp',var_lsp)
+write(*,*) 'Does Ls peri vary ?',var_lsp
 
 ! #---------- Stopping criteria parameters
-  Max_iter_pem=100000000
-  call getin('Max_iter_pem', Max_iter_pem)
+Max_iter_pem = 100000000
+call getin('Max_iter_pem',Max_iter_pem)
 
-  water_ice_criterion=0.2
-  call getin('water_ice_criterion', water_ice_criterion)
+water_ice_criterion = 0.2
+call getin('water_ice_criterion',water_ice_criterion)
 
-  co2_ice_criterion=0.2
-  call getin('co2_ice_criterion', co2_ice_criterion)
+co2_ice_criterion = 0.2
+call getin('co2_ice_criterion',co2_ice_criterion)
 
-  ps_criterion = 0.15
-  call getin('ps_criterion',ps_criterion)
+ps_criterion = 0.15
+call getin('ps_criterion',ps_criterion)
 
-  dt_pem=1
-  call getin('dt_pem', dt_pem)
+dt_pem = 1
+call getin('dt_pem', dt_pem)
 
 ! #---------- Subsurface parameters
-  soil_pem=.true.
-  call getin('soil_pem', soil_pem)
+soil_pem = .true.
+call getin('soil_pem',soil_pem)
 
-  adsorption_pem = .true.
-  call getin('adsorption_pem',adsorption_pem)
+adsorption_pem = .true.
+call getin('adsorption_pem',adsorption_pem)
 
-  co2glaciersflow = .true.
-  call getin('co2glaciersflow', co2glaciersflow)
+co2glaciersflow = .true.
+call getin('co2glaciersflow',co2glaciersflow)
 
-  h2oglaciersflow = .true.
-  call getin('h2oglaciersflow', h2oglaciersflow)
+h2oglaciersflow = .true.
+call getin('h2oglaciersflow',h2oglaciersflow)
 
-  reg_thprop_dependp = .false.
-  call getin('reg_thprop_dependp',reg_thprop_dependp)
-  write(*,*)  'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
+reg_thprop_dependp = .false.
+call getin('reg_thprop_dependp',reg_thprop_dependp)
+write(*,*)  'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
 
 ! #---------- Layering parameters
-  fluxgeo = 0.
-  call getin('fluxgeo',fluxgeo)
-  write(*,*) 'Flux Geothermal is set to',fluxgeo
-   
-  depth_breccia   = 10.
-  call getin('depth_breccia',depth_breccia)
-  write(*,*) 'Depth of breccia is set to',depth_breccia
+fluxgeo = 0.
+call getin('fluxgeo',fluxgeo)
+write(*,*) 'Flux Geothermal is set to',fluxgeo
 
-  depth_bedrock   = 1000.
-  call getin('depth_bedrock',depth_bedrock)
-  write(*,*) 'Depth of bedrock is set to',depth_bedrock
+depth_breccia = 10.
+call getin('depth_breccia',depth_breccia)
+write(*,*) 'Depth of breccia is set to',depth_breccia
 
-   icetable_equilibrium = .true.
-   call getin('icetable_equilibrium',icetable_equilibrium)
-   write(*,*)  'Is the ice table computed at equilibrium?', icetable_equilibrium
+depth_bedrock = 1000.
+call getin('depth_bedrock',depth_bedrock)
+write(*,*) 'Depth of bedrock is set to',depth_bedrock
 
-   icetable_dynamic = .false.
-   call getin('icetable_dynamic',icetable_dynamic)
-   write(*,*)  'Is the ice table computed with the dynamic method?', icetable_dynamic
-  if ((.not.soil_pem).and.((icetable_equilibrium).or.(icetable_dynamic))) then
-       write(*,*) 'Ice table  must be used when soil_pem = T'
-       call abort_physic(modname,"Ice table  must be used when soil_pem = T",1)
-  endif
+icetable_equilibrium = .true.
+call getin('icetable_equilibrium',icetable_equilibrium)
+write(*,*)  'Is the ice table computed at equilibrium?', icetable_equilibrium
 
-  if ((.not.soil_pem).and.adsorption_pem) then
-       write(*,*) 'Adsorption must be used when soil_pem = T'
-       call abort_physic(modname,"Adsorption must be used when soil_pem = T",1)
-  endif
-  
-  if ((.not.soil_pem).and.(fluxgeo.gt.0.)) then
-       write(*,*) 'Soil is not activated but Flux Geo > 0.'
-       call abort_physic(modname,"Soil is not activated but Flux Geo > 0.",1)
-  endif
- 
-  if ((.not.soil_pem).and.reg_thprop_dependp) then
-     write(*,*) 'Regolith properties vary with Ps only when soil is set to true'
-     call abort_physic(modname,'Regolith properties vary with Ps only when soil is set to true',1)
-  endif
+icetable_dynamic = .false.
+call getin('icetable_dynamic',icetable_dynamic)
+write(*,*)  'Is the ice table computed with the dynamic method?', icetable_dynamic
+if ((.not. soil_pem) .and. ((icetable_equilibrium) .or. (icetable_dynamic))) then
+    write(*,*) 'Ice table  must be used when soil_pem = T'
+    call abort_physic(modname,"Ice table  must be used when soil_pem = T",1)
+endif
 
-  if (evol_orbit_pem .and. year_bp_ini == 0.) then
-     write(*,*)  'You want to follow the file obl_ecc_lsp.asc for changing orb parameters,'
-     write(*,*)  'but you did not specify from which year to start.'
-     call abort_physic(modname,"evol_orbit_pem=.true. but year_bp_ini=0",1)
-  endif 
+if ((.not. soil_pem) .and. adsorption_pem) then
+    write(*,*) 'Adsorption must be used when soil_pem = T'
+    call abort_physic(modname,"Adsorption must be used when soil_pem = T",1)
+endif
 
-  water_reservoir_nom = 1e4
-  call getin('water_reservoir_nom',water_reservoir_nom)
+if ((.not. soil_pem) .and. (fluxgeo > 0.)) then
+    write(*,*) 'Soil is not activated but Flux Geo > 0.'
+    call abort_physic(modname,"Soil is not activated but Flux Geo > 0.",1)
+endif
 
-  END SUBROUTINE conf_pem
+if ((.not. soil_pem) .and. reg_thprop_dependp) then
+    write(*,*) 'Regolith properties vary with Ps only when soil is set to true'
+    call abort_physic(modname,'Regolith properties vary with Ps only when soil is set to true',1)
+endif
+
+if (evol_orbit_pem .and. year_bp_ini == 0.) then
+    write(*,*)  'You want to follow the file obl_ecc_lsp.asc for changing orb parameters,'
+    write(*,*)  'but you did not specify from which year to start.'
+    call abort_physic(modname,"evol_orbit_pem=.true. but year_bp_ini=0",1)
+endif
+
+water_reservoir_nom = 1.e4
+call getin('water_reservoir_nom',water_reservoir_nom)
+
+END SUBROUTINE conf_pem
 
 END MODULE conf_pem_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90	(revision 3075)
+++ /trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90	(revision 3076)
@@ -1,46 +1,44 @@
-    MODULE constants_marspem_mod
+MODULE constants_marspem_mod
 
+implicit none
 
-    IMPLICIT NONE
 ! Duration of a year and day 
-      INTEGER,PARAMETER :: sols_per_my =669 ! Number of Sols per year
-      REAL, PARAMETER :: sec_per_sol=88775.  ! Duration of a sol, in seconds
+integer, parameter :: sols_per_my = 669    ! Number of Sols per year
+real,    parameter :: sec_per_sol = 88775. ! Duration of a sol, in seconds
 
 ! Molecular masses for CO2,H2O and non condensible gaz, following Franz et al. 2017    
-      REAL,PARAMETER :: m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
-      REAL,PARAMETER :: m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
-      REAL,PARAMETER :: m_h2o = 18.01528E-3      ! Molecular weight of h2o (kg/mol)
+real, parameter :: m_co2 = 44.01e-3    ! CO2 molecular mass (kg/mol)   
+real, parameter :: m_noco2 = 33.37e-3  ! Non condensible mol mass (kg/mol)   
+real, parameter :: m_h2o = 18.01528e-3 ! Molecular weight of h2o (kg/mol)
 
 !     Coefficient for Clapeyron law for CO2 condensation temperature (Tco2 = beta/(alpha-log(vmr)),following James et al. 1992
-      REAL,PARAMETER :: alpha_clap_co2 = 23.3494  !Uniteless,James et al. 1992
-      REAL,PARAMETER :: beta_clap_co2 = 3182.48   !Kelvin, James et al. 1992
+real, parameter :: alpha_clap_co2 = 23.3494 ! Uniteless, James et al. 1992
+real, parameter :: beta_clap_co2 = 3182.48  ! Kelvin, James et al. 1992
 
 !     Coefficient for Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005
-      REAL,PARAMETER :: alpha_clap_h2o = 28.9074  ! Uniteless, Murphy and Koop 2005
-      REAL,PARAMETER :: beta_clap_h2o = -6143.7   ! Kelvin, Murphy and Koop 2005
+real, parameter :: alpha_clap_h2o = 28.9074 ! Uniteless, Murphy and Koop 2005
+real, parameter :: beta_clap_h2o = -6143.7  ! Kelvin, Murphy and Koop 2005
 
 !     Density of the regolith (Zent et al., 1995, Buhler and Piqueux 2021)     
-      REAL,PARAMETER ::  rho_regolith = 2000.     ! kg/m^3
+real, parameter :: rho_regolith = 2000. ! kg/m^3
 
 !     Average  Thermal inertia of the surface, breccia, bedrock, following Mellon et al., 2000., Wood et al., 2008
-      REAL,PARAMETER :: TI_regolith_avg = 250.        ! Averaged of the observed thermal inertia for regolith following Mellon et al., 2000[SI]
-      REAL,PARAMETER :: TI_breccia = 750.             ! Thermal inertia of Breccia following Wood 2009 [SI]
-      REAL,PARAMETER :: TI_bedrock = 2300.            ! Thermal inertia of Bedrock following Wood 2009 [SI]
+real, parameter :: TI_regolith_avg = 250. ! Averaged of the observed thermal inertia for regolith following Mellon et al., 2000[SI]
+real, parameter :: TI_breccia = 750.      ! Thermal inertia of Breccia following Wood 2009 [SI]
+real, parameter :: TI_bedrock = 2300.     ! Thermal inertia of Bedrock following Wood 2009 [SI]
 
 !     Porosity of the soil
-      REAL,PARAMETER :: porosity = 0.4                ! porosity of the martian soil, correspond to the value for a random loose packing of monodiperse sphere (Scott, 1960)
+real, parameter :: porosity = 0.4 ! porosity of the martian soil, correspond to the value for a random loose packing of monodiperse sphere (Scott, 1960)
 
 !     Stefan Boltzmann constant
-      REAL,PARAMETER :: sigmaB=5.678E-8 
+real, parameter :: sigmaB = 5.678e-8 
 
 !     Latent heat of CO2
-      REAL,PARAMETER ::  Lco2 =  5.71E5                ! Pilorget and Forget 2016
+real, parameter :: Lco2 =  5.71e5 ! Pilorget and Forget 2016
 
 ! Conversion H2O/CO2 frost to perenial frost and vice versa
-      REAL,PARAMETER :: threshold_water_frost2perenial = 1000. !~ 1m
-      REAL,PARAMETER :: threshold_co2_frost2perenial = 3200.   !~ 2m
+real, parameter :: threshold_water_frost2perenial = 1000. !~ 1m
+real, parameter :: threshold_co2_frost2perenial = 3200.   !~ 2m
 
+END MODULE constants_marspem_mod
 
-      
-    END MODULE constants_marspem_mod
-
Index: /trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3075)
+++ /trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3076)
@@ -26,5 +26,4 @@
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
 
 IMPLICIT NONE
@@ -105,5 +104,6 @@
 subroutine compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
 
-     USE comconst_mod, ONLY: pi,g
+use comconst_mod,  only: pi,g
+use abort_pem_mod, only: abort_pem
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -181,5 +181,6 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-      USE comconst_mod, ONLY: pi
+use comconst_mod,  only: pi
+use abort_pem_mod, only: abort_pem
 
 
Index: unk/LMDZ.COMMON/libf/evolution/info_run_PEM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/info_run_PEM.F90	(revision 3075)
+++ 	(revision )
@@ -1,33 +1,0 @@
-SUBROUTINE info_run_PEM(year_iter,criterion_stop,i_myear,n_myear)
-
-!=======================================================================
-!
-! Purpose: Write in a file called info_run_PEM.txt the reason why the PEM did stop and the number of extrapolation year done
-!          Update the file tmp_PEMyears.txt to count the number of simulated Martian years
-!
-! Author: RV, JBC
-!=======================================================================
-
-  use time_evol_mod, only: convert_years
-
-  IMPLICIT NONE
-
-  integer, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
-  integer, intent(in) :: i_myear, n_myear ! Current simulated Martian year and maximum number of Martian years to be simulated
-
-  logical :: ok
-
-inquire(file = 'info_run_PEM.txt', exist = ok)
-if (ok) then
-    open(12,file = 'info_run_PEM.txt',status = "old",position = "append",action = "write")
-else
-    open(12,file = 'info_run_PEM.txt',status = "new",action = "write")
-endif
-write(12,*) year_iter, i_myear, criterion_stop
-close(12)
-
-open(100,file = 'tmp_PEMyears.txt',status = 'replace')
-write(100,*) i_myear, n_myear, convert_years
-close(100)
-
-END SUBROUTINE info_run_PEM
Index: /trunk/LMDZ.COMMON/libf/evolution/info_run_PEM_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/info_run_PEM_mod.F90	(revision 3076)
+++ /trunk/LMDZ.COMMON/libf/evolution/info_run_PEM_mod.F90	(revision 3076)
@@ -0,0 +1,46 @@
+MODULE info_run_PEM_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE info_run_PEM(year_iter,criterion_stop,i_myear,n_myear)
+
+!=======================================================================
+!
+! Purpose: Write in a file called info_run_PEM.txt the reason why the PEM did stop and the number of extrapolation year done
+!          Update the file tmp_PEMyears.txt to count the number of simulated Martian years
+!
+! Author: RV, JBC
+!=======================================================================
+
+use time_evol_mod, only: convert_years
+
+implicit none
+
+!----- Arguments
+integer, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
+integer, intent(in) :: i_myear, n_myear          ! Current simulated Martian year and maximum number of Martian years to be simulated
+
+!----- Local variables
+logical :: ok
+
+!----- Code
+inquire(file = 'info_run_PEM.txt', exist = ok)
+if (ok) then
+    open(12,file = 'info_run_PEM.txt',status = "old",position = "append",action = "write")
+else
+    open(12,file = 'info_run_PEM.txt',status = "new",action = "write")
+endif
+write(12,*) year_iter, i_myear, criterion_stop
+close(12)
+
+open(100,file = 'tmp_PEMyears.txt',status = 'replace')
+write(100,*) i_myear, n_myear, convert_years
+close(100)
+
+END SUBROUTINE info_run_PEM
+
+END MODULE info_run_PEM_mod
Index: unk/LMDZ.COMMON/libf/evolution/interpolate_TIPEM_TIGCM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/interpolate_TIPEM_TIGCM.F90	(revision 3075)
+++ 	(revision )
@@ -1,38 +1,0 @@
-   subroutine interpolate_TIPEM_TIGCM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_PEM,TI_GCM)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!
-!!! Purpose: Transfer the thermal inertia from the PEM vertical  grid to the GCM vertical grid
-!!! 
-!!! 
-!!! Author: LL
-!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-      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)          !  Thermal inertia in the PEM vertical grid [J/m^2/K/s^{1/2}]
-      real,intent(inout) :: TI_GCM(ngrid,nsoil_GCM,nslope)       ! Thermal inertia in the GCM vertical grid [J/m^2/K/s^{1/2}]
-
-!local variable
-      integer :: ig,islope,iloop                                 ! loop variables
-     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/interpolate_TIPEM_TIGCM_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/interpolate_TIPEM_TIGCM_mod.F90	(revision 3076)
+++ /trunk/LMDZ.COMMON/libf/evolution/interpolate_TIPEM_TIGCM_mod.F90	(revision 3076)
@@ -0,0 +1,39 @@
+MODULE interpolate_TIPEM_TIGCM_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE interpolate_TIPEM_TIGCM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_PEM,TI_GCM)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Transfer the thermal inertia from the PEM vertical  grid to the GCM vertical grid
+!!!
+!!!
+!!! Author: LL
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+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, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM    ! Thermal inertia in the PEM vertical grid [J/m^2/K/s^{1/2}]
+
+real, dimension(ngrid,nsoil_GCM,nslope), intent(inout) :: TI_GCM ! Thermal inertia in the GCM vertical grid [J/m^2/K/s^{1/2}]
+
+!----- Code
+TI_GCM(:,:,:) = TI_PEM(:,:nsoil_GCM,:)
+
+END SUBROUTINE interpolate_TIPEM_TIGCM
+
+END MODULE interpolate_TIPEM_TIGCM_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90	(revision 3075)
+++ /trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90	(revision 3076)
@@ -1,7 +1,8 @@
-module lask_param_mod
+MODULE lask_param_mod
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
 !!! Purpose: Define parameters from Laskar et al., 2004 evolution
-!!! 
+!!!
 !!! Author: RV, JBC
 !!!
@@ -10,45 +11,45 @@
 implicit none
 
-  real,save,allocatable :: yearlask(:) ! year before present from Laskar et al. Tab
-  real,save,allocatable :: obllask(:)  ! obliquity    [deg]
-  real,save,allocatable :: ecclask(:)  ! excentricity [deg]
-  real,save,allocatable :: lsplask(:)  ! ls perihelie [deg]
-  integer, save         :: last_ilask  ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
+real, dimension(:), allocatable, save :: yearlask   ! year before present from Laskar et al. Tab
+real, dimension(:), allocatable, save :: obllask    ! obliquity    [deg]
+real, dimension(:), allocatable, save :: ecclask    ! excentricity [deg]
+real, dimension(:), allocatable, save :: lsplask    ! ls perihelie [deg]
+integer,                         save :: last_ilask ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
 
+!=======================================================================
 contains
+!=======================================================================
 
-  subroutine ini_lask_param_mod
-  
-  implicit none
+SUBROUTINE ini_lask_param_mod
 
-  integer :: nlask ! number of lines in Laskar files
-  integer :: ierr
+implicit none
 
-  nlask = 0
-  open(1,file = 'obl_ecc_lsp.asc')
-  do
-      read(1,*,iostat = ierr)
-      if (ierr /= 0) exit
-      nlask = nlask + 1
-  enddo
-  close(1)
-  allocate(yearlask(nlask))
-  allocate(obllask(nlask))
-  allocate(ecclask(nlask))
-  allocate(lsplask(nlask))
-  
-  end subroutine ini_lask_param_mod
+integer :: nlask ! number of lines in Laskar files
+integer :: ierr
 
+nlask = 0
+open(1,file = 'obl_ecc_lsp.asc')
+do
+    read(1,*,iostat = ierr)
+    if (ierr /= 0) exit
+    nlask = nlask + 1
+enddo
+close(1)
+allocate(yearlask(nlask),obllask(nlask),ecclask(nlask),lsplask(nlask))
 
-  subroutine end_lask_param_mod
+END SUBROUTINE ini_lask_param_mod
 
-  implicit none
+!=======================================================================
 
-  if (allocated(yearlask)) deallocate(yearlask)
-  if (allocated(obllask)) deallocate(obllask)
-  if (allocated(ecclask)) deallocate(ecclask)
-  if (allocated(lsplask)) deallocate(lsplask)
+SUBROUTINE end_lask_param_mod
 
-  end subroutine end_lask_param_mod
-  
-end module lask_param_mod
+implicit none
+
+if (allocated(yearlask)) deallocate(yearlask)
+if (allocated(obllask)) deallocate(obllask)
+if (allocated(ecclask)) deallocate(ecclask)
+if (allocated(lsplask)) deallocate(lsplask)
+
+END SUBROUTINE end_lask_param_mod
+
+END MODULE lask_param_mod
Index: unk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90	(revision 3075)
+++ 	(revision )
@@ -1,86 +1,0 @@
-!
-! $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
-
-!=======================================================================
-!
-! Purpose: Read in the data_GCM_Yr*.nc the number of time step
-!
-! Author: RV
-!=======================================================================
-
-  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_counter>"
-        ierr = nf90_inq_varid (fID, "time_counter", vID)
-        IF (ierr .NE. nf90_noerr) THEN
-          write(*,*)"read_data_GCM: Le champ <time_counter> 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("nb_time_step_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 "time_counter" dimension
-        ierr = nf90_inq_dimid(fID,"time_counter",vID)
-        ierr = nf90_inquire_dimension(fID,vID,len=timelen)
-        ENDIF
-      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)
-
-  write(*,*) "The number of timestep of the PCM run data=", timelen
-
-  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/nb_time_step_GCM_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM_mod.F90	(revision 3076)
+++ /trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM_mod.F90	(revision 3076)
@@ -0,0 +1,100 @@
+MODULE nb_time_step_GCM_mod
+
+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
+
+character(256) :: msg, var, modname
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE nb_time_step_GCM(fichnom,timelen)
+
+implicit none
+
+!=======================================================================
+!
+! Purpose: Read in the data_GCM_Yr*.nc the number of time step
+!
+! Author: RV
+!=======================================================================
+
+include "dimensions.h"
+
+!=======================================================================
+! Arguments:
+character(len = *), intent(in) :: fichnom !--- FILE NAME
+!=======================================================================
+!   Local Variables
+integer        :: iq, fID, vID, idecal, ierr
+integer        :: timelen ! number of times stored in the file
+!-----------------------------------------------------------------------
+modname = "nb_time_step_GCM"
+
+!  Open initial state NetCDF file
+var = fichnom
+call error_msg(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
+
+ierr = nf90_inq_varid (fID, "temps", vID)
+if (ierr /= nf90_noerr) then
+    write(*,*)"read_data_GCM: Le champ <temps> est absent"
+    write(*,*)"read_data_GCM: J essaie <time_counter>"
+    ierr = nf90_inq_varid (fID, "time_counter", vID)
+    if (ierr /= nf90_noerr) then
+        write(*,*)"read_data_GCM: Le champ <time_counter> est absent"
+        write(*,*)"read_data_GCM: J essaie <Time>"
+        ierr = nf90_inq_varid (fID, "Time", vID)
+        if (ierr /= nf90_noerr) then
+            write(*,*)"read_data_GCM: Le champ <Time> est absent"
+            write(*,*)trim(nf90_strerror(ierr))
+            call abort_gcm("nb_time_step_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 "time_counter" dimension
+        ierr = nf90_inq_dimid(fID,"time_counter",vID)
+        ierr = nf90_inquire_dimension(fID,vID,len = timelen)
+    endif
+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 error_msg(NF90_CLOSE(fID),"close",fichnom)
+
+write(*,*) "The number of timestep of the PCM run data=", timelen
+
+END SUBROUTINE nb_time_step_GCM
+
+!=======================================================================
+
+SUBROUTINE error_msg(ierr,typ,nam)
+
+implicit none
+
+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)//">"
+    case default
+        write(*,*) 'There is no message for this error.'
+        error stop
+end select
+call abort_gcm(trim(modname),trim(msg),ierr)
+
+END SUBROUTINE error_msg
+
+END MODULE nb_time_step_GCM_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90	(revision 3075)
+++ /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90	(revision 3076)
@@ -3,5 +3,5 @@
 !=======================================================================
 !
-! Purpose: Compute the maximum number of iteration for PEM based on the 
+! Purpose: Compute the maximum number of iteration for PEM based on the
 ! stopping criterion given by the orbital parameters
 !
@@ -9,7 +9,9 @@
 !=======================================================================
 
-IMPLICIT NONE
-
-CONTAINS
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
 
 SUBROUTINE orbit_param_criterion(i_myear,year_iter_max)
@@ -26,33 +28,33 @@
     use planete_mod, only: e_elips, obliquit, lsperi
 #endif
-  use time_evol_mod,  only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
-  use comconst_mod,   only: pi
-  use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, ini_lask_param_mod, last_ilask
-
-  IMPLICIT NONE
+use time_evol_mod,  only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
+use comconst_mod,   only: pi
+use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, ini_lask_param_mod, last_ilask
+
+implicit none
 
 !--------------------------------------------------------
 ! Input Variables
 !--------------------------------------------------------
-  integer, intent(in) :: i_myear ! Martian year of the simulation
+integer, intent(in) :: i_myear ! Martian year of the simulation
 
 !--------------------------------------------------------
 ! Output Variables
 !--------------------------------------------------------
-  integer, intent(out) :: year_iter_max ! Maximum number of iteration of the PEM before orb changes too much
-
-!--------------------------------------------------------
-! Local variables 
-!--------------------------------------------------------
-  integer :: Year                                           ! Year of the simulation
-  real    :: Lsp                                            ! Ls perihelion
-  real    :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change
-  real    :: max_obl, max_ecc, max_lsp                      ! Maximum value of orbit param given the admissible change
-  real    :: min_obl, min_ecc, min_lsp                      ! Minimum value of orbit param given the admissible change
-  integer :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
-  integer :: ilask, iylask                                  ! Loop variables
-  real    :: xa, xb, ya, yb                                 ! Variables for interpolation
-  logical :: found_obl, found_ecc, found_lsp                ! Flag variables for orbital parameters
-  logical :: crossed_modulo_d, crossed_modulo_i             ! Flag variables for modulo of Lsp
+integer, intent(out) :: year_iter_max ! Maximum number of iteration of the PEM before orb changes too much
+
+!--------------------------------------------------------
+! Local variables
+!--------------------------------------------------------
+integer :: Year                                           ! Year of the simulation
+real    :: Lsp                                            ! Ls perihelion
+real    :: max_change_obl, max_change_ecc, max_change_lsp ! Maximum admissible change
+real    :: max_obl, max_ecc, max_lsp                      ! Maximum value of orbit param given the admissible change
+real    :: min_obl, min_ecc, min_lsp                      ! Minimum value of orbit param given the admissible change
+integer :: max_obl_iter, max_ecc_iter, max_lsp_iter       ! Maximum year iteration before reaching an inadmissible value of orbit param
+integer :: ilask, iylask                                  ! Loop variables
+real    :: xa, xb, ya, yb                                 ! Variables for interpolation
+logical :: found_obl, found_ecc, found_lsp                ! Flag variables for orbital parameters
+logical :: crossed_modulo_d, crossed_modulo_i             ! Flag variables for modulo of Lsp
 
 ! **********************************************************************
@@ -84,5 +86,5 @@
 endif
 
-! Constant max change case 
+! Constant max change case
 max_change_obl = 0.6
 max_change_ecc = 2.8e-3
@@ -106,5 +108,5 @@
 write(*,*) 'Maximum Lsp accepted =', max_lsp
 write(*,*) 'Minimum Lsp accepted =', min_lsp
-! End Constant max change case 
+! End Constant max change case
 
 ! If we do not want some orb parameter to change, they should not be a stopping criterion,
@@ -175,5 +177,5 @@
             endif
         else
-            if (crossed_modulo_d) then ! If decreasing Lsp "crossed" the modulo but min_lsp has not been met yet 
+            if (crossed_modulo_d) then ! If decreasing Lsp "crossed" the modulo but min_lsp has not been met yet
                 ya = ya - 360.
                 yb = yb - 360.
@@ -208,5 +210,5 @@
 END SUBROUTINE orbit_param_criterion
 
-!********************************************************************************   
+!***********************************************************************
 
 END MODULE orbit_param_criterion_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3075)
+++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3076)
@@ -27,36 +27,46 @@
 PROGRAM pem
 
-use phyetat0_mod,               only: phyetat0
-use phyredem,                   only: physdem0, physdem1
-use netcdf,                     only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
-use turb_mod,                   only: q2, wstar
-use comslope_mod,               only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
-use logic_mod,                  only: iflag_phys
-use mod_const_mpi,              only: COMM_LMDZ
+use phyetat0_mod,                 only: phyetat0
+use phyredem,                     only: physdem0, physdem1
+use netcdf,                       only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
+use turb_mod,                     only: q2, wstar
+use comslope_mod,                 only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
+use logic_mod,                    only: iflag_phys
+use mod_const_mpi,                only: COMM_LMDZ
 use infotrac
-use geometry_mod,               only: latitude_deg
-use conf_pem_mod,               only: conf_pem
-use pemredem,                   only: pemdem0, pemdem1
-use glaciers_mod,               only: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow
-use criterion_pem_stop_mod,     only: criterion_waterice_stop, criterion_co2_stop
-use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
-                                      m_noco2, threshold_water_frost2perenial, threshold_co2_frost2perenial
-use evol_co2_ice_s_mod,         only: evol_co2_ice_s
-use evol_h2o_ice_s_mod,         only: evol_h2o_ice_s
-use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
-                                      TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
-                                      tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
-                                      fluxgeo,                          & ! Geothermal flux for the PEM and GCM
-                                      water_reservoir                     ! Water ressources
-use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
-                                      ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
-                                      co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
-use time_evol_mod,              only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
-use orbit_param_criterion_mod,  only: orbit_param_criterion
-use recomp_orb_param_mod,       only: recomp_orb_param
-use ice_table_mod,              only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
-                                      ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
-use soil_thermalproperties_mod, only: update_soil_thermalproperties
-use time_phylmdz_mod,           only: daysec, dtphys, day_end
+use geometry_mod,                 only: latitude_deg
+use conf_pem_mod,                 only: conf_pem
+use pemredem,                     only: pemdem0, pemdem1
+use glaciers_mod,                 only: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow
+use criterion_pem_stop_mod,       only: criterion_waterice_stop, criterion_co2_stop
+use constants_marspem_mod,        only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
+                                        m_noco2, threshold_water_frost2perenial, threshold_co2_frost2perenial
+use evol_co2_ice_s_mod,           only: evol_co2_ice_s
+use evol_h2o_ice_s_mod,           only: evol_h2o_ice_s
+use comsoil_h_PEM,                only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
+                                        TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
+                                        tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
+                                        fluxgeo,                          & ! Geothermal flux for the PEM and GCM
+                                        water_reservoir                     ! Water ressources
+use adsorption_mod,               only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
+                                        ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
+                                        co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
+use time_evol_mod,                only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
+use orbit_param_criterion_mod,    only: orbit_param_criterion
+use recomp_orb_param_mod,         only: recomp_orb_param
+use ice_table_mod,                only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
+                                        ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
+use soil_thermalproperties_mod,   only: update_soil_thermalproperties
+use time_phylmdz_mod,             only: daysec, dtphys, day_end
+use abort_pem_mod,                only: abort_pem
+use soil_settings_PEM_mod,        only: soil_settings_PEM
+use compute_tendencies_slope_mod, only: compute_tendencies_slope
+use info_run_PEM_mod,             only: info_run_PEM
+use interpolate_TIPEM_TIGCM_mod,  only: interpolate_TIPEM_TIGCM
+use nb_time_step_GCM_mod,         only: nb_time_step_GCM
+use pemetat0_mod,                 only: pemetat0
+use read_data_GCM_mod,            only: read_data_GCM
+use recomp_tend_co2_slope_mod,    only: recomp_tend_co2_slope
+use soil_pem_compute_mod,         only: soil_pem_compute
 
 #ifndef CPP_STD
@@ -85,13 +95,10 @@
 
 #ifndef CPP_1D
-    use iniphysiq_mod,    only: iniphysiq
-    use control_mod,      only: iphysiq, day_step, nsplit_phys
-    use comconst_mod,     only: rad, g, cpp, pi
+    use iniphysiq_mod,            only: iniphysiq
+    use control_mod,              only: iphysiq, day_step, nsplit_phys
+    use comconst_mod,             only: rad, g, cpp, pi
 #else
-    use time_phylmdz_mod, only: iphysiq, day_step
-    use comcstfi_h,       only: pi, rad, g, cpp
-#endif
-
-#ifdef CPP_1D
+    use time_phylmdz_mod,         only: iphysiq, day_step
+    use comcstfi_h,               only: pi, rad, g, cpp
     use regular_lonlat_mod,       only: init_regular_lonlat
     use physics_distribution_mod, only: init_physics_distribution
@@ -99,7 +106,8 @@
     use init_testphys1d_mod,      only: init_testphys1d
     use comvert_mod,              only: ap, bp
+    use writerestart1D_mod,       only: writerestart1D
 #endif
 
-IMPLICIT NONE
+implicit none
 
 include "dimensions.h"
@@ -858,6 +866,6 @@
                 Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
                 Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
-                call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
-                call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
+                call soil_pem_compute(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
+                call soil_pem_compute(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
                 tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
                 do ig = 1,ngrid
@@ -916,6 +924,6 @@
 !------------------------
     write(*,*) "Adaptation of the new co2 tendencies to the current pressure"
-    call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, &
-                               global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
+    call recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),emis,vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, &
+                               global_ave_press_GCM,global_ave_press_new)
 
 !------------------------
Index: /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3075)
+++ /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3076)
@@ -1,15 +1,26 @@
-SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness, &
-                      tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
-                      global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,  &
-                      m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
-   
-  use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
-  use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, water_reservoir_nom, depth_breccia, depth_bedrock, index_breccia, index_bedrock
-  use comsoil_h,                  only: volcapa,inertiedat 
-  use adsorption_mod,             only: regolith_adsorption, adsorption_pem
-  use ice_table_mod,              only: computeice_table_equilibrium
-  use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
-  use soil_thermalproperties_mod, only: update_soil_thermalproperties
-  use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
+MODULE pemetat0_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness,   &
+                    tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
+                    global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,                &
+                    m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
+
+use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
+use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, water_reservoir_nom, depth_breccia, depth_bedrock, index_breccia, index_bedrock
+use comsoil_h,                  only: volcapa,inertiedat
+use adsorption_mod,             only: regolith_adsorption, adsorption_pem
+use ice_table_mod,              only: computeice_table_equilibrium
+use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
+use soil_thermalproperties_mod, only: update_soil_thermalproperties
+use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
+use abort_pem_mod,              only: abort_pem
+use soil_pem_compute_mod,       only: soil_pem_compute
+use soil_pem_ini_mod,           only: soil_pem_ini
 
 #ifndef CPP_STD
@@ -20,58 +31,57 @@
 #endif
 
-  implicit none
-
-  include "callkeys.h"
-
-  character(len=*), intent(in) :: filename                    ! name of the startfi_PEM.nc
-  integer,intent(in) :: ngrid                                 ! # of physical grid points
-  integer,intent(in) :: nsoil_GCM                             ! # of vertical grid points in the GCM
-  integer,intent(in) :: nsoil_PEM                             ! # of vertical grid points in the PEM
-  integer,intent(in) :: nslope                                ! # of sub-grid slopes
-  integer,intent(in) :: timelen                               ! # time samples
-  real, intent(in) :: tsurf_ave_yr1(ngrid,nslope)             ! surface temperature at the first year of GCM call [K]
-  real,intent(in) :: tsurf_ave_yr2(ngrid,nslope)              ! surface temperature at the second  year of GCM call [K]
-  real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
-  real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
-  real,intent(in) :: ps_inst(ngrid,timelen)                   ! surface pressure [Pa]
-  real,intent(in) :: timestep                                 ! time step [s]
-  real,intent(in) :: tend_h2oglaciers(ngrid,nslope)           ! tendencies on h2o glaciers 
-  real,intent(in) :: tend_co2glaciers(ngrid,nslope)           ! tendencies on co2 glaciers
-  real,intent(in) :: co2ice(ngrid,nslope)                     ! co2 ice amount [kg/m^2]
-  real,intent(in) :: waterice(ngrid,nslope)                   ! water ice amount [kg/m^2]
-  real, intent(in) :: watersurf_ave(ngrid,nslope)             ! surface water ice density, yearly averaged  (kg/m^3)
-  real, intent(inout) :: watersoil_ave(ngrid,nsoil_PEM,nslope)! surface water ice density, yearly averaged (kg/m^3)
-  real, intent(in) :: global_ave_pressure                     ! mean average pressure on the planet [Pa]
+implicit none
+
+include "callkeys.h"
+
+character(len=*),               intent(in) :: filename            ! name of the startfi_PEM.nc
+integer,                        intent(in) :: ngrid               ! # of physical grid points
+integer,                        intent(in) :: nsoil_GCM           ! # of vertical grid points in the GCM
+integer,                        intent(in) :: nsoil_PEM           ! # of vertical grid points in the PEM
+integer,                        intent(in) :: nslope              ! # of sub-grid slopes
+integer,                        intent(in) :: timelen             ! # time samples
+real,                           intent(in) :: timestep            ! time step [s]
+real,                           intent(in) :: global_ave_pressure ! mean average pressure on the planet [Pa]
+real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr1       ! surface temperature at the first year of GCM call [K]
+real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr2       ! surface temperature at the second  year of GCM call [K]
+real, dimension(ngrid,timelen), intent(in) :: q_co2               ! MMR tracer co2 [kg/kg]
+real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
+real, dimension(ngrid,timelen), intent(in) :: ps_inst             ! surface pressure [Pa]
+real, dimension(ngrid,nslope),  intent(in) :: tend_h2oglaciers    ! tendencies on h2o glaciers
+real, dimension(ngrid,nslope),  intent(in) :: tend_co2glaciers    ! tendencies on co2 glaciers
+real, dimension(ngrid,nslope),  intent(in) :: co2ice              ! co2 ice amount [kg/m^2]
+real, dimension(ngrid,nslope),  intent(in) :: waterice            ! water ice amount [kg/m^2]
+real, dimension(ngrid,nslope),  intent(in) :: watersurf_ave       ! surface water ice density, yearly averaged  (kg/m^3)
 ! outputs
-
-  real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope)                ! soil (mid-layer) thermal inertia in the PEM grid [SI]
-  real,intent(inout) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)             ! soil (mid-layer) temperature [K]
-  real,intent(inout) :: ice_table(ngrid,nslope)                       ! Ice table depth [m]
-  real,intent(inout) :: ice_table_thickness(ngrid,nslope)             ! Ice table thickness [m]
-  real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,nslope,timelen)    ! instantaneous soil (mid-layer) temperature [K]
-  real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of co2 adsorbed [kg/m^2]
-  real,intent(out) :: deltam_co2_regolith_phys(ngrid)                 ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
-  real,intent(inout) :: m_h2o_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of h2o adsorbed [kg/m^2]
-  real,intent(out) :: deltam_h2o_regolith_phys(ngrid)                 ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
-  real,intent(inout) :: water_reservoir(ngrid)                        ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
+real, dimension(ngrid),                          intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
+real, dimension(ngrid),                          intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
+real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
+real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
+real, dimension(ngrid,nslope),                   intent(inout) :: ice_table           ! Ice table depth [m]
+real, dimension(ngrid,nslope),                   intent(inout) :: ice_table_thickness ! Ice table thickness [m]
+real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst          ! instantaneous soil (mid-layer) temperature [K]
+real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
+real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
+real, dimension(ngrid),                          intent(inout) :: water_reservoir     ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
+real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_ave       ! surface water ice density, yearly averaged (kg/m^3)
 ! local
-   real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)                     ! soil temperature saved in the start [K]
-   real :: TI_startPEM(ngrid,nsoil_PEM,nslope)                        ! soil thermal inertia saved in the start [SI]
-   LOGICAL :: found                                                   ! check if variables are found in the start
-   LOGICAL :: found2                                                  ! check if variables are found in the start
-   integer :: iloop,ig,islope,it,isoil                                ! index for loops
-   real :: kcond                                                      ! Thermal conductivity, intermediate variable [SI]
-   real :: delta                                                      ! Depth of the interface regolith-breccia, breccia -bedrock [m]
-   CHARACTER*2 :: num                                                 ! intermediate string to read PEM start sloped variables
-   real :: tsoil_saved(ngrid,nsoil_PEM)                               ! saved soil temperature [K]
-   real :: tsoil_tmp_yr1(ngrid,nsoil_PEM,nslope)                      ! intermediate soil temperature during yr1[K]
-   real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope)                      ! intermediate soil temperature during yr 2 [K]
-   real :: alph_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computation []
-   real :: beta_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computatio []                               
-   LOGICAL :: startpem_file                                           ! boolean to check if we read the startfile or not
-#ifdef CPP_STD   
-   logical :: watercaptag(ngrid)
-
-   watercaptag(:)=.false.
+real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM               ! soil temperature saved in the start [K]
+real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM                  ! soil thermal inertia saved in the start [SI]
+logical                                 :: found                        ! check if variables are found in the start
+logical                                 :: found2                       ! check if variables are found in the start
+integer                                 :: iloop, ig, islope, it, isoil ! index for loops
+real                                    :: kcond                        ! Thermal conductivity, intermediate variable [SI]
+real                                    :: delta                        ! Depth of the interface regolith-breccia, breccia -bedrock [m]
+character(2)                            :: num                          ! intermediate string to read PEM start sloped variables
+real, dimension(ngrid,nsoil_PEM)        :: tsoil_saved                  ! saved soil temperature [K]
+real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                ! intermediate soil temperature during yr1[K]
+real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                ! intermediate soil temperature during yr 2 [K]
+real, dimension(ngrid,nsoil_PEM - 1)    :: alph_tmp                     ! Intermediate for tsoil computation []
+real, dimension(ngrid,nsoil_PEM - 1)    :: beta_tmp                     ! Intermediate for tsoil computatio []
+logical                                 :: startpem_file                ! boolean to check if we read the startfile or not
+
+#ifdef CPP_STD
+   logical, dimension(ngrid) :: watercaptag
+   watercaptag(:) = .false.
 #endif
 
@@ -81,5 +91,5 @@
 !!!
 !!! Order: 0. Previous year of the PEM run
-!!!        1. Thermal Inertia 
+!!!        1. Thermal Inertia
 !!!        2. Soil Temperature
 !!!        3. Ice table
@@ -99,53 +109,51 @@
 
 !1. Run
-
 if (startpem_file) then
-   ! open pem initial state file:
-   call open_startphy(filename)
+    ! open pem initial state file:
+    call open_startphy(filename)
 
     if (soil_pem) then
- 
+
 !1. Thermal Inertia
 ! a. General case
-
-DO islope=1,nslope
-  write(num,fmt='(i2.2)') islope
-  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,index_breccia,islope).lt.TI_breccia) then 
-!!! transition
-             delta = depth_breccia
-             TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
-            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
-                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
-            do iloop=index_breccia+2,index_bedrock
-              TI_PEM(ig,iloop,islope) = TI_breccia
-            enddo      
-          else ! we keep the high ti values
-            do iloop=index_breccia+1,index_bedrock
-              TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
-            enddo 
-          endif ! TI PEM and breccia comparison
-!! transition
-          delta = depth_bedrock
-          TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
-               (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
-               ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
-          do iloop=index_bedrock+2,nsoil_PEM
-            TI_PEM(ig,iloop,islope) = TI_bedrock
-          enddo   
-      enddo
-   else ! found
-     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
-
-write(*,*) 'PEMETAT0: THERMAL INERTIA DONE'
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        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,index_breccia,islope) < TI_breccia) then
+                    !!! transition
+                    delta = depth_breccia
+                    TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
+                                                        (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
+                                                        ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
+                    do iloop=index_breccia+2,index_bedrock
+                        TI_PEM(ig,iloop,islope) = TI_breccia
+                    enddo
+                else ! we keep the high ti values
+                    do iloop = index_breccia + 1,index_bedrock
+                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
+                    enddo
+                endif ! TI PEM and breccia comparison
+                !! transition
+                delta = depth_bedrock
+                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
+                                                        (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
+                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
+                do iloop = index_bedrock + 2,nsoil_PEM
+                    TI_PEM(ig,iloop,islope) = TI_bedrock
+                enddo
+            enddo
+        else ! found
+            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
+
+    write(*,*) 'PEMETAT0: THERMAL INERTIA doNE'
 
 ! b. Special case for inertiedat, inertiedat_PEM
@@ -154,14 +162,14 @@
  do iloop = 1,nsoil_GCM
    inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
- enddo 
+ enddo
 !!! zone de transition
 delta = depth_breccia
 do ig = 1,ngrid
-if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then 
+if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
 inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
             (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
                       ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
 
- do iloop = index_breccia+2,index_bedrock 
+ do iloop = index_breccia+2,index_bedrock
        inertiedat_PEM(ig,iloop) = TI_breccia
   enddo
@@ -170,5 +178,5 @@
    do iloop=index_breccia+1,index_bedrock
       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
-   enddo 
+   enddo
 endif ! comparison ti breccia
 enddo!ig
@@ -189,7 +197,7 @@
 endif ! not found
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !2. Soil Temperature
-DO islope=1,nslope
+do islope=1,nslope
   write(num,fmt='(i2.2)') islope
    call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
@@ -199,5 +207,5 @@
 !      do ig = 1,ngrid
 !        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
-!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))   
+!        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
 !       do iloop=index_breccia+2,index_bedrock
 !            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
@@ -205,5 +213,5 @@
 !      enddo
 !        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
-!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))  
+!        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
 !
 !      do iloop=index_bedrock+2,nsoil_PEM
@@ -215,7 +223,7 @@
 
      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-     call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-
-    
+     call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+
+
    else
 ! predictor corrector: restart from year 1 of the GCM and build the evolution of
@@ -223,8 +231,8 @@
 
     tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
-    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
-    call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
-    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
-    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
+    call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
+    call soil_pem_compute(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
+    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
+    call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
 
 
@@ -244,9 +252,9 @@
         enddo
       enddo
-ENDDO ! islope
-
-write(*,*) 'PEMETAT0: SOIL TEMP  DONE'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+enddo ! islope
+
+write(*,*) 'PEMETAT0: SOIL TEMP  doNE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !3. Ice Table
   call get_field("ice_table",ice_table,found)
@@ -261,10 +269,10 @@
    endif
 
-write(*,*) 'PEMETAT0: ICE TABLE  DONE'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+write(*,*) 'PEMETAT0: ICE TABLE  doNE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !4. CO2 & H2O Adsorption
  if(adsorption_pem) then
-  DO islope=1,nslope
+  do islope=1,nslope
    write(num,fmt='(i2.2)') islope
    call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
@@ -273,8 +281,8 @@
        exit
     endif
-    
-  ENDDO
-
-  DO islope=1,nslope
+
+  enddo
+
+  do islope=1,nslope
    write(num,fmt='(i2.2)') islope
    call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
@@ -283,6 +291,6 @@
       exit
     endif
-    
-  ENDDO
+
+  enddo
 
     call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
@@ -293,13 +301,13 @@
     endif
     if((.not.found2)) then
-       deltam_h2o_regolith_phys(:) = 0.  
+       deltam_h2o_regolith_phys(:) = 0.
     endif
-write(*,*) 'PEMETAT0: CO2 & H2O adsorption DONE'
+write(*,*) 'PEMETAT0: CO2 & H2O adsorption doNE'
     endif
 endif ! soil_pem
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !. 5 water reservoir
-#ifndef CPP_STD   
+#ifndef CPP_STD
    call get_field("water_reservoir",water_reservoir,found)
    if(.not.found) then
@@ -326,5 +334,5 @@
       do ig = 1,ngrid
 
-          if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then 
+          if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then
 !!! transition
              delta = depth_breccia
@@ -335,9 +343,9 @@
           do iloop=index_breccia+2,index_bedrock
             TI_PEM(ig,iloop,islope) = TI_breccia
-         enddo      
+         enddo
          else ! we keep the high ti values
            do iloop=index_breccia+1,index_bedrock
                   TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
-           enddo 
+           enddo
          endif
 
@@ -349,5 +357,5 @@
           do iloop=index_bedrock+2,nsoil_PEM
             TI_PEM(ig,iloop,islope) = TI_bedrock
-         enddo   
+         enddo
       enddo
 enddo
@@ -357,7 +365,7 @@
  enddo
 !!! zone de transition
-delta = depth_breccia 
+delta = depth_breccia
 do ig = 1,ngrid
-      if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then 
+      if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
 inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
             (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
@@ -365,11 +373,11 @@
 
 
- do iloop = index_breccia+2,index_bedrock 
+ do iloop = index_breccia+2,index_bedrock
 
        inertiedat_PEM(ig,iloop) = TI_breccia
-   
+
  enddo
-else 
-   do iloop = index_breccia+1,index_bedrock 
+else
+   do iloop = index_breccia+1,index_bedrock
        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
     enddo
@@ -394,8 +402,8 @@
  enddo
 
-write(*,*) 'PEMETAT0: TI DONE'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!b) Soil temperature 
+write(*,*) 'PEMETAT0: TI doNE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!b) Soil temperature
     do islope = 1,nslope
 !     do ig = 1,ngrid
@@ -414,11 +422,11 @@
 !            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
 !      enddo
-     
+
 !       enddo
 
       call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-      call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-
-   
+      call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+
+
     do it = 1,timelen
         do isoil = nsoil_GCM+1,nsoil_PEM
@@ -433,7 +441,7 @@
         enddo
 enddo !islope
-write(*,*) 'PEMETAT0: TSOIL DONE'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+write(*,*) 'PEMETAT0: TSOIL doNE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !c) Ice table
        call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
@@ -442,46 +450,48 @@
          call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
        enddo
-       write(*,*) 'PEMETAT0: Ice table DONE'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       write(*,*) 'PEMETAT0: Ice table doNE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !d) Regolith adsorbed
- if(adsorption_pem) then
+if (adsorption_pem) then
     m_co2_regolith_phys(:,:,:) = 0.
     m_h2o_regolith_phys(:,:,:) = 0.
 
     call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
-                                m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
-   
+                             m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
+
     deltam_co2_regolith_phys(:) = 0.
-    deltam_h2o_regolith_phys(:) = 0.  
- endif
-
-write(*,*) 'PEMETAT0: CO2 adsorption DONE'
+    deltam_h2o_regolith_phys(:) = 0.
+endif
+
+write(*,*) 'PEMETAT0: CO2 adsorption doNE'
 endif !soil_pem
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !. e) water reservoir
-#ifndef CPP_STD   
-      do ig=1,ngrid
-        if(watercaptag(ig)) then
-           water_reservoir(ig)=water_reservoir_nom
+#ifndef CPP_STD
+    do ig = 1,ngrid
+        if (watercaptag(ig)) then
+            water_reservoir(ig) = water_reservoir_nom
         else
-           water_reservoir(ig)=0.
+            water_reservoir(ig) = 0.
         endif
-      enddo
+    enddo
 #endif
 
 endif ! of if (startphy_file)
 
-if(soil_pem) then
-!! Sanity check
-    DO ig = 1,ngrid
-      DO islope = 1,nslope
-        DO iloop = 1,nsoil_PEM
-          if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1) 
-        ENDDO
-      ENDDO
-     ENDDO
-endif!soil_pem
-
-END SUBROUTINE
+if (soil_pem) then ! Sanity check
+    do ig = 1,ngrid
+        do islope = 1,nslope
+            do iloop = 1,nsoil_PEM
+                if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1)
+            enddo
+        enddo
+    enddo
+endif !soil_pem
+
+END SUBROUTINE pemetat0
+
+END MODULE pemetat0_mod
+
Index: unk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 3075)
+++ 	(revision )
@@ -1,354 +1,0 @@
-SUBROUTINE read_data_GCM(fichnom,timelen, iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries, & 
-             min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
-             watersurf_density_ave,watersoil_density)
-
-      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
-      USE comsoil_h_PEM, ONLY: soil_pem
-      use constants_marspem_mod,only: m_co2,m_noco2
-      IMPLICIT NONE
-
-!=======================================================================
-!
-! Purpose: Read initial confitions file from the GCM
-!
-! Authors: RV & LL
-!=======================================================================
-
-  include "dimensions.h"
-
-!===============================================================================
-! Arguments:
-
-  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
-  INTEGER, INTENT(IN) :: timelen                   ! number of times stored in the file
-  INTEGER :: iim_input,jjm_input,ngrid,nslope      ! number of points in the lat x lon dynamical grid, number of subgrid slopes
-! Ouputs
-  REAL, INTENT(OUT) ::  min_co2_ice(ngrid,nslope) ! Minimum of co2 ice  per slope of the year [kg/m^2]
-  REAL, INTENT(OUT) ::  min_h2o_ice(ngrid,nslope) ! Minimum of h2o ice per slope of the year [kg/m^2]
-  REAL, INTENT(OUT)  :: vmr_co2_gcm_phys(ngrid,timelen) ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
-  REAL, INTENT(OUT) ::  ps_timeseries(ngrid,timelen)! Surface Pressure [Pa]
-  REAL, INTENT(OUT) ::  q_co2(ngrid,timelen)        ! CO2 mass mixing ratio in the first layer [kg/m^3]
-  REAL, INTENT(OUT) ::  q_h2o(ngrid,timelen)        ! H2O mass mixing ratio in the first layer [kg/m^3]
-  REAL, INTENT(OUT) ::  co2_ice_slope(ngrid,nslope,timelen) ! co2 ice amount per  slope of the year [kg/m^2]
-!SOIL
-  REAL, INTENT(OUT) ::  tsurf_ave(ngrid,nslope)         ! Average surface temperature of the concatenated file [K]
-  REAL, INTENT(OUT) ::  tsoil_ave(ngrid,nsoilmx,nslope) ! Average soil temperature of the concatenated file [K]
-  REAL ,INTENT(OUT) ::  tsurf_gcm(ngrid,nslope,timelen)                  ! Surface temperature of the concatenated file, time series [K]
-  REAL , INTENT(OUT) ::  tsoil_gcm(ngrid,nsoilmx,nslope,timelen)         ! Soil temperature of the concatenated file, time series [K]
-  REAL , INTENT(OUT) ::  watersurf_density_ave(ngrid,nslope)             ! Water density at the surface [kg/m^3]
-  REAL , INTENT(OUT) ::  watersoil_density(ngrid,nsoilmx,nslope,timelen) ! Water density in the soil layer, time series [kg/m^3]
-!===============================================================================
-!   Local Variables 
-  CHARACTER(LEN=256) :: msg, var, modname               ! for reading
-  INTEGER :: iq, fID, vID, idecal                       ! for reading
-  INTEGER :: ierr                                       ! for reading
-  CHARACTER(len=12) :: start_file_type="earth" ! default start file type
-
-  REAL,ALLOCATABLE :: time(:) ! times stored in start
-  INTEGER :: indextime ! index of selected time
-
-  INTEGER :: edges(4),corner(4)
-  INTEGER :: i,j,l,t                                                   ! loop variables
-  real  ::  A , B, mmean                                           ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
-
-  INTEGER :: islope                                                    ! loop for variables
-  CHARACTER*2 :: num                                                   ! for reading sloped variables
-  REAL ::  h2o_ice_s_dyn(iim_input+1,jjm_input+1,nslope,timelen)       ! h2o ice per slope of the concatenated file [kg/m^2]
-  REAL ::  watercap_slope(iim_input+1,jjm_input+1,nslope,timelen)
-  REAL ::  vmr_co2_gcm(iim_input+1,jjm_input+1,timelen)                ! CO2 volume mixing ratio in the first layer  [mol/m^3]
-  REAL ::  ps_GCM(iim_input+1,jjm_input+1,timelen)                     ! Surface Pressure [Pa]
-  REAL ::  min_co2_ice_dyn(iim_input+1,jjm_input+1,nslope)
-  REAL ::  min_h2o_ice_dyn(iim_input+1,jjm_input+1,nslope)
-  REAL ::  tsurf_ave_dyn(iim_input+1,jjm_input+1,nslope)               ! Average surface temperature of the concatenated file [K]
-  REAL ::  tsoil_ave_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope)       ! Average soil temperature of the concatenated file [K]
-  REAL ::  tsurf_gcm_dyn(iim_input+1,jjm_input+1,nslope,timelen)       ! Surface temperature of the concatenated file, time series [K]
-  REAL ::  tsoil_gcm_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen)! Soil temperature of the concatenated file, time series [K]
-  REAL ::  q_co2_dyn(iim_input+1,jjm_input+1,timelen)                  ! CO2 mass mixing ratio in the first layer [kg/m^3]
-  REAL ::  q_h2o_dyn(iim_input+1,jjm_input+1,timelen)                  ! H2O mass mixing ratio in the first layer [kg/m^3]
-  REAL ::  co2_ice_slope_dyn(iim_input+1,jjm_input+1,nslope,timelen)  ! co2 ice amount per  slope of the year [kg/m^2]
-  REAL ::  watersurf_density_dyn(iim_input+1,jjm_input+1,nslope,timelen)! Water density at the surface, time series [kg/m^3]
-  REAL ::  watersurf_density(ngrid,nslope,timelen)                     ! Water density at the surface, time series [kg/m^3]
-  REAL ::  watersoil_density_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Water density in the soil layer, time series [kg/m^3]
-
-!-----------------------------------------------------------------------
-  modname="read_data_gcm"
-
-      A =(1/m_co2 - 1/m_noco2)
-      B=1/m_noco2
-
-  write(*,*) "Opening ", fichnom, "..."
-
-!  Open initial state NetCDF file
-  var=fichnom
-  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
-
-     write(*,*) "Downloading data for vmr co2..."
-
-  CALL get_var3("co2_layer1"   ,q_co2_dyn)
-
-     write(*,*) "Downloading data for vmr co2 done"
-     write(*,*) "Downloading data for vmr h20..."
-
-  CALL get_var3("h2o_layer1"   ,q_h2o_dyn)
-
-     write(*,*) "Downloading data for vmr h2o done"
-     write(*,*) "Downloading data for surface pressure ..."
-
-  CALL get_var3("ps"   ,ps_GCM)
-
-     write(*,*) "Downloading data for surface pressure done"
-     write(*,*) "nslope=", nslope
-
-if(nslope.gt.1) then
-
-     write(*,*) "Downloading data for co2ice_slope ..."
-
-DO islope=1,nslope
-  write(num,fmt='(i2.2)') islope
-  call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
-ENDDO
-
-     write(*,*) "Downloading data for co2ice_slope done"
-     write(*,*) "Downloading data for h2o_ice_s_slope ..."
-
-DO islope=1,nslope
-  write(num,fmt='(i2.2)') islope
-  call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
-ENDDO
-
-     write(*,*) "Downloading data for h2o_ice_s_slope done"
-
-#ifndef CPP_STD
-
-     write(*,*) "Downloading data for watercap_slope ..."
-DO islope=1,nslope
-       write(num,fmt='(i2.2)') islope
-       call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
-!        watercap_slope(:,:,:,:)= 0.
-ENDDO            
-     write(*,*) "Downloading data for watercap_slope done"
-
-#endif
-    
- write(*,*) "Downloading data for tsurf_slope ..."
-
-DO islope=1,nslope
-  write(num,fmt='(i2.2)') islope
-  call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
-ENDDO
-
-     write(*,*) "Downloading data for tsurf_slope done"
-
-#ifndef CPP_STD
-
-     if(soil_pem) then
-
-     write(*,*) "Downloading data for soiltemp_slope ..."
-
-DO islope=1,nslope
-  write(num,fmt='(i2.2)') islope
-  call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
-ENDDO
-
-     write(*,*) "Downloading data for soiltemp_slope done"
-
-     write(*,*) "Downloading data for watersoil_density ..."
-
-DO islope=1,nslope
-  write(num,fmt='(i2.2)') islope
-  call get_var4("Waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
-ENDDO
-
-     write(*,*) "Downloading data for  watersoil_density  done"
-
-     write(*,*) "Downloading data for  watersurf_density  ..."
-
-DO islope=1,nslope
-  write(num,fmt='(i2.2)') islope
-  call get_var3("Waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
-ENDDO
-
-     write(*,*) "Downloading data for  watersurf_density  done"
-
-  endif !soil_pem
-
-#endif
-
-  else !nslope=1 no slope, we copy all the values
-
-    CALL get_var3("h2o_ice_s", h2o_ice_s_dyn(:,:,1,:))
-    CALL get_var3("co2ice", co2_ice_slope_dyn(:,:,1,:))
-    call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
-#ifndef CPP_STD
-!    call get_var3("watercap", watercap_slope(:,:,1,:))
-  watercap_slope(:,:,1,:)=0.
-  watersurf_density_dyn(:,:,:,:)=0.
-  watersoil_density_dyn(:,:,:,:,:)=0.
-#endif
-
-    if(soil_pem) then
-      call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
-    endif !soil_pem
-  endif !nslope=1
-
-! Compute the minimum over the year for each point
-  write(*,*) "Computing the min of h2o_ice_slope"
-  min_h2o_ice_dyn(:,:,:)=minval(h2o_ice_s_dyn+watercap_slope,4)
-  write(*,*) "Computing the min of co2_ice_slope"
-  min_co2_ice_dyn(:,:,:)=minval(co2_ice_slope_dyn,4)
-
-!Compute averages
-
-    write(*,*) "Computing average of tsurf"
-    tsurf_ave_dyn(:,:,:)=SUM(tsurf_gcm_dyn(:,:,:,:),4)/timelen
-
-#ifndef CPP_1D
-  DO islope = 1,nslope
-    DO t=1,timelen
-      CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
-    ENDDO
-  ENDDO
-#endif
-
-  if(soil_pem) then
-    write(*,*) "Computing average of tsoil"
-    tsoil_ave_dyn(:,:,:,:)=SUM(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
-    write(*,*) "Computing average of watersurf_density"
-    watersurf_density_ave(:,:) = SUM(watersurf_density(:,:,:),3)/timelen
-  endif
-
-! By definition, a density is positive, we get rid of the negative values
-  DO i=1,iim+1
-    DO j = 1, jjm_input+1
-       DO islope=1,nslope
-          if (min_co2_ice_dyn(i,j,islope).LT.0) then
-            min_co2_ice_dyn(i,j,islope)  = 0.
-          endif
-          if (min_h2o_ice_dyn(i,j,islope).LT.0) then
-            min_h2o_ice_dyn(i,j,islope)  = 0.
-          endif
-       ENDDO
-    ENDDO
-  ENDDO
-
-  DO i=1,iim+1
-    DO j = 1, jjm_input+1
-      DO t = 1, timelen
-         if (q_co2_dyn(i,j,t).LT.0) then
-              q_co2_dyn(i,j,t)=1E-10
-         elseif (q_co2_dyn(i,j,t).GT.1) then
-              q_co2_dyn(i,j,t)=1.
-         endif
-         if (q_h2o_dyn(i,j,t).LT.0) then
-              q_h2o_dyn(i,j,t)=1E-30
-         elseif (q_h2o_dyn(i,j,t).GT.1) then
-              q_h2o_dyn(i,j,t)=1.
-         endif
-         mmean=1/(A*q_co2_dyn(i,j,t) +B)
-         vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
-      ENDDO
-    ENDDO
-  ENDDO
-
-#ifndef CPP_1D
-
-     CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
-     call gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,ps_GCM,ps_timeseries)
-     CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,q_co2_dyn,q_co2)
-     CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,q_h2o_dyn,q_h2o)
-
-     DO islope = 1,nslope
-       CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
-       CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
-       if(soil_pem) then
-       DO l=1,nsoilmx
-         CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
-         DO t=1,timelen
-           CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
-           CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
-         ENDDO
-       ENDDO
-       endif !soil_pem
-       DO t=1,timelen
-         CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
-         CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
-       ENDDO
-     ENDDO
-
-     CALL gr_dyn_fi(nslope,iim_input+1,jjm_input+1,ngrid,tsurf_ave_dyn,tsurf_ave)
-
-#else
-
-  vmr_co2_gcm_phys(1,:)=vmr_co2_gcm(1,1,:)
-  ps_timeseries(1,:)=ps_GCM(1,1,:)
-  q_co2(1,:)=q_co2_dyn(1,1,:)
-  q_h2o(1,:)=q_h2o_dyn(1,1,:)
-  min_co2_ice(1,:)=min_co2_ice_dyn(1,1,:)
-  min_h2o_ice(1,:)=min_h2o_ice_dyn(1,1,:)
-  if(soil_pem) then
-    tsoil_ave(1,:,:)=tsoil_ave_dyn(1,1,:,:)
-    tsoil_gcm(1,:,:,:)=tsoil_gcm_dyn(1,1,:,:,:)
-    watersoil_density(1,:,:,:)=watersoil_density_dyn(1,1,:,:,:)
-  endif !soil_pem
-  tsurf_GCM(1,:,:)=tsurf_GCM_dyn(1,1,:,:)
-  co2_ice_slope(1,:,:)=co2_ice_slope_dyn(1,1,:,:)
-  tsurf_ave(1,:)=tsurf_ave_dyn(1,1,:)
-
-#endif
-
-  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_GCM_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM_mod.F90	(revision 3076)
+++ /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM_mod.F90	(revision 3076)
@@ -0,0 +1,368 @@
+MODULE read_data_GCM_mod
+
+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
+
+character(256) :: msg, var, modname ! for reading
+integer        :: fID, vID          ! for reading
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE read_data_GCM(fichnom,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries,           &
+                         min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
+                         watersurf_density_ave,watersoil_density)
+use comsoil_h,             only: nsoilmx
+use comsoil_h_PEM,         only: soil_pem
+use constants_marspem_mod, only: m_co2, m_noco2
+
+implicit none
+
+!=======================================================================
+!
+! Purpose: Read initial confitions file from the GCM
+!
+! Authors: RV & LL
+!=======================================================================
+
+include "dimensions.h"
+
+!=======================================================================
+! Arguments:
+
+character(len = *), intent(in) :: fichnom                             !--- FILE NAME
+integer,            intent(in) :: timelen                             ! number of times stored in the file
+integer                        :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
+! Ouputs
+real, dimension(ngrid,nslope),         intent(out) :: min_co2_ice      ! Minimum of co2 ice  per slope of the year [kg/m^2]
+real, dimension(ngrid,nslope),         intent(out) :: min_h2o_ice      ! Minimum of h2o ice per slope of the year [kg/m^2]
+real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_gcm_phys ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
+real, dimension(ngrid,timelen),        intent(out) :: ps_timeseries    ! Surface Pressure [Pa]
+real, dimension(ngrid,timelen),        intent(out) :: q_co2            ! CO2 mass mixing ratio in the first layer [kg/m^3]
+real, dimension(ngrid,timelen),        intent(out) :: q_h2o            ! H2O mass mixing ratio in the first layer [kg/m^3]
+real, dimension(ngrid,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per  slope of the year [kg/m^2]
+!SOIL
+real, dimension(ngrid,nslope),                 intent(out) :: tsurf_ave             ! Average surface temperature of the concatenated file [K]
+real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_ave             ! Average soil temperature of the concatenated file [K]
+real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_gcm             ! Surface temperature of the concatenated file, time series [K]
+real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_gcm             ! Soil temperature of the concatenated file, time series [K]
+real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_ave ! Water density at the surface [kg/m^3]
+real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
+!===============================================================================
+!   Local Variables
+character(12)                   :: start_file_type = "earth" ! default start file type
+real, dimension(:), allocatable :: time      ! times stored in start
+integer                         :: indextime ! index of selected time
+integer                         :: edges(4), corner(4)
+integer                         :: i, j, l, t          ! loop variables
+real                            ::  A , B, mmean       ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
+integer                         :: islope              ! loop for variables
+character(2)                    :: num                 ! for reading sloped variables
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn         ! h2o ice per slope of the concatenated file [kg/m^2]
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap_slope
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm           ! CO2 volume mixing ratio in the first layer  [mol/m^3]
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                ! Surface Pressure [Pa]
+real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
+real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
+real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_ave_dyn         ! Average surface temperature of the concatenated file [K]
+real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_ave_dyn         ! Average soil temperature of the concatenated file [K]
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_gcm_dyn         ! Surface temperature of the concatenated file, time series [K]
+real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_gcm_dyn         ! Soil temperature of the concatenated file, time series [K]
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn             ! CO2 mass mixing ratio in the first layer [kg/m^3]
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn             ! H2O mass mixing ratio in the first layer [kg/m^3]
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn     ! co2 ice amount per  slope of the year [kg/m^2]
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn ! Water density at the surface, time series [kg/m^3]
+real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn ! Water density in the soil layer, time series [kg/m^3]
+real, dimension(ngrid,nslope,timelen)                               :: watersurf_density     ! Water density at the surface, time series [kg/m^3]
+
+!-----------------------------------------------------------------------
+modname="read_data_gcm"
+
+A = (1/m_co2 - 1/m_noco2)
+B = 1/m_noco2
+
+write(*,*) "Opening ", fichnom, "..."
+
+!  Open initial state NetCDF file
+var = fichnom
+call error_msg(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
+
+write(*,*) "Downloading data for vmr co2..."
+
+call get_var3("co2_layer1"   ,q_co2_dyn)
+
+write(*,*) "Downloading data for vmr co2 done"
+write(*,*) "Downloading data for vmr h20..."
+
+call get_var3("h2o_layer1"   ,q_h2o_dyn)
+
+write(*,*) "Downloading data for vmr h2o done"
+write(*,*) "Downloading data for surface pressure ..."
+
+call get_var3("ps"   ,ps_GCM)
+
+write(*,*) "Downloading data for surface pressure done"
+write(*,*) "nslope=", nslope
+
+if (nslope > 1) then
+    write(*,*) "Downloading data for co2ice_slope ..."
+
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
+    enddo
+    
+    write(*,*) "Downloading data for co2ice_slope done"
+    write(*,*) "Downloading data for h2o_ice_s_slope ..."
+    
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
+    enddo
+
+    write(*,*) "Downloading data for h2o_ice_s_slope done"
+
+#ifndef CPP_STD
+    write(*,*) "Downloading data for watercap_slope ..."
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
+!        watercap_slope(:,:,:,:) = 0.
+    enddo
+    write(*,*) "Downloading data for watercap_slope done"
+#endif
+
+    write(*,*) "Downloading data for tsurf_slope ..."
+
+    do islope=1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
+    enddo
+
+    write(*,*) "Downloading data for tsurf_slope done"
+
+#ifndef CPP_STD
+    if (soil_pem) then
+        write(*,*) "Downloading data for soiltemp_slope ..."
+
+        do islope = 1,nslope
+            write(num,'(i2.2)') islope
+            call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
+        enddo
+
+        write(*,*) "Downloading data for soiltemp_slope done"
+        write(*,*) "Downloading data for watersoil_density ..."
+
+        do islope = 1,nslope
+            write(num,'(i2.2)') islope
+            call get_var4("Waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
+        enddo
+
+        write(*,*) "Downloading data for  watersoil_density  done"
+        write(*,*) "Downloading data for  watersurf_density  ..."
+
+        do islope = 1,nslope
+            write(num,'(i2.2)') islope
+            call get_var3("Waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
+        enddo
+
+        write(*,*) "Downloading data for  watersurf_density  done"
+
+    endif !soil_pem
+#endif
+else !nslope = 1 no slope, we copy all the values
+    call get_var3("h2o_ice_s", h2o_ice_s_dyn(:,:,1,:))
+    call get_var3("co2ice", co2_ice_slope_dyn(:,:,1,:))
+    call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
+    
+#ifndef CPP_STD
+!    call get_var3("watercap", watercap_slope(:,:,1,:))
+    watercap_slope(:,:,1,:) = 0.
+    watersurf_density_dyn(:,:,:,:) = 0.
+    watersoil_density_dyn(:,:,:,:,:) = 0.
+#endif
+
+    if (soil_pem) call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
+endif !nslope=1
+
+! Compute the minimum over the year for each point
+write(*,*) "Computing the min of h2o_ice_slope"
+min_h2o_ice_dyn(:,:,:) = minval(h2o_ice_s_dyn + watercap_slope,4)
+write(*,*) "Computing the min of co2_ice_slope"
+min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn,4)
+
+!Compute averages
+write(*,*) "Computing average of tsurf"
+tsurf_ave_dyn(:,:,:) = sum(tsurf_gcm_dyn(:,:,:,:),4)/timelen
+
+#ifndef CPP_1D
+    do islope = 1,nslope
+        do t = 1,timelen
+            call gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
+        enddo
+    enddo
+#endif
+
+if (soil_pem) then
+    write(*,*) "Computing average of tsoil"
+    tsoil_ave_dyn(:,:,:,:) = sum(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
+    write(*,*) "Computing average of watersurf_density"
+    watersurf_density_ave(:,:) = sum(watersurf_density(:,:,:),3)/timelen
+endif
+
+! By definition, a density is positive, we get rid of the negative values
+where (min_co2_ice_dyn < 0.) min_co2_ice_dyn = 0.
+where (min_h2o_ice_dyn < 0.) min_h2o_ice_dyn = 0.
+
+do i = 1,iim + 1
+    do j = 1,jjm_input + 1
+        do t = 1, timelen
+            if (q_co2_dyn(i,j,t) < 0) then
+                q_co2_dyn(i,j,t) = 1.e-10
+            else if (q_co2_dyn(i,j,t) > 1) then
+                q_co2_dyn(i,j,t) = 1.
+            endif
+            if (q_h2o_dyn(i,j,t) < 0) then
+                q_h2o_dyn(i,j,t) = 1.e-30
+            else if (q_h2o_dyn(i,j,t) > 1) then
+                q_h2o_dyn(i,j,t) = 1.
+            endif
+            mmean = 1/(A*q_co2_dyn(i,j,t) + B)
+            vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
+        enddo
+    enddo
+enddo
+
+#ifndef CPP_1D
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_GCM,ps_timeseries)
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
+
+    do islope = 1,nslope
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
+        if (soil_pem) then
+            do l = 1,nsoilmx
+                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
+                do t=1,timelen
+                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
+                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
+                enddo
+            enddo
+        endif !soil_pem
+        do t = 1,timelen
+            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
+            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
+        enddo
+    enddo
+
+    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_ave_dyn,tsurf_ave)
+#else
+    vmr_co2_gcm_phys(1,:) = vmr_co2_gcm(1,1,:)
+    ps_timeseries(1,:) = ps_GCM(1,1,:)
+    q_co2(1,:) = q_co2_dyn(1,1,:)
+    q_h2o(1,:) = q_h2o_dyn(1,1,:)
+    min_co2_ice(1,:) = min_co2_ice_dyn(1,1,:)
+    min_h2o_ice(1,:) = min_h2o_ice_dyn(1,1,:)
+    if (soil_pem) then
+        tsoil_ave(1,:,:) = tsoil_ave_dyn(1,1,:,:)
+        tsoil_gcm(1,:,:,:) = tsoil_gcm_dyn(1,1,:,:,:)
+        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
+    endif !soil_pem
+    tsurf_GCM(1,:,:) = tsurf_GCM_dyn(1,1,:,:)
+    co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
+    tsurf_ave(1,:) = tsurf_ave_dyn(1,1,:)
+#endif
+
+END SUBROUTINE read_data_gcm
+
+SUBROUTINE check_dim(n1,n2,str1,str2)
+
+implicit none
+
+integer,            intent(in) :: n1, n2
+character(len = *), intent(in) :: str1, str2
+
+character(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)
+endif
+
+END SUBROUTINE check_dim
+
+!=======================================================================
+
+SUBROUTINE get_var1(var,v)
+
+implicit none
+
+character(len = *), intent(in)  :: var
+real, dimension(:), intent(out) :: v
+
+call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
+call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
+
+END SUBROUTINE get_var1
+
+!=======================================================================
+
+SUBROUTINE get_var3(var,v) ! on U grid
+
+implicit none
+
+character(len = *),     intent(in)  :: var
+real, dimension(:,:,:), intent(out) :: v
+
+call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
+call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
+
+END SUBROUTINE get_var3
+
+!=======================================================================
+
+SUBROUTINE get_var4(var,v)
+
+implicit none
+
+character(len = *),       intent(in)  :: var
+real, dimension(:,:,:,:), intent(out) :: v
+
+call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
+call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
+
+END SUBROUTINE get_var4
+
+!=======================================================================
+
+SUBROUTINE error_msg(ierr,typ,nam)
+
+implicit none
+
+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)//">"
+    case default
+        write(*,*) 'There is no message for this error.'
+        error stop
+end select
+call abort_gcm(trim(modname),trim(msg),ierr)
+
+END SUBROUTINE error_msg
+
+END MODULE read_data_GCM_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 3075)
+++ /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 3076)
@@ -6,12 +6,16 @@
 ! Author: RV, JBC
 !=======================================================================
-IMPLICIT NONE
 
-CONTAINS
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
 
 SUBROUTINE recomp_orb_param(i_myear,year_iter)
 
-  use time_evol_mod, only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
-#ifndef CPP_STD      
+use time_evol_mod, only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
+use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, end_lask_param_mod, last_ilask
+#ifndef CPP_STD
     use comconst_mod, only: pi
     use planete_h,    only: e_elips, obliquit, lsperi, periheli, aphelie, p_elips, peri_day, year_day
@@ -21,13 +25,11 @@
 #endif
 
-  use lask_param_mod, only: yearlask, obllask, ecclask, lsplask, end_lask_param_mod, last_ilask
-
-  IMPLICIT NONE
+implicit none
 
 !--------------------------------------------------------
 ! Input Variables
 !--------------------------------------------------------
-  integer, intent(in) :: i_myear   ! Number of simulated Martian years
-  integer, intent(in) :: year_iter ! Number of iterations of the PEM
+integer, intent(in) :: i_myear   ! Number of simulated Martian years
+integer, intent(in) :: year_iter ! Number of iterations of the PEM
 
 !--------------------------------------------------------
@@ -36,23 +38,22 @@
 
 !--------------------------------------------------------
-! Local variables 
+! Local variables
 !--------------------------------------------------------
-  integer :: Year           ! Year of the simulation
-  real    :: Lsp            ! Ls perihelion
-  integer :: ilask          ! Loop variable
-  real    :: halfaxe        ! Million of km
-  real    :: unitastr
-  real    :: xa, xb, ya, yb ! Variables for interpolation
-  logical :: found_year     ! Flag variable
+integer :: Year           ! Year of the simulation
+real    :: Lsp            ! Ls perihelion
+integer :: ilask          ! Loop variable
+real    :: halfaxe        ! Million of km
+real    :: unitastr
+real    :: xa, xb, ya, yb ! Variables for interpolation
+logical :: found_year     ! Flag variable
 
 ! **********************************************************************
 ! 0. Initializations
 ! **********************************************************************
-#ifdef CPP_STD 
-    real aphelie
-    real periheli
+#ifdef CPP_STD
+    real :: aphelie, periheli
 
-    aphelie=apoastr
-    periheli=periastr
+    aphelie = apoastr
+    periheli = periastr
 #endif
 
@@ -111,5 +112,5 @@
 p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr
 
-#ifndef CPP_STD  
+#ifndef CPP_STD
     call call_dayperi(Lsp,e_elips,peri_day,year_day)
 #endif
@@ -122,6 +123,6 @@
 END SUBROUTINE recomp_orb_param
 
-!********************************************************************************   
-      
+!********************************************************************************
+
 END MODULE recomp_orb_param_mod
 
Index: unk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90	(revision 3075)
+++ 	(revision )
@@ -1,56 +1,0 @@
-!
-! $Id $
-!
-SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,co2ice_slope, emissivity_slope, &
-                                 vmr_co2_gcm,vmr_co2_pem,ps_GCM_2,global_ave_press_GCM,global_ave_press_new)
-
-      use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB,Lco2, sols_per_my, sec_per_sol
-
-      IMPLICIT NONE
-
-!=======================================================================
-!
-!  Routine that compute the evolution of the tendencie for co2 ice
-!
-!=======================================================================
-
-!   arguments:
-!   ----------
-
-!   INPUT
-  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) :: global_ave_press_GCM                       ! global averaged pressure at previous timestep
-  REAL, intent(in) :: global_ave_press_new                       ! global averaged pressure at current timestep
-  REAL, intent(in) ::  tendencies_co2_ice_phys_ini(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
-  REAL, intent(in) :: co2ice_slope(ngrid,nslope)                 ! CO2 ice per mesh and sub-grid slope(kg/m^2)
-  REAL, intent(in) :: emissivity_slope(ngrid,nslope)             ! Emissivity per mesh and sub-grid slope(1)
-
-!   OUTPUT
-  REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
-
-!   local:
-!   ----
-
-  INTEGER :: i,t,islope
-  REAL ::  coef, ave
-
-! Evolution of the water ice for each physical point
-  do i=1,ngrid
-    do islope=1,nslope
-      coef=sols_per_my*sec_per_sol*emissivity_slope(i,islope)*sigmaB/Lco2
-      ave=0.
-      if(co2ice_slope(i,islope).gt.1e-4 .and. abs(tendencies_co2_ice_phys(i,islope)).gt.1e-5) then
-        do t=1,timelen
-           ave=ave+(beta_clap_co2/(alpha_clap_co2-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100.)))**4  &
-              -(beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*(global_ave_press_new/global_ave_press_GCM)/100.)))**4
-         enddo
-      endif
-      if(ave.lt.1e-4) ave = 0.
-      tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys_ini(i,islope)-coef*ave/timelen
-    enddo
-  enddo
-
-END SUBROUTINE recomp_tend_co2_slope
Index: /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90	(revision 3076)
+++ /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90	(revision 3076)
@@ -0,0 +1,60 @@
+MODULE recomp_tend_co2_slope_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,co2ice_slope,emissivity_slope, &
+                                 vmr_co2_gcm,vmr_co2_pem,ps_GCM_2,global_ave_press_GCM,global_ave_press_new)
+
+use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB,Lco2, sols_per_my, sec_per_sol
+
+implicit none
+
+!=======================================================================
+!
+!  Routine that compute the evolution of the tendencie for co2 ice
+!
+!=======================================================================
+
+!   arguments:
+!   ----------
+!   INPUT
+integer,                        intent(in) :: timelen, ngrid, nslope
+real, dimension(ngrid,timelen), intent(in) :: vmr_co2_gcm                 ! physical point field : Volume mixing ratio of co2 in the first layer
+real, dimension(ngrid,timelen), intent(in) :: vmr_co2_pem                 ! physical point field : Volume mixing ratio of co2 in the first layer
+real, dimension(ngrid,timelen), intent(in) :: ps_GCM_2                    ! physical point field : Surface pressure in the GCM
+real,                           intent(in) :: global_ave_press_GCM        ! global averaged pressure at previous timestep
+real,                           intent(in) :: global_ave_press_new        ! global averaged pressure at current timestep
+real, dimension(ngrid,nslope),  intent(in) :: tendencies_co2_ice_phys_ini ! physical point field : Evolution of perenial ice over one year
+real, dimension(ngrid,nslope),  intent(in) :: co2ice_slope                ! CO2 ice per mesh and sub-grid slope(kg/m^2)
+real, dimension(ngrid,nslope),  intent(in) :: emissivity_slope            ! Emissivity per mesh and sub-grid slope(1)
+!   OUTPUT
+real, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
+
+!   local:
+!   ----
+integer :: i, t, islope
+real    :: coef, ave
+
+! Evolution of the water ice for each physical point
+do i = 1,ngrid
+    do islope = 1,nslope
+        coef = sols_per_my*sec_per_sol*emissivity_slope(i,islope)*sigmaB/Lco2
+        ave = 0.
+        if (co2ice_slope(i,islope) > 1.e-4 .and. abs(tendencies_co2_ice_phys(i,islope)) > 1.e-5) then
+            do t=1,timelen
+                ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100.)))**4 &
+                      - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*(global_ave_press_new/global_ave_press_GCM)/100.)))**4
+            enddo
+        endif
+        if (ave < 1e-4) ave = 0.
+        tendencies_co2_ice_phys(i,islope) = tendencies_co2_ice_phys_ini(i,islope) - coef*ave/timelen
+    enddo
+enddo
+
+END SUBROUTINE recomp_tend_co2_slope
+
+END MODULE recomp_tend_co2_slope_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90	(revision 3075)
+++ /trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90	(revision 3076)
@@ -1,3 +1,3 @@
-program reshape_XIOS_output
+PROGRAM reshape_XIOS_output
 
 !=======================================================================
@@ -10,224 +10,193 @@
 ! Authors: RV & LL
 !=======================================================================
-    use netcdf
-    implicit none
-    integer :: status, ncid, ncid1, ncid2
-    integer :: nDims, nVars, nGlobalAtts, unlimDimID
-    integer i,j
-
-    integer :: include_parents
-
-    integer, dimension(:),allocatable :: dimids
-    integer, dimension(:),allocatable :: varids
-
-    integer, dimension(:),allocatable :: dimids_2
-    integer, dimension(:),allocatable :: varids_2
-
-    integer, dimension(:),allocatable :: dimid_var
-
-    real, dimension(:), allocatable :: tempvalues_1d
-    real, dimension(:), allocatable :: values_1d
-
-    real, dimension(:,:), allocatable :: tempvalues_2d
-    real, dimension(:,:), allocatable :: values_2d
-
-    real, dimension(:,:,:), allocatable :: tempvalues_3d
-    real, dimension(:,:,:), allocatable :: values_3d
-
-    real, dimension(:,:,:,:), allocatable :: tempvalues_4d
-    real, dimension(:,:,:,:), allocatable :: values_4d
-
-  character*1 str2
-  character*30 :: name_
-  character*30 :: namevar
-  integer  :: xtype_var
-  integer :: len_
-  integer :: len_1,len_2
-  integer :: len_lat, len_lon, len_time, len_soil
-  integer :: dimid_lon, dimid_lat, dimid_time, dimid_soil
-  integer :: dimid_2
-  integer :: numdims
-  integer :: numatts
-  integer :: numyear
-
-DO numyear=1, 2
-write(*,*) 'numyear',numyear
-write(str2(1:1),'(i1.1)') numyear
-!nf90_open                 ! open existing netCDF dataset
-!integer :: ncid, status
-!...
-status = nf90_open(path = "data2reshape"//str2//".nc", mode = nf90_nowrite, ncid = ncid1)
-if(status /= nf90_noerr) call handle_err(status)
-
-status = nf90_create(path = "datareshaped"//str2//".nc", cmode=or(nf90_noclobber,nf90_64bit_offset), ncid = ncid2)
-if(status /= nf90_noerr) call handle_err(status)
-
-status = nf90_inquire(ncid1, ndims, nvars, nglobalatts, unlimdimid)
-if(status /= nf90_noerr) call handle_err(status)
-
-allocate(dimids(ndims))
-allocate(varids(nvars))
-
-allocate(dimids_2(ndims))
-allocate(varids_2(nvars))
-
-status = nf90_inq_dimids(ncid1, ndims, dimids, include_parents)
-if(status /= nf90_noerr) call handle_err(status)
-status = nf90_inq_varids(ncid1, nvars, varids)
-if(status /= nf90_noerr) call handle_err(status)
-
-do i=1,ndims
-  status = nf90_inquire_dimension(ncid1, dimids(i), name_, len_)
-  if(status /= nf90_noerr) call handle_err(status)
-  if(name_.eq."lon" .or. name_.eq."longitude")  then
-     dimid_lon=dimids(i)
-     len_lon=len_
-     len_=len_+1
-  elseif(name_.eq."lat".or. name_.eq."latitude") then
-     dimid_lat=dimids(i)
-     len_lat=len_
-  elseif(name_.eq."time_counter".or. name_.eq. "Time") then
-     dimid_time=dimids(i)
-     len_time=len_
-  elseif(name_.eq."soil_layers".or. name_.eq. "subsurface_layers") then
-     dimid_soil=dimids(i)
-     len_soil=len_
-  endif
-  status = nf90_def_dim(ncid2, name_, len_, dimid_2)
-  if(status /= nf90_noerr) call handle_err(status)
-  dimids_2(i)=dimid_2
+use netcdf
+
+implicit none
+
+integer                               :: state, ncid, ncid1, ncid2, nDims, nVars, nGlobalAtts, unlimDimID
+integer                               :: i, j, include_parents
+integer, dimension(:),    allocatable :: dimids, varids, dimids_2, varids_2, dimid_var
+real, dimension(:),       allocatable :: tempvalues_1d, values_1d
+real, dimension(:,:),     allocatable :: tempvalues_2d, values_2d
+real, dimension(:,:,:),   allocatable :: tempvalues_3d, values_3d
+real, dimension(:,:,:,:), allocatable :: tempvalues_4d, values_4d
+character(1)                          :: str2
+character(30)                         :: name_, namevar
+integer                               :: xtype_var, len_, len_1, len_2, len_lat, len_lon, len_time, len_soil
+integer                               :: dimid_lon, dimid_lat, dimid_time, dimid_soil, dimid_2, numdims, numatts, numyear
+
+do numyear = 1,2
+    write(*,*) 'numyear',numyear
+    write(str2(1:1),'(i1.1)') numyear
+    !nf90_open                 ! open existing netCDF dataset
+    !integer :: ncid, state
+    !...
+    state = nf90_open(path = "data2reshape"//str2//".nc", mode = nf90_nowrite, ncid = ncid1)
+    if (state /= nf90_noerr) call handle_err(state)
+
+    state = nf90_create(path = "datareshaped"//str2//".nc", cmode=or(nf90_noclobber,nf90_64bit_offset), ncid = ncid2)
+    if (state /= nf90_noerr) call handle_err(state)
+
+    state = nf90_inquire(ncid1, ndims, nvars, nglobalatts, unlimdimid)
+    if (state /= nf90_noerr) call handle_err(state)
+
+    allocate(dimids(ndims))
+    allocate(varids(nvars))
+
+    allocate(dimids_2(ndims))
+    allocate(varids_2(nvars))
+
+    state = nf90_inq_dimids(ncid1, ndims, dimids, include_parents)
+    if (state /= nf90_noerr) call handle_err(state)
+    state = nf90_inq_varids(ncid1, nvars, varids)
+    if (state /= nf90_noerr) call handle_err(state)
+
+    do i = 1,ndims
+        state = nf90_inquire_dimension(ncid1, dimids(i), name_, len_)
+        if (state /= nf90_noerr) call handle_err(state)
+        if (name_ == "lon" .or. name_ == "longitude") then
+            dimid_lon = dimids(i)
+            len_lon = len_
+            len_ = len_ + 1
+        elseif (name_ == "lat".or. name_ == "latitude") then
+            dimid_lat=dimids(i)
+            len_lat=len_
+        elseif (name_ == "time_counter".or. name_ ==  "Time") then
+            dimid_time=dimids(i)
+            len_time=len_
+        elseif (name_ == "soil_layers".or. name_ ==  "subsurface_layers") then
+            dimid_soil=dimids(i)
+            len_soil = len_
+        endif
+        state = nf90_def_dim(ncid2, name_,len_,dimid_2)
+        if (state /= nf90_noerr) call handle_err(state)
+        dimids_2(i) = dimid_2
+    enddo
+
+    do i = 1,nvars
+        state = nf90_inquire_variable(ncid1, varids(i),name = namevar,xtype = xtype_var,ndims = numdims,natts = numatts)
+        write(*,*) "namevar00= ", namevar
+        if (state /= nf90_noerr) call handle_err(state)
+        allocate(dimid_var(numdims))
+        state = nf90_inquire_variable(ncid1,varids(i),name = namevar,xtype = xtype_var,ndims = numdims,dimids = dimid_var,natts = numatts)
+        if (state /= nf90_noerr) call handle_err(state)
+        if (numdims == 1) then
+            if (namevar == "lon") then
+                allocate(tempvalues_1d(len_lon))
+                allocate(values_1d(len_lon + 1))
+                state = nf90_get_var(ncid1,varids(i),tempvalues_1d)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_def_var(ncid2,namevar,xtype_var, dimid_var, varids_2(i))
+                if (state /= nf90_noerr) call handle_err(state)
+                values_1d(1:len_lon) = tempvalues_1d(:)
+                values_1d(len_lon + 1) = values_1d(1)
+                state = nf90_enddef(ncid2)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_put_var(ncid2, varids_2(i), values_1d)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_redef(ncid2)
+                if (state /= nf90_noerr) call handle_err(state)
+                deallocate(tempvalues_1d)
+                deallocate(values_1d)
+            else
+                state = nf90_inquire_dimension(ncid1,dimid_var(1),name_,len_)
+                if (state /= nf90_noerr) call handle_err(state)
+                allocate(tempvalues_1d(len_))
+                state = nf90_get_var(ncid1,varids(i),tempvalues_1d)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_def_var(ncid2,namevar,xtype_var,dimid_var,varids_2(i))
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_enddef(ncid2)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_put_var(ncid2, varids_2(i), tempvalues_1d)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_redef(ncid2)
+                if (state /= nf90_noerr) call handle_err(state)
+                deallocate(tempvalues_1d)
+            endif
+        else if (numdims == 2) then
+            if (namevar == "area") then
+                allocate(tempvalues_2d(len_lon,len_lat))
+                allocate(values_2d(len_lon + 1,len_lat))
+                state = nf90_get_var(ncid1,varids(i),tempvalues_2d)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_def_var(ncid2,namevar,xtype_var,dimid_var,varids_2(i))
+                if (state /= nf90_noerr) call handle_err(state)
+                values_2d(1:len_lon,:) = tempvalues_2d(:,:)
+                values_2d(len_lon+1,:) = values_2d(1,:)
+                state = nf90_enddef(ncid2)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_put_var(ncid2,varids_2(i),values_2d)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_redef(ncid2)
+                if (state /= nf90_noerr) call handle_err(state)
+                deallocate(tempvalues_2d)
+                deallocate(values_2d)
+            else
+                state = nf90_inquire_dimension(ncid1,dimid_var(1),name_,len_1)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_inquire_dimension(ncid1,dimid_var(2),name_,len_2)
+                if (state /= nf90_noerr) call handle_err(state)
+                allocate(tempvalues_2d(len_1,len_2))
+                state = nf90_get_var(ncid1, varids(i), tempvalues_2d)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_def_var(ncid2,namevar,xtype_var,dimid_var,varids_2(i))
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_enddef(ncid2)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_put_var(ncid2, varids_2(i), tempvalues_2d)
+                if (state /= nf90_noerr) call handle_err(state)
+                state = nf90_redef(ncid2)
+                if (state /= nf90_noerr) call handle_err(state)
+                deallocate(tempvalues_2d)
+            endif
+        elseif (numdims == 3) then
+            allocate(tempvalues_3d(len_lon,len_lat,len_time))
+            allocate(values_3d(len_lon + 1,len_lat,len_time))
+            state = nf90_get_var(ncid1,varids(i),tempvalues_3d)
+            if (state /= nf90_noerr) call handle_err(state)
+            state = nf90_def_var(ncid2,namevar,xtype_var,dimid_var,varids_2(i))
+            if (state /= nf90_noerr) call handle_err(state)
+            values_3d(1:len_lon,:,:) = tempvalues_3d(:,:,:)
+            values_3d(len_lon+1,:,:) = values_3d(1,:,:)
+            state = nf90_enddef(ncid2)
+            if (state /= nf90_noerr) call handle_err(state)
+            state = nf90_put_var(ncid2, varids_2(i), values_3d)
+            if (state /= nf90_noerr) call handle_err(state)
+            state = nf90_redef(ncid2)
+            if (state /= nf90_noerr) call handle_err(state)
+            deallocate(tempvalues_3d)
+            deallocate(values_3d)
+        else if (numdims == 4) then
+            allocate(tempvalues_4d(len_lon,len_lat,len_soil,len_time))
+            allocate(values_4d(len_lon+1,len_lat,len_soil,len_time))
+            state = nf90_get_var(ncid1, varids(i), tempvalues_4d)
+            if (state /= nf90_noerr) call handle_err(state)
+            state = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
+            if (state /= nf90_noerr) call handle_err(state)
+            state = nf90_enddef(ncid2)
+            values_4d(1:len_lon,:,:,:) = tempvalues_4d(:,:,:,:)
+            values_4d(len_lon+1,:,:,:) = values_4d(1,:,:,:)
+            if (state /= nf90_noerr) call handle_err(state)
+            state = nf90_put_var(ncid2, varids_2(i), values_4d)
+            if (state /= nf90_noerr) call handle_err(state)
+            state = nf90_redef(ncid2)
+            if (state /= nf90_noerr) call handle_err(state)
+            deallocate(tempvalues_4d)
+            deallocate(values_4d)
+        endif
+        deallocate(dimid_var)
+    enddo
+
+    state = nf90_enddef(ncid2)
+    if (state /= nf90_noerr) call handle_err(state)
+    state = nf90_close(ncid1)
+    if (state /= nf90_noerr) call handle_err(state)
+    state = nf90_close(ncid2)
+    if (state /= nf90_noerr) call handle_err(state)
+
+    deallocate(dimids,varids,dimids_2,varids_2)
 enddo
 
-do i=1,nvars
-  status = nf90_inquire_variable(ncid1, varids(i), name=namevar, xtype=xtype_var, ndims = numdims,natts = numatts)
-  write(*,*) "namevar00= ", namevar
-  if(status /= nf90_noerr) call handle_err(status)
-  allocate(dimid_var(numdims))
-  status = nf90_inquire_variable(ncid1, varids(i), name=namevar, xtype=xtype_var, ndims = numdims, dimids=dimid_var, natts = numatts)
-  if(status /= nf90_noerr) call handle_err(status)
-  if(numdims.eq.1) then
-    if(namevar.eq."lon") then
-      allocate(tempvalues_1d(len_lon))
-      allocate(values_1d(len_lon+1))
-      status = nf90_get_var(ncid1, varids(i), tempvalues_1d)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
-      if(status /= nf90_noerr) call handle_err(status)
-      values_1d(1:len_lon)=tempvalues_1d(:)
-      values_1d(len_lon+1)=values_1d(1)
-      status = nf90_enddef(ncid2)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_put_var(ncid2, varids_2(i), values_1d)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_redef(ncid2)
-      if(status /= nf90_noerr) call handle_err(status)
-      deallocate(tempvalues_1d)
-      deallocate(values_1d)  
-    else
-      status = nf90_inquire_dimension(ncid1, dimid_var(1), name_, len_)
-      if(status /= nf90_noerr) call handle_err(status)
-      allocate(tempvalues_1d(len_))
-      status = nf90_get_var(ncid1, varids(i), tempvalues_1d)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_enddef(ncid2)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_put_var(ncid2, varids_2(i), tempvalues_1d)  
-      if(status /= nf90_noerr) call handle_err(status) 
-      status = nf90_redef(ncid2)
-      if(status /= nf90_noerr) call handle_err(status)
-      deallocate(tempvalues_1d)    
-    endif
-  elseif(numdims.eq.2) then
-    if(namevar.eq."area") then
-      allocate(tempvalues_2d(len_lon,len_lat))
-      allocate(values_2d(len_lon+1,len_lat))     
-      status = nf90_get_var(ncid1, varids(i), tempvalues_2d)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
-      if(status /= nf90_noerr) call handle_err(status)
-      values_2d(1:len_lon,:)=tempvalues_2d(:,:)
-      values_2d(len_lon+1,:)=values_2d(1,:)
-      status = nf90_enddef(ncid2)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_put_var(ncid2, varids_2(i), values_2d)   
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_redef(ncid2)
-      if(status /= nf90_noerr) call handle_err(status)
-      deallocate(tempvalues_2d)
-      deallocate(values_2d) 
-    else
-      status = nf90_inquire_dimension(ncid1, dimid_var(1), name_, len_1)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_inquire_dimension(ncid1, dimid_var(2), name_, len_2)
-      if(status /= nf90_noerr) call handle_err(status)
-      allocate(tempvalues_2d(len_1,len_2))
-      status = nf90_get_var(ncid1, varids(i), tempvalues_2d)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_enddef(ncid2)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_put_var(ncid2, varids_2(i), tempvalues_2d) 
-      if(status /= nf90_noerr) call handle_err(status) 
-      status = nf90_redef(ncid2)
-      if(status /= nf90_noerr) call handle_err(status)
-      deallocate(tempvalues_2d)
-    endif
-  elseif(numdims.eq.3) then
-      allocate(tempvalues_3d(len_lon,len_lat,len_time))
-      allocate(values_3d(len_lon+1,len_lat,len_time))
-      status = nf90_get_var(ncid1, varids(i), tempvalues_3d)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
-      if(status /= nf90_noerr) call handle_err(status)
-      values_3d(1:len_lon,:,:)=tempvalues_3d(:,:,:)
-      values_3d(len_lon+1,:,:)=values_3d(1,:,:)
-      status = nf90_enddef(ncid2)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_put_var(ncid2, varids_2(i), values_3d)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_redef(ncid2)
-      if(status /= nf90_noerr) call handle_err(status)
-      deallocate(tempvalues_3d)
-      deallocate(values_3d)
-  elseif(numdims.eq.4) then
-      allocate(tempvalues_4d(len_lon,len_lat,len_soil,len_time))
-      allocate(values_4d(len_lon+1,len_lat,len_soil,len_time))
-      status = nf90_get_var(ncid1, varids(i), tempvalues_4d)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_def_var(ncid2, namevar, xtype_var, dimid_var, varids_2(i))
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_enddef(ncid2)
-      values_4d(1:len_lon,:,:,:)=tempvalues_4d(:,:,:,:)
-      values_4d(len_lon+1,:,:,:)=values_4d(1,:,:,:)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_put_var(ncid2, varids_2(i), values_4d)
-      if(status /= nf90_noerr) call handle_err(status)
-      status = nf90_redef(ncid2)
-      if(status /= nf90_noerr) call handle_err(status)
-      deallocate(tempvalues_4d)
-      deallocate(values_4d)
-  endif
-
-  deallocate(dimid_var)
-enddo
-
-status = nf90_enddef(ncid2)
-if(status /= nf90_noerr) call handle_err(status)
-status = nf90_close(ncid1)
-if(status /= nf90_noerr) call handle_err(status)
-status = nf90_close(ncid2)
-if(status /= nf90_noerr) call handle_err(status)
-
-
-deallocate(dimids)
-deallocate(varids)
-deallocate(dimids_2)
-deallocate(varids_2)
-
-enddo
-
-end program reshape_XIOS_output
-
+END PROGRAM reshape_XIOS_output
+
Index: unk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM.F90	(revision 3075)
+++ 	(revision )
@@ -1,95 +1,0 @@
-      SUBROUTINE soil_TIfeedback_PEM(ngrid,nsoil,icecover,   newtherm_i)
-
-      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.      ! surface water ice thermal inertia [SI]
-      REAL :: rho_ice = 920.            ! density of water ice [kg/m^3]
-      REAL :: prev_thermi(ngrid,nsoil)  ! previous thermal inertia [SI]
-!Inputs
-!------
-
-      REAL ,INTENT(IN):: icecover(ngrid)         ! tracer on the surface (kg.m-2)
-
-!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_TIfeedback_PEM_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM_mod.F90	(revision 3076)
+++ /trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM_mod.F90	(revision 3076)
@@ -0,0 +1,98 @@
+MODULE soil_TIfeedback_PEM_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE soil_TIfeedback_PEM(ngrid,nsoil,icecover,newtherm_i)
+
+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
+!=======================================================================
+
+!Inputs
+!------
+integer,                intent(in) :: ngrid    ! Number of horizontal grid points
+integer,                intent(in) :: nsoil    ! Number of soil layers
+real, dimension(ngrid), intent(in) :: icecover ! tracer on the surface (kg.m-2)
+!Outputs
+!-------
+real,intent(inout) :: newtherm_i(ngrid,nsoil)    ! New soil thermal inertia
+!Local variables
+!---------------
+integer :: ig                       ! Grid point (ngrid)
+integer :: ik                       ! Grid point (nsoil)
+integer :: iref                     ! Ice/Regolith boundary index
+real    :: icedepth                 ! Ice cover thickness (m)
+real    :: inert_h2o_ice = 800.     ! surface water ice thermal inertia [SI]
+real    :: rho_ice = 920.           ! density of water ice [kg/m^3]
+real    :: prev_thermi(ngrid,nsoil) ! previous thermal inertia [SI]
+
+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 >= 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 < 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 >= layer_PEM(ik)) .and. (icedepth < 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 SUBROUTINE soil_TIfeedback_PEM
+
+END MODULE soil_TIfeedback_PEM_mod
Index: unk/LMDZ.COMMON/libf/evolution/soil_pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90	(revision 3075)
+++ 	(revision )
@@ -1,121 +1,0 @@
-      subroutine soil_pem_routine(ngrid,nsoil,firstcall, &
-               therm_i,timestep,tsurf,tsoil)
-
-
-      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
-                          mthermdiff_PEM, thermdiff_PEM, coefq_PEM, &
-                          coefd_PEM, mu_PEM,alph_PEM,beta_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
-!-----------------------------------------------------------------------
-
-#include "dimensions.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 [SI]
-      real,intent(in) :: timestep       ! time step [s]
-      real,intent(in) :: tsurf(ngrid)   ! surface temperature [K]
- 
-! outputs:
-      real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature [K]
-! 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_pem_compute_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_pem_compute_mod.F90	(revision 3076)
+++ /trunk/LMDZ.COMMON/libf/evolution/soil_pem_compute_mod.F90	(revision 3076)
@@ -0,0 +1,115 @@
+MODULE soil_pem_compute_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE soil_pem_compute(ngrid,nsoil,firstcall,therm_i,timestep,tsurf,tsoil)
+
+use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_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
+!-----------------------------------------------------------------------
+
+#include "dimensions.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, dimension(ngrid,nsoil), intent(in) :: therm_i   ! thermal inertia [SI]
+real,                         intent(in) :: timestep  ! time step [s]
+real, dimension(ngrid),       intent(in) :: tsurf     ! surface temperature [K]
+ 
+! outputs:
+real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
+! 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. - 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. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
+    enddo
+enddo
+
+END SUBROUTINE soil_pem_compute
+
+END MODULE soil_pem_compute_mod
Index: unk/LMDZ.COMMON/libf/evolution/soil_pem_ini.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_pem_ini.F90	(revision 3075)
+++ 	(revision )
@@ -1,112 +1,0 @@
-      subroutine soil_pem_ini(ngrid,nsoil,therm_i,tsurf,tsoil)
-
-
-      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
-                          mthermdiff_PEM, thermdiff_PEM, coefq_PEM, &
-                          coefd_PEM, mu_PEM,alph_PEM,beta_PEM,fluxgeo
-      use comsoil_h,only: volcapa
-      implicit none
-
-!-----------------------------------------------------------------------
-!  Author: LL
-!  Purpose: Compute soil temperature using an implict 1st order scheme, stationnary solution
-!  
-!  Note: depths of layers and mid-layers, soil thermal inertia and 
-!        heat capacity are commons in comsoil_PEM.h
-!-----------------------------------------------------------------------
-
-#include "dimensions.h"
-
-!-----------------------------------------------------------------------
-!  arguments
-!  ---------
-!  inputs:
-      integer,intent(in) :: ngrid       ! number of (horizontal) grid-points 
-      integer,intent(in) :: nsoil       ! number of soil layers  
-      real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia [SI]
-      real,intent(in) :: tsurf(ngrid)   ! surface temperature [K]
-
-! outputs:
-      real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature [K]
-! local variables:
-      integer ig,ik,iloop     
-
-! 0. Initialisations and preprocessing step
-    
-! 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.
-	! q_{k+1/2}
-
-      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
-
-
-
-!  1. Compute beta_PEM coefficients 
-! 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
-
-!  2. Compute soil temperatures
-do iloop = 1,10 !just convergence
-! 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
-      
-enddo ! iloop 
-      end
-
Index: /trunk/LMDZ.COMMON/libf/evolution/soil_pem_ini_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_pem_ini_mod.F90	(revision 3076)
+++ /trunk/LMDZ.COMMON/libf/evolution/soil_pem_ini_mod.F90	(revision 3076)
@@ -0,0 +1,108 @@
+MODULE soil_pem_ini_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE soil_pem_ini(ngrid,nsoil,therm_i,tsurf,tsoil)
+
+use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
+use comsoil_h,     only: volcapa
+
+implicit none
+
+!-----------------------------------------------------------------------
+!  Author: LL
+!  Purpose: Compute soil temperature using an implict 1st order scheme, stationnary solution
+!
+!  Note: depths of layers and mid-layers, soil thermal inertia and
+!        heat capacity are commons in comsoil_PEM.h
+!-----------------------------------------------------------------------
+
+#include "dimensions.h"
+
+!-----------------------------------------------------------------------
+!  arguments
+!  ---------
+!  inputs:
+integer,                      intent(in) :: ngrid   ! number of (horizontal) grid-points
+integer,                      intent(in) :: nsoil   ! number of soil layers
+real, dimension(ngrid,nsoil), intent(in) :: therm_i ! thermal inertia [SI]
+real, dimension(ngrid),       intent(in) :: tsurf   ! surface temperature [K]
+
+! outputs:
+real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
+! local variables:
+integer :: ig, ik, iloop
+
+! 0. Initialisations and preprocessing step
+! 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.
+! q_{k+1/2}
+
+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
+
+!  1. Compute beta_PEM coefficients
+! 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. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
+    enddo
+enddo
+
+!  2. Compute soil temperatures
+do iloop = 1,10 !just convergence
+    ! 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. - 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
+
+enddo ! iloop
+
+END SUBROUTINE soil_pem_ini
+
+END MODULE soil_pem_ini_mod
Index: unk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F	(revision 3075)
+++ 	(revision )
@@ -1,113 +1,0 @@
-      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, 
-     &   depth_breccia,depth_bedrock,index_breccia,index_bedrock  
-      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: LL, based on work by  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
-!
-!
-!  The various actions and variable read/initialized are:
-!  1. Read/build layer (and midlayer) depth
-!  2. 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)	!  Thermal inertia  in the GCM [SI]
-      real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope)	! Thermal inertia   in the PEM [SI]
-
-!======================================================================
-! 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
-      ! 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
-
-! 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
-
-
-! 3. Index for subsurface layering
-! ------------------
-      index_breccia= 1
-      do iloop = 1,nsoil_PEM-1
-        if (depth_breccia.ge.layer_PEM(iloop)) then
-          index_breccia=iloop 
-        else
-         exit
-        endif
-      enddo
-
-      index_bedrock= 1
-      do iloop = 1,nsoil_PEM-1
-        if (depth_bedrock.ge.layer_PEM(iloop)) then
-          index_bedrock=iloop 
-        else
-         exit
-        endif
-      enddo
-
-
-      end subroutine soil_settings_PEM
Index: /trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90	(revision 3076)
+++ /trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90	(revision 3076)
@@ -0,0 +1,109 @@
+MODULE soil_settings_PEM_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+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, depth_breccia, depth_bedrock, index_breccia, index_bedrock
+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: LL, based on work by  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
+!
+!  The various actions and variable read/initialized are:
+!  1. Read/build layer (and midlayer) depth
+!  2. 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, dimension(ngrid,nsoil_GCM,nslope), intent(in) :: TI_GCM    ! Thermal inertia  in the GCM [SI]
+
+real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Thermal inertia   in the PEM [SI]
+!=======================================================================
+! 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
+! 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
+
+! 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 > 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
+
+! 3. Index for subsurface layering
+! ------------------
+index_breccia = 1
+do iloop = 1,nsoil_PEM - 1
+    if (depth_breccia >= layer_PEM(iloop)) then
+        index_breccia = iloop
+    else
+        exit
+    endif
+enddo
+
+index_bedrock = 1
+do iloop = 1,nsoil_PEM - 1
+    if (depth_bedrock >= layer_PEM(iloop)) then
+        index_bedrock = iloop
+    else
+        exit
+    endif
+enddo
+
+END SUBROUTINE soil_settings_PEM
+
+END MODULE soil_settings_PEM_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90	(revision 3075)
+++ /trunk/LMDZ.COMMON/libf/evolution/time_evol_mod.F90	(revision 3076)
@@ -1,17 +1,17 @@
 MODULE time_evol_mod
 
-IMPLICIT NONE  
+implicit none
 
-  integer :: year_bp_ini         !     Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
-  integer :: dt_pem              !     Time step used by the PEM in Planetary years
-  real    :: convert_years       !     Conversion ratio from Planetary years to Earth years
-  real    :: water_ice_criterion !     Percentage of change of the surface of water ice sublimating before stopping the PEM 
-  real    :: co2_ice_criterion   !     Percentage of change of the surface of co2 ice sublimating before stopping the PEM 
-  real    :: ps_criterion        !     Percentage of change of averaged surface pressure before stopping the PEM
-  integer :: Max_iter_pem        !     Maximal number of iteration when converging to a steady state, read in evol.def
-  logical :: evol_orbit_pem      !     True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
-  logical :: var_obl             !     True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity 
-  logical :: var_ecc             !     True if we want the PEM to follow obl_ecc_lsp.asc parameters for excenticity
-  logical :: var_lsp             !     True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
+integer :: year_bp_ini         ! Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
+integer :: dt_pem              ! Time step used by the PEM in Planetary years
+real    :: convert_years       ! Conversion ratio from Planetary years to Earth years
+real    :: water_ice_criterion ! Percentage of change of the surface of water ice sublimating before stopping the PEM
+real    :: co2_ice_criterion   ! Percentage of change of the surface of co2 ice sublimating before stopping the PEM
+real    :: ps_criterion        ! Percentage of change of averaged surface pressure before stopping the PEM
+integer :: Max_iter_pem        ! Maximal number of iteration when converging to a steady state, read in evol.def
+logical :: evol_orbit_pem      ! True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
+logical :: var_obl             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity
+logical :: var_ecc             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for excenticity
+logical :: var_lsp             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
 
 END MODULE time_evol_mod
Index: /trunk/LMDZ.MARS/deftank/pem/launch_pem.sh
===================================================================
--- /trunk/LMDZ.MARS/deftank/pem/launch_pem.sh	(revision 3075)
+++ /trunk/LMDZ.MARS/deftank/pem/launch_pem.sh	(revision 3076)
@@ -168,7 +168,7 @@
     mv startfi.nc startfi_evol.nc
     if [ -f "start.nc" ]; then
-        cp start.nc start_evol.nc
+        mv start.nc start_evol.nc
     elif [ -f "start1D.txt" ]; then
-        cp start1D.txt start1D_evol.txt
+        mv start1D.txt start1D_evol.txt
     fi
     ./$exePEM > out_runPEM${iPEM} 2>&1
