Index: trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90 (revision 4170) +++ trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90 (revision 4174) @@ -31,12 +31,12 @@ ! --------- 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,11 +120,11 @@ ! 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,12 +160,12 @@ ! ---- 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,12 +200,14 @@ ! 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,95 +390,3 @@ !======================================================================= -!======================================================================= -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