Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 4170)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 4174)
@@ -33,5 +33,5 @@
 use display,            only: print_ini, print_end, print_msg, LVL_NFO, LVL_WRN
 use evolution,          only: n_yr_run, n_yr_sim, ntot_yr_sim, nmax_yr_run, dt, idt, r_plnt2earth_yr, pem_ini_date
-use geometry,           only: ngrid, nslope, nday, nsoil_PCM, nsoil, cell_area, total_surface
+use geometry,           only: ngrid, nslope, nsoil_PCM, nsoil, cell_area, total_surface
 use glaciers,           only: h2oice_flow, co2ice_flow, flow_co2glaciers, flow_h2oglaciers
 use ice_table,          only: icetable_equilibrium, icetable_dynamic, evolve_ice_table
@@ -48,7 +48,7 @@
 use soil_therm_inertia, only: update_soil_TI
 use sorption,           only: do_sorption, compute_totmass_adsorbed, evolve_regolith_adsorption
-use stopping_crit,      only: stopFlags, stopping_crit_pressure, stopping_crit_h2o, stopping_crit_h2oice, stopping_crit_co2ice
+use stopping_crit,      only: stopFlags, stopping_crit_pressure, stopping_crit_h2oice, stopping_crit_co2ice
 use surface,            only: emissivity_PCM
-use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs
+use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2o_fluxes
 use surf_temp,          only: tsurf_PCM, adapt_tsurf2disappearedice
 use tendencies,         only: compute_tendice, evolve_tend_co2ice, evolve_flux_ssice
@@ -80,6 +80,4 @@
 real(dp), dimension(:,:), allocatable :: zshift_surf   ! Elevation shift for the surface [m]
 real(dp), dimension(:,:), allocatable :: zlag          ! Newly built lag thickness [m]
-! Tendency-related:
-real(dp), dimension(:,:), allocatable :: d_h2oice_new ! Adjusted tendency of perennial H2O ice to keep balance between donor and recipient [kg/m2/y]
 ! Layering-related:
 type(ptrarray), dimension(:,:), allocatable :: current          ! Current active stratum in the layering
@@ -95,5 +93,4 @@
 real(qp) :: totmass_ini                                                ! Initial total CO2 mass [kg]
 real(qp) :: totmass_adsh2o                                             ! Current total adsorbed H2O mass [kg]
-real(qp) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm   ! Variables to balance H2O ice reservoirs
 
 ! CODE
@@ -235,5 +232,4 @@
     allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
     if (do_layering) then
-        allocate(d_h2oice_new(ngrid,nslope))
         h2oice_depth_old(:,:) = h2oice_depth(:,:)
         ! The sublimating tendency coming from subsurface ice is given through the surface ice tendency
@@ -241,5 +237,5 @@
         do islope = 1,nslope
             do i = 1,ngrid
-                call evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice_new(i,islope),new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p)
+                call evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice(i,islope),new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p)
                 call print_layering(layerings_map(i,islope))
             end do
@@ -248,7 +244,10 @@
         call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
         ! Balance H2O ice reservoirs
-        call stopping_crit_h2o(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopcrit)
-        call balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
-        deallocate(d_h2oice_new)
+        call balance_h2o_fluxes(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,stopcrit)
+        if (stopcrit%stop_code() > 0) then
+            deallocate(zshift_surf,zlag)
+            call print_msg(stopcrit%stop_message(),LVL_NFO)
+            exit
+        end if
     else
         zlag(:,:) = 0._dp
@@ -376,5 +375,5 @@
     if (backup_rate > 0) then
         if (mod(idt,backup_rate) == 0) call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, &
-                                                            icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map,idt)
+                                                            icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,h2o_ads_reg,co2_ads_reg,layerings_map,idt)
     end if
 
@@ -408,5 +407,5 @@
 
 call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, &
-                     icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map)
+                     icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,h2o_ads_reg,co2_ads_reg,layerings_map)
 
 ! Update the duration information of the workflow