Changeset 3967

Timestamp:

Nov 20, 2025, 3:01:08 PM (4 weeks ago)

Author:

jbclement

Message:

PEM:
New management of stopping criteria: giving explicit names and collecting everything in the same place.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• evol_ice_mod.F90 (modified) (3 diffs)
• pem.F90 (modified) (7 diffs)
• stopping_crit_mod.F90 (modified) (17 diffs)

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

@@ -794 +794 @@
 == 14/11/2025 == JBC
 Clearing all the tags and code related to the Generic PCM.
+
+== 20/11/2025 == JBC
+New management of stopping criteria: giving explicit names and collecting everything in the same place.

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

@@ -49 +49 @@
 !=======================================================================
 
-SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
+SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopCrit)
 
 use time_evol_mod,     only: dt
 use glaciers_mod,      only: rho_h2oice
-use stopping_crit_mod, only: stopping_crit_h2o
+use stopping_crit_mod, only: stopping_crit_h2o, stopFlags
 
 implicit none
@@ -73 +73 @@
 real, dimension(ngrid,nslope), intent(inout) :: h2o_ice  ! H2O ice (kg/m^2)
 real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Tendency of H2O ice (kg/m^2/year)
-integer,                       intent(inout) :: stopPEM  ! Stopping criterion code
+type(stopFlags),               intent(inout) :: stopCrit ! Stopping criterion
 real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]
 
@@ -88 +88 @@
     zshift_surf = d_h2oice*dt/rho_h2oice
 else ! In 3D
-    call stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopPEM)
-    if (stopPEM > 0) return
+    call stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopCrit)
+    if (stopCrit%stop_code() > 0) return
 
     call balance_h2oice_reservoirs(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)

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

@@ -41 +41 @@
 use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
                                       metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice, computeTcondCO2
-use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2, stopping_crit_h2o
+use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2, stopping_crit_h2o, stopFlags
 use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
 use evol_ice_mod,               only: evol_co2_ice, evol_h2o_ice, balance_h2oice_reservoirs
@@ -167 +167 @@
 real,   dimension(:,:,:), allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
 real,   dimension(:,:,:), allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step
-integer                               :: stopPEM              ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu
+type(stopFlags)                       :: stopCrit             ! Stopping criteria
 real                                  :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
 real,   dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Grid points x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
@@ -612 +612 @@
 write(*,*) '********* PEM cycle *********'
 i_myear_leg = 0
-stopPEM = 0
 if (layering_algo) then
     allocate(h2o_ice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
@@ -700 +699 @@
 
         allocate(d_h2oice_new(ngrid,nslope))
-        call stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopPEM)
+        call stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopCrit)
         call balance_h2oice_reservoirs(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
     
@@ -722 +721 @@
     else
         zlag = 0.
-        call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
+        call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopCrit)
         call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
     endif
@@ -928 +927 @@
 !------------------------
     write(*,*) "> Checking the stopping criteria"
-    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
-    call stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopPEM,ngrid,ps_avg_global_ini,ps_avg_global_new,nslope)
+    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopCrit,ngrid)
+    call stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopCrit,ngrid,ps_avg_global_ini,ps_avg_global_new,nslope)
     i_myear_leg = i_myear_leg + dt
     i_myear = i_myear + dt
-    if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
-    if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
+    if (stopCrit%stop_code() == 0 .and. i_myear_leg >= n_myear_leg) stopCrit%max_iter_reached = .true.
+    if (stopCrit%stop_code() == 0 .and. i_myear >= n_myear) stopCrit%max_years_reached = .true.
     call system_clock(c2)
-    if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
-    if (stopPEM > 0) then
-        select case (stopPEM)
-            case(1)
-                write(*,'(a,i0,a)') " **** STOPPING because surface of h2o ice sublimating is too low: ", stopPEM, ". See message above."
-            case(2)
-                write(*,'(a,i0,a)') " **** STOPPING because tendencies on h2o ice = 0: ", stopPEM, ". See message above."
-            case(3)
-                write(*,'(a,i0,a)') " **** STOPPING because surface of co2 ice sublimating is too low: ", stopPEM, ". See message above."
-            case(4)
-                write(*,'(a,i0,a)') " **** STOPPING because surface global pressure changed too much: ", stopPEM, ". See message above."
-            case(5)
-                write(*,'(a,i0)') " **** STOPPING because maximum number of iterations is reached (possibly due to orbital parameters): ", stopPEM
-            case(6)
-                write(*,'(a,i0)') " **** STOPPING because maximum number of Martian years to be simulated is reached: ", stopPEM
-            case(7)
-                write(*,'(a,i0)') " **** STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM
-            case default
-                write(*,'(a,i0)') " **** STOPPING with unknown stopping criterion code: ", stopPEM
-        end select
+    if (stopCrit%stop_code() == 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopCrit%time_limit_reached = .true.
+    if (stopCrit%is_any_set()) then
+        write(*,*) stopCrit%stop_message()
         exit
     else
@@ -1168 +1150 @@
              co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,layerings_map)
 
-call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
+call info_PEM(i_myear_leg,stopCrit%stop_code(),i_myear,n_myear)
 
 write(*,*)

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

@@ -3 +3 @@
 implicit none
 
-!=======================================================================
+type :: stopFlags
+    logical :: surface_h2oice_change = .false.
+    logical :: zero_dh2oice          = .false.
+    logical :: surface_co2ice_change = .false.
+    logical :: pressure_change       = .false.
+    logical :: max_iter_reached      = .false.
+    logical :: max_years_reached     = .false.
+    logical :: time_limit_reached    = .false.
 contains
-!=======================================================================
+    procedure :: is_any_set
+    procedure :: stop_code
+    procedure :: stop_message
+end type
+
+!=======================================================================
+contains
+!=======================================================================
+! To detect if any flag is true
+FUNCTION is_any_set(this) RESULT(stopPEM)
+
+class(stopFlags), intent(in) :: this
+logical :: stopPEM
+
+stopPEM = this%surface_h2oice_change .or. &
+          this%zero_dh2oice          .or. &
+          this%surface_co2ice_change .or. &
+          this%pressure_change       .or. &
+          this%max_iter_reached      .or. &
+          this%max_years_reached     .or. &
+          this%time_limit_reached
+
+END FUNCTION is_any_set
+
+!=======================================================================
+! To give the digit code corresponding to the stopping flag
+FUNCTION stop_code(this) result(code)
+    
+class(StopFlags), intent(in) :: this
+integer :: code
+
+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%max_iter_reached)      then; code = 5
+elseif (this%max_years_reached)     then; code = 6
+elseif (this%time_limit_reached)    then; code = 7
+endif
+
+END FUNCTION stop_code
+
+!=======================================================================
+! To give the message corresponding to the stopping flag
+FUNCTION stop_message(this) result(msg)
+    
+class(StopFlags), intent(in) :: this
+character(:), allocatable :: msg
+integer :: stopCode
+
+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%max_iter_reached)      then; msg = " **** STOPPING because maximum number of iterations is reached (possibly due to orbital parameters)."
+elseif (this%max_years_reached)     then; msg = " **** STOPPING because maximum number of Martian years to be simulated is reached."
+elseif (this%time_limit_reached)    then; msg = " **** STOPPING because the job time limit is going to be reached."
+endif
+
+END FUNCTION stop_message
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -14 +81 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-SUBROUTINE stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
+SUBROUTINE stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopCrit,ngrid)
 
 use time_evol_mod, only: h2o_ice_crit
@@ -35 +102 @@
 ! Outputs
 !--------
-integer, intent(inout) :: stopPEM ! Stopping criterion code
+type(stopFlags), intent(inout) :: stopCrit ! Stopping criterion
 ! Locals
 ! ------
@@ -42 +109 @@
 
 !=======================================================================
-if (stopPEM > 0) return
+if (stopCrit%stop_code() > 0) return
 
 ! Computation of the present surface of h2o ice still sublimating
@@ -54 +121 @@
 ! Check of the criterion
 if (h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit)) then
-    stopPEM = 1
+    stopCrit%surface_h2oice_change = .true.
     write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
     write(*,*) "h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit)", h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit)
@@ -61 +128 @@
     write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
 else if (h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit)) then
-    stopPEM = 1
+    stopCrit%surface_h2oice_change = .true.
     write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
     write(*,*) "h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit)", h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit)
@@ -73 +140 @@
 !=======================================================================
 
-SUBROUTINE stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopPEM,ngrid,ps_avg_global_ini,ps_avg_global,nslope)
+SUBROUTINE stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopCrit,ngrid,ps_avg_global_ini,ps_avg_global,nslope)
 
 use time_evol_mod, only: co2_ice_crit, ps_criterion
@@ -96 +163 @@
 ! Outputs
 !--------
-integer, intent(inout) :: stopPEM ! Stopping criterion code
+type(stopFlags), intent(inout) :: stopCrit ! Stopping criterion
 
 ! Locals
@@ -104 +171 @@
 
 !=======================================================================
-if (stopPEM > 0) return
+if (stopCrit%stop_code() > 0) return
 
 ! Computation of the present surface of co2 ice still sublimating
@@ -116 +183 @@
 ! Check of the criterion
 if (co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2_ice_crit)) then
-    stopPEM = 3
+    stopCrit%surface_co2ice_change = .true.
     write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
     write(*,*) "co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2_ice_crit)", co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2_ice_crit)
@@ -123 +190 @@
     write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
 else if (co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2_ice_crit)) then
-    stopPEM = 3
+    stopCrit%surface_co2ice_change = .true.
     write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
     write(*,*) "co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2_ice_crit)", co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2_ice_crit)
@@ -132 +199 @@
 
 if (ps_avg_global < ps_avg_global_ini*(1. - ps_criterion)) then
-    stopPEM = 4
+    stopCrit%pressure_change = .true.
     write(*,*) "Reason of stopping: the global pressure reaches the threshold"
     write(*,*) "ps_avg_global < ps_avg_global_ini*(1. - ps_criterion)", ps_avg_global < ps_avg_global_ini*(1. - ps_criterion)
@@ -139 +206 @@
     write(*,*) "Percentage of change accepted =", ps_criterion*100.
 else if (ps_avg_global > ps_avg_global_ini*(1. + ps_criterion)) then
-    stopPEM = 4
+    stopCrit%pressure_change = .true.
     write(*,*) "Reason of stopping: the global pressure reaches the threshold"
     write(*,*) "ps_avg_global > ps_avg_global_ini*(1. + ps_criterion)", ps_avg_global > ps_avg_global_ini*(1. + ps_criterion)
@@ -151 +218 @@
 !=======================================================================
 
-SUBROUTINE stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopPEM)
+SUBROUTINE stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopCrit)
 
 use time_evol_mod, only: dt
@@ -174 +241 @@
 ! Outputs
 !--------
-integer, intent(inout) :: stopPEM ! Stopping criterion code
-real,    intent(out) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to conserve H2O
+type(stopFlags), intent(inout) :: stopCrit ! Stopping criterion
+real, intent(out) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to conserve H2O
 
 ! Locals
@@ -182 +249 @@
 
 !=======================================================================
-if (stopPEM > 0) return
+if (stopCrit%stop_code() > 0) return
 
 ! We compute the amount of water going out from the atmosphere (S_atm_2_h2o) and going into the atmophere (S_h2o_2_atm)
@@ -218 +285 @@
     write(*,*) "Amount of condensing ice  =", S_atm_2_h2oice
     write(*,*) "Amount of sublimating ice =", S_h2oice_2_atm
-    stopPEM = 2
+    stopCrit%zero_dh2oice = .true.
 endif