Changeset 4174 for trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90

Timestamp:

Apr 8, 2026, 11:36:16 AM (2 days ago)

Author:

jbclement

Message:

PEM:
Refactor of H2O flux balancing:

• centralizes flux balancing logic;
• removes intermediate variables, especially related to of H2O mass bookkeeping;
• adds new stopping criterion when H2O flux balance fails8 (sanity check);
• introduces configurable weighting strategies for H2O flux balancing.

JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90 (modified) (5 diffs)

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

@@ -31 +31 @@
 ! ---------
 type :: stopFlags
-    logical(k4) :: surface_h2oice_change  = .false.
-    logical(k4) :: zero_dh2oice           = .false.
-    logical(k4) :: surface_co2ice_change  = .false.
-    logical(k4) :: pressure_change        = .false.
-    logical(k4) :: nmax_yr_run_reached    = .false.
-    logical(k4) :: nmax_yr_runorb_reached = .false.
-    logical(k4) :: nmax_yr_sim_reached    = .false.
-    logical(k4) :: time_limit_reached     = .false.
+    logical(k4) :: surface_h2oice_change     = .false.
+    logical(k4) :: h2o_flux_balance_unclosed = .false.
+    logical(k4) :: surface_co2ice_change     = .false.
+    logical(k4) :: pressure_change           = .false.
+    logical(k4) :: nmax_yr_run_reached       = .false.
+    logical(k4) :: nmax_yr_runorb_reached    = .false.
+    logical(k4) :: nmax_yr_sim_reached       = .false.
+    logical(k4) :: time_limit_reached        = .false.
 contains
     procedure :: is_any_set
@@ -120 +120 @@
 ! CODE
 ! ----
-stopPEM = this%surface_h2oice_change  .or. &
-          this%zero_dh2oice           .or. &
-          this%surface_co2ice_change  .or. &
-          this%pressure_change        .or. &
-          this%nmax_yr_run_reached    .or. &
-          this%nmax_yr_runorb_reached .or. &
-          this%nmax_yr_sim_reached    .or. &
+stopPEM = this%surface_h2oice_change     .or. &
+          this%h2o_flux_balance_unclosed .or. &
+          this%surface_co2ice_change     .or. &
+          this%pressure_change           .or. &
+          this%nmax_yr_run_reached       .or. &
+          this%nmax_yr_runorb_reached    .or. &
+          this%nmax_yr_sim_reached       .or. &
           this%time_limit_reached
 
@@ -160 +160 @@
 ! ----
 code = 0
-if     (this%surface_h2oice_change)  then; code = 1
-elseif (this%zero_dh2oice)           then; code = 2
-elseif (this%surface_co2ice_change)  then; code = 3
-elseif (this%pressure_change)        then; code = 4
-elseif (this%nmax_yr_run_reached)    then; code = 5
-elseif (this%nmax_yr_runorb_reached) then; code = 6
-elseif (this%nmax_yr_sim_reached)    then; code = 7
-elseif (this%time_limit_reached)     then; code = 8
+if     (this%surface_h2oice_change)     then; code = 1
+elseif (this%h2o_flux_balance_unclosed) then; code = 2
+elseif (this%surface_co2ice_change)     then; code = 3
+elseif (this%pressure_change)           then; code = 4
+elseif (this%nmax_yr_run_reached)       then; code = 5
+elseif (this%nmax_yr_runorb_reached)    then; code = 6
+elseif (this%nmax_yr_sim_reached)       then; code = 7
+elseif (this%time_limit_reached)        then; code = 8
 end if
 
@@ -200 +200 @@
 ! CODE
 ! ----
-if     (this%surface_h2oice_change)  then; msg = "**** STOPPING because surface of H2O ice has changed too much (see message above)."
-elseif (this%zero_dh2oice)           then; msg = "**** STOPPING because tendencies on H2O ice = 0 (see message above)."
-elseif (this%surface_co2ice_change)  then; msg = "**** STOPPING because surface of CO2 ice has changed too much (see message above)."
-elseif (this%pressure_change)        then; msg = "**** STOPPING because surface global pressure has changed too much (see message above)."
-elseif (this%nmax_yr_run_reached)    then; msg = "**** STOPPING because maximum number of years is reached."
-elseif (this%nmax_yr_runorb_reached) then; msg = "**** STOPPING because maximum number of years due to orbital parameters is reached."
-elseif (this%nmax_yr_sim_reached)    then; msg = "**** STOPPING because the number of years to be simulated is reached."
-elseif (this%time_limit_reached)     then; msg = "**** STOPPING because the job time limit is going to be reached."
+if     (this%surface_h2oice_change)     then; msg = "**** STOPPING because surface of H2O ice has changed too much (see message above)."
+elseif (this%h2o_flux_balance_unclosed) then; msg = "**** STOPPING because H2O flux imbalance cannot be closed under physical bounds (see message above)."
+elseif (this%surface_co2ice_change)     then; msg = "**** STOPPING because surface of CO2 ice has changed too much (see message above)."
+elseif (this%pressure_change)           then; msg = "**** STOPPING because surface global pressure has changed too much (see message above)."
+elseif (this%nmax_yr_run_reached)       then; msg = "**** STOPPING because maximum number of years is reached."
+elseif (this%nmax_yr_runorb_reached)    then; msg = "**** STOPPING because maximum number of years due to orbital parameters is reached."
+elseif (this%nmax_yr_sim_reached)       then; msg = "**** STOPPING because the number of years to be simulated is reached."
+elseif (this%time_limit_reached)        then; msg = "**** STOPPING because the job time limit is going to be reached."
+else
+    msg = "**** STOPPING for unknown reason (this should not happen!)."
 end if
 
@@ -388 +390 @@
 !=======================================================================
 
-!=======================================================================
-SUBROUTINE 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)
-!-----------------------------------------------------------------------
-! NAME
-!     stopping_crit_h2o
-!
-! DESCRIPTION
-!     Check if PEM oscillates infinitely with H2O (no net exchange).
-!
-! AUTHORS & DATE
-!     L. Lange, 2024
-!     JB Clement, 2025
-!
-! NOTES
-!
-!-----------------------------------------------------------------------
-
-! DEPENDENCIES
-! ------------
-use geometry,  only: ngrid, nslope, cell_area
-use evolution, only: dt
-use slopes,    only: subslope_dist, def_slope_mean
-use maths,     only: pi
-use display,   only: print_msg, LVL_WRN
-use utility,   only: real2str
-
-! DECLARATION
-! -----------
-implicit none
-
-! ARGUMENTS
-! ---------
-real(dp), dimension(:),   intent(in)    :: delta_h2o_ads            ! Mass of H2O adsorbed/desorbded in the soil (kg/m^2)
-real(dp), dimension(:),   intent(in)    :: delta_icetable ! Mass of H2O that condensed/sublimated at the ice table
-real(dp), dimension(:,:), intent(in)    :: h2o_ice, d_h2oice
-real(qp),                 intent(out)   :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm
-type(stopFlags),          intent(inout) :: stopcrit
-
-! LOCAL VARIABLES
-! ---------------
-integer(di) :: i, islope
-
-! CODE
-! ----
-! Check to escape the subroutine (not relevant in 1D or if the PEM should stop already)
-if (ngrid == 1 .or. stopcrit%stop_code() > 0) then
-    S_atm_2_h2o = 1._qp
-    S_h2o_2_atm = S_atm_2_h2o
-    S_atm_2_h2oice = 1._qp
-    S_h2oice_2_atm = S_atm_2_h2oice
-    return
-end if
-
-! We compute the amount of water going out from the atmosphere (S_atm_2_h2o) and going into the atmophere (S_h2o_2_atm)
-S_atm_2_h2o = 0._qp
-S_h2o_2_atm = 0._qp
-S_atm_2_h2oice = 0._qp
-S_h2oice_2_atm = 0._qp
-do i = 1,ngrid
-    if (delta_h2o_ads(i) > 0._dp) then
-        S_atm_2_h2o = S_atm_2_h2o + delta_h2o_ads(i)*cell_area(i)
-    else
-        S_h2o_2_atm = S_h2o_2_atm - delta_h2o_ads(i)*cell_area(i)
-    end if
-    if (delta_icetable(i) > 0._dp) then
-        S_atm_2_h2o = S_atm_2_h2o + delta_icetable(i)*cell_area(i)
-    else
-        S_h2o_2_atm = S_h2o_2_atm - delta_icetable(i)*cell_area(i)
-    end if
-    do islope = 1,nslope
-        if (d_h2oice(i,islope) > 0._dp) then
-            S_atm_2_h2o = S_atm_2_h2o + d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180._dp)
-            S_atm_2_h2oice = S_atm_2_h2oice + d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180._dp)
-        else if (d_h2oice(i,islope) < 0._dp .and. h2o_ice(i,islope) > 0._dp) then
-            S_h2o_2_atm = S_h2o_2_atm - d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180._dp)
-            S_h2oice_2_atm = S_h2oice_2_atm - d_h2oice(i,islope)*dt*cell_area(i)*subslope_dist(i,islope)/cos(def_slope_mean(islope)*pi/180._dp)
-        end if
-    end do
-end do
-
-! Since relative atmospheric water is kept constant, we need to equate condensing reservoirs to the sublimating ones.
-! It is not possible if one of them is missing.
-if ((abs(S_atm_2_h2oice) < minieps_qp .or. abs(S_h2oice_2_atm) < minieps_qp)) then
-    call print_msg("There is no sublimating or condensing H2O ice! This can be due to the absence of H2O ice in the PCM run.",LVL_WRN)
-    call print_msg("Amount of condensing ice  [kg] = "//real2str(S_atm_2_h2oice),LVL_WRN)
-    call print_msg("Amount of sublimating ice [kg] = "//real2str(S_h2oice_2_atm),LVL_WRN)
-    stopcrit%zero_dh2oice = .true.
-end if
-
-END SUBROUTINE stopping_crit_h2o
-!=======================================================================
-
 END MODULE stopping_crit