Index: trunk/LMDZ.COMMON/libf/evolution/backup.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/backup.F90 (revision 4173) +++ trunk/LMDZ.COMMON/libf/evolution/backup.F90 (revision 4174) @@ -70,5 +70,5 @@ !======================================================================= SUBROUTINE save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, & - icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map,backup_idt) + icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,h2o_ads_reg,co2_ads_reg,layerings_map,backup_idt) !----------------------------------------------------------------------- ! NAME @@ -107,5 +107,5 @@ real(dp), dimension(:), intent(in) :: ps_avg, ps_dev real(dp), intent(in) :: ps_avg_glob, ps_avg_glob_ini -real(dp), dimension(:,:), intent(in) :: tsurf_avg, tsurf_dev, icetable_depth, icetable_thickness, h2oice_depth, flux_ssice_avg +real(dp), dimension(:,:), intent(in) :: tsurf_avg, tsurf_dev, icetable_depth, icetable_thickness, h2oice_depth real(dp), dimension(:,:,:), intent(in) :: tsoil_avg, tsoil_dev, ice_porefilling, h2o_ads_reg, co2_ads_reg type(layering), dimension(:,:), intent(in) :: layerings_map @@ -115,5 +115,5 @@ ! LOCAL VARIABLES ! --------------- -real(dp), dimension(ngrid,nslope) :: h2o_ice4PCM, co2_ice4PCM, tsurf4PCM, flux_geo4PCM, albedo4PCM, emissivity4PCM, flux_ssice4PCM, h2oice_depth4PCM +real(dp), dimension(ngrid,nslope) :: h2o_ice4PCM, co2_ice4PCM, tsurf4PCM, flux_geo4PCM, albedo4PCM, emissivity4PCM, h2oice_depth4PCM real(dp), dimension(ngrid,nlayer) :: teta4PCM, air_mass4PCM real(dp), dimension(ngrid,nsoil_PCM,nslope) :: tsoil4PCM, inertiesoil4PCM Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 4173) +++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 4174) @@ -956,2 +956,9 @@ - Keeping a clear separation between the subsurface ice flux and the surface ice tendency. - Making sure that subsurface ice depth is well given to the PCM at the end of the PEM (ice table VS layering). + +== 08/04/2026 == JBC +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. Index: trunk/LMDZ.COMMON/libf/evolution/layered_deposits.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/layered_deposits.F90 (revision 4173) +++ trunk/LMDZ.COMMON/libf/evolution/layered_deposits.F90 (revision 4174) @@ -1255,9 +1255,9 @@ ! ARGUMENTS ! --------- -type(layering), intent(inout) :: this -type(stratum), pointer, intent(inout) :: current -real(dp), intent(inout) :: d_co2ice, d_h2oice -logical(k4), intent(inout) :: new_str, new_lag -real(dp), intent(out) :: zshift_surf, zlag +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 Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 4173) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 4174) @@ -33,5 +33,5 @@ use display, only: print_ini, print_end, print_msg, LVL_NFO, LVL_WRN 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, nday, nsoil_PCM, nsoil, cell_area, total_surface +use geometry, only: ngrid, nslope, nsoil_PCM, nsoil, cell_area, total_surface use glaciers, only: h2oice_flow, co2ice_flow, flow_co2glaciers, flow_h2oglaciers use ice_table, only: icetable_equilibrium, icetable_dynamic, evolve_ice_table @@ -48,7 +48,7 @@ use soil_therm_inertia, only: update_soil_TI use sorption, only: do_sorption, compute_totmass_adsorbed, evolve_regolith_adsorption -use stopping_crit, only: stopFlags, stopping_crit_pressure, stopping_crit_h2o, stopping_crit_h2oice, stopping_crit_co2ice +use stopping_crit, only: stopFlags, stopping_crit_pressure, stopping_crit_h2oice, stopping_crit_co2ice use surface, only: emissivity_PCM -use surf_ice, only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs +use surf_ice, only: evolve_co2ice, evolve_h2oice, balance_h2o_fluxes use surf_temp, only: tsurf_PCM, adapt_tsurf2disappearedice use tendencies, only: compute_tendice, evolve_tend_co2ice, evolve_flux_ssice @@ -80,6 +80,4 @@ real(dp), dimension(:,:), allocatable :: zshift_surf ! Elevation shift for the surface [m] real(dp), dimension(:,:), allocatable :: zlag ! Newly built lag thickness [m] -! Tendency-related: -real(dp), dimension(:,:), allocatable :: d_h2oice_new ! Adjusted tendency of perennial H2O ice to keep balance between donor and recipient [kg/m2/y] ! Layering-related: type(ptrarray), dimension(:,:), allocatable :: current ! Current active stratum in the layering @@ -95,5 +93,4 @@ real(qp) :: totmass_ini ! Initial total CO2 mass [kg] real(qp) :: totmass_adsh2o ! Current total adsorbed H2O mass [kg] -real(qp) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to balance H2O ice reservoirs ! CODE @@ -235,5 +232,4 @@ allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope)) if (do_layering) then - allocate(d_h2oice_new(ngrid,nslope)) h2oice_depth_old(:,:) = h2oice_depth(:,:) ! The sublimating tendency coming from subsurface ice is given through the surface ice tendency @@ -241,5 +237,5 @@ do islope = 1,nslope do i = 1,ngrid - call evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice_new(i,islope),new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p) + 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)) end do @@ -248,7 +244,10 @@ call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth) ! Balance H2O ice reservoirs - call 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) - call balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new) - deallocate(d_h2oice_new) + call balance_h2o_fluxes(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,stopcrit) + if (stopcrit%stop_code() > 0) then + deallocate(zshift_surf,zlag) + call print_msg(stopcrit%stop_message(),LVL_NFO) + exit + end if else zlag(:,:) = 0._dp @@ -376,5 +375,5 @@ if (backup_rate > 0) then if (mod(idt,backup_rate) == 0) call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, & - icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map,idt) + icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,h2o_ads_reg,co2_ads_reg,layerings_map,idt) end if @@ -408,5 +407,5 @@ call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, & - icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map) + icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,h2o_ads_reg,co2_ads_reg,layerings_map) ! Update the duration information of the workflow Index: trunk/LMDZ.COMMON/libf/evolution/sorption.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/sorption.F90 (revision 4173) +++ trunk/LMDZ.COMMON/libf/evolution/sorption.F90 (revision 4174) @@ -142,5 +142,5 @@ use slopes, only: subslope_dist, def_slope_mean use atmosphere, only: ap, bp -use physics, only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2, A, B +use physics, only: alpha_clap_h2o, beta_clap_h2o, m_h2o, A, B use maths, only: pi Index: trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90 (revision 4173) +++ 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 Index: trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90 (revision 4173) +++ trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90 (revision 4174) @@ -16,5 +16,5 @@ ! DEPENDENCIES ! ------------ -use numerics, only: dp, qp, di, k4 +use numerics, only: dp, qp, di, k4, minieps ! DECLARATION @@ -24,4 +24,11 @@ ! PARAMETERS ! ---------- +integer(di), parameter :: WEIGHT_UNIFORM = 1_di ! Equal weight for all adjustable points +integer(di), parameter :: WEIGHT_ABSFLUX = 2_di ! Weight by absolute local H2O surface flux +integer(di), parameter :: WEIGHT_CAPACITY = 3_di ! Weight by available local adjustment range +integer(di), parameter :: WEIGHT_COMBINED = 4_di ! Capacity weighted by current absolute flux +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 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] @@ -338,7 +345,7 @@ ! DEPENDENCIES ! ------------ -use geometry, only: ngrid, nslope +use geometry, only: ngrid use evolution, only: dt -use stopping_crit, only: stopping_crit_h2o, stopFlags +use stopping_crit, only: stopFlags use display, only: print_msg, LVL_NFO @@ -354,9 +361,4 @@ real(dp), dimension(:,:), intent(out) :: zshift_surf -! LOCAL VARIABLES -! --------------- -real(qp) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! H2O balance variables -real(dp), dimension(ngrid,nslope) :: d_h2oice_new ! Tendencies computed to keep balance - ! CODE ! ---- @@ -366,18 +368,13 @@ where (d_h2oice(:,:) < 0._dp .and. abs(d_h2oice(:,:)*dt) > h2o_ice(:,:)) ! Not enough ice to sublimate everything ! We sublimate what we can - d_h2oice_new(:,:) = h2o_ice(:,:)/dt - ! It means the tendency is zero next time - d_h2oice(:,:) = 0._dp - else where - d_h2oice_new(:,:) = d_h2oice(:,:) + d_h2oice(:,:) = -h2o_ice(:,:)/dt end where else ! In 3D - call 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) + call balance_h2o_fluxes(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,stopcrit) if (stopcrit%stop_code() > 0) return - call balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new) end if -h2o_ice(:,:) = h2o_ice(:,:) + d_h2oice_new(:,:)*dt -zshift_surf(:,:) = zshift_surf(:,:) + d_h2oice_new(:,:)*dt/rho_h2oice +h2o_ice(:,:) = h2o_ice(:,:) + d_h2oice(:,:)*dt +zshift_surf(:,:) = zshift_surf(:,:) + d_h2oice(:,:)*dt/rho_h2oice END SUBROUTINE evolve_h2oice @@ -385,76 +382,198 @@ !======================================================================= -SUBROUTINE balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_balanced) -!----------------------------------------------------------------------- -! NAME -! balance_h2oice_reservoirs -! -! DESCRIPTION -! Balance H2O flux from and into atmosphere across reservoirs. -! -! AUTHORS & DATE -! JB Clement, 2025 -! -! NOTES -! -!----------------------------------------------------------------------- - -! DEPENDENCIES -! ------------ -use evolution, only: dt -use geometry, only: ngrid, nslope - -! DECLARATION -! ----------- -implicit none - -! ARGUMENTS -! --------- -real(qp), intent(in) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm -real(dp), dimension(:,:), intent(in) :: h2o_ice -real(dp), dimension(:,:), intent(inout) :: d_h2oice -real(dp), dimension(:,:), intent(out) :: d_balanced +SUBROUTINE balance_h2o_fluxes(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,stopcrit,weight_type) +!----------------------------------------------------------------------- +! NAME +! balance_h2o_fluxes +! +! DESCRIPTION +! Balance H2O fluxes from and into atmosphere across reservoirs. +! +! AUTHORS & DATE +! JB Clement, 04/2025 +! +! NOTES +! +!----------------------------------------------------------------------- + +! DEPENDENCIES +! ------------ +use evolution, only: dt +use geometry, only: ngrid, nslope, cell_area +use stopping_crit, only: stopFlags +use utility, only: real2str, int2str +use display, only: print_msg, LVL_WRN + +! DECLARATION +! ----------- +implicit none + +! ARGUMENTS +! --------- +real(dp), dimension(:), intent(in) :: delta_h2o_ads, delta_icetable +real(dp), dimension(:,:), intent(in) :: h2o_ice +integer(di), intent(in), optional :: weight_type +real(dp), dimension(:,:), intent(inout) :: d_h2oice +type(stopFlags), intent(inout) :: stopcrit ! LOCAL VARIABLES ! --------------- -integer(di) :: i, islope -real(qp) :: S_corr_subl, S_corr_cond, S_target_subl, S_target_cond, S_ghostice ! Balance variables -real(dp) :: d_target - -! CODE -! ---- -! Compute global targets -S_corr_cond = (S_h2o_2_atm - S_atm_2_h2o)/2._qp ! Correction = deviation to the mean -S_corr_subl = (S_atm_2_h2o - S_h2o_2_atm)/2._qp ! Correction = deviation to the mean -S_target_cond = abs(S_atm_2_h2oice + S_corr_cond) -S_target_subl = abs(S_h2oice_2_atm + S_corr_subl) - -! Compute balanced tendencies with positivity limiter -d_balanced(:,:) = 0._dp -S_ghostice = 0._qp +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 +real(dp) :: flux_scale, tol_flux +real(qp) :: B_flux, B_flux_surfice, B_flux_nonsurfice, Delta, Delta_used, wsum +real(dp), dimension(:,:), allocatable :: w, d_corr, flux_new +real(dp), dimension(:,:), allocatable :: flux_min, flux_max +real(dp), dimension(:,:), allocatable :: capacity, capacity_up, capacity_down +logical(k4), dimension(:,:), allocatable :: active, adjustable +logical(k4), dimension(:,:), allocatable :: clipped_high, clipped_low + +! CODE +! ---- +! Initialization +allocate(w(ngrid,nslope),d_corr(ngrid,nslope),flux_new(ngrid,nslope)) +allocate(flux_min(ngrid,nslope),flux_max(ngrid,nslope)) +allocate(capacity(ngrid,nslope),capacity_up(ngrid,nslope),capacity_down(ngrid,nslope)) +allocate(active(ngrid,nslope),adjustable(ngrid,nslope)) +allocate(clipped_high(ngrid,nslope),clipped_low(ngrid,nslope)) + +! Build physical bounds and geometry weights do i = 1,ngrid do islope = 1,nslope - if (d_h2oice(i,islope) > 0._dp) then ! Condensation - d_balanced(i,islope) = d_h2oice(i,islope)*real(S_target_cond/S_atm_2_h2oice,dp) - else if (d_h2oice(i,islope) < 0._dp) then ! Sublimation - d_target = d_h2oice(i,islope)*real(S_target_subl/S_h2oice_2_atm,dp) - if (abs(d_target*dt) <= h2o_ice(i,islope)) then ! Enough ice to sublimate - d_balanced(i,islope) = d_target - else ! Not enough ice to sublimate - ! Sublimate all the available ice - d_balanced(i,islope) = -h2o_ice(i,islope)/dt - ! If fully depleted, zero tendency for next time - d_h2oice(i,islope) = 0._dp - ! Compute the amount of H2O ice unable to sublimate - S_ghostice = S_ghostice + real(abs(d_target*dt),qp) - real(h2o_ice(i,islope),qp) - end if + flux_min(i,islope) = -h2o_ice(i,islope)/dt + flux_max(i,islope) = huge(1._dp) + if (d_h2oice(i,islope) > 0._dp) then + flux_min(i,islope) = 0._dp + else if (d_h2oice(i,islope) < 0._dp) then + flux_max(i,islope) = 0._dp end if + capacity(i,islope) = flux_max(i,islope) - flux_min(i,islope) end do end do - -! Enforce conservation removing the ghost ice from places where ice condensed -where (d_balanced(:,:) > 0._dp) d_balanced(:,:) = d_balanced(:,:)*real((S_target_cond - S_ghostice)/S_target_cond,dp) - -END SUBROUTINE balance_h2oice_reservoirs +active(:,:) = .true. + +! Choose weighting method +method_used = weight_method +if (present(weight_type)) then + if (weight_type >= WEIGHT_UNIFORM .and. weight_type <= WEIGHT_DIRECTION) then + weight_method = weight_type + method_used = weight_type + else + call print_msg('Unknown H2O flux weighting method. Falling back to WEIGHT_DIRECTION.',LVL_WRN) + method_used = WEIGHT_DIRECTION + end if +end if + +! Compute initial imbalance +B_flux_nonsurfice = 0._qp +B_flux_surfice = 0._qp +do i = 1,ngrid + B_flux_nonsurfice = B_flux_nonsurfice + real(delta_h2o_ads(i) + delta_icetable(i),qp)/real(cell_area(i),qp)/real(dt,qp) + do islope = 1,nslope + B_flux_surfice = B_flux_surfice + real(d_h2oice(i,islope),qp) + end do +end do +B_flux = B_flux_nonsurfice + B_flux_surfice +Delta = -B_flux +flux_scale = max(1._dp,abs(real(B_flux,dp))) ! Scale flux imbalance to get a relevant tolerance for the balancing procedure (e.g., if imbalance is very small, we don't want to require an even smaller imbalance to stop iterating) +tol_flux = max(1.e-12_dp,1.e-10_dp*flux_scale) ! Tolerance for closing the flux balance, scaled to the initial imbalance + +! Iterative correction +do iter = 1,max_iter + ! Stopping criterion + if (abs(Delta) <= real(tol_flux,qp)) exit + + ! Prepare correction + w(:,:) = 0._dp + ! Compute directional capacities + capacity_up(:,:) = flux_max(:,:) - d_h2oice(:,:) + capacity_down(:,:) = d_h2oice(:,:) - flux_min(:,:) + + ! Compute adjustable mask + adjustable(:,:) = active(:,:) + if (Delta > 0._qp) then + where (capacity_up <= eps) adjustable = .false. + else + where (capacity_down <= eps) adjustable = .false. + end if + if (.not. any(adjustable)) exit + + ! Compute weights + select case (method_used) + case (WEIGHT_UNIFORM) ! Flat redistribution across adjustable cells + where (adjustable) w = 1._dp + case (WEIGHT_ABSFLUX) ! Prioritize cells already carrying strong fluxes + where (adjustable) w = abs(d_h2oice) + case (WEIGHT_CAPACITY) ! Prioritize cells with large admissible adjustment room + where (adjustable) w = capacity + case (WEIGHT_COMBINED) ! Mix current activity (abs flux) and admissible room + where (adjustable) w = capacity*max(abs(d_h2oice),eps) + case (WEIGHT_RESERVOIR) ! Favor points with larger local ice reservoirs + where (adjustable) w = max(h2o_ice,eps) + case (WEIGHT_DIRECTION) ! Follow directional room needed by the sign of Delta + if (Delta > 0._qp) then + where (adjustable) w = max(capacity_up,eps) + else + where (adjustable) w = max(capacity_down,eps) + end if + case default + where (adjustable) w = 1._dp + end select + + wsum = sum(real(w,qp),mask = adjustable) + if (wsum <= 0._qp) exit ! Fallback all weights are zero + + ! Compute correction + d_corr(:,:) = 0._dp + where (adjustable) d_corr = real(Delta*real(w,qp)/wsum,dp) + flux_new(:,:) = d_h2oice(:,:) + d_corr(:,:) + + ! Apply bounds (clipping) + clipped_high(:,:) = .false. + clipped_low(:,:) = .false. + where (adjustable .and. flux_new > flux_max) + d_corr = flux_max - d_h2oice + d_h2oice = flux_max + clipped_high = .true. + end where + where (adjustable .and. flux_new < flux_min) + d_corr = flux_min - d_h2oice + d_h2oice = flux_min + clipped_low = .true. + end where + where (adjustable .and. .not. (clipped_high .or. clipped_low)) + d_h2oice = flux_new + end where + where (clipped_high .or. clipped_low) + active = .false. + end where + + ! Update imbalance + Delta_used = sum(real(d_corr,qp)) + Delta = Delta - Delta_used + + if (abs(Delta_used) < tiny_corr) exit +end do + +! Diagnostics +if (abs(Delta) > real(tol_flux,qp)) then + call print_msg('Unable to close H2O flux balance on surface-ice reservoirs.',LVL_WRN) + call print_msg('Weighting method used = '//int2str(method_used),LVL_WRN) + call print_msg('Residual imbalance [kg/m2/y] = '//real2str(Delta),LVL_WRN) + call print_msg('Initial imbalance [kg/m2/y] = '//real2str(B_flux),LVL_WRN) + call print_msg('Initial surface ice flux [kg/m2/y] = '//real2str(B_flux_surfice),LVL_WRN) + call print_msg('Initial non-surface ice flux [kg/m2/y] = '//real2str(B_flux_nonsurfice),LVL_WRN) + stopcrit%h2o_flux_balance_unclosed = .true. +end if + +! Cleanup +deallocate(w,d_corr,flux_new,flux_min,flux_max) +deallocate(capacity,capacity_up,capacity_down) +deallocate(active,adjustable,clipped_high,clipped_low) + +END SUBROUTINE balance_h2o_fluxes !=======================================================================