Index: /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 2848)
+++ /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 2849)
@@ -32,5 +32,5 @@
  REAL ::  rho_regolith = 2000.    ! density of the reoglith, Buhler & Piqueux 2021
  real :: alpha_clapeyron = -6143.7! eq. (2) in Murphy & Koop 2005
- real :: beta_clapeyron = 29.9074 ! eq. (2) in Murphy & Koop 2005
+ real :: beta_clapeyron = 28.9074 ! eq. (2) in Murphy & Koop 2005
  real :: mi = 2.84e-7             ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
  real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
@@ -88,6 +88,7 @@
   do islope = 1,nslope
     do iloop = 1,n_1km
+       K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
      if(TI_PEM(ig,iloop,islope).lt.inertie_thresold)  then
-       K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
+
        theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1+K*pvapor_avg(ig)))**mu
        m_h2o_adsorbed(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*mi*rho_regolith
Index: /trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90	(revision 2848)
+++ /trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90	(revision 2849)
@@ -1,3 +1,3 @@
-   SUBROUTINE computeice_table(timelen,ngrid,nslope,nsoil_GCM,nsoil_PEM,tsoil,tsurf,q_co2,q_h2o,ps,ice_table)
+   SUBROUTINE computeice_table(ngrid,nslope,nsoil_PEM,rhowatersurf_ave,rhowatersoil_ave,ice_table)
 #ifndef CPP_STD
     USE comsoil_h, only:  inertiedat, volcapa
@@ -6,127 +6,29 @@
     implicit none 
 
-    integer,intent(in) :: timelen,ngrid,nslope,nsoil_PEM,nsoil_GCM
-    real,intent(in) :: tsoil(ngrid,nsoil_PEM,nslope,timelen)    ! soil temperature, timeseries [K]
-    real,intent(in) :: tsurf(ngrid,nslope,timelen)              ! surface temperature [K]
-    real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
-    real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
-    real,intent(in) :: ps(ngrid,timelen)                        ! surface pressure [Pa]
-    real,intent(out) :: ice_table(ngrid,nslope)                 ! ice table [m]
+    integer,intent(in) :: ngrid,nslope,nsoil_PEM
+    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)    ! surface temperature, timeseries [K]
+    real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope)    ! soil water density, yearly average [kg/m^3]
+    real,intent(inout) :: ice_table(ngrid,nslope)                 ! ice table [m]
+    real :: z1,z2
+    integer ig, islope,isoil
+    real :: diff_rho(nsoil_PEM)                    ! difference of vapor content
 
-    real :: m_h2o = 18.01528E-3
-    real :: m_co2 = 44.01E-3  
-    real :: m_noco2 = 33.37E-3  
-    real :: A,B,z1,z2
-    real :: alpha = -6143.7 
-    real :: beta = 29.9074
-
-    integer ig, islope,isoil,it
-    real,allocatable :: mass_mean(:,:)                            ! mean mass above the surface
-    real,allocatable :: zplev_mean(:,:)                           ! pressure above the surface
-    real,allocatable :: pvapor(:,:)                               ! partial pressure above the surface
-
-    real,allocatable :: rhovapor(:,:,:)
-    real,allocatable :: rhovapor_avg(:,:)                          ! mean vapor_density at the surface yearly averaged
-
-    real :: psv_surf
-    real,allocatable :: rho_soil(:,:,:,:)            ! water vapor in the soil
-    real,allocatable :: rho_soil_avg(:,:,:)                ! water vapor in the soil yearly averaged
-
-    real,allocatable :: diff_rho(:,:,:)                    ! difference of vapor content
-
-     allocate(rhovapor(ngrid,nslope,timelen))
-     allocate(rhovapor_avg(ngrid,nslope))
-     allocate(pvapor(ngrid,timelen))
-
-     allocate(mass_mean(ngrid,timelen))
-     allocate(zplev_mean(ngrid,timelen))
-
-! 0. Some initializations
-
-      A =(1/m_co2 - 1/m_noco2)
-      B=1/m_noco2
-! 1. Compute rho surface yearly averaged
-
-!   1.1 Compute the partial pressure of vapor
-!a. the molecular mass into the column
-     do ig = 1,ngrid
-       mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
-     enddo
-
-! b. pressure level
-     do it = 1,timelen
-       do ig = 1,ngrid
-         zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
-       enddo
-     enddo
-
-! c. Vapor pressure
-     pvapor(:,:) = mass_mean(:,:)/m_h2o*q_h2o(:,:)*zplev_mean(:,:)
-    
-     deallocate(mass_mean)
-     deallocate(zplev_mean)
-  
-!    1.2   Check if there is frost at the surface and then compute the density
-!    at the surface
-     do ig = 1,ngrid
-       do islope = 1,nslope
-         do it = 1,timelen
-           psv_surf = exp(alpha/tsurf(ig,islope,it) +beta)
-          if ((isnan(psv_surf)).or.(isnan(pvapor(ig,it)))) then
-          write(*,*) 'pb vapor',ig,islope,it
-         stop
-         endif       
-         rhovapor(ig,islope,it) = min(psv_surf,pvapor(ig,it))/tsurf(ig,islope,it)
-         enddo
-       enddo
-     enddo
-     deallocate(pvapor)
-
-!    1.3 Density at the surface yearly averaged
-     rhovapor_avg(:,:) = SUM(rhovapor(:,:,:),3)/timelen
-
-     deallocate(rhovapor)
-
-! 2. Compute rho soil vapor
-   
-     allocate(rho_soil_avg(ngrid,nslope,nsoil_PEM))
-     allocate(rho_soil(ngrid,nslope,nsoil_PEM,timelen))
 
      do ig = 1,ngrid
-       do islope = 1,nslope
-         do isoil = 1,nsoil_PEM
-            do it = 1,timelen 
-              rho_soil(ig,islope,isoil,it) = exp(alpha/tsoil(ig,isoil,islope,it) +beta)/tsoil(ig,isoil,islope,it)        
-             enddo
-           enddo
-       enddo
-     enddo
-
-    rho_soil_avg(:,:,:) = SUM( rho_soil(:,:,:,:),4)/timelen
-    deallocate(rho_soil)
-
-! 3. Computing ice table
-   
-    ice_table (:,:) = -1e4
-
-         allocate(diff_rho(ngrid,nslope,nsoil_PEM))
+      do islope = 1,nslope
+           ice_table(ig,islope) = -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) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
+
          enddo
-    
-  deallocate(rhovapor_avg)
-  deallocate(rho_soil_avg)
 
-     do ig = 1,ngrid
-       do islope = 1,nslope
-         if(diff_rho(ig,islope,1) > 0) then
+         if(diff_rho(1) > 0) then
            ice_table(ig,islope) = 0.
          else
            do isoil = 1,nsoil_PEM -1
-             if((diff_rho(ig,islope,isoil).lt.0).and.(diff_rho(ig,islope,isoil+1).gt.0.)) then
-               z1 = (diff_rho(ig,islope,isoil) - diff_rho(ig,islope,isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
-               z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(ig,islope,isoil+1)
+             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,islope) = -z2/z1
                exit
@@ -136,6 +38,6 @@
         enddo
       enddo
+
     
-    deallocate(diff_rho)
 
 !=======================================================================
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 2848)
+++ /trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90	(revision 2849)
@@ -37,4 +37,6 @@
 
 !=======================================================================
+
+  STOPPING=.false.
 
   pos_tend=0.
Index: /trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 2848)
+++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 2849)
@@ -78,5 +78,5 @@
                            co2iceflow, beta_slope, capcal_slope,&
                            albedo_slope,emiss_slope,qsurf_slope,&
-                           iflat,ini_comslope_h
+                           iflat,major_slope,ini_comslope_h
       use time_phylmdz_mod, only: daysec,dtphys
       USE comconst_mod, ONLY: rad,g,r,cpp,pi
@@ -272,8 +272,15 @@
      REAL,ALLOCATABLE  :: alph_locslope(:,:)
      REAL,ALLOCATABLE  :: beta_locslope(:,:)   
-
+     REAL,ALLOCATABLE  :: watersurf_density_timeseries(:,:,:,:)
+     REAL,ALLOCATABLE  :: watersoil_density_timeseries(:,:,:,:,:)
+     REAL,ALLOCATABLE  :: watersurf_density_phys_timeseries(:,:,:)
+     REAL,ALLOCATABLE  :: watersurf_density_phys_ave(:,:)
+     REAL,ALLOCATABLE  :: watersoil_density_phys_PEM_timeseries(:,:,:,:)
+     REAL,ALLOCATABLE  :: watersoil_density_phys_PEM_ave(:,:,:)
      REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)
      REAL, ALLOCATABLE :: delta_co2_adsorbded(:)
      REAL :: totmass_adsorbded
+   real :: alpha_clap_h2o = -6143.7
+   real :: beta_clap_h2o = 28.9074
 
 #ifdef CPP_STD
@@ -289,5 +296,5 @@
 ! Loop variable
      LOGICAL :: bool_sublim
-     INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop
+     INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil
      REAL :: beta = 3182.48
      REAL :: alpha = 23.3494
@@ -406,5 +413,5 @@
               watercap,inertiesoil,nslope,tsurf_slope,           &
               tsoil_slope,co2ice_slope,def_slope,def_slope_mean, &
-              subslope_dist,albedo_slope,emiss_slope, TI_GCM_start,     &
+              subslope_dist,major_slope,albedo_slope,emiss_slope, TI_GCM_start,     &
               qsurf_slope,watercap_slope)
 
@@ -551,9 +558,15 @@
      allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
      allocate(delta_co2_adsorbded(ngrid))
-
+     allocate(watersurf_density_timeseries(iim+1,jjm+1,nslope,timelen))
+     allocate(watersoil_density_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
+     allocate(watersurf_density_phys_timeseries(ngrid,nslope,timelen))
+     allocate(watersurf_density_phys_ave(ngrid,nslope))
+     allocate(watersoil_density_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
+     allocate(watersoil_density_phys_PEM_ave(ngrid,nsoilmx_PEM,nslope))
      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)
+                       nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &      
+                       watersurf_density_timeseries,watersoil_density_timeseries)
 
 ! 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
@@ -569,5 +582,6 @@
 
      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)
+                  nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &      
+                  watersurf_density_timeseries,watersoil_density_timeseries)
 
      print *, "Downloading data Y2 done"
@@ -622,4 +636,5 @@
 
       allocate(ice_depth(ngrid,nslope))
+      ice_depth(:,:) = 0.
       allocate(TI_GCM_phys(ngrid,nsoilmx,nslope))
 
@@ -636,11 +651,13 @@
       deallocate(co2_ice_GCM_slope)
       deallocate(TI_GCM)
-
-  if(soil_pem) then
-
       deallocate(tsurf_GCM_timeseries)
+
+
+   if(soil_pem) then
       call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM)
-
       DO islope = 1,nslope
+        DO t=1,timelen
+          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersurf_density_timeseries(:,:,islope,t),watersurf_density_phys_timeseries(:,islope,t))
+        ENDDO
         DO l=1,nsoilmx
            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope))
@@ -648,12 +665,21 @@
            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))
+           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersoil_density_timeseries(:,:,l,islope,t),watersoil_density_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))
            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope))
+           DO t=1,timelen
+           watersoil_density_phys_PEM_timeseries(:,l,islope,t) =  watersoil_density_phys_PEM_timeseries(:,nsoilmx,islope,t) 
+           ENDDO
         ENDDO
       ENDDO
-
+       watersoil_density_phys_PEM_ave(:,:,:) = SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
+      watersurf_density_phys_ave(:,:) = SUM(watersurf_density_phys_timeseries(:,:,:),3)/timelen
+      deallocate(watersurf_density_timeseries)
+      deallocate(watersurf_density_phys_timeseries)
+      deallocate(watersoil_density_timeseries)
       deallocate(tsoil_ave_yr1)
       deallocate(tsoil_ave)
@@ -789,7 +815,8 @@
 !---------------------------- 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, &
+      call pemetat0(.true.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
       tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsurf_phys_GCM_timeseries,tsoil_phys_PEM_timeseries,&
-      tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), co2_adsorbded_phys,delta_co2_adsorbded)
+      tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM, co2_adsorbded_phys,delta_co2_adsorbded,&
+      watersurf_density_phys_ave,watersoil_density_phys_PEM_ave)
 
     if(soil_pem) then
@@ -848,4 +875,5 @@
      year_iter=0
      print *, "Max_iter_pem", Max_iter_pem
+     print *, 'year_iter_max',year_iter_max
      do while (year_iter.LT.year_iter_max)
 
@@ -856,8 +884,15 @@
            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
+     enddo
+       print *, 'Global average pressure old time step',global_ave_press_old
+       print *, 'Global average pressure new time step',global_ave_press_new
+
+     do i=1,ngrid
        if(soil_pem) then
          global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
        endif
      enddo
+       print *, 'Global average pressure old time step',global_ave_press_old
+       print *, 'Global average pressure new time step',global_ave_press_new
 
 ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
@@ -974,5 +1009,5 @@
               flag_co2flow_mesh(ig) = 1.
             ENDIF ! co2ice > hmax
-          ENDIF ! iflattsoil_lope
+          ENDIF ! iflat
         ENDDO !islope 
        ENDDO !ig
@@ -1032,8 +1067,8 @@
               emiss_slope(ig,islope) = emissiv         
             endif
-          elseif( co2ice_slope(ig,islope).GT. 1E-5) THEN !Put tsurf as tcond co2
+          elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice_phys_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
+              if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-3) then
                 ave = ave + beta/(alpha-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
               else
@@ -1056,4 +1091,5 @@
       enddo
 
+
     if(soil_pem) then
 
@@ -1071,39 +1107,43 @@
   do islope = 1,nslope
      TI_locslope(:,:) = TI_PEM(:,:,islope)
-     Tsoil_locslope(:,:) = tsoil_PEM(:,:,islope)
-     Tsurf_locslope(:) = tsurf_ave_phys(:,islope)
      alph_locslope(:,:) = alph_PEM(:,:,islope)
      beta_locslope(:,:) = beta_PEM(:,:,islope)
-     call soil_pem_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))
-     enddo
-     TI_PEM(:,:,islope) = TI_locslope(:,:)
-     tsoil_PEM(:,:,islope) = Tsoil_locslope(:,:)
-     tsurf_ave_phys(:,islope)  =  Tsurf_locslope(:)
-     alph_PEM(:,:,islope) = alph_locslope(:,:)
-     beta_PEM(:,:,islope) = beta_locslope(:,:)
-  enddo
-
-  print *, "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(*,*) "Problem : There is nan tsoil", ig, iloop, islope, tsoil_PEM(ig,iloop,islope)
-           stop 
-           endif
+       Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
+       Tsurf_locslope(:) =  tsurf_phys_GCM_timeseries(:,islope,t)
+
+       call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
+       call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
+
+
+       tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
+       do ig = 1,ngrid
+         do isoil = 1,nsoilmx_PEM
+          watersoil_density_phys_PEM_timeseries(ig,l,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
+          if(isnan(Tsoil_locslope(ig,isoil))) then
+             write(*,*)'Tsoil=',ig,isoil,Tsoil_locslope(ig,isoil)
+          endif
          enddo
        enddo
      enddo
+     alph_PEM(:,:,islope) = alph_locslope(:,:) 
+     beta_PEM(:,:,islope) = beta_locslope(:,:) 
+  enddo
+
+ 
+     tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
+     watersoil_density_phys_PEM_ave(:,:,:)= SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
+  print *, "Update of soil temperature done"
+
+  deallocate(TI_locslope)
+  deallocate(Tsoil_locslope)
+  deallocate(Tsurf_locslope)
+  deallocate(alph_locslope)
+  deallocate(beta_locslope)
 
   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)
+     call computeice_table(ngrid,nslope,nsoilmx_PEM,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth)
        
       print *, "Update soil propreties"
@@ -1153,4 +1193,7 @@
       if (STOPPING_water) then
         print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water 
+      endif
+      if (year_iter.ge.year_iter_max) then
+        print *, "STOPPING because maximum number of iterations reached"
       endif
       if (STOPPING_1_water) then
@@ -1210,5 +1253,5 @@
       ENDDO ! of DO ig=1,ngrid
 
-! III_a.2 Tsurf and Tsoil update (for startfi)
+! III_a.2  Tsoil update (for startfi)
 
    if(soil_pem) then
@@ -1219,11 +1262,4 @@
    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
-      ENDDO ! of DO ig=1,ngrid
 
 #ifndef CPP_STD
@@ -1260,4 +1296,15 @@
       ENDDO
     ENDDO
+
+
+      DO ig = 1,ngrid
+        tsurf(ig) = 0.
+        DO islope = 1,nslope
+          tsurf(ig) = tsurf(ig) + (emiss_slope(ig,islope)*tsurf_slope(ig,islope))**4 &
+                     * subslope_dist(ig,islope)      
+        ENDDO
+        tsurf(ig) = tsurf(ig)**(0.25)/emis(ig)
+      ENDDO ! of DO ig=1,ngrid
+
 #endif
 
Index: /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2848)
+++ /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2849)
@@ -1,9 +1,9 @@
    SUBROUTINE pemetat0(startpem_file,filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
-                      tsurf_ave_yr1,tsurf_ave,q_co2,q_h2o,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, 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, &
+                      global_ave_pressure, m_co2_regolith_phys,deltam_co2_regolith_phys, watersurf_ave,watersoil_ave)
    
    use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
    use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
    use comsoil_h, only:  volcapa,inertiedat 
-   use ini_soil_mod, only:  ini_icetable
    use soil_evolution_mod, only: soil_pem,soil_pem_CN
    use adsorption_mod, only : regolith_co2adsorption
@@ -32,5 +32,7 @@
   real,intent(in) :: waterice(ngrid,nslope) 
   real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope) 
-
+  real, intent(in) :: watersurf_ave(ngrid,nslope)     ! surface water ice density, yearly averaged  (kg/m^3)
+  real, intent(inout) :: watersoil_ave(ngrid,nsoil_PEM,nslope)    ! surface water ice density, yearly averaged (kg/m^3)
+  real, intent(in) :: global_ave_pressure                       ! mean average pressure on the planet
 ! outputs
 
@@ -59,4 +61,6 @@
    real :: co2_ads_prev(ngrid)
    real :: year_PEM_read
+   real :: alpha_clap_h2o = -6143.7
+   real :: beta_clap_h2o = 28.9074
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -221,4 +225,10 @@
         enddo
      enddo
+      do isoil = nsoil_GCM+1,nsoil_PEM
+        do ig = 1,ngrid
+        watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope)
+        enddo
+      enddo
+
     
 ENDDO
@@ -230,16 +240,9 @@
 !3. Ice Table
 
-   call get_field("ice_table",ice_table,found)
-   if(.not.found) then
-      write(*,*)'PEM settings: failed loading <Ice Table>'
-      write(*,*)'will reconstruct the values of ice table'
-
-      call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:),  tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
-
-   else 
-   ! update ice table
-     call computeice_table(timelen,ngrid,nslope,nsoil_PEM,tsoil_inst,tsurf_inst,q_co2,q_h2o, ps_inst, ice_table)
-
-   endif
+
+
+     call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table)
+     call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM)
+
 
 print *,'PEMETAT0: ICE TABLE  DONE'
@@ -259,17 +262,11 @@
    deltam_co2_regolith_phys(:) = 0.
    exit
+   else
+   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
+
    endif
 ENDDO
 
-if (found) then
-   DO iloop = 1,nsoil_GCM
-      tsoil_tmp_yr1(:,iloop,:) = tsoil_PEM_yr1(:,iloop,:)
-
-   ENDDO  
-
-   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
-
-
-endif
+
 print *,'PEMETAT0: CO2 adsorption done  '
 
@@ -364,6 +361,4 @@
 !b) Soil temperature 
     do islope = 1,nslope
-
-     write(*,*) "islope=",islope
      do ig = 1,ngrid
         kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
@@ -381,5 +376,4 @@
             tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
       enddo
-!     write(*,*) "ig, islope, T=", ig,islope,tsoil_PEM(ig,:,islope)
      enddo
     
@@ -389,12 +383,17 @@
         enddo
      enddo
+      do isoil = nsoil_GCM+1,nsoil_PEM
+        do ig = 1,ngrid
+        watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope)
+        enddo
+      enddo
 enddo
 print *,'PEMETAT0: TSOIL DONE  '
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !c) Ice table    
-      do islope = 1,nslope
-      call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:), &  
-       tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope))
-      enddo
+     call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table)
+     call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM)
+
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
Index: /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2848)
+++ /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2849)
@@ -3,5 +3,6 @@
 !
 SUBROUTINE read_data_GCM(fichnom,min_h2o_ice_s,min_co2_ice_s,iim_input,jjm_input,nlayer,vmr_co2_gcm,ps_GCM,timelen, & 
-             min_co2_ice_slope,min_h2o_ice_slope,nslope,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2_GCM,q_h2o_GCM,co2_ice_slope)
+             min_co2_ice_slope,min_h2o_ice_slope,nslope,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2_GCM,q_h2o_GCM,co2_ice_slope, &
+             watersurf_density,watersoil_density)
 
       use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
@@ -39,5 +40,4 @@
   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) ::  ps_GCM(iim_input+1,jjm_input+1,timelen)
 
@@ -48,4 +48,6 @@
   REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen) ! Surface temperature of the concatenated file
   REAL , INTENT(OUT) ::  tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file
+  REAL , INTENT(OUT) ::  watersurf_density(iim_input+1,jjm_input+1,nslope,timelen) ! Soil temperature of the concatenated file
+  REAL , INTENT(OUT) ::  watersoil_density(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file
 
   REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Thermal Inertia  of the concatenated file
@@ -79,5 +81,4 @@
 
       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))
@@ -161,4 +162,22 @@
 
      print *, "Downloading data for inertiesoil_slope done"
+
+     print *, "Downloading data for watersoil_density ..."
+
+DO islope=1,nslope
+  write(num,fmt='(i2.2)') islope
+  call get_var4("Waterdensity_soil_slope"//num,watersoil_density(:,:,:,islope,:))
+ENDDO
+
+     print *, "Downloading data for  watersoil_density  done"
+
+     print *, "Downloading data for  watersurf_density  ..."
+
+DO islope=1,nslope
+  write(num,fmt='(i2.2)') islope
+  call get_var3("Waterdensity_surface"//num,watersurf_density(:,:,islope,:))
+ENDDO
+
+     print *, "Downloading data for  watersurf_density  done"
 
   endif !soil_pem
@@ -234,4 +253,10 @@
   ENDDO
 
+
+      deallocate(co2_ice_s)
+      deallocate(h2o_ice_s_slope)
+      deallocate(h2o_ice_s)
+
+
   CONTAINS
 
Index: /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 2848)
+++ /trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90	(revision 2849)
@@ -8,10 +8,13 @@
 
       USE temps_mod_evol, ONLY: year_bp_ini, year_PEM
-#ifndef CPP_STD
+#ifndef CPP_STD      
+      USE comconst_mod, ONLY: pi
       USE planete_h, ONLY: e_elips, obliquit, timeperi
 #else
       use planete_mod, only: e_elips, obliquit, timeperi
+      USE comcstfi_mod, only: pi
+
 #endif
-      USE comcstfi_mod, only: pi
+
       USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, &
    				end_lask_param_mod,last_ilask
Index: /trunk/LMDZ.COMMON/libf/evolution/update_soil.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/update_soil.F90	(revision 2848)
+++ /trunk/LMDZ.COMMON/libf/evolution/update_soil.F90	(revision 2849)
@@ -26,8 +26,6 @@
  REAL ::  L =  51058.
  REAL ::  inertie_thresold = 800. ! look for ice
- REAL ::  inertie_co2glaciers = 2120 ! Mellon et al. 2000
  REAL ::  inertie_averaged = 250 ! Mellon et al. 2000
- REAL ::  ice_inertia = 2120  ! Inertia of ice
- REAL ::  surfaceice_inertia = 800  ! Inertia of ice
+ REAL ::  ice_inertia = 1200  ! Inertia of ice
  REAL ::  P610 = 610.0
  REAL ::  m_h2o = 18.01528E-3
@@ -47,29 +45,27 @@
 ! 0. Initialisation
 
-  do ig = 1,ngrid
-    do islope = 1,nslope
-     if((abs(tend_h2oglaciers(ig,islope)).lt.1e-5).and.(abs(waterice(ig,islope)).gt.0)) then
-        ispermanent_h2oglaciers(ig,islope) = 1
-     else
-        ispermanent_h2oglaciers(ig,islope) = 0
-     endif 
+!  do ig = 1,ngrid
+!    do islope = 1,nslope
+!     if((abs(tend_h2oglaciers(ig,islope)).lt.1e-5).and.(abs(waterice(ig,islope)).gt.(1.e-4))) then
+!        ispermanent_h2oglaciers(ig,islope) = 1
+!     else
+!        ispermanent_h2oglaciers(ig,islope) = 0
+!     endif 
+!
+!     if((abs(tend_co2glaciers(ig,islope)).lt.1e-5).and.(abs(co2ice(ig,islope)).gt.(1.e-4))) then
+!        ispermanent_co2glaciers(ig,islope) = 1
+!     else
+!        ispermanent_co2glaciers(ig,islope) = 0
+!     endif
+!    enddo
+!  enddo
 
-     if((abs(tend_co2glaciers(ig,islope)).lt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then
-        ispermanent_co2glaciers(ig,islope) = 1
-     else
-        ispermanent_co2glaciers(ig,islope) = 0
-     endif
-    enddo
 
-  enddo
-
- ispermanent_co2glaciers(:,:) = 0
- ispermanent_h2oglaciers(:,:) = 0
 
 !  1.Ice TI feedback
 
-    do islope = 1,nslope
-     call soil_TIfeedback_PEM(ngrid,nsoil_PEM,waterice(:,islope),   TI_PEM(:,:,islope))
-    enddo
+!    do islope = 1,nslope
+!     call soil_TIfeedback_PEM(ngrid,nsoil_PEM,waterice(:,islope),   TI_PEM(:,:,islope))
+!    enddo
 
 ! 2. Modification of the regolith thermal inertia.
