Index: /trunk/LMDZ.COMMON/libf/evolution/ComputeCO2Adsorption.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/ComputeCO2Adsorption.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/ComputeCO2Adsorption.F90	(revision 2835)
@@ -16,6 +16,4 @@
  REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)
  REAL,INTENT(INOUT) :: delta_mreg(ngrid)
-
-
   
 ! Constants:
@@ -71,5 +69,4 @@
 enddo
 
-
 ! 2. Check the exchange between the atmosphere and the regolith
 
Index: unk/LMDZ.COMMON/libf/evolution/Example_ice_evol.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/Example_ice_evol.F90	(revision 2834)
+++ 	(revision )
@@ -1,27 +1,0 @@
-         ! Update surface ice distribution to iterate to steady state if requested
-         if(ice_update)then
-
-            do ig=1,ngrid
-
-               delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
-               
-               ! add multiple years of evolution
-               qsurf_hist(ig,igcm_h2o_ice) = qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep 
-
-               ! if ice has gone -ve, set to zero
-               if(qsurf_hist(ig,igcm_h2o_ice).lt.0.0)then
-                  qsurf_hist(ig,igcm_h2o_ice) = 0.0 
-               endif
-
-               ! if ice is seasonal, set to zero (NEW)
-               if(ice_min(ig).lt.0.01)then
-                  qsurf_hist(ig,igcm_h2o_ice) = 0.0 
-               endif
-
-            enddo
-
-            ! enforce ice conservation
-            ice_tot= SUM(qsurf_hist(:,igcm_h2o_ice)*cell_area(:) )/SUM(cell_area(:))
-            qsurf_hist(:,igcm_h2o_ice) = qsurf_hist(:,igcm_h2o_ice)*(icesrf/ice_tot)
-
-         endif
Index: /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 2835)
@@ -3,6 +3,4 @@
   contains
 
-
-
   subroutine regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen, ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, theta_h2o_adsorbded,m_h2o_adsorbed)
 
@@ -11,5 +9,4 @@
 
  implicit none
-
 
  INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen ! size dimension
@@ -20,9 +17,7 @@
  REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)   ! Soil temperature (K)
 
-
 ! output
  REAL,INTENT(OUT) ::  m_h2o_adsorbed(ngrid,nsoil_PEM,nslope)     ! Density of h2o adsorbed (kg/m^3)
  REAL,INTENT(OUT) :: theta_h2o_adsorbded(ngrid,nsoil_PEM,nslope) ! Fraction of the pores occupied by H2O molecules
-
 
 ! constant
@@ -40,8 +35,4 @@
  real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
 
-
-
-
-
 ! local variable
  real :: A,B                                                   ! Used to compute the mean mass above the surface
@@ -54,10 +45,6 @@
  real, allocatable :: pvapor_avg(:)                            ! yearly average vapor pressure above the surface
 
-
-
 ! 0. Some initializations
 
-
-   
      allocate(mass_mean(ngrid,timelen))
      allocate(zplev_mean(ngrid,timelen))
@@ -77,5 +64,4 @@
      enddo
 
-
 ! b. pressure level
      do it = 1,timelen
@@ -92,8 +78,4 @@
      deallocate(zplev_mean)
      deallocate(mass_mean)
-
-
-
-
 
 ! 2. we compute the mass of co2 adsorded in each layer of the meshes  
@@ -110,6 +92,5 @@
         theta_h2o_adsorbded(ig,iloop,islope) = (K*p_sat/(1+K*p_sat))**mu
         m_h2o_adsorbed(ig,iloop,islope) =as*theta_h2o_adsorbded(ig,iloop,islope)*mi*rho_regolith
-      endif
- 
+     endif
    enddo
   enddo
@@ -118,5 +99,4 @@
  RETURN
   end subroutine
-
 
   SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_completesoil,delta_mreg)
@@ -140,6 +120,4 @@
  REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)   ! Density of co2 adsorbed (kg/m^3)
  REAL,INTENT(INOUT) :: delta_mreg(ngrid)                            ! Difference density of co2 adsorbed (kg/m^3)
-
-
   
 ! Constants:
@@ -170,9 +148,5 @@
  real, allocatable :: pco2_avg(:)                              ! yearly averaged
 
-
-
 ! 0. Some initializations
-
-
 
      allocate(mass_mean(ngrid,timelen))
@@ -188,6 +162,4 @@
      delta_mreg(:) = 0.
 
-
-
 !0.1 Look at perenial ice
   do ig = 1,ngrid
@@ -207,5 +179,4 @@
    enddo
 
-
 !   0.2  Compute the partial pressure of CO2
 !a. the molecular mass into the column
@@ -213,5 +184,4 @@
        mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
      enddo
-
 
 ! b. pressure level
@@ -226,5 +196,4 @@
      pco2_avg(:) = sum(pco2(:,:),2)/timelen
 
- 
      deallocate(zplev_mean)
      deallocate(mass_mean)
@@ -234,6 +203,4 @@
 ! 1. Compute the fraction of the pores occupied by H2O
  call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen, ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,theta_h2o_adsorbed, m_h2o_adsorbed)
-
-
 
 ! 2.  we compute the mass of co2 adsorded in each layer of the meshes  
@@ -256,5 +223,4 @@
  enddo
 enddo
-
 
 ! 2. Check the exchange between the atmosphere and the regolith
@@ -278,5 +244,4 @@
    m_co2_completesoil(:,:,:) = dm_co2_regolith_slope(:,:,:)
 
-
 !=======================================================================
       RETURN
Index: /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod.F90	(revision 2835)
@@ -47,5 +47,4 @@
   ENDDO
 
-
 !  We reorganise the difference on the physical grid
   tendencies_h2o_ice_phys(1)=tendencies_h2o_ice(1,1)
@@ -61,5 +60,4 @@
   tendencies_h2o_ice_phys(ig0) = tendencies_h2o_ice(1,jjm_input+1)
 
-
 END SUBROUTINE compute_tendencies
 
Index: /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod_slope.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod_slope.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod_slope.F90	(revision 2835)
@@ -6,5 +6,4 @@
 
       IMPLICIT NONE
-
 
 !=======================================================================
@@ -34,7 +33,5 @@
 !=======================================================================
 
-
 !  We compute the difference
-!  tendencies_h2o_ice(:,:,:)=min_h2o_ice_Y2(:,:,:)-min_h2o_ice_Y1(:,:,:)
 
   DO j=1,jjm_input+1
@@ -46,42 +43,19 @@
   ENDDO
 
-     print *, "jjm_input+1", jjm_input+1
-     print *, "iim_input+1", iim_input+1
-     print *, "nslope+1", nslope+1
-
 !  If the difference is too small; there is no evolution
   DO j=1,jjm_input+1
     DO i = 1, iim_input
        DO islope = 1, nslope
-!         print *, "tendencies_h2o_ice(i,j,islope)LAAA", tendencies_h2o_ice(i,j,islope)
          if(abs(tendencies_h2o_ice(i,j,islope)).LT.1.0E-10) then
             tendencies_h2o_ice(i,j,islope)=0.
          endif
-!         print *, "tendencies_h2o_ice(i,j,islope)HERE", tendencies_h2o_ice(i,j,islope)
        enddo
     ENDDO
   ENDDO
 
-
-!  We reorganise the difference on the physical grid
-  tendencies_h2o_ice_phys(1,:)=tendencies_h2o_ice(1,1,:)
-
-  ig0 = 2
-  DO j=2,jjm_input
-    DO i = 1, iim_input
-       tendencies_h2o_ice_phys(ig0,:)  =tendencies_h2o_ice(i,j,:)
-       ig0= ig0 + 1
-    ENDDO
+  DO islope = 1,nslope
+    CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tendencies_h2o_ice(:,:,islope),tendencies_h2o_ice_phys(:,islope))
   ENDDO
-
-  tendencies_h2o_ice_phys(ig0,:) = tendencies_h2o_ice(1,jjm_input+1,:)
-
-!  print *, "tendencies_h2o_ice_physze", tendencies_h2o_ice_phys(:,:)
-
 
 END SUBROUTINE compute_tendencies_slope
 
-
-
-
-
Index: /trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90	(revision 2835)
@@ -4,5 +4,4 @@
     USE comsoil_h_PEM, only: layer_PEM 
     implicit none 
-
 
     integer,intent(in) :: timelen,ngrid,nslope,nsoil_PEM,nsoil_GCM
@@ -14,5 +13,4 @@
     real,intent(out) :: ice_table(ngrid,nslope)                 ! ice table [m]
 
-
     real :: m_h2o = 18.01528E-3
     real :: m_co2 = 44.01E-3  
@@ -21,5 +19,4 @@
     real :: alpha = -6143.7 
     real :: beta = 29.9074
-
 
     integer ig, islope,isoil,it
@@ -37,6 +34,4 @@
     real,allocatable :: diff_rho(:,:,:)                    ! difference of vapor content
 
-
-
      allocate(rhovapor(ngrid,nslope,timelen))
      allocate(rhovapor_avg(ngrid,nslope))
@@ -45,7 +40,4 @@
      allocate(mass_mean(ngrid,timelen))
      allocate(zplev_mean(ngrid,timelen))
-
-
-
 
 ! 0. Some initializations
@@ -60,5 +52,4 @@
        mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
      enddo
-
 
 ! b. pressure level
@@ -91,15 +82,8 @@
      deallocate(pvapor)
 
-
-
-
 !    1.3 Density at the surface yearly averaged
      rhovapor_avg(:,:) = SUM(rhovapor(:,:,:),3)/timelen
 
-
      deallocate(rhovapor)
-
-
-
 
 ! 2. Compute rho soil vapor
@@ -118,10 +102,6 @@
      enddo
 
-
-
     rho_soil_avg(:,:,:) = SUM( rho_soil(:,:,:,:),4)/timelen
     deallocate(rho_soil)
-
-
 
 ! 3. Computing ice table
@@ -132,19 +112,10 @@
 
          do isoil = 1,nsoil_PEM
-
-
-             diff_rho(:,:,isoil) = rhovapor_avg(:,:) - rho_soil_avg(:,:,isoil)
-!             write(*,*) 'diff =',ig,islope,isoil,diff_rho(ig,islope,isoil),rhovapor_avg(ig,islope) ,rho_soil_avg(ig,islope,isoil)
+           diff_rho(:,:,isoil) = rhovapor_avg(:,:) - rho_soil_avg(:,:,isoil)
+!          write(*,*) 'diff =',ig,islope,isoil,diff_rho(ig,islope,isoil),rhovapor_avg(ig,islope) ,rho_soil_avg(ig,islope,isoil)
          enddo
     
-
   deallocate(rhovapor_avg)
-      deallocate(rho_soil_avg)
-
-
-
-
-
-
+  deallocate(rho_soil_avg)
 
      do ig = 1,ngrid
@@ -164,12 +135,6 @@
         enddo
       enddo
-
     
     deallocate(diff_rho)
-
-
-
-
-
 
 !=======================================================================
@@ -177,6 +142,3 @@
 
 
-
-
-
       END
Index: unk/LMDZ.COMMON/libf/evolution/conf_evol.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/conf_evol.F90	(revision 2834)
+++ 	(revision )
@@ -1,54 +1,0 @@
-!
-! $Id: conf_evol.F $
-
-SUBROUTINE conf_evol( tapedef, etatinit )
-
-#ifdef CPP_IOIPSL
-  use IOIPSL
-#else
-  ! if not using IOIPSL, we still need to use (a local version of) getin
-  use ioipsl_getincom
-#endif
-
-  USE temps_mod_evol, ONLY: nyear, dt_pem, alpha_criterion
-
-  IMPLICIT NONE
-!-----------------------------------------------------------------------
-! Read the run_pem.def file
-!
-!     Arguments :
-!
-!     nyear   : Maximun number of year over which the PEM can interpolate
-!     dt_pem  : in years, the time step used by the PEM    
-!     alpha_criterion : percentage of change before stopping the PEM 
-!
-  LOGICAL,INTENT(IN) :: etatinit
-  INTEGER,INTENT(IN) :: tapedef
-
-!   Declarations :
-!   --------------
-  include "dimensions.h"
-  include "paramet.h"
-  include "comdissnew.h"
-  include "iniprint.h"
-
-!   local:
-!   ------
-
-  REAL clonn,clatt,grossismxx,grossismyy
-  REAL dzoomxx,dzoomyy, tauxx,tauyy
-  LOGICAL  fxyhypbb, ysinuss
-  LOGICAL use_filtre_fft
-
-
-  nyear=100
-  CALL getin('nyear', nyear)
-
-  dt_pem=1
-  CALL getin('dt_pem', dt_pem)
-
-  alpha_criterion=0.2
-  CALL getin('alpha_criterion', alpha_criterion)
-
-
-END SUBROUTINE conf_evol
Index: /trunk/LMDZ.COMMON/libf/evolution/criterion_co2_ice_stop_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/criterion_co2_ice_stop_mod.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/criterion_co2_ice_stop_mod.F90	(revision 2835)
@@ -23,6 +23,4 @@
   REAL,    intent(in) ::  latitude(ngrid)          ! physical point field : Latitude
   REAL,    intent(in) ::  initial_co2_ice(n_band_lat)  ! Initial/Actual surface of water ice
-
-
 
 !   OUTPUT
@@ -59,7 +57,4 @@
        present_co2(j).GT.initial_co2_ice(j)*(1+alpha_criterion)) then
          STOPPING=.TRUE. 
-         print *, "j", j
-         print *, "present_co2(j)", present_co2(j)
-         print *, "initial_co2_ice(j)", initial_co2_ice(j)
     endif
   enddo
Index: /trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_slope.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_slope.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_slope.F90	(revision 2835)
@@ -20,5 +20,5 @@
 !   INPUT
   INTEGER, intent(in) :: ngrid,nslope                  ! # of grid physical grid points 
-  REAL,    intent(in) :: cell_area(ngrid)       ! physical point field : Area of the cells
+  REAL,    intent(in) :: cell_area(ngrid)              ! physical point field : Area of the cells
   REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
   REAL,    intent(in) :: ini_surf
@@ -27,5 +27,4 @@
   REAL,    intent(in) :: global_ave_press_new
 
-
 !   OUTPUT
   LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
@@ -33,5 +32,5 @@
 !   local:
 !   -----
-  INTEGER :: i,islope                    ! Loop
+  INTEGER :: i,islope   ! Loop
   REAL :: present_surf  ! Initial/Actual surface of water ice
 
@@ -46,25 +45,18 @@
    do islope=1,nslope
       if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
-         print *, "i", i
-         print *, "initial_h2o_ice(i,islope)", initial_h2o_ice(i,islope)
-         print *, "qsurf(i,:)", qsurf(i,:)
-         print *, "cell_area(i)", cell_area(i)
-         print *, "present_surf",present_surf
          present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
       endif
    enddo
   enddo
-
-!  print *, "initial_h2o_ice", initial_h2o_ice
-!  print *, "qsurf", qsurf
-
-  print *, "present_surf", present_surf
-  print *, "ini_surf", ini_surf
-  print *, "ini_surf*0.8", ini_surf*(1-alpha_criterion)
   
 !   check of the criterion
   if(present_surf.LT.ini_surf*(1-alpha_criterion) .OR. &
      present_surf.GT.ini_surf*(1+alpha_criterion)) then
-  STOPPING=.TRUE. 
+    STOPPING=.TRUE. 
+    print *, "Reason of stopping : The surface of co2 ice sublimating reach the threshold:"
+    print *, "Current surface of co2 ice sublimating=", present_surf
+    print *, "Initial surface of co2 ice sublimating=", ini_surf
+    print *, "Percentage of change accepted=", alpha_criterion*100
+    print *, "present_surf<ini_surf*(1-alpha_criterion)", (present_surf.LT.ini_surf*(1-alpha_criterion))
   endif
 
@@ -73,12 +65,12 @@
   endif
 
-!  if(global_ave_press_GCM.LT.global_ave_press_new*(1-alpha_criterion) .OR. &
-!     global_ave_press_GCM.GT.global_ave_press_new*(1+alpha_criterion)) then
-!  STOPPING=.TRUE. 
-!  endif
-
   if(global_ave_press_new.LT.global_ave_press_GCM*(0.9) .OR. &
      global_ave_press_new.GT.global_ave_press_GCM*(1.1)) then
-  STOPPING=.TRUE. 
+    STOPPING=.TRUE. 
+    print *, "Reason of stopping : The global pressure reach the threshold:"
+    print *, "Current global pressure=", global_ave_press_new
+    print *, "GCM global pressure=", global_ave_press_GCM
+    print *, "Percentage of change accepted=", 0.1*100
+    print *, "global_ave_press_new<global_ave_press_GCM*(0.9)", (global_ave_press_new.LT.global_ave_press_GCM*(0.9))
   endif
 
@@ -86,5 +78,2 @@
 
 
-
-
-
Index: /trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_water_slope.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_water_slope.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/criterion_ice_stop_mod_water_slope.F90	(revision 2835)
@@ -25,5 +25,4 @@
   REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
 
-
 !   OUTPUT
   LOGICAL, intent(out) :: STOPPING              ! Logical : is the criterion reached?
@@ -39,34 +38,27 @@
     STOPPING=.FALSE.
 
-!   computation of the actual surface
+!   computation of the present surface of water ice sublimating
   present_surf=0.
   do i=1,ngrid
     do islope=1, nslope
       if (initial_h2o_ice(i,islope).GT.0.5 .and. qsurf(i,islope).GT.0.) then
-!         print *, "i", i
-!         print *, "initial_h2o_ice(i,islope)", initial_h2o_ice(i,islope)
-!         print *, "qsurf(i,islope)", qsurf(i,islope)
-!         print *, "cell_area(i)", cell_area(i)
-!         print *, "present_surf",present_surf
          present_surf=present_surf+cell_area(i)*subslope_dist(i,islope)
       endif
     enddo
   enddo
-
-!  print *, "initial_h2o_ice", initial_h2o_ice
-!  print *, "qsurf", qsurf
-
-!  print *, "present_surf", present_surf
-!  print *, "ini_surf", ini_surf
-!  print *, "ini_surf*0.8", ini_surf*(1-alpha_criterion)
   
 !   check of the criterion
   if(present_surf.LT.ini_surf*(1-alpha_criterion) .OR. &
      present_surf.GT.ini_surf*(1+alpha_criterion)) then
-  STOPPING=.TRUE. 
+    STOPPING=.TRUE. 
+    print *, "Reason of stopping : The surface of water ice sublimating reach the threshold:"
+    print *, "Current surface of water ice sublimating=", present_surf
+    print *, "Initial surface of water ice sublimating=", ini_surf
+    print *, "Percentage of change accepted=", alpha_criterion*100
+    print *, "present_surf<ini_surf*(1-alpha_criterion)", (present_surf.LT.ini_surf*(1-alpha_criterion))
   endif
 
   if (ini_surf.LT. 1E-5 .and. ini_surf.GT. -1E-5) then
-       STOPPING=.FALSE.
+    STOPPING=.FALSE.
   endif
 
Index: /trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod_slope.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod_slope.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/evol_co2_ice_s_mod_slope.F90	(revision 2835)
@@ -21,5 +21,4 @@
 
   INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope   ! # of grid points along longitude/latitude/ total
-!  REAL, intent(in) ::  tendencies_h2o_ice_phys(ngrid) ! physical point field : Evolution of perenial ice over one year
   REAL, intent(in) ::  cell_area(ngrid)
 
@@ -28,5 +27,4 @@
   LOGICAL :: STOPPING
   REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
-
 
 !   local:
Index: /trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90	(revision 2835)
@@ -6,5 +6,5 @@
 
   USE temps_mod_evol, ONLY: dt_pem
-  	  use comslope_mod, ONLY: subslope_dist
+  use comslope_mod, ONLY: subslope_dist
 
       IMPLICIT NONE
@@ -22,5 +22,4 @@
 
   INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope   ! # of grid points along longitude/latitude/ total
-!  REAL, intent(in) ::  tendencies_h2o_ice_phys(ngrid) ! physical point field : Evolution of perenial ice over one year
   REAL, intent(in) ::  cell_area(ngrid)
 
@@ -30,17 +29,12 @@
   REAL, intent(inout) ::  tendencies_h2o_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
 
-
 !   local:
 !   ----
 
   INTEGER :: i,j,ig0,islope                                  ! loop variable
-!  REAL :: not_budget, budget
   REAL :: pos_tend, neg_tend, real_coefficient,negative_part
   REAL ::  new_tendencies(ngrid,nslope)
 
-
 !=======================================================================
-
-!  budget=sum(qsurf(:))
 
   pos_tend=0.
@@ -59,12 +53,8 @@
   enddo
 
-!  print *, "pos_tend", pos_tend
-!  print *, "neg_tend", neg_tend
-
   if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then
      do i=1,ngrid
        do islope=1,nslope
        if(tendencies_h2o_ice_phys(i,islope).LT.0) then
-!          print *, "pos_tend/neg_tend", pos_tend/neg_tend
           new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)*(pos_tend/neg_tend)
        else
@@ -74,5 +64,4 @@
      enddo
   elseif(neg_tend.LT.pos_tend .and. neg_tend.GT.0) then
-!          print *, "neg_tend/pos_tend", neg_tend/pos_tend
      do i=1,ngrid
        do islope=1,nslope
@@ -85,6 +74,8 @@
      enddo
   elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then 
-
-!      call criterion_ice_stop(cell_area,1,qsurf*0.,STOPPING,ngrid,cell_area)
+    print *, "Reason of stopping : There is either no water ice sublimating or no water ice increasing !!"
+    print *, "Tendencies on ice sublimating=", neg_tend
+    print *, "Tendencies on ice increasing=", pos_tend
+    print *, "This can be due to the absence of water ice in the PCM run!!"
       call criterion_ice_stop_water_slope(cell_area,1,qsurf(:,:)*0.,STOPPING,ngrid,cell_area)
       do i=1,ngrid
@@ -95,38 +86,19 @@
   endif
 
- 
-  
   negative_part = 0.
-
 
 ! Evolution of the water ice for each physical point
   do i=1,ngrid
     do islope=1, nslope
-!      qsurf(i)=qsurf(i)+tendencies_h2o_ice_phys(i)*dt_pem
       qsurf(i,islope)=qsurf(i,islope)+new_tendencies(i,islope)*dt_pem
-!      budget=budget+tendencies_h2o_ice_phys(i)*dt_pem
       if (qsurf(i,islope).lt.0) then
-!        not_budget=not_budget+qsurf(i)
-!       print *, "NNqsurf(i,islope)", qsurf(i,islope)
-!       print *, "NNnew_tendencies(i,islope)", new_tendencies(i,islope)
-!       print *, "NNtendencies_h2o_ice_phys(i,islope)", tendencies_h2o_ice_phys(i,islope)
         negative_part=negative_part-qsurf(i,islope)*cell_area(i)*subslope_dist(i,islope)
         qsurf(i,islope)=0.
         tendencies_h2o_ice_phys(i,islope)=0.
-!        print *, "NNineg", i
-      endif
-      if(qsurf(i,islope).NE.qsurf(i,islope)) then
-!          print *, "qsurf(i,islope)",qsurf(i,islope)
-!          print *, "new_tendencies",new_tendencies(i,islope)
-!          print *, "tendencies_h2o_ice_phys",tendencies_h2o_ice_phys(i,islope)
-!          print *, "i", i
-!          print *,"islope",islope
       endif
     enddo
   enddo
 
-!  print *, "negative_part", negative_part
   real_coefficient=negative_part/pos_tend
-!  print *, "real_coefficient", real_coefficient
 
   do i=1,ngrid
@@ -140,7 +112,3 @@
 
 
-! Conservation of water ice
-!  qsurf(:)=qsurf(:)*budget/sum(qsurf(:))
-
-
 END SUBROUTINE evol_h2o_ice_s_slope
Index: /trunk/LMDZ.COMMON/libf/evolution/ini_soil_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/ini_soil_mod.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/ini_soil_mod.F90	(revision 2835)
@@ -1,4 +1,3 @@
       MODULE ini_soil_mod
-
 
       IMPLICIT NONE
@@ -6,9 +5,6 @@
       CONTAINS    
 
-
      subroutine ini_icetable(timelen,ngrid,nsoil_PEM, &
                therm_i, timestep,tsurf_ave,tsoil_ave,tsurf_inst, tsoil_inst,q_co2,q_h2o,ps,ice_table)
-
-
 
       use vertical_layers_mod, only: ap,bp
@@ -28,8 +24,4 @@
 
 #include "dimensions.h"
-!#include "dimphys.h"
-
-!#include"comsoil.h"
-
 
 !-----------------------------------------------------------------------
@@ -41,5 +33,4 @@
       integer,intent(in) :: nsoil_PEM	! number of soil layers  in the PEM
 
-
       real,intent(in) :: timestep	    ! time step
       real,intent(in) :: tsurf_ave(ngrid)   ! surface temperature
@@ -49,5 +40,4 @@
       real,intent(in) :: ps(ngrid,timelen)                        ! surface pressure [Pa]
       real,intent(in) :: tsurf_inst(ngrid,timelen) ! soil (mid-layer) temperature
-
 
 ! outputs:
@@ -57,6 +47,4 @@
       real,intent(inout) :: therm_i(ngrid,nsoil_PEM) ! thermal inertia
 
-
-      
 ! local variables:
       integer ig,isoil,it,k,iref
@@ -99,46 +87,39 @@
     real,allocatable :: diff_rho(:)                    ! difference of vapor content
 
-
-      A =(1/m_co2 - 1/m_noco2)
-      B=1/m_noco2
-
-
- ice_table(:) = 1.
-do ig = 1,ngrid
- 
- error_depth = 1.
- countloop = 0
-
-
-    
- do while(( error_depth.gt.tol_error).and.(countloop.lt.countmax).and.(ice_table(ig).gt.-1e-20))
-    
-    countloop = countloop +1
-    Tcol_saved(:) = tsoil_ave(ig,:)
+    A =(1/m_co2 - 1/m_noco2)
+    B=1/m_noco2
+
+    ice_table(:) = 1.
+    do ig = 1,ngrid
+      error_depth = 1.
+      countloop = 0
+    
+      do while(( error_depth.gt.tol_error).and.(countloop.lt.countmax).and.(ice_table(ig).gt.-1e-20))
+    
+        countloop = countloop +1
+        Tcol_saved(:) = tsoil_ave(ig,:)
 
 !2.  Compute ice table
 
 ! 2.1 Compute  water density at the surface, yearly averaged
-    allocate(mass_mean(timelen))
+        allocate(mass_mean(timelen))
 !   1.1 Compute the partial pressure of vapor
 ! a. the molecular mass into the column
        mass_mean(:) = 1/(A*q_co2(ig,:) +B)
 ! b. pressure level
-     allocate(zplev(timelen))
-     do it = 1,timelen
+        allocate(zplev(timelen))
+        do it = 1,timelen
          zplev(it) = ap(1) + bp(1)*ps(ig,it)
-     enddo
+        enddo
 
 ! c. Vapor pressure
-     allocate(pvapor(timelen))
-     pvapor(:) = mass_mean(:)/m_h2o*q_h2o(ig,:)*zplev(:)
-     deallocate(zplev) 
-     deallocate(mass_mean)
-
-
+        allocate(pvapor(timelen))
+        pvapor(:) = mass_mean(:)/m_h2o*q_h2o(ig,:)*zplev(:)
+        deallocate(zplev) 
+        deallocate(mass_mean)
   
 !  d!  Check if there is frost at the surface and then compute the density
 !    at the surface
-     allocate(rhovapor(timelen))
+        allocate(rhovapor(timelen))
 
          do it = 1,timelen
@@ -146,132 +127,109 @@
            rhovapor(it) = min(psv_surf,pvapor(it))/tsurf_inst(ig,it)
          enddo
-     deallocate(pvapor)
-     rhovapor_avg =     SUM(rhovapor(:),1)/timelen                     
-     deallocate(rhovapor)
+        deallocate(pvapor)
+        rhovapor_avg =SUM(rhovapor(:),1)/timelen                     
+        deallocate(rhovapor)
 
 ! 2.2 Compute  water density at the soil layer, yearly averaged
 
-     allocate(rho_soil(timelen))
-     allocate(rho_soil_avg(nsoil_PEM))
-
-
-
-      do isoil = 1,nsoil_PEM
-        do it = 1,timelen 
-              rho_soil(it) = exp(alpha/tsoil_inst(ig,isoil,it) +beta)/tsoil_inst(ig,isoil,it)        
-       enddo
-       rho_soil_avg(isoil) = SUM(rho_soil(:),1)/timelen     
-      enddo
-     deallocate(rho_soil)
+        allocate(rho_soil(timelen))
+        allocate(rho_soil_avg(nsoil_PEM))
+
+        do isoil = 1,nsoil_PEM
+          do it = 1,timelen 
+            rho_soil(it) = exp(alpha/tsoil_inst(ig,isoil,it) +beta)/tsoil_inst(ig,isoil,it)        
+          enddo
+         rho_soil_avg(isoil) = SUM(rho_soil(:),1)/timelen     
+        enddo
+        deallocate(rho_soil)
          
-
 !2.3 Final: compute ice table
-         icedepth_prev = ice_table(ig)
-         ice_table(ig) = -1
-         allocate(diff_rho(nsoil_PEM))
-         do isoil = 1,nsoil_PEM
-             diff_rho(isoil) = rhovapor_avg -  rho_soil_avg(isoil)
+        icedepth_prev = ice_table(ig)
+        ice_table(ig) = -1
+        allocate(diff_rho(nsoil_PEM))
+        do isoil = 1,nsoil_PEM
+          diff_rho(isoil) = rhovapor_avg -  rho_soil_avg(isoil)
+        enddo
+        deallocate(rho_soil_avg)
+
+        if(diff_rho(1) > 0) then
+          ice_table(ig) = 0.
+        else
+          do isoil = 1,nsoil_PEM -1
+            if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
+              z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
+              z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
+              ice_table(ig) = -z2/z1
+              exit
+            endif
+          enddo
+        endif
+    
+        deallocate(diff_rho)
+
+!3. Update Soil Thermal Inertia
+
+        if (ice_table(ig).gt. 0.) then
+          if (ice_table(ig).lt. 1e-10) then
+            do isoil = 1,nsoil_PEM
+              therm_i(ig,isoil)=ice_inertia
+            enddo
+          else 
+      ! 4.1 find the index of the mixed layer
+          iref=0 ! initialize iref
+          do k=1,nsoil_PEM ! loop on layers
+            if (ice_table(ig).ge.layer_PEM(k)) then
+              iref=k ! pure regolith layer up to here
+            else
+              ! correct iref was obtained in previous cycle
+            exit
+            endif
+          enddo
+
+      ! 4.2 Build the new ti
+         do isoil=1,iref
+           therm_i(ig,isoil) =inertiedat_PEM(ig,isoil) 
          enddo
-         deallocate(rho_soil_avg)
-
-
-         if(diff_rho(1) > 0) then
-           ice_table(ig) = 0.
-         else
-           do isoil = 1,nsoil_PEM -1
-             if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
-               z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
-               z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
-               ice_table(ig) = -z2/z1
-               exit
-             endif
-           enddo
-          endif
-
-
-    
-    deallocate(diff_rho)
-
-
-!3. Update Soil Thermal Inertia
-
-  if (ice_table(ig).gt. 0.) then
-  if (ice_table(ig).lt. 1e-10) then
-      do isoil = 1,nsoil_PEM
-      therm_i(ig,isoil)=ice_inertia
-      enddo
-  else 
-      ! 4.1 find the index of the mixed layer
-      iref=0 ! initialize iref
-      do k=1,nsoil_PEM ! loop on layers
-        if (ice_table(ig).ge.layer_PEM(k)) then
-          iref=k ! pure regolith layer up to here
-        else
-         ! correct iref was obtained in previous cycle
-         exit
-        endif
-       
-      enddo
-
-
-
-      ! 4.2 Build the new ti
-      do isoil=1,iref
-         therm_i(ig,isoil) =inertiedat_PEM(ig,isoil) 
-      enddo
-
-
-      if (iref.lt.nsoil_PEM) then
-        if (iref.ne.0) then
-          ! mixed layer
-           therm_i(ig,iref+1)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
-            (((ice_table(ig)-layer_PEM(iref))/(inertiedat_PEM(ig,iref)**2))+ &
+
+         if (iref.lt.nsoil_PEM) then
+           if (iref.ne.0) then
+             ! mixed layer
+             therm_i(ig,iref+1)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
+              (((ice_table(ig)-layer_PEM(iref))/(inertiedat_PEM(ig,iref)**2))+ &
                       ((layer_PEM(iref+1)-ice_table(ig))/(ice_inertia**2))))
             
-
-
-        else ! first layer is already a mixed layer
-          ! (ie: take layer(iref=0)=0)
-          therm_i(ig,1)=sqrt((layer_PEM(1))/ &
+           else ! first layer is already a mixed layer
+            ! (ie: take layer(iref=0)=0)
+            therm_i(ig,1)=sqrt((layer_PEM(1))/ &
                           (((ice_table(ig))/(inertiedat_PEM(ig,1)**2))+ &
                            ((layer_PEM(1)-ice_table(ig))/(ice_inertia**2))))
-        endif ! of if (iref.ne.0)        
-        ! lower layers of pure ice
-        do isoil=iref+2,nsoil_PEM
-          therm_i(ig,isoil)=ice_inertia
-        enddo
-      endif ! of if (iref.lt.(nsoilmx))
-      endif ! permanent glaciers
-
-
+           endif ! of if (iref.ne.0)        
+           ! lower layers of pure ice
+           do isoil=iref+2,nsoil_PEM
+             therm_i(ig,isoil)=ice_inertia
+           enddo
+         endif ! of if (iref.lt.(nsoilmx))
+       endif ! permanent glaciers
 
        call soil_pem_1D(nsoil_PEM,.true.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM) 
-
        call soil_pem_1D(nsoil_PEM,.false.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM) 
 
-
-      do it = 1,timelen 
-        tsoil_inst(ig,:,it) = tsoil_inst(ig,:,it) - (Tcol_saved(:) - tsoil_ave(ig,:))
-      enddo
-
-
-
-      error_depth = abs(icedepth_prev - ice_table(ig))
-      
-     endif
-
-
-enddo
-
- error_depth = 1.
- countloop = 0
-enddo
+       do it = 1,timelen 
+         tsoil_inst(ig,:,it) = tsoil_inst(ig,:,it) - (Tcol_saved(:) - tsoil_ave(ig,:))
+       enddo
+
+       error_depth = abs(icedepth_prev - ice_table(ig))      
+       endif
+
+     enddo
+
+     error_depth = 1.
+     countloop = 0
+   enddo
  
-
       END SUBROUTINE ini_icetable
       subroutine soil_pem_1D(nsoil,firstcall, &
                therm_i,                          &
                timestep,tsurf,tsoil,alph,beta)
-
 
       use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,  &
@@ -309,8 +267,4 @@
       real,intent(inout) :: alph(nsoil-1)
       real,intent(inout) :: beta(nsoil-1)
-
-
-
-
       
 ! local variables:
@@ -330,5 +284,4 @@
 	enddo
 
-
 ! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
         do ik=1,nsoil-1
@@ -336,5 +289,4 @@
                      +(mlayer_PEM(ik)-layer_PEM(ik))*mthermdiff_PEM(ik-1))  &
                          /(mlayer_PEM(ik)-mlayer_PEM(ik-1))
-
 	enddo
 
@@ -365,11 +317,5 @@
 	enddo
 
-
-
-      
-  
 endif ! of if (firstcall)
-
-
 
       IF (.not.firstcall) THEN
@@ -387,7 +333,4 @@
       
           ENDIF
-
-
-
 
 !  2. Compute beta_PEM coefficients (preprocessing for next time step)
@@ -405,13 +348,5 @@
       enddo
 
-
-
       end
 
-
-
-
-      END MODULE ini_soil_mod
-
-
-      
+      END MODULE ini_soil_mod      
Index: /trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90	(revision 2835)
+++ /trunk/LMDZ.COMMON/libf/evolution/lask_param_mod.F90	(revision 2835)
@@ -0,0 +1,37 @@
+module lask_param_mod
+
+implicit none
+
+  real,save,allocatable :: yearlask(:)
+  real,save,allocatable :: oblask(:)
+  real,save,allocatable :: exlask(:)
+  real,save,allocatable :: lsplask(:)
+  integer, save :: last_ilask
+
+contains
+
+  subroutine ini_lask_param_mod(nlask)
+  
+  implicit none
+  integer,intent(in) :: nlask ! number of line
+  
+  allocate(yearlask(nlask))
+  allocate(oblask(nlask))
+  allocate(exlask(nlask))
+  allocate(lsplask(nlask))
+  
+  end subroutine ini_lask_param_mod
+
+
+  subroutine end_lask_param_mod
+
+  implicit none
+
+  if (allocated(yearlask)) deallocate(yearlask)
+  if (allocated(oblask)) deallocate(oblask)
+  if (allocated(exlask)) deallocate(exlask)
+  if (allocated(lsplask)) deallocate(lsplask)
+
+  end subroutine end_lask_param_mod
+  
+end module lask_param_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM.F90	(revision 2835)
@@ -56,4 +56,6 @@
   CALL err(NF90_CLOSE(fID),"close",fichnom)
 
+  print *, "The number of timestep of the PCM run data=", timelen
+
   CONTAINS
 
Index: unk/LMDZ.COMMON/libf/evolution/orbit_param.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/orbit_param.F90	(revision 2834)
+++ 	(revision )
@@ -1,68 +1,0 @@
-!
-! $Id $
-!
-SUBROUTINE orbit_param(year)
-
-      IMPLICIT NONE
-
-
-!=======================================================================
-!
-!  Compute the tendencies of the evolution of water ice over the years
-!
-!=======================================================================
-
-!   arguments:
-!   ----------
-
-!   INPUT
-
-     INTEGER, intent(in) :: iim_input,jjm_input,ngrid                             ! # of grid points along longitude/latitude/ total
-     REAL, intent(in) , dimension(iim_input+1,jjm_input+1):: min_h2o_ice_Y1       ! LON x LAT field : minimum of water ice at each point for the first year
-     REAL, intent(in) , dimension(iim_input+1,jjm_input+1):: min_h2o_ice_Y2       ! LON x LAT field : minimum of water ice at each point for the second year
-
-!   OUTPUT
-     REAL, intent(out) , dimension(iim_input+1,jjm_input+1) :: tendencies_h2o_ice ! LON x LAT field : difference between the minima = evolution of perenial ice
-     REAL, intent(out) , dimension(ngrid)   :: tendencies_h2o_ice_phys            ! physical point field : difference between the minima = evolution of perenial ice
-
-!   local:
-!   ------
-
-     INTEGER :: i,j,ig0                                                           ! loop variable
-
-!=======================================================================
-
-
-!  We compute the difference
-  tendencies_h2o_ice(:,:)=min_h2o_ice_Y2(:,:)-min_h2o_ice_Y1(:,:)
-
-!  If the difference is too small; there is no evolution
-  DO j=1,jjm_input+1
-    DO i = 1, iim_input
-       if(abs(tendencies_h2o_ice(i,j)).LT.1.0E-10) then
-          tendencies_h2o_ice(i,j)=0.
-       endif
-    ENDDO
-  ENDDO
-
-
-!  We reorganise the difference on the physical grid
-  tendencies_h2o_ice_phys(1)=tendencies_h2o_ice(1,1)
-
-  ig0 = 2
-  DO j=2,jjm_input
-    DO i = 1, iim_input
-       tendencies_h2o_ice_phys(ig0)  =tendencies_h2o_ice(i,j)
-       ig0= ig0 + 1
-    ENDDO
-  ENDDO
-
-  tendencies_h2o_ice_phys(ig0) = tendencies_h2o_ice(1,jjm_input+1)
-
-
-END SUBROUTINE orbit_param
-
-
-
-
-
Index: /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90	(revision 2835)
+++ /trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90	(revision 2835)
@@ -0,0 +1,167 @@
+      MODULE orbit_param_criterion_mod
+
+      IMPLICIT NONE
+
+      CONTAINS
+
+      SUBROUTINE orbit_param_criterion(year_iter_max)
+
+      USE temps_mod_evol, ONLY: year_bp_ini, year_PEM
+      USE planete_h, ONLY: e_elips, obliquit, timeperi
+      USE comcstfi_h, only: pi
+      USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
+   				ini_lask_param_mod,last_ilask
+
+      IMPLICIT NONE
+
+!--------------------------------------------------------
+! Input Variables
+!--------------------------------------------------------
+
+!--------------------------------------------------------
+! Output Variables
+!--------------------------------------------------------
+
+      integer,intent(out) :: year_iter_max	! Maximum number of iteration of the PEM before orb changes too much
+
+!--------------------------------------------------------
+! Local variables 
+!--------------------------------------------------------
+
+      real :: Year	                ! Year of the simulation
+      real :: timeperi_ls	                ! Year of the simulation
+      integer nlask,ilask!,last_ilask !Loop variable
+      parameter (nlask = 20001)
+
+      real max_change_obl,max_change_ex,max_change_lsp ! Percentage of change that is considered to be acceptible
+
+      real max_obl,max_ex,max_lsp !Maximum value of orbit param given the acceptable percentage 
+      real min_obl,min_ex,min_lsp !Maximum value of orbit param given the acceptable percentage
+
+      real max_obl_iter,max_ex_iter,max_lsp_iter !Maximum year iteration before reaching an unacceptable value
+
+      logical obl_not_found, ex_not_found,lsp_not_found !Loop variable (first call)
+
+      ! **********************************************************************
+      ! 0. Initializations
+      ! **********************************************************************
+
+          Year=year_bp_ini+year_PEM
+          timeperi_ls=timeperi*360/2/pi
+
+	  call ini_lask_param_mod(nlask)
+
+          print *, "orbit_param_criterion, Year in pem.def=", year_bp_ini
+          print *, "orbit_param_criterion, Year in the startpem.nc =", year_PEM
+          print *, "orbit_param_criterion, Current year=", Year
+          print *, "orbit_param_criterion, Current obl=", obliquit
+          print *, "orbit_param_criterion, Current ex=", e_elips
+          print *, "orbit_param_criterion, Current lsp=", timeperi_ls
+
+          open(73,file='ob_ex_lsp.asc')
+          do ilask=1,nlask
+            read(73,*) yearlask(ilask),oblask(ilask),      &
+             exlask(ilask),lsplask(ilask)
+            yearlask(ilask)=yearlask(ilask)*1000
+
+            if(yearlask(ilask).GT.Year) then
+		last_ilask=ilask+1
+            endif
+          end do
+          close(73)
+
+       print *, "Coresponding line in the ob_ex_lsp.asc file=", last_ilask
+
+! 5% max change case
+
+        max_change_obl=0.05
+        max_change_ex=0.05
+        max_change_lsp=0.05
+
+	max_obl=obliquit*(1+max_change_obl)
+	min_obl=obliquit*(1-max_change_obl)
+
+	max_ex=e_elips*(1+max_change_ex)
+	min_ex=e_elips*(1-max_change_ex)
+
+	max_lsp=timeperi_ls*(1+max_change_lsp)
+	min_lsp=timeperi_ls*(1-max_change_lsp)
+
+!End of 5% max change case
+
+!Constant max change case 
+
+        max_change_obl=0.01
+        max_change_ex=0.01
+        max_change_lsp=20.
+
+	max_obl=obliquit+max_change_obl
+	min_obl=obliquit-max_change_obl
+
+	max_ex=e_elips+max_change_ex
+	min_ex=e_elips-max_change_ex
+
+	max_lsp=timeperi_ls+max_change_lsp
+	min_lsp=timeperi_ls-max_change_lsp
+
+!End Constant max change case 
+
+        obl_not_found=.TRUE.
+        ex_not_found=.TRUE.
+        lsp_not_found=.TRUE.
+
+        max_obl_iter=999999999999
+        max_ex_iter =999999999999
+        max_lsp_iter=999999999999
+
+        do ilask=last_ilask+1,1,-1
+           if((oblask(ilask).GT.max_obl).and. obl_not_found ) then
+              max_obl_iter=(max_obl-oblask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
+                               / (oblask(ilask-1)-oblask(ilask))
+              obl_not_found=.FALSE.
+           elseif((oblask(ilask).LT.min_obl).and. obl_not_found ) then
+              max_obl_iter=(min_obl-oblask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
+                               / (oblask(ilask-1)-oblask(ilask))
+              obl_not_found=.FALSE.
+           endif
+           if((exlask(ilask).GT.max_ex).and. ex_not_found ) then
+              max_ex_iter=(max_ex-exlask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
+                               / (exlask(ilask-1)-exlask(ilask))
+              ex_not_found=.FALSE.
+           elseif((exlask(ilask).LT.min_ex ).and. ex_not_found ) then
+              max_ex_iter=(min_ex-exlask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
+                               / (exlask(ilask-1)-exlask(ilask))
+              ex_not_found=.FALSE.
+           endif
+           if((lsplask(ilask).GT.max_lsp).and. lsp_not_found ) then
+              max_lsp_iter=(max_lsp-lsplask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
+                               / (lsplask(ilask-1)-lsplask(ilask))
+              lsp_not_found=.FALSE.
+           elseif((lsplask(ilask).LT.min_lsp ).and. lsp_not_found ) then
+              max_lsp_iter=(min_lsp-lsplask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) &
+                               / (lsplask(ilask-1)-lsplask(ilask))
+              lsp_not_found=.FALSE.
+           endif
+        enddo
+
+      print *, "Maximum obliquity accepted=", max_obl
+      print *, "Minimum obliquity accepted=", min_obl
+      print *, "Maximum number of iteration for the obl. parameter=", max_obl_iter
+
+      print *, "Maximum excentricity accepted=", max_ex
+      print *, "Minimum excentricity accepted=", min_ex
+      print *, "Maximum number of iteration for the ex. parameter=", max_ex_iter
+
+      print *, "Maximum lsp accepted=", max_obl
+      print *, "Minimum lsp accepted=", min_obl
+      print *, "Maximum number of iteration for the lsp. parameter=", max_lsp_iter
+
+      year_iter_max=min(max_obl_iter,max_ex_iter,max_lsp_iter)
+
+      print *, "So the max. number of iteration (year) for the orbital parameter=", year_iter_max
+
+        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 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 2835)
@@ -3,22 +3,28 @@
 
 ! I   Initialisation
-!    I_a READ run.def , run_pem.def
-!    I_b READ starfi_0.nc
-!    I_c READ GCM data and convert to the physical grid 
-!    I_d Compute tendencies & Save initial situation
+!    I_a READ run.def 
+!    I_b READ of start_evol.nc and starfi_evol.nc
+!    I_c Subslope parametrisation
+!    I_d READ GCM data and convert to the physical grid 
+!    I_e Initialisation of the PEM variable and soil 
+!    I_f Compute tendencies & Save initial situation
+!    I_g Save initial PCM situation
+!    I_h Read the PEMstart
+!    I_i Compute orbit criterion
 
 ! II  Run
 !    II_a update pressure, ice and tracers 
-!    II_b CO2 glaciers flows 
-!    II_c Update surface and soil temperatures
-!    II_d Update the tendencies and stopping criterion
+!    II_b Evolution of the ice
+!    II_c CO2 glaciers flows 
+!    II_d Update surface and soil temperatures
+!    II_e Update the tendencies 
+!    II_f Checking the stopping criterion
 
 ! III Output
-!    III_a Recomp tendencies for the start and startfi
+!    III_a Update surface value for the PCM start files
 !    III_b Write start and starfi.nc
 !    III_c Write start_pem
 
 !------------------------
-
 
 PROGRAM pem
@@ -32,15 +38,12 @@
                            hmons, summit, base,albedo_h2o_frost, &
                           frost_albedo_threshold,emissiv
-      use dimradmars_mod, only: totcloudfrac, albedo, & 
-                                ini_dimradmars_mod
-      use turb_mod, only: q2, wstar, ini_turb_mod
-      use dust_param_mod, only: tauscaling, ini_dust_param_mod
-!      use co2cloud_mod, only: mem_Mccn_co2, mem_Mh2o_co2,&
-!      &                        mem_Nccn_co2,ini_co2cloud
+      use dimradmars_mod, only: totcloudfrac, albedo
+      use turb_mod, only: q2, wstar
+      use dust_param_mod, only: tauscaling
       use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
                         nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
                         nf90_inquire_dimension,nf90_close
       use phyredem, only: physdem0, physdem1
-      use tracer_mod, only: noms,nqmx,igcm_h2o_ice ! tracer names
+      use tracer_mod, only: noms,igcm_h2o_ice ! tracer names
 
 ! For phyredem :
@@ -56,6 +59,4 @@
       USE iniphysiq_mod, ONLY: iniphysiq
       USE infotrac
-      USE temps_mod_evol, ONLY: nyear, dt_pem
-!     USE vertical_layers_mod, ONLY: ap,bp
       USE comcstfi_h, only: pi 
 
@@ -68,17 +69,24 @@
                            iflat
 
-      USE geometry_mod, only: longitude_deg,latitude_deg 
+      USE geometry_mod, only: latitude_deg 
 
       use pemredem, only:  pemdem1
-      USE soil_evolution_mod, ONLY: soil_pem_CN
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
-
+      USE soil_evolution_mod, ONLY: soil_pem, soil_pem_CN
       use comsoil_h_PEM, only: ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
-                              TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM, &        ! soil thermal inertia          
-                              tsoil_PEM ,       &        !number of subsurface layers
-                              mlayer_PEM,layer_PEM, &       ! soil mid layer depths
-                              fluxgeo, &  ! geothermal flux
-                              co2_adsorbded_phys ! mass of co2 in the regolith
+                              TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM,         & ! soil thermal inertia          
+                              tsoil_PEM, mlayer_PEM,layer_PEM,                  & !number of subsurface layers, soil mid layer depths
+                              fluxgeo, co2_adsorbded_phys                         ! geothermal flux, mass of co2 in the regolith
       use adsorption_mod, only : regolith_co2adsorption
+
+!!! For orbit parameters
+      USE temps_mod_evol, ONLY: dt_pem, evol_orbit_pem, Max_iter_pem
+      use planete_h, only: aphelie, periheli, year_day, peri_day, &
+                          obliquit
+      use orbit_param_criterion_mod, only : orbit_param_criterion
+      use recomp_orb_param_mod, only: recomp_orb_param
+
+
   IMPLICIT NONE
 
@@ -104,5 +112,5 @@
 
 ! Variable for reading start.nc
-      character (len = *), parameter :: FILE_NAME_start = "start_0.nc" !Name of the file used for initialsing the PEM
+      character (len = *), parameter :: FILE_NAME_start = "start_evol.nc" !Name of the file used for initialsing the PEM
   !   variables dynamiques
   REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
@@ -120,5 +128,5 @@
 ! Variable for reading starfi.nc
 
-      character (len = *), parameter :: FILE_NAME = "startfi_0.nc" !Name of the file used for initialsing the PEM
+      character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM
       integer :: ncid, varid,status                                !Variable for handling opening of files
       integer :: phydimid, subdimid, nlayerdimid, nqdimid          !Variable ID for Netcdf files
@@ -152,12 +160,7 @@
       REAL , dimension(:,:), allocatable :: min_h2o_ice_s_1        ! LON x LAT field : minimum of water ice at each point for the first year
       REAL , dimension(:,:), allocatable :: min_h2o_ice_s_2        ! LON x LAT field : minimum of water ice at each point for the second year
-      REAL , dimension(:,:,:), allocatable :: h2o_ice_first_last_day     ! LON x LAT x 2 field : First and Last value of a GCM year simulation of water ice
 
       REAL , dimension(:,:), allocatable :: min_co2_ice_s_1        ! LON x LAT field : minimum of water ice at each point for the first year
       REAL , dimension(:,:), allocatable :: min_co2_ice_s_2        ! LON x LAT field : minimum of water ice at each point for the second year
-      REAL , dimension(:,:,:), allocatable :: co2_ice_first_last_day     ! LON x LAT x 2 field : First and Last value of a GCM year simulation of water ice
-
-      REAL, dimension(:),allocatable  :: local_old_press           ! physical point field : Local pressure of initial/previous time step
-      REAL, dimension(:),allocatable  :: local_new_press           ! physical point field : Local pressure of current time step
 
       REAL :: global_ave_press_GCM
@@ -166,19 +169,12 @@
 
       REAL , dimension(:,:), allocatable ::  zplev_new
-      REAL , dimension(:,:), allocatable :: zplev_old
+      REAL , dimension(:,:), allocatable ::  zplev_gcm
       REAL , dimension(:,:,:), allocatable ::  zplev_new_timeseries  ! same but with the time series
-      REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries! same but with the time series
-
-
-
-
-      REAL :: tot_co2_atm,tot_var_co2_atm
-
-
-      INTEGER :: year_iter  !Counter for the number of PEM iteration
+      REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries   ! same but with the time series
+
       LOGICAL :: STOPPING_water   ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
       LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate?
-      LOGICAL :: STOPPING_co2   ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
-      LOGICAL :: STOPPING_1_co2 ! Logical : is there still water ice to sublimate?
+      LOGICAL :: STOPPING_co2     ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
+      LOGICAL :: STOPPING_1_co2   ! Logical : is there still water ice to sublimate?
 
       REAL, dimension(:,:,:),allocatable  :: q_co2_GCM ! Initial amount of co2 in the first layer
@@ -195,19 +191,13 @@
       REAL, ALLOCATABLE ::  q_h2o_GCM(:,:,:)
       REAL ,allocatable ::  q_h2o_PEM_phys(:,:)
-      real ,allocatable :: q_phys(:,:,:) 
       integer :: timelen
       REAL :: ave
 
       REAL, ALLOCATABLE :: p(:,:)  !(ngrid,llmp1)
-      REAL :: extra_mass ! Extra mass of a tracer if it is greater than 1
-
-      REAL :: deficit_mass ! Extra mass of a tracer if it is lower than 0
-
-
-
+      REAL :: extra_mass           ! Extra mass of a tracer if it is greater than 1
 
 !!!!!!!!!!!!!!!!!!!!!!!! SLOPE
       character*2 :: str2
-      REAL ,allocatable :: watercap_slope(:,:)    !(ngrid,nslope)
+      REAL ,allocatable :: watercap_slope(:,:)                           !(ngrid,nslope)
       REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1        ! LON x LAT field : minimum of water ice at each point for the first year
       REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2        ! LON x LAT field : minimum of water ice at each point for the second year
@@ -217,8 +207,8 @@
 
       REAL , dimension(:,:,:,:), allocatable :: co2_ice_GCM_slope        ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM
-      REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope        ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM
-
-
-      REAL, dimension(:,:),allocatable  :: initial_co2_ice_sublim_slope           ! physical point field : Logical array indicating sublimating point of co2 ice
+      REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope     ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM
+
+
+      REAL, dimension(:,:),allocatable  :: initial_co2_ice_sublim_slope    ! physical point field : Logical array indicating sublimating point of co2 ice
       REAL, dimension(:,:),allocatable  :: initial_h2o_ice_slope           ! physical point field : Logical array indicating sublimating point of h2o ice
       REAL, dimension(:,:),allocatable  :: initial_co2_ice_slope           ! physical point field : Logical array indicating sublimating point of co2 ice
@@ -230,5 +220,5 @@
       REAL, dimension(:,:),allocatable  :: tendencies_h2o_ice_phys_slope   ! physical point field : Tendency of evolution of perenial co2 ice over a year
 
-      REAL,SAVE,ALLOCATABLE,DIMENSION(:) ::  co2_hmax                          ! Maximum height  for CO2 deposit on slopes (m)
+      REAL,SAVE,ALLOCATABLE,DIMENSION(:) ::  co2_hmax                      ! Maximum height  for CO2 deposit on slopes (m)
 
       REAL, PARAMETER :: rho_co2 = 1600           ! CO2 ice density (kg/m^3)
@@ -237,59 +227,43 @@
       REAL , dimension(:), allocatable :: flag_co2flow_mesh(:)  !(ngrid)          ! To flag where there is a glacier flow
 
-
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL 
 
-
-     REAL, ALLOCATABLE :: tsurf_ave(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature [K]
-     REAL, ALLOCATABLE  :: tsurf_ave_phys(:,:) ! IG x LAT x SLOPE field : Averaged Surface Temperature [K]
-     REAL, ALLOCATABLE :: tsoil_ave(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
-     REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
+     REAL, ALLOCATABLE :: tsurf_ave(:,:,:)          ! LON x LAT x SLOPE field : Averaged Surface Temperature [K]
+     REAL, ALLOCATABLE  :: tsurf_ave_phys(:,:)      ! IG x LAT x SLOPE field : Averaged Surface Temperature [K]
+     REAL, ALLOCATABLE :: tsoil_ave(:,:,:,:)        ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
+     REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:)    ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
      REAL, ALLOCATABLE :: tsoil_ave_phys_yr1(:,:,:) !IG x SLOPE field : Averaged Soil Temperature [K]
 
-     REAL, ALLOCATABLE :: TI_GCM(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Thermal Inertia  [SI]
-     REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:) !LONX LAT x SLOPE XTULES field : NOn averaged Surf Temperature [K]
+     REAL, ALLOCATABLE :: TI_GCM(:,:,:,:)                  ! LON x LAT x SLOPE field : Averaged Thermal Inertia  [SI]
+     REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:)    !LONX LAT x SLOPE XTULES field : NOn averaged Surf Temperature [K]
      REAL, ALLOCATABLE :: tsurf_phys_GCM_timeseries(:,:,:) !IG x SLOPE XTULES field : NOn averaged Surf Temperature [K]
 
      REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
-     REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
-
-     REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature of the first year of the gcm [K]
+     REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:,:)    !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
+
+     REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:,:)     ! LON x LAT x SLOPE field : Averaged Surface Temperature of the first year of the gcm [K]
      REAL, ALLOCATABLE  :: tsurf_ave_phys_yr1(:,:) ! IG x LAT x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
 
-
-
-     REAL :: tol_soil = 1.e-6
-     INTEGER :: countmax = 15
-     INTEGER :: countloop
-     REAL, ALLOCATABLE :: tsoil_phys_PEM_saved(:,:,:)
-     REAL :: error_tsoil
-
-
-     LOGICAL :: firstcall
-!    REAL, PARAMETER :: daysec=88775.              ! duree du sol (s)  ~88775 s
-     REAL, PARAMETER :: year_day = 669
      REAL, PARAMETER :: year_step = 1
+     INTEGER :: year_iter        !Counter for the number of PEM iteration
+     INTEGER :: year_iter_max
      REAL :: timestep
 
-     REAL, ALLOCATABLE :: inertiesoil(:,:) ! Thermal inertia of the mesh for restart
-
-     REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:) ! Averaged GCM Thermal Inertia  [SI]
+     REAL, ALLOCATABLE :: inertiesoil(:,:)    ! Thermal inertia of the mesh for restart
+
+     REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:)  ! Averaged GCM Thermal Inertia  [SI]
      REAL, ALLOCATABLE :: TI_GCM_start(:,:,:) ! Averaged GCM Thermal Inertia  [SI]
-
-     REAL,ALLOCATABLE  :: q_h2o_PEM_phys_ave(:) ! averaged water vapor content
-     REAL,ALLOCATABLE  :: interp_coef(:)
 
      REAL,ALLOCATABLE  :: ice_depth(:,:)
      REAL,ALLOCATABLE  :: TI_locslope(:,:)
      REAL,ALLOCATABLE  :: Tsoil_locslope(:,:)
-     REAL,ALLOCATABLE  :: Tsoil_locslope_CN(:,:)
      REAL,ALLOCATABLE  :: Tsurf_locslope(:)
-     REAL,ALLOCATABLE  :: Tsurf_locslope_prev(:)
      REAL,ALLOCATABLE  :: alph_locslope(:,:)
      REAL,ALLOCATABLE  :: beta_locslope(:,:)   
-     REAL :: kcond
+
      REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)
      REAL, ALLOCATABLE :: delta_co2_adsorbded(:)
      REAL :: totmass_adsorbded
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
@@ -310,5 +284,8 @@
       time_phys=0. !test
 
-!      n_band_lat=18
+! Some constants
+
+      ngrid=ngridmx
+      nlayer=llm
 
       m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
@@ -317,5 +294,8 @@
       B=1/m_noco2
 
-
+      year_day=669
+      daysec=88775.
+      dtphys=0
+      timestep=year_day*daysec/year_step
 
 !------------------------
@@ -324,11 +304,11 @@
 !    I_a READ run.def 
 
-
 !------------------------
 
 !----------------------------READ run.def ---------------------
       CALL conf_gcm( 99, .TRUE. )
-      CALL conf_evol( 99, .TRUE. )
-
+
+      call infotrac_init
+      nq=nqtot
 
 !------------------------
@@ -336,19 +316,13 @@
 ! I   Initialisation
 !    I_a READ run.def 
-!    I_b READ starfi_0.nc
-
-!----------------------------READ startfi.nc --------------------- 
-
-!----------------------------initialisation --------------------- 
-
+!    I_b READ of start_evol.nc and starfi_evol.nc
+
+!------------------------
+
+!----------------------------Initialisation : READ some constant of startfi_evol.nc --------------------- 
+
+! In the gcm, these values are given to the physic by the dynamic.
+! Here we simply read them in the startfi_evol.nc file
       status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
-      status =nf90_inq_dimid(ncid, "physical_points", phydimid)
-      status =nf90_inquire_dimension(ncid, phydimid, len = ngrid)
-
-      status =nf90_inq_dimid(ncid, "nlayer", nlayerdimid)
-      status =nf90_inquire_dimension(ncid, nlayerdimid, len = nlayer)
-
-      status =nf90_inq_dimid(ncid,"number_of_advected_fields",nqdimid)
-      status =nf90_inquire_dimension(ncid, nqdimid, len = nq)
 
       allocate(longitude(ngrid))
@@ -372,9 +346,5 @@
       status =nf90_close(ncid)
 
-
-
 !----------------------------READ start.nc --------------------- 
-    
-     call infotrac_init
 
      allocate(q(ip1jmp1,llm,nqtot))
@@ -392,35 +362,24 @@
      status = nf90_get_var(ncid, bpvarid, bp)
      status =nf90_close(ncid)
-     
-        
-   
-    daysec=88775.
-    dtphys=0
-
- 
-
-    CALL iniphysiq(iim,jjm,llm, &
+
+     CALL iniphysiq(iim,jjm,llm, &
           (jjm-1)*iim+2,comm_lmdz, &
           daysec,day_ini,dtphys/nsplit_phys, &
           rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
           iflag_phys)
-   DO nnq=1,nqtot
-    if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
-   ENDDO 
-
- 
-!----------------------------reading ---------------------
-
+
+     DO nnq=1,nqtot
+       if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
+     ENDDO 
+
+!----------------------------READ startfi.nc --------------------- 
 
 ! First we read the initial state (starfi.nc)
 
-
-       allocate(watercap_slope(ngrid,nslope))
-       allocate(TI_GCM_start(ngrid,nsoilmx,nslope))
-       allocate(q_h2o_PEM_phys_ave(ngrid))
+      allocate(watercap_slope(ngrid,nslope))
+      allocate(TI_GCM_start(ngrid,nsoilmx,nslope))
       allocate(inertiesoil(ngrid,nsoilmx))
 
-
-         CALL phyetat0 (FILE_NAME,0,0, &
+      CALL phyetat0 (FILE_NAME,0,0, &
               nsoilmx,ngrid,nlayer,nq,   &
               day_ini,time_phys,         &
@@ -432,40 +391,48 @@
               qsurf_slope,watercap_slope)
 
-
-
-     deallocate(TI_GCM_start) !not used then
-
+     if(soil_pem) then
+       deallocate(TI_GCM_start) !not used then
+     endif
+
+! Remove unphysical values of surface tracer
      DO i=1,ngrid
-         DO nnq=1,nqtot
-           DO islope=1,nslope
-             if(qsurf_slope(i,nnq,islope).LT.0) then
-                qsurf_slope(i,nnq,islope)=0.
-             endif
-           enddo
+       DO nnq=1,nqtot
+         DO islope=1,nslope
+           if(qsurf_slope(i,nnq,islope).LT.0) then
+             qsurf_slope(i,nnq,islope)=0.
+           endif
          enddo
        enddo
+     enddo
+
+!------------------------
+
+! I   Initialisation
+!    I_a READ run.def 
+!    I_b READ of start_evol.nc and starfi_evol.nc
+!    I_c Subslope parametrisation
+
+!------------------------
+
+!----------------------------Subslope parametrisation definition --------------------- 
 
 !     Define some slope statistics
-
-
-      
-
-       iflat=1
-       DO islope=2,nslope
-         IF(abs(def_slope_mean(islope)).lt. &
-           abs(def_slope_mean(iflat)))THEN
-           iflat = islope
-         ENDIF     
-
-       ENDDO
-       PRINT*,'Flat slope for islope = ',iflat
-       PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
+     iflat=1
+     DO islope=2,nslope
+       IF(abs(def_slope_mean(islope)).lt. &
+         abs(def_slope_mean(iflat))) THEN
+         iflat = islope
+       ENDIF     
+     ENDDO
+
+     PRINT*,'Flat slope for islope = ',iflat
+     PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
+
 ! CO2 max thickness (for glaciers flows)
-
-       allocate(co2_hmax(nslope)
-       if(nslope.eq.7) ! ugly way to implement that ...
+       allocate(co2_hmax(nslope))
+       if(nslope.eq.7) then ! ugly way to implement that ...
 !        CF documentation that explain how the values are computed
          co2_hmax(1) = 1.5
-         co2_hmax(7) = co2_max(1)
+         co2_hmax(7) = co2_hmax(1)
          co2_hmax(2) = 2.4
          co2_hmax(6) = co2_hmax(2)
@@ -475,17 +442,17 @@
        endif
 
-
-
-
+     allocate(flag_co2flow(ngrid,nslope))
+     allocate(flag_co2flow_mesh(ngrid))
 
        flag_co2flow(:,:) = 0.     
        flag_co2flow_mesh(:) = 0.
 
-!------------------------
+!---------------------------- READ GCM data --------------------- 
 
 ! I   Initialisation
-!    I_a READ run.def , run_pem.def
-!    I_b READ starfi_0.nc
-!    I_c READ GCM data and convert to the physical grid 
+!    I_a READ run.def 
+!    I_b READ of start_evol.nc and starfi_evol.nc
+!    I_c Subslope parametrisation
+!    I_d READ GCM data and convert to the physical grid 
 
 !------------------------
@@ -493,6 +460,6 @@
 ! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the GCM run, saving only the minimum value
 
-
-     call nb_time_step_GCM("concat_year_one.nc",timelen)
+     call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
+
      allocate(min_h2o_ice_s_1(iim+1,jjm+1))
      allocate(min_co2_ice_s_1(iim+1,jjm+1))
@@ -506,13 +473,6 @@
      allocate(tsurf_ave(iim+1,jjm+1,nslope))
      allocate(tsurf_ave_yr1(iim+1,jjm+1,nslope))
-     allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope))
-     allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope))
-     allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope))
      allocate(tsurf_ave_phys(ngrid,nslope))
      allocate(tsurf_ave_phys_yr1(ngrid,nslope))
-     allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope))
-     allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
-     allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
-
      allocate(tsurf_GCM_timeseries(iim+1,jjm+1,nslope,timelen))
      allocate(tsurf_phys_GCM_timeseries(ngrid,nslope,timelen))
@@ -520,15 +480,20 @@
      allocate(co2_ice_GCM_slope(iim+1,jjm+1,nslope,timelen))
      allocate(Tsurfave_before_saved(ngrid,nslope))
-     allocate(tsoil_phys_PEM_saved(ngrid,nsoilmx_PEM,nslope))
+     allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope))
+     allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope))
+     allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope))
+     allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope))
+     allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
+     allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
      allocate(delta_co2_adsorbded(ngrid))
-     allocate(q_phys(ngrid,llm,nqtot))
-
-
-
-     call read_data_GCM("concat_year_one.nc", min_h2o_ice_s_1,min_co2_ice_s_1,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_yr1,timelen,min_co2_ice_slope_1,min_h2o_ice_slope_1,&   
+
+     print *, "Downloading data Y1..."
+
+     call read_data_GCM("data_GCM_Y1.nc", min_h2o_ice_s_1,min_co2_ice_s_1,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_yr1,timelen,min_co2_ice_slope_1,min_h2o_ice_slope_1,&   
                        nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope)
 
-
 ! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the GCM run, saving only the minimum value
+
+     print *, "Downloading data Y1 done"
 
      allocate(min_h2o_ice_s_2(iim+1,jjm+1))
@@ -537,12 +502,12 @@
      allocate(min_h2o_ice_slope_2(iim+1,jjm+1,nslope))
 
-     call read_data_GCM("concat_year_two.nc", min_h2o_ice_s_2,min_co2_ice_s_2,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM,timelen,min_co2_ice_slope_2,min_h2o_ice_slope_2, &
+     print *, "Downloading data Y2"
+
+     call read_data_GCM("data_GCM_Y2.nc", min_h2o_ice_s_2,min_co2_ice_s_2,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM,timelen,min_co2_ice_slope_2,min_h2o_ice_slope_2, &
                   nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope)
 
-
-
+     print *, "Downloading data Y2 done"
 
 ! The variables in the dynamic grid are transfered to the physical grid
-
 
      allocate(vmr_co2_gcm_phys(ngrid,timelen))
@@ -552,52 +517,43 @@
      allocate(q_co2_GCM_phys(ngrid,timelen))
      allocate(q_co2_PEM_phys(ngrid,timelen))
-
+     allocate(ps_phys(ngrid))
+     allocate(ps_phys_timeseries(ngrid,timelen))
+     allocate(ps_phys_timeseries_yr1(ngrid,timelen))
 
      CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,vmr_co2_gcm,vmr_co2_gcm_phys)
      CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_h2o_GCM,q_h2o_GCM_phys)
      CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_co2_GCM,q_co2_GCM_phys)
-
-     do nnq=1,nqtot
-
-        call gr_dyn_fi(llm,iip1,jjp1,ngridmx,q(:,:,nnq),q_phys(:,:,nnq))
-     enddo
-
+     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys)
+     call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries)
+     call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1)
+     CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys)
+     CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1)
 
      deallocate(vmr_co2_gcm)
      deallocate(q_h2o_GCM)
      deallocate(q_co2_GCM)
-
-
-    q_co2_PEM_phys(:,:)=  q_co2_GCM_phys(:,:)
-    q_h2o_PEM_phys(:,:)=  q_h2o_GCM_phys(:,:)
-
-
-     allocate(ps_phys(ngrid))
-     allocate(ps_phys_timeseries(ngrid,timelen))
-     allocate(ps_phys_timeseries_yr1(ngrid,timelen))
-
-
-
-
-
-     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys)
-
-
-! for the time series:
-
-
-
-    call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries)
-    call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1)
-    deallocate(ps_GCM)
-    deallocate(ps_GCM_yr1)
-
-
-!!!!!!!!!!!!!!!!!!!!! Initialisation for soil_PEM
+     deallocate(ps_GCM)
+     deallocate(ps_GCM_yr1)
+     deallocate(tsurf_ave)
+     deallocate(tsurf_ave_yr1)
+
+     q_co2_PEM_phys(:,:)=  q_co2_GCM_phys(:,:)
+     q_h2o_PEM_phys(:,:)=  q_h2o_GCM_phys(:,:)
+
+!------------------------
+
+! I   Initialisation
+!    I_a READ run.def 
+!    I_b READ of start_evol.nc and starfi_evol.nc
+!    I_c Subslope parametrisation
+!    I_d READ GCM data and convert to the physical grid 
+!    I_e Initialisation of the PEM variable and soil 
+
+!------------------------
+
+!---------------------------- Initialisation of the PEM soil and values --------------------- 
 
       call end_comsoil_h_PEM
       call ini_comsoil_h_PEM(ngrid,nslope)
-      timestep=year_day*daysec/year_step
-
 
       allocate(ice_depth(ngrid,nslope))
@@ -605,129 +561,63 @@
 
       DO islope = 1,nslope
+        if(soil_pem) then
           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,TI_GCM(:,:,:,islope),TI_GCM_phys(:,:,islope))
+        endif !soil_pem
+        DO t=1,timelen
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t))
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t))
+        enddo
       ENDDO
 
-
+      deallocate(co2_ice_GCM_slope)
+      deallocate(TI_GCM)
+
+  if(soil_pem) then
+
+      deallocate(tsurf_GCM_timeseries)
       call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM)
 
-      deallocate(TI_GCM)
-
-
       DO islope = 1,nslope
-
         DO l=1,nsoilmx
            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope))
+           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope))
+           DO t=1,timelen
+             CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
+           ENDDO
         ENDDO
         DO l=nsoilmx+1,nsoilmx_PEM
            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope))
-        ENDDO
-      ENDDO
-
-
-
-       deallocate(tsoil_ave_yr1)
-
-
-
-
-
-
-
-
-
-
-      DO islope = 1,nslope
-
-        DO l=1,nsoilmx
-           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope))
-        ENDDO
-        DO l=nsoilmx+1,nsoilmx_PEM
            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope))
         ENDDO
       ENDDO
-     
-
-
-       deallocate(tsoil_ave)
-
-
-      DO islope = 1,nslope
-       
-        DO l=1,nsoilmx
-          DO t=1,timelen
-           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
-          ENDDO
-        ENDDO
-      ENDDO
-
+
+      deallocate(tsoil_ave_yr1)
+      deallocate(tsoil_ave)
       deallocate(tsoil_GCM_timeseries)
-     
- 
-
-
-       CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys)
-       
-       CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1)
-
-       deallocate(tsurf_ave)
-       deallocate(tsurf_ave_yr1)
-
-   
-
-
-      DO islope = 1,nslope
-
-        DO t=1,timelen
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t))
-        enddo
-      enddo
-      deallocate(tsurf_GCM_timeseries)
-
-
-  DO islope = 1,nslope
-
-        DO t=1,timelen
-          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t))
-        enddo
-      enddo
-
-
-
-      deallocate(co2_ice_GCM_slope)
-
-
+  endif !soil_pem
 
 !------------------------
 
 ! I   Initialisation
-!    I_a READ run.def , run_pem.def
-!    I_b READ starfi_0.nc
-!    I_c READ GCM data and convert to the physical grid 
-!    I_d Compute tendencies & Save initial situation
-
-!------------------------
-
-
-
-
-
-!----- Compute tendencies
-
+!    I_a READ run.def 
+!    I_b READ of start_evol.nc and starfi_evol.nc
+!    I_c Subslope parametrisation
+!    I_d READ GCM data and convert to the physical grid 
+!    I_e Initialisation of the PEM variable and soil 
+!    I_f Compute tendencies & Save initial situation
+
+!----- Compute tendencies from the PCM run
 
      allocate(tendencies_h2o_ice(iim+1,jjm+1))
      allocate(tendencies_h2o_ice_phys(ngrid))
-
      allocate(tendencies_co2_ice(iim+1,jjm+1))
      allocate(tendencies_co2_ice_phys(ngrid))
-
      allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope))
      allocate(tendencies_co2_ice_phys_slope(ngrid,nslope))
      allocate(tendencies_co2_ice_phys_slope_ini(ngrid,nslope))
-
-
      allocate(tendencies_h2o_ice_slope(iim+1,jjm+1,nslope))
      allocate(tendencies_h2o_ice_phys_slope(ngrid,nslope))
 
-!  Compute the tendencies of the evolution of water ice over the years
+!  Compute the tendencies of the evolution of ice over the years
 
       call compute_tendencies(tendencies_h2o_ice,min_h2o_ice_s_1,&
@@ -740,5 +630,4 @@
              min_co2_ice_slope_2,iim,jjm,ngrid,tendencies_co2_ice_phys_slope,nslope)
 
-
       tendencies_co2_ice_phys_slope_ini(:,:)=tendencies_co2_ice_phys_slope(:,:)
 
@@ -746,6 +635,16 @@
              min_h2o_ice_slope_2,iim,jjm,ngrid,tendencies_h2o_ice_phys_slope,nslope)
 
-     
-!----- Save initial situation
+!------------------------
+
+! I   Initialisation
+!    I_a READ run.def 
+!    I_b READ of start_evol.nc and starfi_evol.nc
+!    I_c Subslope parametrisation
+!    I_d READ GCM data and convert to the physical grid 
+!    I_e Initialisation of the PEM variable and soil 
+!    I_f Compute tendencies & Save initial situation
+!    I_g Save initial PCM situation
+
+!---------------------------- Save initial PCM situation --------------------- 
 
      allocate(initial_h2o_ice(ngrid))
@@ -754,17 +653,23 @@
      allocate(initial_co2_ice_slope(ngrid,nslope))
      allocate(initial_h2o_ice_slope(ngrid,nslope))
-     year_iter=0
 
 ! We save the places where water ice is sublimating
+! We compute the surface of water ice sublimating
+  ini_surf=0.
+  ini_surf_co2=0.
+  ini_surf_h2o=0.
+  Total_surface=0.
   do i=1,ngrid
-      if (tendencies_h2o_ice_phys(i).LT.0) then
+    Total_surface=Total_surface+cell_area(i)
+    if (tendencies_h2o_ice_phys(i).LT.0) then
          initial_h2o_ice(i)=1.
-      else
+         ini_surf=ini_surf+cell_area(i)
+    else
          initial_h2o_ice(i)=0.         
-      endif
+    endif
     do islope=1,nslope
-
       if (tendencies_co2_ice_phys_slope(i,islope).LT.0) then
          initial_co2_ice_sublim_slope(i,islope)=1.
+         ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
       else
          initial_co2_ice_sublim_slope(i,islope)=0.         
@@ -777,4 +682,5 @@
       if (tendencies_h2o_ice_phys_slope(i,islope).LT.0) then
          initial_h2o_ice_slope(i,islope)=1.
+         ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
       else
          initial_h2o_ice_slope(i,islope)=0.         
@@ -783,35 +689,16 @@
   enddo
 
-! We compute the surface of water ice sublimating
-  ini_surf=0.
-  ini_surf_co2=0.
-  ini_surf_h2o=0.
-  Total_surface=0.
-  do i=1,ngrid
-    if (initial_h2o_ice(i).GT.0.5) then
-       ini_surf=ini_surf+cell_area(i)
-    endif
-    do islope=1,nslope
-      if (initial_co2_ice_sublim_slope(i,islope).GT.0.5) then
-         ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
-      endif
-      if (initial_h2o_ice_slope(i,islope).GT.0.5) then
-         ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
-      endif
-    enddo
-    Total_surface=Total_surface+cell_area(i)
-  enddo
-
-     print *, "ini_surf_co2=", ini_surf_co2
-     print *, "ini_surf=", ini_surf
-     print *, "ini_surf_h2o=", ini_surf_h2o
-
-
-     allocate(local_old_press(ngrid))
-     allocate(local_new_press(ngrid))
+     print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
+     print *, "Total initial surface of h2o ice sublimating=", ini_surf
+     print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
+     print *, "Total surface of the planet=", Total_surface
     
-     allocate(zplev_new(ngrid,nlayer+1))
-     allocate(zplev_old(ngrid,nlayer+1))
-
+     allocate(zplev_gcm(ngrid,nlayer+1))
+
+     DO l=1,nlayer+1
+       DO ig=1,ngrid
+         zplev_gcm(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
+       ENDDO
+     ENDDO
 
      global_ave_press_old=0.
@@ -821,66 +708,69 @@
 
      global_ave_press_GCM=global_ave_press_old
-     print *, "global_ave_press_old", global_ave_press_old
-     
-
-
-!----- Time loop
-     print *, "Before timeloop"
-     allocate(flag_co2flow(ngrid,nslope))
-     allocate(flag_co2flow_mesh(ngrid))
-
-     firstcall=.TRUE.
-
-
-
-
+     global_ave_press_new=global_ave_press_old
+     print *, "Initial global average pressure=", global_ave_press_GCM
 
 !------------------------
 
 ! I   Initialisation
-!    I_a READ run.def , run_pem.def
-!    I_b READ starfi_0.nc
-!    I_c READ GCM data and convert to the physical grid 
-!    I_d Compute tendencies & Save initial situation
-!    I_e Read startfi_PEM
-
-!------------------------
-
-
-! now we complete with the PEMstart
+!    I_a READ run.def 
+!    I_b READ of start_evol.nc and starfi_evol.nc
+!    I_c Subslope parametrisation
+!    I_d READ GCM data and convert to the physical grid 
+!    I_e Initialisation of the PEM variable and soil 
+!    I_f Compute tendencies & Save initial situation
+!    I_g Save initial PCM situation
+!    I_h Read the PEMstart
+
+!---------------------------- Read the PEMstart --------------------- 
 
       call pemetat0(.false.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
-      tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries_yr1,ps_phys_timeseries,tsurf_phys_GCM_timeseries,tsoil_phys_PEM_timeseries,&
+      tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsurf_phys_GCM_timeseries,tsoil_phys_PEM_timeseries,&
       tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), co2_adsorbded_phys,delta_co2_adsorbded)
-totmass_adsorbded = 0.
+
+
+    if(soil_pem) then
+     totmass_adsorbded = 0.
+
      do ig = 1,ngrid
       do islope =1, nslope
        do l = 1,nsoilmx_PEM - 1
-
           totmass_adsorbded = totmass_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
           subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
           cell_area(ig)
        enddo
-
       enddo
-
      enddo
 
-     write(*,*) "tot mass in the regolith", totmass_adsorbded
-     stop
+     write(*,*) "Tot mass in the regolith=", totmass_adsorbded
+     deallocate(tsoil_ave_phys_yr1)  
+    endif !soil_pem
      deallocate(tsurf_ave_phys_yr1)
-     deallocate(ps_phys_timeseries_yr1)
-     deallocate(tsoil_ave_phys_yr1)    
-
+     deallocate(ps_phys_timeseries_yr1)  
+
+!------------------------
+
+! I   Initialisation
+!    I_a READ run.def 
+!    I_b READ of start_evol.nc and starfi_evol.nc
+!    I_c Subslope parametrisation
+!    I_d READ GCM data and convert to the physical grid 
+!    I_e Initialisation of the PEM variable and soil 
+!    I_f Compute tendencies & Save initial situation
+!    I_g Save initial PCM situation
+!    I_h Read the PEMstar
+!    I_i Compute orbit criterion
+
+     CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
+
+     if(evol_orbit_pem) then
+       call orbit_param_criterion(year_iter_max)
+     else
+       year_iter_max=Max_iter_pem
+     endif
 
 !--------------------------- END INITIALISATION ---------------------
 
-
-
-
 !---------------------------- RUN ---------------------
-
-
-
 
 !------------------------
@@ -890,92 +780,49 @@
 
 !------------------------
-
-    
-     do while (.true.)
-! II.a.1. global pressure
-     global_ave_press_new=global_ave_press_old
-
+     year_iter=0
+     print *, "Max_iter_pem", Max_iter_pem
+     do while (year_iter.LT.year_iter_max)
+
+! II.a.1. Compute updated global pressure
+     print *, "Recomputing the new pressure..."
      do i=1,ngrid
        do islope=1,nslope
            global_ave_press_new=global_ave_press_new-g*cell_area(i)*tendencies_co2_ice_phys_slope(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
        enddo
-     
-      global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
+       if(soil_pem) then
+         global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
+       endif
      enddo
-! II.a.2. tracers
+
+! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
      allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
-        DO l=1,nlayer+1
-         DO ig=1,ngrid
-          zplev_old(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
-          zplev_old_timeseries(ig,l,:) = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
-         ENDDO
-        ENDDO
-
-     do i=1,ip1jmp1
-
-       ps(i)=ps(i)*global_ave_press_new/global_ave_press_old
-
-     enddo
+     print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..."
+     DO l=1,nlayer+1
+       DO ig=1,ngrid
+         zplev_old_timeseries(ig,l,:) = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
+       ENDDO
+     ENDDO
+
+! II.a.3. Surface pressure timeseries
+     print *, "Recomputing the surface pressure timeserie adapted to the new pressure..."
      do i = 1,ngrid
-       ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_old
        ps_phys_timeseries(i,:) = ps_phys_timeseries(i,:)*global_ave_press_new/global_ave_press_old
      enddo
 
+! II.a.4. New pressure levels timeseries
      allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
-
-        do l=1,nlayer+1
-         do ig=1,ngrid
-          zplev_new(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
-          zplev_new_timeseries(ig,l,:)  = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
-         enddo
-        enddo
-
-
-      print *, "1 is okay"
-
-        DO nnq=1,nqtot
-        if (noms(nnq).NE."co2") then
-          DO l=1,llm-1
-            DO ig=1,ngrid
-               q(ig,l,nnq)=q(ig,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))
-            ENDDO
-            q(:,llm,nnq)=q(:,llm-1,nnq)
-          ENDDO
-        else
-          DO l=1,llm-1
-            DO ig=1,ngrid
-                q(ig,l,nnq)=q(ig,l,nnq)*(zplev_old(ig,l)-zplev_old(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))  &
-                                + (   (zplev_new(ig,l)-zplev_new(ig,l+1))  -       &
-                                      (zplev_old(ig,l)-zplev_old(ig,l+1))     )  / &
-                                      (zplev_new(ig,l)-zplev_new(ig,l+1))
-            ENDDO
-           q(:,llm,nnq)=q(:,llm-1,nnq)
-          ENDDO
-        endif
-        ENDDO
-
-
- DO nnq=1,nqtot
-       DO ig=1,ngrid
-         DO l=1,llm-1
-          if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number" .or. noms(nnq).NE."ccn_number") ) then
-              extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
-              q(ig,l,nnq)=1.
-              q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
-          endif
-          if(q(ig,l,nnq).LT.0) then
-              q(ig,l,nnq)=1E-30
-          endif
-         ENDDO
+     print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..."
+     do l=1,nlayer+1
+       do ig=1,ngrid
+         zplev_new_timeseries(ig,l,:)  = ap(l)  + bp(l)*ps_phys_timeseries(ig,:)
        enddo
-   enddo
-
-
-
-
-  l=1
-  DO ig=1,ngrid
-      DO t = 1, timelen
-
+     enddo
+
+! II.a.5. Tracers timeseries
+      print *, "Recomputing of tracer VMR timeseries for the new pressure..."
+
+     l=1
+     DO ig=1,ngrid
+       DO t = 1, timelen
          q_h2o_PEM_phys(ig,t)=q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
          if(q_h2o_PEM_phys(ig,t).LT.0) then
@@ -985,12 +832,9 @@
            q_h2o_PEM_phys(ig,t)=1.
          endif
-
-   enddo
-  enddo
-
-
-
-  DO ig=1,ngrid
-      DO t = 1, timelen
+       enddo
+     enddo
+
+     DO ig=1,ngrid
+       DO t = 1, timelen
          q_co2_PEM_phys(ig,t)=q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))  &
                                 + (   (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))  -       &
@@ -1004,32 +848,31 @@
          mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
          vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
-      ENDDO
-  ENDDO
+       ENDDO
+     ENDDO
    
      deallocate(zplev_new_timeseries)
      deallocate(zplev_old_timeseries)
 
-
-! II.a.3. Evolution of the ice
-
-
-     call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
-      print *, "4 is okay"
-
-     call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
-      print *, "5 is okay"
-
-
-!------------------------
-
 ! II  Run
 !    II_a update pressure, ice and tracers 
-!    II_b CO2 glaciers flows 
-
-!------------------------
-
-    
-
-
+!    II_b Evolution of the ice
+
+! II.b. Evolution of the ice
+      print *, "Evolution of h2o ice"
+     call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
+
+      print *, "Evolution of co2 ice"
+     call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
+
+!------------------------
+
+! II  Run
+!    II_a update pressure, ice and tracers 
+!    II_b Evolution of the ice
+!    II_c CO2 glaciers flows 
+
+!------------------------
+
+      print *, "Co2 glacier flow"
        DO ig = 1,ngrid
         DO islope = 1,nslope
@@ -1038,7 +881,5 @@
             IF(co2ice_slope(ig,islope).ge.rho_co2*co2_hmax(islope) * &
                   cos(pi*def_slope_mean(islope)/180.)) THEN 
-
 ! Second: determine the flatest slopes possible:
-       
                 IF(islope.gt.iflat) THEN
                   iaval=islope-1
@@ -1048,5 +889,4 @@
                 do while ((iaval.ne.iflat).and.  &
                     (subslope_dist(ig,iaval).eq.0))
-
                   IF(iaval.gt.iflat) THEN
                      iaval=iaval-1
@@ -1054,9 +894,6 @@
                      iaval=iaval+1
                   ENDIF
-             
                 enddo
-
-
-                co2ice_slope(ig,iaval) = co2ice_slope(ig,iaval) + &
+              co2ice_slope(ig,iaval) = co2ice_slope(ig,iaval) + &
                (co2ice_slope(ig,islope) - rho_co2* co2_hmax(islope) *     &
                cos(pi*def_slope_mean(islope)/180.)) *             &
@@ -1065,9 +902,9 @@
                cos(pi*def_slope_mean(islope)/180.)                
 
-                co2ice_slope(ig,islope)=rho_co2*co2_hmax(islope) *        &
-                  cos(pi*def_slope_mean(islope)/180.) 
-
-             flag_co2flow(ig,islope) = 1.
-             flag_co2flow_mesh(ig) = 1.
+              co2ice_slope(ig,islope)=rho_co2*co2_hmax(islope) *        &
+               cos(pi*def_slope_mean(islope)/180.) 
+
+              flag_co2flow(ig,islope) = 1.
+              flag_co2flow_mesh(ig) = 1.
             ENDIF ! co2ice > hmax
           ENDIF ! iflattsoil_lope
@@ -1075,107 +912,93 @@
        ENDDO !ig
 
-      print *, "6 is okay"
-
 !------------------------
 
 ! II  Run
 !    II_a update pressure, ice and tracers 
-!    II_b CO2 glaciers flows 
-!    II_c Update surface and soil temperatures
-
-!------------------------
-
-
-! II_c.1 Update Tsurf
+!    II_b Evolution of the ice
+!    II_c CO2 glaciers flows 
+!    II_d Update surface and soil temperatures
+
+!------------------------
+
+! II_d.1 Update Tsurf
+
+      print *, "Updating the new Tsurf"
       bool_sublim=0
       Tsurfave_before_saved(:,:) = tsurf_ave_phys(:,:)
-       DO ig = 1,ngrid
+      DO ig = 1,ngrid
         DO islope = 1,nslope
-           if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-5) THEN !co2ice disappeared, look for closest point without co2ice
-
-
-              if(latitude_deg(ig).gt.0) then
-                do ig_loop=ig,ngrid
-                  DO islope_loop=islope,iflat,-1
-                     if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
-                            tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
-                            bool_sublim=1
-                            exit
-                     endif
-                  enddo
-                  if (bool_sublim.eq.1) then
-                   exit
+          if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-5) THEN !co2ice disappeared, look for closest point without co2ice
+            if(latitude_deg(ig).gt.0) then
+              do ig_loop=ig,ngrid
+                DO islope_loop=islope,iflat,-1
+                  if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
+                    tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
+                    bool_sublim=1
+                    exit
                   endif
                 enddo
-              else
-               do ig_loop=ig,1,-1
-                  DO islope_loop=islope,iflat
-                     if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
-                            tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
-                            bool_sublim=1
-                            exit
-                     endif
-                  enddo
-                  if (bool_sublim.eq.1) then
-                   exit
+                if (bool_sublim.eq.1) then
+                  exit
+                endif
+              enddo
+            else
+              do ig_loop=ig,1,-1
+                DO islope_loop=islope,iflat
+                  if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
+                    tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
+                    bool_sublim=1
+                    exit
                   endif
                 enddo
-
+                if (bool_sublim.eq.1) then
+                  exit
+                endif
+              enddo
+            endif
+            initial_co2_ice_slope(ig,islope)=0
+            if ((co2ice_slope(ig,islope).lt.1e-5).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then   
+              albedo_slope(ig,1,islope) = albedo_h2o_frost
+              albedo_slope(ig,2,islope) = albedo_h2o_frost
+            else
+              albedo_slope(ig,1,islope) = albedodat(ig)
+              albedo_slope(ig,2,islope) = albedodat(ig) 
+              emiss_slope(ig,islope) = emissiv         
+            endif
+          elseif( co2ice_slope(ig,islope).GT. 1E-5) THEN !Put tsurf as tcond co2
+            ave=0.
+            do t=1,timelen
+              if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-5) then
+                ave = ave + beta/(alpha-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
+              else
+                ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
               endif
-              initial_co2_ice_slope(ig,islope)=0
-
-              if ((co2ice_slope(ig,islope).lt.1e-5).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then
-                 
-               albedo_slope(ig,1,islope) = albedo_h2o_frost
-               albedo_slope(ig,2,islope) = albedo_h2o_frost
-
-              else
-                            
-               albedo_slope(ig,1,islope) = albedodat(ig)
-               albedo_slope(ig,2,islope) = albedodat(ig) 
-               emiss_slope(ig,islope) = emissiv         
-
-              endif
-
-           elseif( co2ice_slope(ig,islope).GT. 1E-5) THEN !Put tsurf as tcond co2
-  
-             ave=0.
-             do t=1,timelen
-              if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-5) then
-                 ave = ave + beta/(alpha-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
-              else
-                 ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
-              endif
-             enddo
-
-             tsurf_ave_phys(ig,islope)=ave/timelen
-           endif
+            enddo
+            tsurf_ave_phys(ig,islope)=ave/timelen
+          endif
         enddo
-       enddo
-
-      print *, "7 is okay"
-   
-     do t = 1,timelen
-     tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:)) 
-     enddo
-! for the start
-     do ig = 1,ngrid
-      do islope = 1,nslope
-     tsurf_slope(ig,islope) =  tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope))
-     enddo
-     enddo
-
-! II_c.2 Update soil temperature
-
-allocate(TI_locslope(ngrid,nsoilmx_PEM))
-allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
-allocate(Tsurf_locslope(ngrid))
-allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
-allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
-
-
-
-
-print *,"8 is ok"
+      enddo
+
+      do t = 1,timelen
+        tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:)) 
+      enddo
+      ! for the start
+      do ig = 1,ngrid
+        do islope = 1,nslope
+          tsurf_slope(ig,islope) =  tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope))
+        enddo
+      enddo
+
+    if(soil_pem) then
+
+! II_d.2 Update soil temperature
+
+  allocate(TI_locslope(ngrid,nsoilmx_PEM))
+  allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
+  allocate(Tsurf_locslope(ngrid))
+  allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
+  allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
+
+  print *,"Updating soil temperature"
 
 ! Soil averaged
@@ -1186,6 +1009,6 @@
      alph_locslope(:,:) = alph_PEM(:,:,islope)
      beta_locslope(:,:) = beta_PEM(:,:,islope)
-     call soil_pem(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
-     call soil_pem(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
+     call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
+     call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
      do t = 1,timelen
        tsoil_phys_PEM_timeseries(:,:,islope,t) =  tsoil_phys_PEM_timeseries(:,:,islope,t) + (Tsoil_locslope(:,:)- tsoil_PEM(:,:,islope))
@@ -1198,41 +1021,34 @@
   enddo
 
-
-      print *, "9 is okay"
-
-
-!!!!!! Check for the time series
-
-   
-      do ig = 1,ngrid
+  print *, "Update of soil temperature done"
+
+! Check Nan in the time series
+     do ig = 1,ngrid
        do islope = 1,nslope
          do iloop = 1,nsoilmx_PEM
- 
-         if(isnan(tsoil_PEM(ig,iloop,islope))) then
-          write(*,*) "is nan tsoil", ig, iloop, islope, tsoil_PEM(ig,iloop,islope)
+           if(isnan(tsoil_PEM(ig,iloop,islope))) then
+            write(*,*) "Problem : There is nan tsoil", ig, iloop, islope, tsoil_PEM(ig,iloop,islope)
            stop 
-         endif
+           endif
          enddo
-        enddo
-      enddo
-
-
-write(*,*) "before ice table, no stop"
-
-! II_c.3 Update the ice table
+       enddo
+     enddo
+
+  write(*,*) "Compute ice table"
+
+! II_d.3 Update the ice table
      call computeice_table(timelen,ngrid,nslope,nsoilmx,nsoilmx_PEM, tsoil_phys_PEM_timeseries,tsurf_phys_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & 
                            ps_phys_timeseries, ice_depth)
        
-      print *, "10 is okay"
-! II_c.3 Update the soil thermal properties
-      call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),ps_phys, &
-        cell_area,ice_depth,TI_PEM)
-
-! II_c.4 Update the mass of the regolith adsorbded
-   
-  
+      print *, "Update soil propreties"
+! II_d.4 Update the soil thermal properties
+      call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
+        ice_depth,TI_PEM)
+
+! II_d.5 Update the mass of the regolith adsorbded
    call regolith_co2adsorption(ngrid,nslope,nsoilmx_PEM,timelen,ps_phys_timeseries,tsoil_PEM,TI_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope, & 
                                co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), q_co2_PEM_phys,q_h2o_PEM_phys,co2_adsorbded_phys,delta_co2_adsorbded)
 
+  endif !soil_pem
 
 !------------------------
@@ -1240,113 +1056,120 @@
 ! II  Run
 !    II_a update pressure, ice and tracers 
-!    II_b CO2 glaciers flows 
-!    II_c Update surface and soil temperatures
-!    II_d Update the tendencies and stopping criterion
-
-!------------------------
-
-
+!    II_b Evolution of the ice
+!    II_c CO2 glaciers flows 
+!    II_d Update surface and soil temperatures
+!    II_e Update the tendencies 
+
+!------------------------
+
+      print *, "Adaptation of the new co2 tendencies to the current pressure"
       call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,tendencies_co2_ice_phys_slope_ini,vmr_co2_gcm_phys,vmr_co2_pem_phys,ps_phys_timeseries,&     
                                global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
 
-          print *, "12 is okay"
+!------------------------
+
+! II  Run
+!    II_a update pressure, ice and tracers 
+!    II_b Evolution of the ice
+!    II_c CO2 glaciers flows 
+!    II_d Update surface and soil temperatures
+!    II_e Update the tendencies 
+!    II_f Checking the stopping criterion
+
+!------------------------
+      call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
+
+      call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
+
       year_iter=year_iter+dt_pem
 
-!We check with we should stop
-
-
-      call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
-
-
-          print *, "13 is okay"
-      call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
-
- print *, "criterion ok"
-
-
-      print *, "STOPPING_water", STOPPING_water, "STOPPING1_water", STOPPING_1_water
-      print *, "STOPPING_co2", STOPPING_co2, "STOPPING1_co2", STOPPING_1_co2
-      print *, "year_iter=", year_iter
+      print *, "Checking all the stopping criterion."
+      if (STOPPING_water) then
+        print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water 
+      endif
+      if (STOPPING_1_water) then
+        print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
+      endif
+      if (STOPPING_co2) then
+        print *, "STOPPING because surface of co2 ice sublimating is too low or global pressure changed too much, see message above", STOPPING_co2
+      endif
+      if (STOPPING_1_co2) then
+        print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2
+      endif
+
       if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2)  then
         exit
+      else
+        print *, "We continue!"
+        print *, "Number of iteration of the PEM : year_iter=", year_iter  
       endif
 
       global_ave_press_old=global_ave_press_new
-      zplev_old=zplev_new
-
-
-      enddo  ! big loop of the pem
+
+      enddo  ! big time iteration loop of the pem
 
 
 !---------------------------- END RUN PEM ---------------------
 
-
-
-
 !---------------------------- OUTPUT ---------------------
 
-
-
-
 !------------------------
 
 ! III Output
-!    III_a Recomp tendencies for the start and startfi
-
-!------------------------
-
-
-! III_a.1 Ice update
+!    III_a Update surface value for the PCM start files
+
+!------------------------
+
+! III_a.1 Ice update (for startfi)
+! Co2 ice
       DO ig = 1,ngrid
-          co2ice(ig) = 0.
-             DO islope = 1,nslope
-                 co2ice(ig) = co2ice(ig) + co2ice_slope(ig,islope) &
-                            * subslope_dist(ig,islope) /          &
-                           cos(pi*def_slope_mean(islope)/180.)
-             ENDDO
+        co2ice(ig) = 0.
+        DO islope = 1,nslope
+          co2ice(ig) = co2ice(ig) + co2ice_slope(ig,islope) &
+                      * subslope_dist(ig,islope) /          &
+                      cos(pi*def_slope_mean(islope)/180.)
+        ENDDO
       ENDDO ! of DO ig=1,ngrid
-
-      DO ig = 1,ngrid !!!!! TO DOC GHANGE IF QSURF ?ü!?!?!?
-          qsurf(ig,igcm_h2o_ice) = 0.
-             DO islope = 1,nslope
-                 qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
-                            * subslope_dist(ig,islope) /          &
-                           cos(pi*def_slope_mean(islope)/180.)
-             ENDDO
+! H2o ice
+      DO ig = 1,ngrid
+        qsurf(ig,igcm_h2o_ice) = 0.
+        DO islope = 1,nslope
+          qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
+                                   * subslope_dist(ig,islope) /          &
+                                  cos(pi*def_slope_mean(islope)/180.)
+        ENDDO
       ENDDO ! of DO ig=1,ngrid
 
-
-
-
-! III_a.2 Tsurf and Tsoil update
-
-
-   call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
-   tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
-
+! III_a.2 Tsurf and Tsoil update (for startfi)
+
+   if(soil_pem) then
+     call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
+     tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
+   else
+      TI_GCM_phys(:,:,:)=TI_GCM_start(:,:,:)
+   endif !soil_pem
 
       DO ig = 1,ngrid
-          tsurf(ig) = 0.
-          
-             DO islope = 1,nslope
-                 tsurf(ig) = tsurf(ig) + tsurf_slope(ig,islope) &
-                            * subslope_dist(ig,islope)      
-             ENDDO
+        tsurf(ig) = 0.
+        DO islope = 1,nslope
+          tsurf(ig) = tsurf(ig) + tsurf_slope(ig,islope) &
+                     * subslope_dist(ig,islope)      
+        ENDDO
       ENDDO ! of DO ig=1,ngrid
 
       DO ig = 1,ngrid
-         DO iloop = 1,nsoilmx
+        DO iloop = 1,nsoilmx
           tsoil(ig,iloop) = 0.
           inertiesoil(ig,iloop) = 0.
-             DO islope = 1,nslope
-                  tsoil(ig,iloop) =  tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
+          DO islope = 1,nslope
+            tsoil(ig,iloop) =  tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
                             * subslope_dist(ig,islope)   
-                  inertiesoil(ig,iloop) =  inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) &
+            inertiesoil(ig,iloop) =  inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) &
                             * subslope_dist(ig,islope)         
-             ENDDO
           ENDDO
+        ENDDO
       ENDDO ! of DO ig=1,ngrid
 
-! III_a.3 Surface optical properties
+! III_a.3 Surface optical properties (for startfi)
 
     DO ig = 1,ngrid
@@ -1357,23 +1180,77 @@
                          *subslope_dist(ig,islope) 
         ENDDO
-       ENDDO
       ENDDO
-
+    ENDDO
 
     DO ig = 1,ngrid
-
-        emis(ig) =0. 
-        DO islope = 1,nslope
-          emis(ig)= emis(ig)+emiss_slope(ig,islope) &
+      emis(ig) =0. 
+      DO islope = 1,nslope
+        emis(ig)= emis(ig)+emiss_slope(ig,islope) &
                          *subslope_dist(ig,islope) 
-        ENDDO
-       ENDDO
-      
-
-
-!------------------------
+      ENDDO
+    ENDDO
+
+! III_a.4 Pressure (for start)
+     do i=1,ip1jmp1
+       ps(i)=ps(i)*global_ave_press_new/global_ave_press_GCM
+     enddo
+
+     do i = 1,ngrid
+       ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_GCM
+     enddo
+
+! III_a.5 Tracer (for start)
+     allocate(zplev_new(ngrid,nlayer+1))
+
+     do l=1,nlayer+1
+       do ig=1,ngrid
+         zplev_new(ig,l) = ap(l)  + bp(l)*ps_phys(ig)
+       enddo
+     enddo
+
+     DO nnq=1,nqtot
+       if (noms(nnq).NE."co2") then
+         DO l=1,llm-1
+           DO ig=1,ngrid
+             q(ig,l,nnq)=q(ig,l,nnq)*(zplev_gcm(ig,l)-zplev_gcm(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))
+           ENDDO
+           q(:,llm,nnq)=q(:,llm-1,nnq)
+         ENDDO
+       else
+         DO l=1,llm-1
+           DO ig=1,ngrid
+             q(ig,l,nnq)=q(ig,l,nnq)*(zplev_gcm(ig,l)-zplev_gcm(ig,l+1))/(zplev_new(ig,l)-zplev_new(ig,l+1))  &
+                             + (   (zplev_new(ig,l)-zplev_new(ig,l+1))  -       &
+                                   (zplev_gcm(ig,l)-zplev_gcm(ig,l+1))     )  / &
+                                   (zplev_new(ig,l)-zplev_new(ig,l+1))
+           ENDDO
+           q(:,llm,nnq)=q(:,llm-1,nnq)
+         ENDDO
+       endif
+     ENDDO
+
+! Conserving the tracers's mass for GCM start files
+     DO nnq=1,nqtot
+       DO ig=1,ngrid
+         DO l=1,llm-1
+           if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number" .or. noms(nnq).NE."ccn_number") ) then
+             extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
+             q(ig,l,nnq)=1.
+             q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
+           endif
+           if(q(ig,l,nnq).LT.0) then
+             q(ig,l,nnq)=1E-30
+           endif
+         ENDDO
+       enddo
+     enddo
+     
+!------------------------
+     if(evol_orbit_pem) then
+      call recomp_orb_param(year_iter)
+     endif
 
 ! III Output
-!    III_a Recomp tendencies for the start and startfi
+!    III_a Update surface value for the PCM start files
 !    III_b Write start and starfi.nc
 
@@ -1386,17 +1263,7 @@
       ztime_fin=0.
 
-         print *, "RVip1jmp1", ip1jmp1
      allocate(p(ip1jmp1,nlayer+1))
-         print *, "ngrid", ngrid
-         print *, "nlayer", nlayer
-          CALL pression (ip1jmp1,ap,bp,ps,p)
-     print *, "masse 1", masse(1,:)
-          CALL massdair(p,masse)
-     print *, "masse 2", masse(1,:)
-
-     do nnq=1,nqtot
-        call gr_fi_dyn(llm,ngridmx,iip1,jjp1,q_phys(:,:,nnq),q(:,:,nnq))
-     enddo
-
+     CALL pression (ip1jmp1,ap,bp,ps,p)
+     CALL massdair(p,masse)
 
       CALL dynredem0("restart_evol.nc", day_ini, phis)
@@ -1404,10 +1271,9 @@
       CALL dynredem1("restart_evol.nc",   &
       		time_0,vcov,ucov,teta,q,masse,ps)
-
+      print *, "restart_evol.nc has been written"
 
 ! III_b.2 WRITE restartfi.nc 
-
       
-             call physdem0("restartfi_evol.nc",longitude,latitude, &
+      call physdem0("restartfi_evol.nc",longitude,latitude, &
                         nsoilmx,ngrid,nlayer,nq,              &
                         ptimestep,pday,0.,cell_area,          &
@@ -1416,6 +1282,5 @@
                         subslope_dist)
 
-
-          call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,  &
+     call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,  &
                      ptimestep,ztime_fin,                        &
                      tsurf,tsoil,co2ice,albedo,emis,             &
@@ -1425,9 +1290,9 @@
                      emiss_slope,qsurf_slope,watercap_slope, TI_GCM_phys)
 
-
+      print *, "restartfi_evol.nc has been written"
 !------------------------
 
 ! III Output
-!    III_a Recomp tendencies for the start and startfi
+!    III_a Update surface value for the PCM start files
 !    III_b Write start and starfi.nc
 !    III_c Write start_pem
@@ -1435,19 +1300,10 @@
 !------------------------
 
-  
-
-         call pemdem1("restartfi_PEM.nc",ztime_fin,nsoilmx_PEM,ngrid,nslope , &
+         call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
                       tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys)
 
-
- 
-
+      print *, "restartfi_PEM.nc has been written"
+      print *, "The PEM had run for ", year_iter, " years."
       print *, "LL & RV : So far so good"
 
-
-
-
-
-
 END PROGRAM pem
-
Index: /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2835)
@@ -1,13 +1,14 @@
    SUBROUTINE pemetat0(startpem_file,filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
-                      tsurf_ave_yr1,tsurf_ave,q_co2,q_h2o,ps_yr1_inst,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, m_co2_regolith_phys,deltam_co2_regolith_phys)
+                      tsurf_ave_yr1,tsurf_ave,q_co2,q_h2o,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, m_co2_regolith_phys,deltam_co2_regolith_phys)
    
-   use iostart_PEM, only:  open_startphy, close_startphy, get_field 
+   use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
    use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
    use comsoil_h, only:  volcapa,inertiedat 
    use ini_soil_mod, only:  ini_icetable
-   use soil_evolution_mod, only: soil_pem_CN
+   use soil_evolution_mod, only: soil_pem,soil_pem_CN
    use adsorption_mod, only : regolith_co2adsorption
+   USE temps_mod_evol, ONLY: year_PEM
+
   implicit none
-
  
   character(len=*), intent(in) :: filename
@@ -23,5 +24,4 @@
   real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
   real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
-  real,intent(in) :: ps_yr1_inst(ngrid,timelen)
   real,intent(in) :: ps_inst(ngrid,timelen)                        ! surface pressure [Pa]
   real,intent(in) :: tsurf_inst(ngrid,nslope,timelen) ! soil (mid-layer) temperature
@@ -34,8 +34,4 @@
 
 ! outputs
-
-
-
-
 
   real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) thermal inertia (SI)
@@ -59,9 +55,9 @@
    real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope)
    real :: tsoil_tmp(ngrid,nsoil_PEM,nslope)
-   real :: alph_tmp(nsoil_PEM-1)
-   real :: beta_tmp(nsoil_PEM-1)
-   real :: co2_ads_prev(ngrid) 
-   real :: ps_ave_yr1(ngrid)
-   real :: ps_ave_yr2(ngrid)
+   real :: alph_tmp(ngrid,nsoil_PEM-1)
+   real :: beta_tmp(ngrid,nsoil_PEM-1)
+   real :: co2_ads_prev(ngrid)
+   real :: year_PEM_read
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -78,17 +74,19 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-
-! 0. Some initializations
- ps_ave_yr1(:) = sum(ps_yr1_inst(:,:),2)/timelen
- ps_ave_yr2(:) = sum(ps_inst(:,:),2)/timelen
-
 !1. Run
-
-
-
 
 if (startpem_file) then
    ! open pem initial state file:
    call open_startphy(filename)
+
+    if(soil_pem) then
+
+   call get_var("Time",year_PEM_read,found)
+   year_PEM=INT(year_PEM_read)
+   if(.not.found) then
+      write(*,*)'PEMetat0: failed loading year_PEM; take default=0'
+   else
+      write(*,*)'year_PEM of startpem=', year_PEM
+   endif
  
 !1. Thermal Inertia
@@ -103,5 +101,4 @@
 
       do ig = 1,ngrid
-
           if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then 
 !!! transition
@@ -110,26 +107,21 @@
             (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
                       ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
-
-          do iloop=nsoil_GCM+2,n_1km
-            TI_PEM(ig,iloop,islope) = TI_breccia
-         enddo      
-         else ! we keep the high ti values
-           do iloop=nsoil_GCM+1,n_1km
-                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
-           enddo 
-
-         endif ! TI PEM and breccia comparison
-
-
+            do iloop=nsoil_GCM+2,n_1km
+              TI_PEM(ig,iloop,islope) = TI_breccia
+            enddo      
+          else ! we keep the high ti values
+            do iloop=nsoil_GCM+1,n_1km
+              TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
+            enddo 
+          endif ! TI PEM and breccia comparison
 !! transition
-             delta = 1000.
-             TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
-            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
-                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
+          delta = 1000.
+          TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
+               (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
+               ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
           do iloop=n_1km+2,nsoil_PEM
             TI_PEM(ig,iloop,islope) = TI_bedrock
-         enddo   
-      enddo
-
+          enddo   
+      enddo
    else
      do iloop = nsoil_GCM+1,nsoil_PEM
@@ -141,6 +133,4 @@
 print *,'PEMETAT0: THERMAL INERTIA DONE'
 
-
-
 ! b. Special case for inertiedat, inertiedat_PEM
 call get_field("inertiedat_PEM",inertiedat_PEM,found)
@@ -148,7 +138,5 @@
  do iloop = 1,nsoil_GCM
    inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
- enddo
-
-  
+ enddo 
 !!! zone de transition
 delta = 50.
@@ -159,5 +147,4 @@
                       ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
 
-
  do iloop = nsoil_GCM+2,n_1km 
        inertiedat_PEM(ig,iloop) = TI_breccia
@@ -179,6 +166,4 @@
 enddo
 
-
-
  do iloop = n_1km+2, nsoil_PEM
    do ig = 1,ngrid
@@ -187,5 +172,4 @@
  enddo
 endif ! not found
-
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -221,8 +205,8 @@
 
     tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
-    call soil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
-    call soil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
+    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
+    call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
     tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
-    call soil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
+    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)
 
 
@@ -237,9 +221,7 @@
         enddo
      enddo
-
-
-
     
 ENDDO
+
 print *,'PEMETAT0: SOIL TEMP  DONE'
 
@@ -248,5 +230,4 @@
 !3. Ice Table
 
- 
    call get_field("ice_table",ice_table,found)
    if(.not.found) then
@@ -254,7 +235,5 @@
       write(*,*)'will reconstruct the values of ice table'
 
-
       call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:),  tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
-
 
    else 
@@ -262,13 +241,9 @@
      call computeice_table(timelen,ngrid,nslope,nsoil_PEM,tsoil_inst,tsurf_inst,q_co2,q_h2o, ps_inst, ice_table)
 
-
    endif
 
 print *,'PEMETAT0: ICE TABLE  DONE'
 
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 !4. CO2 Adsorption
@@ -299,12 +274,13 @@
 print *,'PEMETAT0: CO2 adsorption done  '
 
-
+endif ! soil_pem
 
 call close_startphy
 
-
 else !No startfi, let's build all by hand
 
-
+    year_PEM=0
+
+    if(soil_pem) then
 
 !a) Thermal inertia
@@ -326,7 +302,5 @@
                   TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
            enddo 
-
          endif
-
 
 !! transition
@@ -339,7 +313,5 @@
          enddo   
       enddo
-
 enddo
-
 
  do iloop = 1,nsoil_GCM
@@ -440,16 +412,12 @@
 print *,'PEMETAT0: CO2 adsorption done  '
 
-
-
-
-
-
-
-
-
+endif !soil_pem
 
 endif ! of if (startphy_file)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    
+     if(soil_pem) then
+
     DO ig = 1,ngrid
       DO islope = 1,nslope
@@ -459,6 +427,4 @@
       ENDDO
     ENDDO
-
-
 
 !! small test
@@ -474,4 +440,6 @@
      ENDDO
 
+     endif!soil_pem
+
       write(*,*) "construction ok, no nan" 
      
Index: /trunk/LMDZ.COMMON/libf/evolution/pemredem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/pemredem.F90	(revision 2835)
@@ -8,5 +8,4 @@
                          day_ini,time,nslope,def_slope,subslope_dist)
 ! create physics restart file and write time-independent variables
-!  use comsoil_h, only: inertiedat, volcapa, mlayer
   use comsoil_h_PEM, only: mlayer_PEM,fluxgeo,inertiedat_PEM
   use iostart_PEM, only : open_restartphy, close_restartphy, & 
@@ -30,6 +29,4 @@
   real, intent(in) :: def_slope(nslope+1) !boundaries for bining of the slopes
   real, intent(in) :: subslope_dist(ngrid,nslope) !undermesh statistics
-
-
   
   ! Create physics start file
@@ -59,5 +56,5 @@
 end subroutine pemdem0
 
-subroutine pemdem1(filename,time,nsoil_PEM,ngrid,nslope, &
+subroutine pemdem1(filename,year_iter,nsoil_PEM,ngrid,nslope, &
                     tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table,m_co2_regolith)
   ! write time-dependent variable to restart file
@@ -66,4 +63,6 @@
 
   use comsoil_h_PEM, only: mlayer_PEM,fluxgeo,inertiedat_PEM
+  USE temps_mod_evol, ONLY: year_PEM
+  use soil_evolution_mod, only: soil_pem
                 
   implicit none
@@ -74,5 +73,4 @@
   integer,intent(in) :: nsoil_PEM
   integer,intent(in) :: ngrid
-  real,intent(in) :: time
   real,intent(IN) :: tsoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
   real,intent(IN) :: inertiesoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
@@ -82,5 +80,7 @@
   integer :: islope
   CHARACTER*2 :: num  
+  integer, intent(IN) :: year_iter
 
+  real :: year_tot
  
   integer :: iq
@@ -88,5 +88,6 @@
   ! indexes of water vapour & water ice tracers (if any):
 
-  
+  year_tot=real(year_iter+year_PEM)
+  print *, "Year of simulation=", year_tot  
 
   ! Open file
@@ -95,7 +96,7 @@
   ! First variable to write must be "Time", in order to correctly
   ! set time counter in file
-  call put_var("Time","Temps de simulation",time)
-  
- 
+  call put_var("Time","Year of simulation",year_tot)
+
+    if(soil_pem) then
   
   ! Multidimensionnal variables (undermesh slope statistics)
@@ -103,19 +104,20 @@
   write(num,fmt='(i2.2)') islope
   call put_field("tsoil_PEM_slope"//num,"Soil temperature by slope type", &
-                 tsoil_slope_PEM(:,:,islope),time)
+                 tsoil_slope_PEM(:,:,islope),year_tot)
   call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type", &
-                 inertiesoil_slope_PEM(:,:,islope),time)
+                 inertiesoil_slope_PEM(:,:,islope),year_tot)
 
   call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith", &
-                    m_co2_regolith(:,:,islope), time)
+                    m_co2_regolith(:,:,islope), year_tot)
 
 ENDDO
 
   call put_field("ice_table","Depth of ice table", &
-                 ice_table,time)
+                 ice_table,year_tot)
 
   call put_field("inertiedat_PEM","Thermal inertie of PEM ", &
-                 inertiedat_PEM,time)
+                 inertiedat_PEM,year_tot)
 
+   endif !soil_pem
 
   ! Close file
Index: /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2835)
@@ -9,4 +9,5 @@
                         nf90_inquire_dimension,nf90_close
       use comsoil_h, only: nsoilmx
+      USE soil_evolution_mod, ONLY: soil_pem
 
       IMPLICIT NONE
@@ -23,4 +24,5 @@
 ! Arguments:
   CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
+  INTEGER, INTENT(IN) :: timelen                   ! number of times stored in the file
 
   INTEGER :: iim_input,jjm_input,nlayer,nslope
@@ -34,15 +36,9 @@
   REAL, INTENT(OUT) ::  min_co2_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
   REAL, INTENT(OUT) ::  min_h2o_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
-!  REAL, ALLOCATABLE ::  vmr_co2_gcm(:,:,:)                     !!!!vmr_co2_phys_gcm(iim_input+1,jjm_input+1,timelen)
-  REAL, INTENT(OUT) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,2676)                     !!!!vmr_co2_phys_gcm(iim_input+1,jjm_input+1,timelen)
-!  REAL, ALLOCATABLE ::  q_h2o_GCM(:,:,:)
-  REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,2676)
-  REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,2676)
-!  REAL, ALLOCATABLE ::  q_co2_GCM(:,:,:)
+  REAL, INTENT(OUT) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,timelen)      !!!!vmr_co2_phys_gcm(iim_input+1,jjm_input+1,timelen)
+  REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,timelen)
+  REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,timelen)
   REAL, ALLOCATABLE ::  q1_co2_GCM(:,:,:)
-!  real, INTENT(OUT) ::  vmr_co2_phys_gcm(:,:)                  !!!!vmr_co2_gcm(ngrid,timelen)
-!  REAL, ALLOCATABLE ::  ps_GCM(:,:,:)
-  REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,2676)
-
+  REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,timelen)
 
 !SOIL
@@ -50,10 +46,10 @@
   REAL, INTENT(OUT) ::  tsoil_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average soil temperature of the concatenated file
 
-  REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,2676) ! Surface temperature of the concatenated file
-  REAL , INTENT(OUT) ::  tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,2676) ! Soil temperature of the concatenated file
-
-  REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,2676) ! Thermal Inertia  of the concatenated file
+  REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen) ! Surface temperature of the concatenated file
+  REAL , INTENT(OUT) ::  tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file
+
+  REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Thermal Inertia  of the concatenated file
   REAL, INTENT(OUT) ::  TI_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average Thermal Inertia  of the concatenated file
-  REAL, INTENT(OUT) ::  co2_ice_slope(iim_input+1,jjm_input+1,nslope,2676) ! Minimum of co2_ice slope of the year
+  REAL, INTENT(OUT) ::  co2_ice_slope(iim_input+1,jjm_input+1,nslope,timelen) ! Minimum of co2_ice slope of the year
 !===============================================================================
 !   Local Variables 
@@ -65,5 +61,4 @@
 
   REAL,ALLOCATABLE :: time(:) ! times stored in start
-  INTEGER :: timelen ! number of times stored in the file
   INTEGER :: indextime ! index of selected time
 
@@ -74,5 +69,4 @@
   INTEGER :: islope
   CHARACTER*2 :: num
-
 
 !-----------------------------------------------------------------------
@@ -84,70 +78,43 @@
       B=1/m_noco2
 
+      allocate(co2_ice_s(iim+1,jjm+1,timelen))
+      allocate(q1_co2_GCM(iim+1,jjm+1,timelen))
+      allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
+      allocate(h2o_ice_s(iim+1,jjm+1,timelen))
+
+  print *, "Opening ", fichnom, "..."
+
 !  Open initial state NetCDF file
   var=fichnom
   CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
 
-      ierr = nf90_inq_varid (fID, "temps", vID)
-      IF (ierr .NE. nf90_noerr) THEN
-        write(*,*)"pemetat0: Le champ <temps> est absent"
-        write(*,*)"pemetat0: J essaie <Time>"
-        ierr = nf90_inq_varid (fID, "Time", vID)
-        IF (ierr .NE. nf90_noerr) THEN
-           write(*,*)"pemetat0: Le champ <Time> est absent"
-           write(*,*)trim(nf90_strerror(ierr))
-     print *, "ICIIII0"
-           CALL ABORT_gcm("pemetat0", "", 1)
-        ENDIF
-        ! Get the length of the "Time" dimension
-     print *, "ICIIIITIME"
-        ierr = nf90_inq_dimid(fID,"Time",vID)
-        ierr = nf90_inquire_dimension(fID,vID,len=timelen)
-      ELSE   
-        ! Get the length of the "temps" dimension
-     print *, "ICIIIITEMPS"
-        ierr = nf90_inq_dimid(fID,"temps",vID)
-        ierr = nf90_inquire_dimension(fID,vID,len=timelen)
-      ENDIF
-
-     print *, "ICIIII"
-
-      allocate(co2_ice_s(iim+1,jjm+1,timelen))
-
-     print *, "ICIIIIAAAA"
-
-!      allocate(q_co2_GCM(iim+1,jjm+1,timelen))
-
-     print *, "ICIIIIBBBB"
-
-!      allocate(q_h2o_GCM(iim+1,jjm+1,timelen))
-
-     print *, "ICIIIICCCC"
-
-      allocate(q1_co2_GCM(iim+1,jjm+1,timelen))
-
-     print *, "ICIIII2"
-
-
-
-      allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
-      allocate(h2o_ice_s(iim+1,jjm+1,timelen))
-          print *, "ICIIII3"
+     print *, "Downloading data for h2oice ..."
 
 ! Get h2o_ice_s of the concatenated file
   CALL get_var3("h2o_ice_s"   ,h2o_ice_s)
 
-  print *, "A"
+     print *, "Downloading data for h2oice done"
+     print *, "Downloading data for co2ice ..."
 
   CALL get_var3("co2ice"   ,co2_ice_s)
+
+     print *, "Downloading data for co2ice done"
+     print *, "Downloading data for vmr co2..."
+
   CALL get_var3("co2_cropped"   ,q_co2_GCM)
+
+     print *, "Downloading data for vmr co2 done"
+     print *, "Downloading data for vmr h20..."
+
   CALL get_var3("h2o_cropped"   ,q_h2o_GCM)
 
-  print *, "B"
+     print *, "Downloading data for vmr h2o done"
+     print *, "Downloading data for surface pressure ..."
 
   CALL get_var3("ps"   ,ps_GCM)
 
-  print *, "C"
-
-  print *, "nslope", nslope
+     print *, "Downloading data for surface pressure done"
+     print *, "nslope=", nslope
+     print *, "Downloading data for co2ice_slope ..."
 
 DO islope=1,nslope
@@ -156,5 +123,6 @@
 ENDDO
 
-  print *, "co2ice_slope"
+     print *, "Downloading data for co2ice_slope done"
+     print *, "Downloading data for h2o_ice_s_slope ..."
 
 DO islope=1,nslope
@@ -163,5 +131,6 @@
 ENDDO
 
-  print *, "h2o_ice_s_slope"
+     print *, "Downloading data for h2o_ice_s_slope done"
+     print *, "Downloading data for tsurf_slope ..."
 
 DO islope=1,nslope
@@ -170,5 +139,9 @@
 ENDDO
 
-  print *, "tsurf_slope"
+     print *, "Downloading data for tsurf_slope done"
+
+     if(soil_pem) then
+
+     print *, "Downloading data for tsoil_slope ..."
 
 DO islope=1,nslope
@@ -177,5 +150,6 @@
 ENDDO
 
-  print *, "tsoil_slope"
+     print *, "Downloading data for tsoil_slope done"
+     print *, "Downloading data for inertiesoil_slope ..."
 
 DO islope=1,nslope
@@ -184,27 +158,30 @@
 ENDDO
 
-  print *, "inertiesoil_slope"
-
-
-
-
+     print *, "Downloading data for inertiesoil_slope done"
+
+  endif
 
 ! Compute the minimum over the year for each point
+  print *, "Computing the min of h2o_ice"
   min_h2o_ice_s(:,:)=minval(h2o_ice_s,3)
+  print *, "Computing the min of co2_ice"
   min_co2_ice_s(:,:)=minval(co2_ice_s,3)
 
+  print *, "Computing the min of h2o_ice_slope"
+  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
+  print *, "Computing the min of co2_ice_slope"
   min_co2_ice_slope(:,:,:)=minval(co2_ice_slope,4)
-  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
 
 !Compute averages
 
-!  DO i=1,timelen
+    print *, "Computing average of tsurf"
     tsurf_ave(:,:,:)=SUM(tsurf_gcm(:,:,:,:),4)/timelen
+
+  if(soil_pem) then
+    print *, "Computing average of tsoil"
     tsoil_ave(:,:,:,:)=SUM(tsoil_gcm(:,:,:,:,:),5)/timelen
+    print *, "Computing average of TI"
     TI_ave(:,:,:,:)=SUM(TI_gcm(:,:,:,:,:),5)/timelen
-!  ENDDO
-
-
-
+  endif
 
 ! By definition, a density is positive, we get rid of the negative values
@@ -229,7 +206,4 @@
     DO j = 1, jjm+1
       DO t = 1, timelen
-!         q1_co2_GCM(i,j,t)=q_co2_GCM(i,j,t)
-!         ps_GCM(i,j,t)=ps(i)
-!c             Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
          if (q_co2_GCM(i,j,t).LT.0) then
               q_co2_GCM(i,j,t)=1E-10
@@ -248,9 +222,4 @@
   ENDDO
 
-
-
-
-
-
   CONTAINS
 
Index: unk/LMDZ.COMMON/libf/evolution/read_data_humidity.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/read_data_humidity.F90	(revision 2834)
+++ 	(revision )
@@ -1,139 +1,0 @@
-!
-! $Id $
-!
-SUBROUTINE read_data_humidity(fichnom,iim_input,jjm_input,nlayer,timelen, nslope,q_h2o_GCM,q_co2_GCM,ps_GCM,tsurf_gcm)
-
-      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
-                        nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
-                        nf90_inquire_dimension,nf90_close
-      use comsoil_h, only: nsoilmx
-
-      IMPLICIT NONE
-
-!=======================================================================
-!
-! Read initial confitions file to compute the ice table
-!
-!=======================================================================
-
-  include "dimensions.h"
-
-!===============================================================================
-! Arguments:
-  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
-
-  INTEGER, INTENT(IN) :: iim_input,jjm_input,nlayer,nslope,timelen
-
-! atmosphere
-
-  REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,timelen)
-  REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,timelen)
-
-  REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,timelen)
-
-
-! Surface
-  REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen) ! Surface temperature of the concatenated file
-
-
-
-!===============================================================================
-!   Local Variables 
-  CHARACTER(LEN=256) :: msg, var, modname
-  INTEGER,PARAMETER :: length=100
-  INTEGER :: iq, fID, vID, idecal
-  INTEGER :: ierr
-  CHARACTER(len=12) :: start_file_type="earth" ! default start file type
-
-  REAL,ALLOCATABLE :: time(:) ! times stored in start
-  INTEGER :: indextime ! index of selected time
-
-  INTEGER :: edges(4),corner(4)
-  INTEGER :: i,j,t
-  real,save :: m_co2, m_noco2, A , B, mmean
-
-  INTEGER :: islope
-  CHARACTER*2 :: num
-
-
-!-----------------------------------------------------------------------
-  modname="pemetat0"
-
-      m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
-      m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
-      A =(1/m_co2 - 1/m_noco2)
-      B=1/m_noco2
-
-!  Open initial state NetCDF file
-  var=fichnom
-  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
-
-  CALL get_var3("co2_cropped" ,q_co2_GCM)
-  CALL get_var3("h2o_cropped" ,q_h2o_GCM)
-  CALL get_var3("ps"   ,ps_GCM)
-
-
-  DO islope=1,nslope
-    write(num,fmt='(i2.2)') islope
-    call get_var3("tsurf_slope"//num,tsurf_gcm(:,:,islope,:))
-  ENDDO
-
-  print *, "tsurf_slope"
-
-
-
-
-
-  CONTAINS
-
-SUBROUTINE check_dim(n1,n2,str1,str2)
-  INTEGER,          INTENT(IN) :: n1, n2
-  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
-  CHARACTER(LEN=256) :: s1, s2
-  IF(n1/=n2) THEN
-    s1='value of '//TRIM(str1)//' ='
-    s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
-    WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(s1),n1,TRIM(s2),n2
-    CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
-  END IF
-END SUBROUTINE check_dim
-
-
-SUBROUTINE get_var1(var,v)
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-END SUBROUTINE get_var1
-
-
-SUBROUTINE get_var3(var,v) ! on U grid
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:,:,:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-
-END SUBROUTINE get_var3
-
-SUBROUTINE get_var4(var,v) 
-  CHARACTER(LEN=*), INTENT(IN)  :: var
-  REAL,             INTENT(OUT) :: v(:,:,:,:)
-  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
-  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
-END SUBROUTINE get_var4
-
-SUBROUTINE err(ierr,typ,nam)
-  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
-  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
-  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
-  IF(ierr==NF90_NoERR) RETURN
-  SELECT CASE(typ)
-    CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
-    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
-    CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
-    CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
-  END SELECT
-  CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
-END SUBROUTINE err
-
-END SUBROUTINE read_data_humidity
Index: /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 2835)
+++ /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 2835)
@@ -0,0 +1,72 @@
+      MODULE recomp_orb_param_mod
+
+      IMPLICIT NONE
+
+      CONTAINS
+
+      SUBROUTINE recomp_orb_param(final_iter)
+
+      USE temps_mod_evol, ONLY: year_bp_ini, year_PEM
+      USE planete_h, ONLY: e_elips, obliquit, timeperi
+      USE comcstfi_h, only: pi
+      USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
+   				end_lask_param_mod,last_ilask
+
+      IMPLICIT NONE
+
+!--------------------------------------------------------
+! Input Variables
+!--------------------------------------------------------
+
+      integer,intent(in) :: final_iter ! Number of year iteration of the PEM 
+
+!--------------------------------------------------------
+! Output Variables
+!--------------------------------------------------------
+
+!--------------------------------------------------------
+! Local variables 
+!--------------------------------------------------------
+
+      real :: Year	                ! Year of the simulation
+      real :: timeperi_ls	        ! time of peri in ls
+      integer ilask			! Loop variable
+
+
+      ! **********************************************************************
+      ! 0. Initializations
+      ! **********************************************************************
+
+
+          Year=year_bp_ini+year_PEM+final_iter
+
+          timeperi_ls=timeperi*360/2/pi
+
+!          Year=-953200
+          print *, "recomp_orb_param, Old year=", year_bp_ini+year_PEM
+          print *, "recomp_orb_param, New year=", Year
+          print *, "recomp_orb_param, Old obl=", obliquit
+          print *, "recomp_orb_param, Old ex=", e_elips
+          print *, "recomp_orb_param, Old lsp=", timeperi_ls
+          print *, "recomp_orb_param, Old timeperi=", timeperi
+
+        do ilask=last_ilask,1,-1
+           if(yearlask(ilask) .GT.Year) then
+              obliquit=oblask(ilask+1)+(oblask(ilask)-oblask(ilask+1))*(Year-yearlask(ilask+1))/(yearlask(ilask)-yearlask(ilask+1))
+              e_elips=exlask(ilask+1)+(exlask(ilask)-exlask(ilask+1))*(Year-yearlask(ilask+1))/(yearlask(ilask)-yearlask(ilask+1))
+              timeperi_ls=lsplask(ilask+1)+(lsplask(ilask)-lsplask(ilask+1))*(Year-yearlask(ilask+1))/(yearlask(ilask)-yearlask(ilask+1))
+              exit
+           endif
+        enddo
+        timeperi=timeperi_ls*2*pi/360
+
+       print *, "recomp_orb_param, Final year of the PEM run:", Year
+       print *, "recomp_orb_param, New obl=", obliquit
+       print *, "recomp_orb_param, New ex=", e_elips
+       print *, "recomp_orb_param, New timeperi=", timeperi
+
+      END SUBROUTINE recomp_orb_param
+
+!********************************************************************************   
+      
+      END MODULE recomp_orb_param_mod
Index: /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90	(revision 2835)
@@ -20,14 +20,10 @@
   REAL, INTENT(in) ::  vmr_co2_pem(ngrid,timelen)                ! physical point field : Volume mixing ratio of co2 in the first layer
   REAL, intent(in) :: ps_GCM_2(ngrid,timelen)                 ! physical point field : Surface pressure in the GCM
-!  REAL, INTENT(in) ::  q_co2_GCM(ngrid)                ! physical point field : Density of co2 in the first layer
-!  REAL, intent(in) :: ps_GCM(ngrid)                 ! physical point field : Density of co2 in the first layer
   REAL, intent(in) :: global_ave_press_GCM
   REAL, intent(in) :: global_ave_press_new
   REAL, intent(in) ::  tendencies_co2_ice_phys_ini(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
 
-
 !   OUTPUT
   REAL, intent(inout) ::  tendencies_co2_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
-
 
 !   local:
@@ -45,29 +41,14 @@
   coef=669*24*3600*eps*sigma/L
 
-  print *, "coef", coef
-  print *, "global_ave_press_GCM", global_ave_press_GCM
-  print *, "global_ave_press_new", global_ave_press_new
-
 ! Evolution of the water ice for each physical point
   do i=1,ngrid
     do islope=1,nslope
-       ave=0.
-    do t=1,timelen
-
-!       write(*,*)'i,t=',i,t,islope, alpha,beta,ave,vmr_co2_gcm(i,t),vmr_co2_pem(i,t),ps_GCM_2(i,t),global_ave_press_GCM,global_ave_press_new
-       ave=ave+(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100)))**4  &
+      ave=0.
+      do t=1,timelen
+        ave=ave+(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100)))**4  &
               -(beta/(alpha-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*global_ave_press_GCM/global_ave_press_new/100)))**4
+      enddo
+      tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys_ini(i,islope)-coef*ave/timelen
     enddo
-!    print *, "i", i
-!    print *, "tendencies_co2_ice_phys_ini bef", tendencies_co2_ice_phys_ini(i,islope)
-!    tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys_ini(i,islope)+coef*ave/timelen
-    tendencies_co2_ice_phys(i,islope)=tendencies_co2_ice_phys_ini(i,islope)-coef*ave/timelen
-
-!    print *, "tendencies after", tendencies_co2_ice_phys(i,islope)
-!    print *, "ave", ave
-!    print *, "timelen", timelen
-!    print *, "vmr_co2_pem(i,t)*ps_GCM_2(i,t)", vmr_co2_pem(i,t)*ps_GCM_2(i,t)
-!    print *, "-------------------"
-  enddo
   enddo
 
Index: /trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/soil_TIfeedback_PEM.F90	(revision 2835)
@@ -1,5 +1,4 @@
       SUBROUTINE soil_TIfeedback_PEM(ngrid,nsoil,icecover,   newtherm_i)
 
-!      use tracer_mod, only:  rho_ice
       use comsoil_h_PEM, only: layer_PEM, inertiedat_PEM
 
@@ -24,5 +23,4 @@
 !  Author: Adapted from J.-B. Madeleine Mars 2008 ( Updated November 2012) by LL, 2022
 !=======================================================================
- 
 
 !Local variables
@@ -42,11 +40,9 @@
 
       REAL ,INTENT(IN):: icecover(ngrid)         ! tracer on the surface (kg.m-2)
-                                       
-                                        ! water ice
+
 !Outputs
 !-------
 
       REAL,INTENT(INOUT) :: newtherm_i(ngrid,nsoil)    ! New soil thermal inertia
-
 
       prev_thermi(:,:) = newtherm_i(:,:)
Index: /trunk/LMDZ.COMMON/libf/evolution/soil_evolution_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_evolution_mod.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/soil_evolution_mod.F90	(revision 2835)
@@ -1,9 +1,9 @@
       MODULE soil_evolution_mod
 
+      IMPLICIT NONE
 
-      IMPLICIT NONE
+      LOGICAL soil_pem  ! True by default, to run with the subsurface physic. Read in pem.def
       
       CONTAINS    
-
 
      subroutine soil_pem_CN(ngrid,nsoil, &
@@ -25,8 +25,4 @@
 
 #include "dimensions.h"
-!#include "dimphys.h"
-
-!#include"comsoil.h"
-
 
 !-----------------------------------------------------------------------
@@ -44,6 +40,4 @@
       real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
 
-
-      
 ! local variables:
       integer ig,ik,isoil
@@ -58,10 +52,8 @@
       real :: Tcol(nsoil)
 
-
 do ig = 1,ngrid
 
    Tcol(:) = tsoil(ig,:)
 !0. Build thermal properties of the ground
-
 
   do isoil = 1,nsoil
@@ -87,6 +79,4 @@
 beta_pem(1) = k_soil(1)*timestep/(rhoc(1)*layer_PEM(2)*layer_PEM(1))
 
-
-
 !1.c At the bottom: dT/dz = -Fgeo/k.
 beta_pem(nsoil) = timestep*k_soil(nsoil)/(2.*rhoc(nsoil)*(layer_PEM(isoil) - layer_PEM(isoil-1))**2)
@@ -107,12 +97,7 @@
   timestep/rhoc(nsoil)*fluxgeo/(layer_PEM(nsoil)-layer_PEM(nsoil-1)) 
 
-
-
 !2. Update by tridiagonal inversion 
   call tridag(d1,d2,d3,re,Tcol,nsoil) 
-  
-
-
-      tsoil(ig,:) = Tcol(:)
+  tsoil(ig,:) = Tcol(:)
   enddo
 
Index: /trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/soil_pem.F90	(revision 2835)
@@ -1,3 +1,3 @@
-      subroutine soil_pem(ngrid,nsoil,firstcall, &
+      subroutine soil_pem_routine(ngrid,nsoil,firstcall, &
                therm_i,                          &
                timestep,tsurf,tsoil,alph_PEM,beta_PEM)
@@ -20,8 +20,4 @@
 
 #include "dimensions.h"
-!#include "dimphys.h"
-
-!#include"comsoil.h"
-
 
 !-----------------------------------------------------------------------
@@ -40,14 +36,7 @@
       real,intent(inout) :: alph_PEM(ngrid,nsoil-1)
       real,intent(inout) :: beta_PEM(ngrid,nsoil-1)
-
-
-
-
       
 ! local variables:
-      integer ig,ik
-
-
-     
+      integer ig,ik    
 
 ! 0. Initialisations and preprocessing step
@@ -57,9 +46,7 @@
       do ig=1,ngrid
         do ik=0,nsoil-1
-	  mthermdiff_PEM(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
-    
+	  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
@@ -69,5 +56,4 @@
                      +(mlayer_PEM(ik)-layer_PEM(ik))*mthermdiff_PEM(ig,ik-1))  &
                          /(mlayer_PEM(ik)-mlayer_PEM(ik-1))
-
 	enddo
       enddo
@@ -100,43 +86,23 @@
 	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)* &  
+        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
-
-
-
+          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)
@@ -157,6 +123,4 @@
       enddo
 
-
-
       end
 
Index: /trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F	(revision 2835)
@@ -56,6 +56,4 @@
       real alpha,lay1 ! coefficients for building layers
       real xmin,xmax ! to display min and max of a field
-
-
       
 !======================================================================
@@ -65,5 +63,4 @@
 
 ! 1.4 Build mlayer(), if necessary
-!      if (interpol) then
       ! default mlayer distribution, following a power law:
       !  mlayer(k)=lay1*alpha**(k-1/2)
@@ -73,5 +70,4 @@
 	  mlayer_PEM(iloop)=lay1*(alpha**(iloop-0.5))
 	enddo
-!      endif
 
 ! 1.5 Build layer(); following the same law as mlayer()
@@ -84,10 +80,6 @@
       enddo
 
-
-
-
 ! 2. Thermal inertia (note: it is declared in comsoil_h)
 ! ------------------
-
 
       do ig = 1,ngrid
@@ -103,3 +95,4 @@
         enddo
       enddo
+
       end subroutine soil_settings_PEM
Index: /trunk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90	(revision 2835)
@@ -3,14 +3,10 @@
 IMPLICIT NONE  
 
-  INTEGER   nyear           !     nyear   : Maximun number of year over which the PEM can interpolate
+  INTEGER   year_bp_ini     !     year_bp_ini : Initial year of the simulation of the PEM (in evol.def)
   INTEGER   dt_pem          !     dt_pem  : in years, the time step used by the PEM    
   REAL      alpha_criterion !     alpha_criterion : percentage of change before stopping the PEM 
-
-!$OMP THREADPRIVATE(nyear)        
-
-!WARNING: when adding a threadprivate variable in this module
-!        do not forget to add it to the copyin clause when opening an OpenMP
-!        parallel section. e.g. in gcm before call leapfrog_loc and/or
-!        possibly in iniphysiq
+  INTEGER   year_PEM        !     year written in startfiPEM.nc
+  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 ob_ex_lsp.asc, read in evol.def
 
 END MODULE temps_mod_evol
Index: /trunk/LMDZ.COMMON/libf/evolution/update_soil.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/update_soil.F90	(revision 2834)
+++ /trunk/LMDZ.COMMON/libf/evolution/update_soil.F90	(revision 2835)
@@ -1,5 +1,4 @@
-   SUBROUTINE update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,ps_new,&
-                          cellarea,ice_depth,TI_PEM)
-
+   SUBROUTINE update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,p_avg_new,&
+                          ice_depth,TI_PEM)
 
  USE comsoil_h, only:  inertiedat, volcapa
@@ -10,10 +9,8 @@
  INTEGER,INTENT(IN) ::  ngrid, nslope, nsoil_PEM
  REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope)
- REAL,INTENT(IN) :: ps_new(ngrid)
- REAL,INTENT(IN) :: cellarea(ngrid)
+ REAL,INTENT(IN) :: p_avg_new
  REAL,INTENT(IN) :: co2ice(ngrid,nslope)
  REAL,INTENT(IN) :: waterice(ngrid,nslope)
  REAL,INTENT(in) :: ice_depth(ngrid,nslope)
- 
 
 ! Outputs:
@@ -44,22 +41,9 @@
  REAL :: regolith_inertia(ngrid,nslope) ! TI of the regolith 
  REAL :: d(ngrid,nsoil_PEM,nslope)
- REAL :: p_avg_new
  REAL :: Total_surface
  INTEGER :: ispermanent_co2glaciers(ngrid,nslope)
  INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)
 
-
 ! 0. Initialisation
-
-
-
-
-  Total_surface = sum(cellarea)
-  p_avg_new = 0.
-  do ig = 1,ngrid
-    p_avg_new = p_avg_new + ps_new(ig)*cellarea(ig)/Total_surface
-  enddo
-
-
 
   do ig = 1,ngrid
@@ -80,20 +64,9 @@
   enddo
 
-
  ispermanent_co2glaciers(:,:) = 0
  ispermanent_h2oglaciers(:,:) = 0
 
-
 !  1.Ice TI feedback
 
-!  do ig=1,ngrid
-!    do islope=1,nslope
-!      if ((ispermanent_co2glaciers(ig,islope).eq.1).or.(ispermanent_h2oglaciers(ig,islope).eq.1)) then
-!       do iloop = 1,n_1km
-!         TI_PEM(ig,iloop,islope) = surfaceice_inertia
-!       enddo
-!      endif
-!    enddo
-!   enddo
     do islope = 1,nslope
      call soil_TIfeedback_PEM(ngrid,nsoil_PEM,waterice(:,islope),   TI_PEM(:,:,islope))
@@ -101,6 +74,4 @@
 
 ! 2. Modification of the regolith thermal inertia.
-
-
 
   do ig=1,ngrid
@@ -115,9 +86,4 @@
    enddo
 enddo
-
-
-
-
-
 
 !  3. Build new TI for the PEM
@@ -142,8 +108,5 @@
          exit
         endif
-       
       enddo
-
-
 
       ! 4.2 Build the new ti
@@ -173,10 +136,4 @@
     enddo !ig
 
-
-
-
-
-
-
 !=======================================================================
       RETURN
