Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 4180) @@ -963,2 +963,9 @@ - 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'. Index: trunk/LMDZ.COMMON/libf/evolution/config.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/config.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/config.F90 (revision 4180) @@ -17,5 +17,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, di, k4, minieps +use numerics, only: dp, di, k4, eps ! DECLARATION @@ -240,11 +240,13 @@ ! 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,13 +257,15 @@ ! ---- ! 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,5 +462,5 @@ 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 Index: trunk/LMDZ.COMMON/libf/evolution/display.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/display.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/display.F90 (revision 4180) @@ -19,5 +19,5 @@ 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,4 +277,40 @@ 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,5 +501,5 @@ 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) Index: trunk/LMDZ.COMMON/libf/evolution/glaciers.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/glaciers.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/glaciers.F90 (revision 4180) @@ -18,5 +18,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, di, k4, minieps +use numerics, only: dp, di, k4, eps ! DECLARATION @@ -301,5 +301,5 @@ 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 Index: trunk/LMDZ.COMMON/libf/evolution/io_netcdf.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/io_netcdf.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/io_netcdf.F90 (revision 4180) @@ -16,5 +16,5 @@ ! 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,5 +989,5 @@ 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,5 +1074,5 @@ 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,5 +1159,5 @@ 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,5 +1244,5 @@ 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,5 +1329,5 @@ 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 Index: trunk/LMDZ.COMMON/libf/evolution/layered_deposits.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/layered_deposits.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/layered_deposits.F90 (revision 4180) @@ -17,5 +17,5 @@ ! 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,4 +505,6 @@ end do this%nb_str = 0 +nullify(this%top) +nullify(this%bottom) END SUBROUTINE del_layering @@ -699,5 +701,8 @@ 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,4 +747,7 @@ ! CODE ! ---- +h2o_ice(:,:) = 0._dp +co2_ice(:,:) = 0._dp +h2oice_depth(:,:) = 0._dp do ig = 1,ngrid do islope = 1,nslope @@ -1070,17 +1078,13 @@ ! 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,6 +1095,6 @@ 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,6 +1132,6 @@ ! 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,4 +1192,8 @@ 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,5 +1232,140 @@ !======================================================================= -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,5 +1382,5 @@ ! ! NOTES -! Creates dust lag layers when ice sublimation occurs. +! !----------------------------------------------------------------------- @@ -1255,21 +1398,24 @@ ! 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,9 +1444,9 @@ 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,22 +1476,6 @@ ! 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,4 +1518,5 @@ 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,22 +1540,6 @@ ! 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,4 +1576,5 @@ 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,22 +1615,6 @@ ! 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,4 +1651,5 @@ 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,4 +1668,8 @@ 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 Index: trunk/LMDZ.COMMON/libf/evolution/maths.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/maths.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/maths.F90 (revision 4180) @@ -17,5 +17,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, di, k4, minieps +use numerics, only: dp, di, k4 ! DECLARATION Index: trunk/LMDZ.COMMON/libf/evolution/numerics.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/numerics.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/numerics.F90 (revision 4180) @@ -74,5 +74,7 @@ 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 Index: trunk/LMDZ.COMMON/libf/evolution/orbit.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/orbit.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/orbit.F90 (revision 4180) @@ -15,5 +15,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, di, k4, largest_nb, minieps +use numerics, only: dp, di, k4, largest_nb, eps ! DECLARATION @@ -477,5 +477,5 @@ 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) Index: trunk/LMDZ.COMMON/libf/evolution/output.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/output.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/output.F90 (revision 4180) @@ -16,5 +16,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, di, k4, minieps +use numerics, only: dp, di, k4, eps ! DECLARATION @@ -469,5 +469,5 @@ ! 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 Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 4180) @@ -31,5 +31,5 @@ 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,5 +38,5 @@ 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,10 +233,9 @@ 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,16 +254,15 @@ 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,5 +335,5 @@ ! 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,13 +344,13 @@ 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,13 +410,6 @@ 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() Index: trunk/LMDZ.COMMON/libf/evolution/planet.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/planet.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/planet.F90 (revision 4180) @@ -21,5 +21,5 @@ ! ------------ use numerics, only: dp, qp, k4 -use layered_deposits, only: layering +use layered_deposits, only: layering, del_layering ! DECLARATION @@ -231,4 +231,8 @@ implicit none +! LOCAL VARIABLES +! --------------- +integer :: i, islope + ! CODE ! ---- @@ -253,4 +257,10 @@ 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,7 +396,15 @@ !----------------------------------------------------------------------- -! DECLARATION -! ----------- -implicit none +! DEPENDENCIES +! ------------ +use geometry, only: ngrid, nslope + +! DECLARATION +! ----------- +implicit none + +! LOCAL VARIABLES +! --------------- +integer :: i, islope ! CODE @@ -407,5 +425,12 @@ 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,5 +449,4 @@ 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) Index: trunk/LMDZ.COMMON/libf/evolution/soil_therm_inertia.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/soil_therm_inertia.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/soil_therm_inertia.F90 (revision 4180) @@ -18,5 +18,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, di, k4, minieps +use numerics, only: dp, di, k4, eps ! DECLARATION @@ -143,5 +143,5 @@ 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,5 +191,5 @@ 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 Index: trunk/LMDZ.COMMON/libf/evolution/sorption.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/sorption.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/sorption.F90 (revision 4180) @@ -18,5 +18,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, qp, di, k4, minieps +use numerics, only: dp, qp, di, k4, eps ! DECLARATION @@ -230,5 +230,5 @@ 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,9 +347,9 @@ 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,5 +368,5 @@ 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)) Index: trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90 (revision 4180) @@ -16,5 +16,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, qp, di, k4, minieps_qp +use numerics, only: dp, di, k4 ! DECLARATION Index: trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90 (revision 4180) @@ -16,5 +16,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, qp, di, k4, minieps +use numerics, only: dp, qp, di, k4, eps, tol ! DECLARATION @@ -30,5 +30,5 @@ 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,6 +420,4 @@ ! --------------- 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,5 +554,6 @@ 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 Index: trunk/LMDZ.COMMON/libf/evolution/surf_temp.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/surf_temp.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/surf_temp.F90 (revision 4180) @@ -16,5 +16,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, di, k4, minieps +use numerics, only: dp, di, k4, eps ! DECLARATION @@ -182,5 +182,5 @@ 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 Index: trunk/LMDZ.COMMON/libf/evolution/tendencies.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/tendencies.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/tendencies.F90 (revision 4180) @@ -18,5 +18,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, di, k4, minieps, tol +use numerics, only: dp, di, k4, eps ! DECLARATION @@ -62,8 +62,8 @@ ! 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,10 +170,10 @@ ! 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,7 +190,7 @@ 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,5 +203,5 @@ ! 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 Index: trunk/LMDZ.COMMON/libf/evolution/utility.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/utility.F90 (revision 4174) +++ trunk/LMDZ.COMMON/libf/evolution/utility.F90 (revision 4180) @@ -16,5 +16,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, qp, di, li, k4, minieps +use numerics, only: dp, qp, di, li, k4, eps ! DECLARATION @@ -302,5 +302,5 @@ 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