Changeset 4180 for trunk/LMDZ.COMMON/libf/evolution

Timestamp:

Apr 10, 2026, 7:17:55 PM (4 hours ago)

Author:

jbclement

Message:

PEM:

• Rework layering-related logic, especially clarify interactions between surface and subsurface water tendencies and disable CO2 lag resistance (inconsistent without updating pressure and mass balance + PCM).
• Prevent simultaneous activation of layering and ice flows.
• Add warning when flux_geo /= 0 while soil is disabled.
• Add new utility function "is_lvl_enabled" for displaying.
• Replace deprecated 'minieps' with 'eps'/'tol'.

JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• config.F90 (modified) (4 diffs)
• display.F90 (modified) (3 diffs)
• glaciers.F90 (modified) (2 diffs)
• io_netcdf.F90 (modified) (6 diffs)
• layered_deposits.F90 (modified) (19 diffs)
• maths.F90 (modified) (1 diff)
• numerics.F90 (modified) (1 diff)
• orbit.F90 (modified) (2 diffs)
• output.F90 (modified) (2 diffs)
• pem.F90 (modified) (7 diffs)
• planet.F90 (modified) (6 diffs)
• soil_therm_inertia.F90 (modified) (3 diffs)
• sorption.F90 (modified) (4 diffs)
• stopping_crit.F90 (modified) (1 diff)
• surf_ice.F90 (modified) (4 diffs)
• surf_temp.F90 (modified) (2 diffs)
• tendencies.F90 (modified) (5 diffs)
• utility.F90 (modified) (2 diffs)

trunk/LMDZ.COMMON/libf/evolution/changelog.txt

@@ -963 +963 @@
     - adds new stopping criterion when H2O flux balance fails8 (sanity check);
     - introduces configurable weighting strategies for H2O flux balancing.
+
+== 10/04/2026 == JBC
+- Rework layering-related logic, especially clarify interactions between surface and subsurface water tendencies and disable CO2 lag resistance (inconsistent without updating pressure and mass balance + PCM).
+- Prevent simultaneous activation of layering and ice flows.
+- Add warning when flux_geo /= 0 while soil is disabled.
+- Add new utility function "is_lvl_enabled" for displaying.
+- Replace deprecated 'minieps' with 'eps'/'tol'.

trunk/LMDZ.COMMON/libf/evolution/config.F90

@@ -17 +17 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, di, k4, minieps
+use numerics, only: dp, di, k4, eps
 
 ! DECLARATION
@@ -240 +240 @@
 ! DEPENDENCIES
 ! ------------
-use stoppage,  only: stop_clean
-use soil,      only: do_soil, reg_thprop_dependp, flux_geo
-use sorption,  only: do_sorption
-use orbit,     only: evo_orbit
-use evolution, only: pem_ini_date
-use ice_table, only: icetable_equilibrium, icetable_dynamic
-use display,   only: print_msg, LVL_WRN
+use stoppage,         only: stop_clean
+use soil,             only: do_soil, reg_thprop_dependp, flux_geo
+use sorption,         only: do_sorption
+use orbit,            only: evo_orbit
+use evolution,        only: pem_ini_date
+use ice_table,        only: icetable_equilibrium, icetable_dynamic
+use glaciers,         only: h2oice_flow, co2ice_flow
+use layered_deposits, only: do_layering
+use display,          only: print_msg, LVL_WRN
 
 ! DECLARATION
@@ -255 +257 @@
 ! ----
 ! Warnings (possible incompatibilities)
-if (evo_orbit .and. abs(pem_ini_date) < minieps) call print_msg('''evo_orbit = .true.'' but the initial date of the PEM is set to 0!',LVL_WRN)
+if (evo_orbit .and. abs(pem_ini_date) < eps) call print_msg('''evo_orbit = .true.'' but the initial date of the PEM is set to 0!',LVL_WRN)
+if (abs(flux_geo) > eps .and. .not. do_soil) call print_msg('soil is not activated but flux_geo /= 0!',LVL_WRN)
 
 ! Errors (true incompatibilities)
 if (.not. do_soil) then
     if (icetable_equilibrium .or. icetable_dynamic) call stop_clean(__FILE__,__LINE__,'ice table must be used when do_soil = true!',1)
-    if (abs(flux_geo) > minieps) call stop_clean(__FILE__,__LINE__,'soil is not activated but flux_geo /= 0!',1)
     if (reg_thprop_dependp) call stop_clean(__FILE__,__LINE__,'regolith properties vary according to Ps only when soil = true!',1)
     if (do_sorption) call stop_clean(__FILE__,__LINE__,'do_soil must be true when do_sorption = true!',1)
 end if
+if (do_layering .and. h2oice_flow) call stop_clean(__FILE__,__LINE__,'layering and H2O ice flow cannot be activated at the same time!',1)
+if (do_layering .and. co2ice_flow) call stop_clean(__FILE__,__LINE__,'layering and CO2 ice flow cannot be activated at the same time!',1)
 
 END SUBROUTINE check_config_incompatibility
@@ -458 +462 @@
 call print_msg('> Initializing soil parameters',LVL_NFO)
 volcapa = controle(35)
-if (abs(volcapa) < minieps) call stop_clean(__FILE__,__LINE__,'volcapa is 0 in "'//startfi_name//'"!',1)
+if (abs(volcapa) < eps) call stop_clean(__FILE__,__LINE__,'volcapa is 0 in "'//startfi_name//'"!',1)
 
 ! Initialize orbital data

trunk/LMDZ.COMMON/libf/evolution/display.F90

@@ -19 +19 @@
 use, intrinsic :: iso_fortran_env, only: stderr => error_unit
 use, intrinsic :: iso_fortran_env, only: stdin => input_unit
-use numerics,                      only: dp, di, li, k4, minieps
+use numerics,                      only: dp, di, li, k4, eps
 
 ! DECLARATION
@@ -277 +277 @@
 
 END SUBROUTINE print_msg
+!=======================================================================
+
+!=======================================================================
+FUNCTION is_lvl_enabled(lvl) RESULT(enabled)
+!-----------------------------------------------------------------------
+! NAME
+!     is_lvl_enabled
+!
+! DESCRIPTION
+!     Return true if a message at level 'lvl' should be displayed.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+integer(di), intent(in) :: lvl
+
+! LOCAL VARIABLES
+! ---------------
+logical(k4) :: enabled
+
+! CODE
+! ----
+enabled = .false.
+if (lvl <= verbosity_lvl) enabled = .true.
+
+END FUNCTION is_lvl_enabled
 !=======================================================================
 
@@ -465 +501 @@
 if (gottacatch_emall_()) then
     call print_msg('',LVL_NFO)
-    if (abs(n_yr_run - first_gen) < minieps) then
+    if (abs(n_yr_run - first_gen) < eps) then
         do i = 1,size(surprise)
             call print_msg(trim(surprise(i)),LVL_NFO)

trunk/LMDZ.COMMON/libf/evolution/glaciers.F90

@@ -18 +18 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, di, k4, minieps
+use numerics, only: dp, di, k4, eps
 
 ! DECLARATION
@@ -301 +301 @@
                     iaval = islope + 1
                 end if
-                do while (iaval /= iflat .and. abs(subslope_dist(ig,iaval)) < minieps)
+                do while (iaval /= iflat .and. abs(subslope_dist(ig,iaval)) < eps)
                     if (iaval > iflat) then
                         iaval = iaval - 1

trunk/LMDZ.COMMON/libf/evolution/io_netcdf.F90

@@ -16 +16 @@
 ! DEPEDENCIES
 ! -----------
-use numerics, only: dp, di, k4, minieps
+use numerics, only: dp, di, k4, eps
 use netcdf,   only: nf90_double, nf90_noerr, nf90_strerror, nf90_write, nf90_nowrite,                &
                     nf90_open, nf90_close, nf90_redef, nf90_enddef, nf90_inquire, nf90_max_var_dims, &
@@ -989 +989 @@
 if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
 if (has_fill) then
-    if (abs(var - fill_value) < minieps) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
+    if (abs(var - fill_value) < eps) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
 end if
 
@@ -1074 +1074 @@
 if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
 if (has_fill) then
-    if (any(abs(var - fill_value) < minieps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
+    if (any(abs(var - fill_value) < eps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
 end if
 
@@ -1159 +1159 @@
 if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
 if (has_fill) then
-    if (any(abs(var - fill_value) < minieps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
+    if (any(abs(var - fill_value) < eps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
 end if
 
@@ -1244 +1244 @@
 if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
 if (has_fill) then
-    if (any(abs(var - fill_value) < minieps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
+    if (any(abs(var - fill_value) < eps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
 end if
 
@@ -1329 +1329 @@
 if (.not. has_fill) call check_nc(nf90_get_att(ncid,varid,"missing_value",fill_value),'getting missing value',has_fill)
 if (has_fill) then
-    if (any(abs(var - fill_value) < minieps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
+    if (any(abs(var - fill_value) < eps)) call stop_clean(__FILE__,__LINE__,'Fill values ('//real2str(fill_value)//') detected in '//var_name//'!',1)
 end if
 

trunk/LMDZ.COMMON/libf/evolution/layered_deposits.F90

@@ -17 +17 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, di, k4, tol
+use numerics, only: dp, di, k4, tol, eps
 use geometry, only: ngrid, nslope
 use surf_ice, only: rho_co2ice, rho_h2oice
@@ -505 +505 @@
 end do
 this%nb_str = 0
+nullify(this%top)
+nullify(this%bottom)
 
 END SUBROUTINE del_layering
@@ -699 +701 @@
     do ig = 1,ngrid
         do k = 1,size(layerings_array,3)
-            call add_stratum(layerings_map(ig,islope),layerings_array(ig,islope,k,1),layerings_array(ig,islope,k,2),layerings_array(ig,islope,k,3),layerings_array(ig,islope,k,4),layerings_array(ig,islope,k,5),layerings_array(ig,islope,k,6))
+            call add_stratum(layerings_map(ig,islope),layerings_array(ig,islope,k,1),       &
+                             layerings_array(ig,islope,k,2),layerings_array(ig,islope,k,3), &
+                             layerings_array(ig,islope,k,4),layerings_array(ig,islope,k,5), &
+                             layerings_array(ig,islope,k,6))
         end do
     end do
@@ -742 +747 @@
 ! CODE
 ! ----
+h2o_ice(:,:) = 0._dp
+co2_ice(:,:) = 0._dp
+h2oice_depth(:,:) = 0._dp
 do ig = 1,ngrid
     do islope = 1,nslope
@@ -1070 +1078 @@
 ! Is the considered stratum a dust lag?
 if (.not. is_dust_lag(current)) return
-do while (is_dust_lag(current))
+do
     h_toplag = h_toplag + thickness(current)
-    if (associated(current%down)) then
-        current => current%down
-    else
-        exit
+    if (.not. associated(current%down)) then ! Bottom of the layering and only dust lag strata
+        h_toplag = 0._dp
+        return
     end if
+    current => current%down
+    if (.not. is_dust_lag(current)) exit
 end do
-
-if (.not. associated(current%down)) then ! Bottom of the layering -> no ice below
-    h_toplag = 0._dp
-    return
-end if
 
 ! Is the stratum below the dust lag made of ice?
@@ -1091 +1095 @@
     end if
 else if (current%h_co2ice > 0._dp .and. current%h_co2ice > current%h_h2oice) then ! No, there is more CO2 ice than H2O ice
+    if (h_toplag < h_patchy_ice) co2_ice = current%h_co2ice*rho_co2ice ! The dust lag is too thin to considered ice as subsurface ice
     h_toplag = 0._dp ! Because it matters only for H2O ice
-    if (h_toplag < h_patchy_ice) co2_ice = current%h_co2ice*rho_co2ice ! But the dust lag is too thin to considered ice as subsurface ice
 end if
 
@@ -1128 +1132 @@
 ! Update of properties in the eroding dust lag
 str%top_elevation = str%top_elevation - h2lift*(1._dp + str%h_pore/str%h_dust)
-str%h_pore        = str%h_pore        - h2lift*str%h_pore/str%h_dust
-str%h_dust        = str%h_dust        - h2lift
+str%h_pore        = max(0._dp,str%h_pore - h2lift*str%h_pore/str%h_dust)
+str%h_dust        = max(0._dp,str%h_dust - h2lift)
 
 ! Remove the eroding dust lag if there is no dust anymore
@@ -1188 +1192 @@
 str%h_dust        = str%h_dust        - hlag_dust
 str%h_pore        = str%h_pore        - h2subl*h_pore/h_ice
+str%h_co2ice      = max(0._dp,str%h_co2ice)
+str%h_h2oice      = max(0._dp,str%h_h2oice)
+str%h_dust        = max(0._dp,str%h_dust)
+str%h_pore        = max(0._dp,str%h_pore)
 
 ! Correct the value of top_elevation for the upper strata
@@ -1224 +1232 @@
 
 !=======================================================================
-SUBROUTINE evolve_layering(this,d_co2ice,d_h2oice,new_str,zshift_surf,new_lag,zlag,current)
+SUBROUTINE get_co2_subl_resistance(current,R_subl)
+!-----------------------------------------------------------------------
+! NAME
+!     get_co2_subl_resistance
+!
+! DESCRIPTION
+!     Compute CO2 sublimation inhibition factor due to upper lag strata.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+type(stratum), pointer, intent(in)  :: current
+real(dp),               intent(out) :: R_subl
+
+! LOCAL VARIABLES
+! ---------------
+real(dp)               :: h_toplag
+type(stratum), pointer :: currentloc
+
+! CODE
+! ----
+R_subl = 1._dp
+if (.not. associated(current%up)) return ! If there is no upper stratum
+
+h_toplag = 0._dp
+currentloc => current%up
+do while (is_dust_lag(currentloc) .and. associated(currentloc%up))
+    h_toplag = h_toplag + thickness(currentloc)
+    currentloc => currentloc%up
+end do
+
+if (currentloc%h_h2oice >= h_patchy_ice) then
+    ! There is a thick H2O ice layer above the dust lag, so CO2 sublimation is fully inhibited
+    R_subl = 0._dp
+else
+    h_toplag = h_toplag + thickness(currentloc)
+    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
+end if
+
+END SUBROUTINE get_co2_subl_resistance
+!=======================================================================
+
+!=======================================================================
+SUBROUTINE settle_h2o_surf_subsurf_tend(h2oice_depth,flux_ssice,d_h2oice,icell,jcell)
+!-----------------------------------------------------------------------
+! NAME
+!     settle_h2o_surf_subsurf_tend
+!
+! DESCRIPTION
+!     Reconcile H2O surface and subsurface tendencies for one grid cell.
+!
+! AUTHORS & DATE
+!     JB Clement, 04/2026
+!
+! NOTES
+!
+!-----------------------------------------------------------------------
+
+! DEPENDENCIES
+! ------------
+use display, only: print_msg, LVL_ERR, LVL_WRN, LVL_DBG, is_lvl_enabled
+
+! DECLARATION
+! -----------
+implicit none
+
+! ARGUMENTS
+! ---------
+real(dp),               intent(in)           :: h2oice_depth
+real(dp),               intent(inout)        :: flux_ssice, d_h2oice
+integer(di),            intent(in), optional :: icell, jcell
+
+! LOCAL VARIABLES
+! ---------------
+logical(k4)   :: has_ssice, ss_sublim, ss_cond, surf_sublim, surf_cond
+character(64) :: cell_tag
+
+! CODE
+! ----
+cell_tag = ''
+if (present(icell) .and. present(jcell)) then
+    if (is_lvl_enabled(LVL_DBG)) write(cell_tag,'(a,i0,a,i0,a)') ' [i=',icell,',islope=',jcell,']'
+end if
+
+has_ssice = h2oice_depth > 0._dp
+ss_sublim = flux_ssice < -eps
+ss_cond = flux_ssice > eps
+surf_sublim = d_h2oice < -eps
+surf_cond = d_h2oice > eps
+
+if (has_ssice) then
+    if (ss_sublim) then
+        if (surf_sublim) then
+            call print_msg('Sublimating subsurface ice but sublimating surface ice is also present!'//trim(cell_tag),LVL_ERR)
+        else if (surf_cond) then
+            call print_msg('Sublimating subsurface ice but condensing surface ice is also present!'//trim(cell_tag),LVL_WRN)
+            call print_msg('This may happen because yearly minimum of frost is growing!',LVL_WRN)
+            call print_msg('Choice: condensing surface ice overrules sublimating subsurface ice!',LVL_WRN)
+            flux_ssice = 0._dp
+        else
+            ! Needed so the layering algorithm receives a retreat tendency
+            d_h2oice = flux_ssice
+        end if
+    else if (ss_cond) then
+        if (surf_cond) then
+            call print_msg('Condensing subsurface ice but condensing surface ice is also present!'//trim(cell_tag),LVL_WRN)
+            call print_msg('This may happen because yearly minimum of frost is growing!',LVL_WRN)
+            call print_msg('Choice: condensing surface ice overrules condensing subsurface ice!',LVL_WRN)
+            flux_ssice = 0._dp
+        else if (surf_sublim) then
+            call print_msg('Condensing subsurface ice but sublimating surface ice is also present!'//trim(cell_tag),LVL_ERR)
+        else
+            call print_msg('Condensing subsurface ice cannot be handled by the current layering algorithm!'//trim(cell_tag),LVL_ERR)
+        end if
+    else
+        flux_ssice = 0._dp ! Very small subsurface flux is set to zero to avoid numerical issues
+    end if
+else
+    if (abs(flux_ssice) > eps) call print_msg('Subsurface ice flux is non-zero while there is no subsurface ice!'//trim(cell_tag),LVL_ERR)
+end if
+
+END SUBROUTINE settle_h2o_surf_subsurf_tend
+!=======================================================================
+
+!=======================================================================
+SUBROUTINE evolve_layering(this,d_co2ice,d_h2oice,h2oice_depth,flux_ssice,new_str,zshift_surf,new_lag,zlag,current,icell,jcell)
 !-----------------------------------------------------------------------
 ! NAME
@@ -1239 +1382 @@
 !
 ! NOTES
-!     Creates dust lag layers when ice sublimation occurs.
+!
 !-----------------------------------------------------------------------
 
@@ -1255 +1398 @@
 ! ARGUMENTS
 ! ---------
-real(dp),               intent(in)    :: d_co2ice, d_h2oice
-type(layering),         intent(inout) :: this
-type(stratum), pointer, intent(inout) :: current
-logical(k4),            intent(inout) :: new_str, new_lag
-real(dp),               intent(out)   :: zshift_surf, zlag
-
-! LOCAL VARIABLES
-! ---------------
-real(dp)               :: dh_co2ice, dh_h2oice, dh_dust, dh_dust_co2, dh_dust_h2o
-real(dp)               :: h_co2ice_old, h_h2oice_old
-real(dp)               :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
-real(dp)               :: h2lift, h2lift_tot, h_pore, h_tot, h_toplag, R_subl
-logical(k4)            :: unexpected
-type(stratum), pointer :: currentloc
-
-! CODE
-! ----
+real(dp),               intent(inout)        :: d_co2ice, d_h2oice, flux_ssice
+real(dp),               intent(in)           :: h2oice_depth
+type(layering),         intent(inout)        :: this
+type(stratum), pointer, intent(inout)        :: current
+logical(k4),            intent(inout)        :: new_str, new_lag
+real(dp),               intent(out)          :: zshift_surf, zlag
+integer(di),            intent(in), optional :: icell, jcell
+
+! LOCAL VARIABLES
+! ---------------
+real(dp)    :: dh_co2ice, dh_h2oice, dh_dust, dh_dust_co2, dh_dust_h2o
+real(dp)    :: h_co2ice_old, h_h2oice_old
+real(dp)    :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot
+real(dp)    :: h2lift, h2lift_tot, h_pore, h_tot!, R_subl
+logical(k4) :: unexpected
+
+! CODE
+! ----
+call settle_h2o_surf_subsurf_tend(h2oice_depth,flux_ssice,d_h2oice,icell,jcell)
+
 dh_co2ice = d_co2ice/rho_co2ice
 dh_h2oice = d_h2oice/rho_h2oice
@@ -1298 +1444 @@
 else if (dh_co2ice >= 0._dp .and. dh_h2oice >= 0._dp .and. dh_dust >= 0._dp) then
     ! Which porosity is considered?
-    if (dh_co2ice > dh_h2oice .and. dh_co2ice > dh_dust) then ! CO2 ice is dominant
+    if (dh_co2ice >= dh_h2oice .and. dh_co2ice >= dh_dust) then ! CO2 ice is dominant
         h_pore = dh_co2ice*co2ice_porosity/(1._dp - co2ice_porosity)
-    else if (dh_h2oice > dh_co2ice .and. dh_h2oice > dh_dust) then ! H2O ice is dominant
+    else if (dh_h2oice >= dh_co2ice .and. dh_h2oice >= dh_dust) then ! H2O ice is dominant
         h_pore = dh_h2oice*h2oice_porosity/(1._dp - h2oice_porosity)
-    else if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice) then ! Dust is dominant
+    else if (dh_dust >= dh_co2ice .and. dh_dust >= dh_h2oice) then ! Dust is dominant
         h_pore = dh_dust*regolith_porosity/(1._dp - regolith_porosity)
     else
@@ -1330 +1476 @@
 
             ! How much is CO2 ice sublimation inhibited by the top dust lag?
-            R_subl = 1.
-            if (associated(current%up)) then ! If there is an upper stratum
-                h_toplag = 0._dp
-                currentloc => current%up
-                do while (is_dust_lag(currentloc) .and. associated(currentloc%up))
-                    h_toplag = h_toplag + thickness(currentloc)
-                    currentloc => currentloc%up
-                end do
-                if (currentloc%h_h2oice >= h_patchy_ice) then
-                    R_subl = 0._dp
-                else if (0._dp < currentloc%h_h2oice .and. currentloc%h_h2oice < h_patchy_ice) then
-                    h_toplag = h_toplag + thickness(currentloc)
-                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
-                else
-                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
-                end if
-            end if
-           h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl
+            !call get_co2_subl_resistance(current,R_subl)
+            !h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl
 
             ! Is there CO2 ice to sublimate?
@@ -1388 +1518 @@
             end if
         end do
+        ! In case the retreat loop exhausted all strata while residual demand remains then the tendency is set to 0
         if (h2subl_co2ice_tot > 0._dp) dh_co2ice = 0._dp ! No CO2 ice available anymore so the tendency is set to 0
         if (h2subl_h2oice_tot > 0._dp) dh_h2oice = 0._dp ! No H2O ice available anymore so the tendency is set to 0
@@ -1409 +1540 @@
 
             ! How much is CO2 ice sublimation inhibited by the top dust lag?
-            R_subl = 1.
-            if (associated(current%up)) then ! If there is an upper stratum
-                h_toplag = 0._dp
-                currentloc => current%up
-                do while (is_dust_lag(currentloc) .and. associated(currentloc%up))
-                    h_toplag = h_toplag + thickness(currentloc)
-                    currentloc => currentloc%up
-                end do
-                if (currentloc%h_h2oice >= h_patchy_ice) then
-                    R_subl = 0._dp
-                else if (0._dp < currentloc%h_h2oice .and. currentloc%h_h2oice < h_patchy_ice) then
-                    h_toplag = h_toplag + thickness(currentloc)
-                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
-                else
-                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
-                end if
-            end if
-           h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl
+            !call get_co2_subl_resistance(current,R_subl)
+            !h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl
 
             ! Is there CO2 ice to sublimate?
@@ -1461 +1576 @@
             end if
         end do
+        ! In case the retreat loop exhausted all strata while residual demand remains then the tendency is set to 0
         if (h2subl_co2ice_tot > 0._dp) dh_co2ice = 0._dp ! No CO2 ice available anymore so the tendency is set to 0
         if (h2lift_tot > 0._dp) dh_dust = 0._dp ! No dust available anymore so the tendency is set to 0
@@ -1499 +1615 @@
 
             ! How much is CO2 ice sublimation inhibited by the top dust lag?
-            R_subl = 1._dp
-            if (associated(current%up)) then ! If there is an upper stratum
-                h_toplag = 0._dp
-                currentloc => current%up
-                do while (is_dust_lag(currentloc) .and. associated(currentloc%up))
-                    h_toplag = h_toplag + thickness(currentloc)
-                    currentloc => currentloc%up
-                end do
-                if (currentloc%h_h2oice >= h_patchy_ice) then
-                    R_subl = 0._dp
-                else if (0._dp < currentloc%h_h2oice .and. currentloc%h_h2oice < h_patchy_ice) then
-                    h_toplag = h_toplag + thickness(currentloc)
-                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
-                else
-                    R_subl = fred_sublco2**(h_toplag/hmin_lag_co2)
-                end if
-            end if
-           h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl
+            !call get_co2_subl_resistance(current,R_subl)
+            !h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl
 
             ! Is there CO2 ice to sublimate?
@@ -1551 +1651 @@
             end if
         end do
+        ! In case the retreat loop exhausted all strata while residual demand remains then the tendency is set to 0
         if (h2subl_co2ice_tot > 0._dp) dh_co2ice = 0._dp ! No CO2 ice available anymore so the tendency is set to 0
         if (h2lift_tot > 0._dp) dh_dust = 0._dp ! No dust available anymore so the tendency is set to 0
@@ -1567 +1668 @@
 end if
 
+! New tendencies for ice
+d_co2ice = dh_co2ice*rho_co2ice
+d_h2oice = dh_h2oice*rho_h2oice
+
 zshift_surf = this%top%top_elevation - zshift_surf
 

trunk/LMDZ.COMMON/libf/evolution/maths.F90

@@ -17 +17 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, di, k4, minieps
+use numerics, only: dp, di, k4
 
 ! DECLARATION

trunk/LMDZ.COMMON/libf/evolution/numerics.F90

@@ -74 +74 @@
 real(qp), parameter :: minieps_qp = epsilon(1._qp)
 real(wp), parameter :: minieps = epsilon(1._wp)
-real(wp), parameter :: tol = 100._wp*minieps
+real(wp), parameter :: eps    = 100._wp*minieps    ! ~ 10^(-14)
+real(wp), parameter :: eps_qp = 100._wp*minieps_qp ! ~ 10^(-14)
+real(wp), parameter :: tol    = 100._wp*eps        ! ~ 10^(-12)
 
 ! Minimum number of significant decimal digits

trunk/LMDZ.COMMON/libf/evolution/orbit.F90

@@ -15 +15 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, di, k4, largest_nb, minieps
+use numerics, only: dp, di, k4, largest_nb, eps
 
 ! DECLARATION
@@ -477 +477 @@
 call print_msg('Obl. (current|accepted min|accepted max): '//real2str(obliquity)//' | '//real2str(min_obl)//real2str(max_obl),LVL_NFO)
 
-max_ecc = min(eccentricity + max_change_ecc,1. - minieps) ! ecc < 1.
+max_ecc = min(eccentricity + max_change_ecc,1. - eps) ! ecc < 1.
 min_ecc = max(eccentricity - max_change_ecc,0._dp) ! ecc >= 0.
 call print_msg('Ecc. (current|accepted min|accepted max): '//real2str(eccentricity)//' | '//real2str(min_ecc)//' | '//real2str(max_ecc),LVL_NFO)

trunk/LMDZ.COMMON/libf/evolution/output.F90

@@ -16 +16 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, di, k4, minieps
+use numerics, only: dp, di, k4, eps
 
 ! DECLARATION
@@ -469 +469 @@
 
 ! If output timing not met, return
-if (abs(mod(idt,output_rate)) > minieps) return
+if (abs(mod(idt,output_rate)) > eps) return
 
 ! If it is the first writing of the current timestep

trunk/LMDZ.COMMON/libf/evolution/pem.F90

@@ -31 +31 @@
 use clim_state_init,    only: read_start, read_startfi, read_startevo
 use config,             only: read_rundef, read_display_config
-use display,            only: print_ini, print_end, print_msg, LVL_NFO, LVL_WRN
+use display,            only: print_ini, print_end, print_msg, is_lvl_enabled, LVL_NFO, LVL_WRN, LVL_DBG
 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, nsoil_PCM, nsoil, cell_area, total_surface
@@ -38 +38 @@
 use layered_deposits,   only: do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering, print_layering
 use maths,              only: pi
-use numerics,           only: dp, qp, di, li, k4, minieps, minieps_qp
+use numerics,           only: dp, qp, di, li, k4, eps, eps_qp
 use orbit,              only: evo_orbit, read_orbitpm, compute_maxyr_orbit, update_orbit
 use output,             only: write_diagevo, dim_ngrid, dim_nsoil
@@ -233 +233 @@
     if (do_layering) then
         h2oice_depth_old(:,:) = h2oice_depth(:,:)
-        ! The sublimating tendency coming from subsurface ice is given through the surface ice tendency
-        where (h2oice_depth(:,:) > 0. .and. flux_ssice_avg(:,:) < 0._dp) d_h2oice(:,:) = flux_ssice_avg(:,:)
         do islope = 1,nslope
             do i = 1,ngrid
-                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))
+                call evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice(i,islope),h2oice_depth(i,islope),flux_ssice_avg(i,islope), &
+                                     new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p,i,islope)
+                if (is_lvl_enabled(LVL_DBG)) call print_layering(layerings_map(i,islope))
             end do
         end do
@@ -255 +254 @@
         call evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit)
         call evolve_co2ice(co2_ice,d_co2ice,zshift_surf)
-    end if
-
-    ! Flow glaciers according to surface ice
-    if (co2ice_flow) then
-        allocate(is_co2ice_flow(ngrid,nslope))
-        call flow_co2glaciers(q_co2_ts,ps_ts,ps_avg_glob_old,ps_avg_glob,co2_ice,is_co2ice_flow)
-    end if
-    if (h2oice_flow) then
-        allocate(is_h2oice_flow(ngrid,nslope))
-        call flow_h2oglaciers(tsurf_avg,h2o_ice,is_h2oice_flow)
-    end if
-    if (do_layering) call surfice2layering(h2o_ice,co2_ice,layerings_map)
+
+        ! Flow glaciers according to surface ice
+        if (co2ice_flow) then
+            allocate(is_co2ice_flow(ngrid,nslope))
+            call flow_co2glaciers(q_co2_ts,ps_ts,ps_avg_glob_old,ps_avg_glob,co2_ice,is_co2ice_flow)
+        end if
+        if (h2oice_flow) then
+            allocate(is_h2oice_flow(ngrid,nslope))
+            call flow_h2oglaciers(tsurf_avg,h2o_ice,is_h2oice_flow)
+        end if
+    end if ! do_layering
 
     ! Adapt surface temperature if surface ice disappeared
@@ -337 +335 @@
 
     ! Checking mass balance for CO2
-    if (abs(CO2cond_ps_PCM - 1._dp) < minieps) then
+    if (abs(CO2cond_ps_PCM - 1._dp) < eps) then
         totmass_co2ice = 0._dp
         totmass_atmco2 = 0._dp
@@ -346 +344 @@
             end do
         end do
-        totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,minieps_qp) ! To avoid division by 0
+        totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,eps_qp) ! To avoid division by 0
         call print_msg('> Relative total CO2 mass balance = '//real2str(100._qp*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini)//' %',LVL_NFO)
         if (abs((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini) > 0.01_qp) then
             call print_msg('Mass balance is not conserved!',LVL_WRN)
-            totmass_ini = max(totmass_atmco2_ini,minieps_qp) ! To avoid division by 0
+            totmass_ini = max(totmass_atmco2_ini,eps_qp) ! To avoid division by 0
             call print_msg('    Atmospheric CO2 mass balance = '//real2str(100._qp*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini)//' %',LVL_WRN)
-            totmass_ini = max(totmass_co2ice_ini,minieps_qp) ! To avoid division by 0
+            totmass_ini = max(totmass_co2ice_ini,eps_qp) ! To avoid division by 0
             call print_msg('    CO2 ice mass balance         = '//real2str(100._qp*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini)//' %',LVL_WRN)
-            totmass_ini = max(totmass_adsco2_ini,minieps_qp) ! To avoid division by 0
+            totmass_ini = max(totmass_adsco2_ini,eps_qp) ! To avoid division by 0
             call print_msg('    Adsorbed CO2 mass balance    = '//real2str(100._qp*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini)//' %',LVL_WRN)
         end if
@@ -412 +410 @@
 call update_workflow_status(n_yr_run,stopcrit%stop_code(),n_yr_sim,ntot_yr_sim)
 
-! Deallocation
-if (do_layering) then
-    deallocate(h2oice_depth,h2oice_depth_old,new_str,new_lag,current)
-    do islope = 1,nslope
-        do i = 1,ngrid
-            call del_layering(layerings_map(i,islope))
-        end do
-    end do
-end if
+! Clean-up
+if (do_layering) deallocate(h2oice_depth,h2oice_depth_old,new_str,new_lag,current)
 call end_allocation()
 

trunk/LMDZ.COMMON/libf/evolution/planet.F90

@@ -21 +21 @@
 ! ------------
 use numerics,         only: dp, qp, k4
-use layered_deposits, only: layering
+use layered_deposits, only: layering, del_layering
 
 ! DECLARATION
@@ -231 +231 @@
 implicit none
 
+! LOCAL VARIABLES
+! ---------------
+integer :: i, islope
+
 ! CODE
 ! ----
@@ -253 +257 @@
 delta_h2o_ads(:) = 0._dp
 delta_co2_ads(:) = 0._dp
+do islope = 1,nslope
+    do i = 1,ngrid
+        layerings_map(i,islope)%nb_str = 0
+        nullify(layerings_map(i,islope)%bottom,layerings_map(i,islope)%top)
+    end do
+end do
 
 END SUBROUTINE allocate_startevo_state
@@ -386 +396 @@
 !-----------------------------------------------------------------------
 
-! DECLARATION
-! -----------
-implicit none
+! DEPENDENCIES
+! ------------
+use geometry, only: ngrid, nslope
+
+! DECLARATION
+! -----------
+implicit none
+
+! LOCAL VARIABLES
+! ---------------
+integer :: i, islope
 
 ! CODE
@@ -407 +425 @@
 if (allocated(tsoil_ts)) deallocate(tsoil_ts)
 if (allocated(tsoil_ts_old)) deallocate(tsoil_ts_old)
-if (allocated(layerings_map)) deallocate(layerings_map)
+if (allocated(layerings_map)) then
+    do islope = 1,nslope
+        do i = 1,ngrid
+            call del_layering(layerings_map(i,islope))
+        end do
+    end do
+    deallocate(layerings_map)
+end if
 if (allocated(h2o_soildensity_avg)) deallocate(h2o_soildensity_avg)
 if (allocated(delta_co2_ads)) deallocate(delta_co2_ads)
@@ -424 +449 @@
 if (allocated(d_co2ice_ini)) deallocate(d_co2ice_ini)
 if (allocated(d_h2oice)) deallocate(d_h2oice)
-if (allocated(layerings_map)) deallocate(layerings_map)
 if (allocated(flux_ssice_avg)) deallocate(flux_ssice_avg)
 

trunk/LMDZ.COMMON/libf/evolution/soil_therm_inertia.F90

@@ -18 +18 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, di, k4, minieps
+use numerics, only: dp, di, k4, eps
 
 ! DECLARATION
@@ -143 +143 @@
     regolith_inertia(:,islope) = inertiedat(:,1)
     do ig = 1,ngrid
-        if (abs(h2o_ice(ig,islope)) < minieps) regolith_inertia(ig,islope) = TI_regolith_avg
+        if (abs(h2o_ice(ig,islope)) < eps) regolith_inertia(ig,islope) = TI_regolith_avg
         if (reg_thprop_dependp) then
             if (TI(ig,1,islope) < reg_inertie_thresold) then
@@ -191 +191 @@
         if (icetable_depth(ig,islope) > -1.e-6_dp) then
         ! 3.0 Case where it is perennial ice
-            if (icetable_depth(ig,islope) < minieps) then
+            if (icetable_depth(ig,islope) < eps) then
                 call get_ice_TI(.true.,1._dp,0._dp,ice_inertia)
                 do isoil = 1,nsoil

trunk/LMDZ.COMMON/libf/evolution/sorption.F90

@@ -18 +18 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, qp, di, k4, minieps
+use numerics, only: dp, qp, di, k4, eps
 
 ! DECLARATION
@@ -230 +230 @@
         deltam_reg_slope(ig,islope) = 0._dp
         do iloop = 1,index_breccia
-            if (TI(ig,iloop,islope) < inertie_thresold .and. abs(h2o_ice(ig,islope)) < minieps .and. abs(co2_ice(ig,islope)) < minieps) then
+            if (TI(ig,iloop,islope) < inertie_thresold .and. abs(h2o_ice(ig,islope)) < eps .and. abs(co2_ice(ig,islope)) < eps) then
                 if (iloop == 1) then
                     deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - h2o_ads_reg(ig,iloop,islope))*(layer(iloop))
@@ -347 +347 @@
     do islope = 1,nslope
         do iloop = 1,index_breccia
-            if (TI(ig,iloop,islope) < inertie_thresold .and. abs(h2o_ice(ig,islope)) < minieps .and. abs(co2_ice(ig,islope)) < minieps) then
+            if (TI(ig,iloop,islope) < inertie_thresold .and. abs(h2o_ice(ig,islope)) < eps .and. abs(co2_ice(ig,islope)) < eps) then
                 dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
                                                          (alpha*pco2_avg(ig) + sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope)))
             else
-                if (abs(co2_ads_reg(ig,iloop,islope)) < minieps) then !!! we are at first call
+                if (abs(co2_ads_reg(ig,iloop,islope)) < eps) then !!! we are at first call
                     dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig) &
                                                              /(alpha*pco2_avg(ig)+sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope)))
@@ -368 +368 @@
         deltam_reg_slope(ig,islope) = 0._dp
         do iloop = 1,index_breccia
-            if (TI(ig,iloop,islope) < inertie_thresold .and. abs(h2o_ice(ig,islope)) < minieps .and. abs(co2_ice(ig,islope)) < minieps) then
+            if (TI(ig,iloop,islope) < inertie_thresold .and. abs(h2o_ice(ig,islope)) < eps .and. abs(co2_ice(ig,islope)) < eps) then
                 if (iloop == 1) then
                     deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - co2_ads_reg(ig,iloop,islope))*(layer(iloop))

trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90

@@ -16 +16 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, qp, di, k4, minieps_qp
+use numerics, only: dp, di, k4
 
 ! DECLARATION

trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90

@@ -16 +16 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, qp, di, k4, minieps
+use numerics, only: dp, qp, di, k4, eps, tol
 
 ! DECLARATION
@@ -30 +30 @@
 integer(di),                              parameter :: WEIGHT_RESERVOIR = 5_di ! Weight by locally available H2O ice mass
 integer(di),                              parameter :: WEIGHT_DIRECTION = 6_di ! Weight by directional capacity (up/down)
-integer(di),                              private   :: weight_method = WEIGHT_DIRECTION ! Default method for balancing
+integer(di),                              private   :: weight_method = WEIGHT_COMBINED ! Default method for balancing
 real(dp),                                 parameter :: rho_co2ice = 1650._dp ! Density of CO2 ice [kg.m-3]
 real(dp),                                 parameter :: rho_h2oice = 920._dp  ! Density of H2O ice [kg.m-3]
@@ -420 +420 @@
 ! ---------------
 integer(di), parameter :: max_iter = 50_di ! Maximum number of iterations for the balancing procedure
-real(dp),    parameter :: eps = 1.e-12_dp ! Small number to prevent division by zero in weights normalization
-real(dp),    parameter :: tiny_corr = minieps ! Minimum correction to consider that progress is made in the balancing procedure
 integer(di)                              :: i, islope, iter
 integer(di)                              :: method_used
@@ -556 +554 @@
     Delta = Delta - Delta_used
 
-    if (abs(Delta_used) < tiny_corr) exit
+    ! Test if the balancing procedure progresses
+    if (abs(Delta_used) < tol) exit
 end do
 

trunk/LMDZ.COMMON/libf/evolution/surf_temp.F90

@@ -16 +16 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, di, k4, minieps
+use numerics, only: dp, di, k4, eps
 
 ! DECLARATION
@@ -182 +182 @@
     do i = 1,nlon
         do islope = 1,nslope
-            if (mask_co2ice_ini(i,j,islope) > 0.5_dp .and. co2ice_ll(i,j,islope) < minieps .and. mask_co2ice_disappeared(i,j,islope) < 0.5_dp) then
+            if (mask_co2ice_ini(i,j,islope) > 0.5_dp .and. co2ice_ll(i,j,islope) < eps .and. mask_co2ice_disappeared(i,j,islope) < 0.5_dp) then
                 found = .false.
                 mask_co2ice_disappeared(i,j,islope) = 1._dp

trunk/LMDZ.COMMON/libf/evolution/tendencies.F90

@@ -18 +18 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, di, k4, minieps, tol
+use numerics, only: dp, di, k4, eps
 
 ! DECLARATION
@@ -62 +62 @@
 
 ! If the difference is too small, then there is no evolution
-where (abs(d_ice) < minieps) d_ice = 0._dp
+where (abs(d_ice) < eps) d_ice = 0._dp
 
 ! If the tendency is negative but there is no ice reservoir for the PEM
-where (d_ice(:,:) < 0._dp .and. abs(perice(:,:)) < minieps) d_ice(:,:) = 0._dp
+where (d_ice(:,:) < 0._dp .and. abs(perice(:,:)) < eps) d_ice(:,:) = 0._dp
 
 END SUBROUTINE compute_tendice
@@ -170 +170 @@
 ! ARGUMENTS
 ! ---------
-real(dp), dimension(:,:),     intent(in) :: h2oice_depth_old ! Old H2O ice depth
-real(dp), dimension(:,:),     intent(in) :: h2oice_depth_new ! New H2O ice depth
-real(dp), dimension(:,:),     intent(in) :: tsurf            ! Surface temperature
-real(dp), dimension(:,:,:,:), intent(in) :: tsoil_ts_old     ! Old soil temperature time series
-real(dp), dimension(:,:,:,:), intent(in) :: tsoil_ts_new     ! New soil temperature time series
-real(dp), dimension(:,:),     intent(inout) :: flux_ssice_avg       ! Tendency of sub-surface ice
+real(dp), dimension(:,:),     intent(in)    :: h2oice_depth_old ! Old H2O ice depth
+real(dp), dimension(:,:),     intent(in)    :: h2oice_depth_new ! New H2O ice depth
+real(dp), dimension(:,:),     intent(in)    :: tsurf            ! Surface temperature
+real(dp), dimension(:,:,:,:), intent(in)    :: tsoil_ts_old     ! Old soil temperature time series
+real(dp), dimension(:,:,:,:), intent(in)    :: tsoil_ts_new     ! New soil temperature time series
+real(dp), dimension(:,:),     intent(inout) :: flux_ssice_avg   ! Tendency of sub-surface ice
 
 ! LOCAL VARIABLES
@@ -190 +190 @@
     do islope = 1,nslope
         ! Higher resistance due to growing lag layer (higher depth)
-        Rz_old = h2oice_depth_old(i,islope)*zcdv/coef_ssdif_PCM(i,islope) ! Old resistance from PCM
-        Rz_new = h2oice_depth_new(i,islope)*zcdv/coef_ssdif_PCM(i,islope) ! New resistance based on new depth
-        R_dec = Rz_old/Rz_new ! Decrease because of resistance
+        Rz_old = h2oice_depth_old(i,islope)*zcdv/max(abs(coef_ssdif_PCM(i,islope)),eps) ! Old resistance from PCM
+        Rz_new = h2oice_depth_new(i,islope)*zcdv/max(abs(coef_ssdif_PCM(i,islope)),eps) ! New resistance based on new depth
+        R_dec = Rz_old/max(Rz_new,eps) ! Decrease because of resistance
 
         ! The maxmimum of the daily averages over one year for the saturation vapor pressure at the ice table location
@@ -203 +203 @@
 
         ! Lower humidity due to growing lag layer (higher depth)
-        if (abs(psv_max_old) < tol) then
+        if (abs(psv_max_old) < eps) then
             hum_dec = 1._dp
         else

trunk/LMDZ.COMMON/libf/evolution/utility.F90

@@ -16 +16 @@
 ! DEPENDENCIES
 ! ------------
-use numerics, only: dp, qp, di, li, k4, minieps
+use numerics, only: dp, qp, di, li, k4, eps
 
 ! DECLARATION
@@ -302 +302 @@
 idigits = 1
 ! If x /= 0 then:
-if (abs(x) >= minieps) idigits = int(log10(abs(x))) + 1
+if (abs(x) >= eps) idigits = int(log10(abs(x))) + 1
 
 END FUNCTION nb_digits