Changeset 4174 for trunk


Ignore:
Timestamp:
Apr 8, 2026, 11:36:16 AM (4 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

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/backup.F90

    r4170 r4174  
    7070!=======================================================================
    7171SUBROUTINE 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, &
    72                            icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map,backup_idt)
     72                           icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,h2o_ads_reg,co2_ads_reg,layerings_map,backup_idt)
    7373!-----------------------------------------------------------------------
    7474! NAME
    … …  
    107107real(dp),       dimension(:),     intent(in)           :: ps_avg, ps_dev
    108108real(dp),                         intent(in)           :: ps_avg_glob, ps_avg_glob_ini
    109 real(dp),       dimension(:,:),   intent(in)           :: tsurf_avg, tsurf_dev, icetable_depth, icetable_thickness, h2oice_depth, flux_ssice_avg
     109real(dp),       dimension(:,:),   intent(in)           :: tsurf_avg, tsurf_dev, icetable_depth, icetable_thickness, h2oice_depth
    110110real(dp),       dimension(:,:,:), intent(in)           :: tsoil_avg, tsoil_dev, ice_porefilling, h2o_ads_reg, co2_ads_reg
    111111type(layering), dimension(:,:),   intent(in)           :: layerings_map
    … …  
    115115! LOCAL VARIABLES
    116116! ---------------
    117 real(dp),    dimension(ngrid,nslope)           :: h2o_ice4PCM, co2_ice4PCM, tsurf4PCM, flux_geo4PCM, albedo4PCM, emissivity4PCM, flux_ssice4PCM, h2oice_depth4PCM
     117real(dp),    dimension(ngrid,nslope)           :: h2o_ice4PCM, co2_ice4PCM, tsurf4PCM, flux_geo4PCM, albedo4PCM, emissivity4PCM, h2oice_depth4PCM
    118118real(dp),    dimension(ngrid,nlayer)           :: teta4PCM, air_mass4PCM
    119119real(dp),    dimension(ngrid,nsoil_PCM,nslope) :: tsoil4PCM, inertiesoil4PCM

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

    r4170 r4174  
    956956- Keeping a clear separation between the subsurface ice flux and the surface ice tendency.
    957957- Making sure that subsurface ice depth is well given to the PCM at the end of the PEM (ice table VS layering).
     958
     959== 08/04/2026 == JBC
     960Refactor of H2O flux balancing:
     961    - centralizes flux balancing logic;
     962    - removes intermediate variables, especially related to of H2O mass bookkeeping;
     963    - adds new stopping criterion when H2O flux balance fails8 (sanity check);
     964    - introduces configurable weighting strategies for H2O flux balancing.

  • trunk/LMDZ.COMMON/libf/evolution/layered_deposits.F90

    r4170 r4174  
    12551255! ARGUMENTS
    12561256! ---------
    1257 type(layering),         intent(inout) ::  this
    1258 type(stratum), pointer, intent(inout) ::  current
    1259 real(dp),               intent(inout) ::  d_co2ice, d_h2oice
    1260 logical(k4),            intent(inout) ::  new_str, new_lag
    1261 real(dp),               intent(out)   ::  zshift_surf, zlag
     1257real(dp),               intent(in)    :: d_co2ice, d_h2oice
     1258type(layering),         intent(inout) :: this
     1259type(stratum), pointer, intent(inout) :: current
     1260logical(k4),            intent(inout) :: new_str, new_lag
     1261real(dp),               intent(out)   :: zshift_surf, zlag
    12621262
    12631263! LOCAL VARIABLES

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

    r4170 r4174  
    3333use display,            only: print_ini, print_end, print_msg, LVL_NFO, LVL_WRN
    3434use evolution,          only: n_yr_run, n_yr_sim, ntot_yr_sim, nmax_yr_run, dt, idt, r_plnt2earth_yr, pem_ini_date
    35 use geometry,           only: ngrid, nslope, nday, nsoil_PCM, nsoil, cell_area, total_surface
     35use geometry,           only: ngrid, nslope, nsoil_PCM, nsoil, cell_area, total_surface
    3636use glaciers,           only: h2oice_flow, co2ice_flow, flow_co2glaciers, flow_h2oglaciers
    3737use ice_table,          only: icetable_equilibrium, icetable_dynamic, evolve_ice_table
    … …  
    4848use soil_therm_inertia, only: update_soil_TI
    4949use sorption,           only: do_sorption, compute_totmass_adsorbed, evolve_regolith_adsorption
    50 use stopping_crit,      only: stopFlags, stopping_crit_pressure, stopping_crit_h2o, stopping_crit_h2oice, stopping_crit_co2ice
     50use stopping_crit,      only: stopFlags, stopping_crit_pressure, stopping_crit_h2oice, stopping_crit_co2ice
    5151use surface,            only: emissivity_PCM
    52 use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs
     52use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2o_fluxes
    5353use surf_temp,          only: tsurf_PCM, adapt_tsurf2disappearedice
    5454use tendencies,         only: compute_tendice, evolve_tend_co2ice, evolve_flux_ssice
    … …  
    8080real(dp), dimension(:,:), allocatable :: zshift_surf   ! Elevation shift for the surface [m]
    8181real(dp), dimension(:,:), allocatable :: zlag          ! Newly built lag thickness [m]
    82 ! Tendency-related:
    83 real(dp), dimension(:,:), allocatable :: d_h2oice_new ! Adjusted tendency of perennial H2O ice to keep balance between donor and recipient [kg/m2/y]
    8482! Layering-related:
    8583type(ptrarray), dimension(:,:), allocatable :: current          ! Current active stratum in the layering
    … …  
    9593real(qp) :: totmass_ini                                                ! Initial total CO2 mass [kg]
    9694real(qp) :: totmass_adsh2o                                             ! Current total adsorbed H2O mass [kg]
    97 real(qp) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm   ! Variables to balance H2O ice reservoirs
    9895
    9996! CODE
    … …  
    235232    allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
    236233    if (do_layering) then
    237         allocate(d_h2oice_new(ngrid,nslope))
    238234        h2oice_depth_old(:,:) = h2oice_depth(:,:)
    239235        ! The sublimating tendency coming from subsurface ice is given through the surface ice tendency
    … …  
    241237        do islope = 1,nslope
    242238            do i = 1,ngrid
    243                 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)
     239                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)
    244240                call print_layering(layerings_map(i,islope))
    245241            end do
    … …  
    248244        call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
    249245        ! Balance H2O ice reservoirs
    250         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)
    251         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)
    252         deallocate(d_h2oice_new)
     246        call balance_h2o_fluxes(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,stopcrit)
     247        if (stopcrit%stop_code() > 0) then
     248            deallocate(zshift_surf,zlag)
     249            call print_msg(stopcrit%stop_message(),LVL_NFO)
     250            exit
     251        end if
    253252    else
    254253        zlag(:,:) = 0._dp
    … …  
    376375    if (backup_rate > 0) then
    377376        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, &
    378                                                             icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map,idt)
     377                                                            icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,h2o_ads_reg,co2_ads_reg,layerings_map,idt)
    379378    end if
    380379
    … …  
    408407
    409408call 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, &
    410                      icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map)
     409                     icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,h2o_ads_reg,co2_ads_reg,layerings_map)
    411410
    412411! Update the duration information of the workflow

  • trunk/LMDZ.COMMON/libf/evolution/sorption.F90

    r4170 r4174  
    142142use slopes,     only: subslope_dist, def_slope_mean
    143143use atmosphere, only: ap, bp
    144 use physics,    only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2, A, B
     144use physics,    only: alpha_clap_h2o, beta_clap_h2o, m_h2o, A, B
    145145use maths,      only: pi
    146146

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

    r4170 r4174  
    3131! ---------
    3232type :: stopFlags
    33     logical(k4) :: surface_h2oice_change  = .false.
    34     logical(k4) :: zero_dh2oice          = .false.
    35     logical(k4) :: surface_co2ice_change  = .false.
    36     logical(k4) :: pressure_change        = .false.
    37     logical(k4) :: nmax_yr_run_reached    = .false.
    38     logical(k4) :: nmax_yr_runorb_reached = .false.
    39     logical(k4) :: nmax_yr_sim_reached    = .false.
    40     logical(k4) :: time_limit_reached     = .false.
     33    logical(k4) :: surface_h2oice_change     = .false.
     34    logical(k4) :: h2o_flux_balance_unclosed = .false.
     35    logical(k4) :: surface_co2ice_change     = .false.
     36    logical(k4) :: pressure_change           = .false.
     37    logical(k4) :: nmax_yr_run_reached       = .false.
     38    logical(k4) :: nmax_yr_runorb_reached    = .false.
     39    logical(k4) :: nmax_yr_sim_reached       = .false.
     40    logical(k4) :: time_limit_reached        = .false.
    4141contains
    4242    procedure :: is_any_set
    … …  
    120120! CODE
    121121! ----
    122 stopPEM = this%surface_h2oice_change  .or. &
    123           this%zero_dh2oice          .or. &
    124           this%surface_co2ice_change  .or. &
    125           this%pressure_change        .or. &
    126           this%nmax_yr_run_reached    .or. &
    127           this%nmax_yr_runorb_reached .or. &
    128           this%nmax_yr_sim_reached    .or. &
     122stopPEM = this%surface_h2oice_change     .or. &
     123          this%h2o_flux_balance_unclosed .or. &
     124          this%surface_co2ice_change     .or. &
     125          this%pressure_change           .or. &
     126          this%nmax_yr_run_reached       .or. &
     127          this%nmax_yr_runorb_reached    .or. &
     128          this%nmax_yr_sim_reached       .or. &
    129129          this%time_limit_reached
    130130
    … …  
    160160! ----
    161161code = 0
    162 if     (this%surface_h2oice_change)  then; code = 1
    163 elseif (this%zero_dh2oice)          then; code = 2
    164 elseif (this%surface_co2ice_change)  then; code = 3
    165 elseif (this%pressure_change)        then; code = 4
    166 elseif (this%nmax_yr_run_reached)    then; code = 5
    167 elseif (this%nmax_yr_runorb_reached) then; code = 6
    168 elseif (this%nmax_yr_sim_reached)    then; code = 7
    169 elseif (this%time_limit_reached)     then; code = 8
     162if     (this%surface_h2oice_change)     then; code = 1
     163elseif (this%h2o_flux_balance_unclosed) then; code = 2
     164elseif (this%surface_co2ice_change)     then; code = 3
     165elseif (this%pressure_change)           then; code = 4
     166elseif (this%nmax_yr_run_reached)       then; code = 5
     167elseif (this%nmax_yr_runorb_reached)    then; code = 6
     168elseif (this%nmax_yr_sim_reached)       then; code = 7
     169elseif (this%time_limit_reached)        then; code = 8
    170170end if
    171171
    … …  
    200200! CODE
    201201! ----
    202 if     (this%surface_h2oice_change)  then; msg = "**** STOPPING because surface of H2O ice has changed too much (see message above)."
    203 elseif (this%zero_dh2oice)           then; msg = "**** STOPPING because tendencies on H2O ice = 0 (see message above)."
    204 elseif (this%surface_co2ice_change)  then; msg = "**** STOPPING because surface of CO2 ice has changed too much (see message above)."
    205 elseif (this%pressure_change)        then; msg = "**** STOPPING because surface global pressure has changed too much (see message above)."
    206 elseif (this%nmax_yr_run_reached)    then; msg = "**** STOPPING because maximum number of years is reached."
    207 elseif (this%nmax_yr_runorb_reached) then; msg = "**** STOPPING because maximum number of years due to orbital parameters is reached."
    208 elseif (this%nmax_yr_sim_reached)    then; msg = "**** STOPPING because the number of years to be simulated is reached."
    209 elseif (this%time_limit_reached)     then; msg = "**** STOPPING because the job time limit is going to be reached."
     202if     (this%surface_h2oice_change)     then; msg = "**** STOPPING because surface of H2O ice has changed too much (see message above)."
     203elseif (this%h2o_flux_balance_unclosed) then; msg = "**** STOPPING because H2O flux imbalance cannot be closed under physical bounds (see message above)."
     204elseif (this%surface_co2ice_change)     then; msg = "**** STOPPING because surface of CO2 ice has changed too much (see message above)."
     205elseif (this%pressure_change)           then; msg = "**** STOPPING because surface global pressure has changed too much (see message above)."
     206elseif (this%nmax_yr_run_reached)       then; msg = "**** STOPPING because maximum number of years is reached."
     207elseif (this%nmax_yr_runorb_reached)    then; msg = "**** STOPPING because maximum number of years due to orbital parameters is reached."
     208elseif (this%nmax_yr_sim_reached)       then; msg = "**** STOPPING because the number of years to be simulated is reached."
     209elseif (this%time_limit_reached)        then; msg = "**** STOPPING because the job time limit is going to be reached."
     210else
     211    msg = "**** STOPPING for unknown reason (this should not happen!)."
    210212end if
    211213
    … …  
    388390!=======================================================================
    389391
    390 !=======================================================================
    391 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)
    392 !-----------------------------------------------------------------------
    393 ! NAME
    394 !     stopping_crit_h2o
    395 !
    396 ! DESCRIPTION
    397 !     Check if PEM oscillates infinitely with H2O (no net exchange).
    398 !
    399 ! AUTHORS & DATE
    400 !     L. Lange, 2024
    401 !     JB Clement, 2025
    402 !
    403 ! NOTES
    404 !
    405 !-----------------------------------------------------------------------
    406 
    407 ! DEPENDENCIES
    408 ! ------------
    409 use geometry,  only: ngrid, nslope, cell_area
    410 use evolution, only: dt
    411 use slopes,    only: subslope_dist, def_slope_mean
    412 use maths,     only: pi
    413 use display,   only: print_msg, LVL_WRN
    414 use utility,   only: real2str
    415 
    416 ! DECLARATION
    417 ! -----------
    418 implicit none
    419 
    420 ! ARGUMENTS
    421 ! ---------
    422 real(dp), dimension(:),   intent(in)    :: delta_h2o_ads            ! Mass of H2O adsorbed/desorbded in the soil (kg/m^2)
    423 real(dp), dimension(:),   intent(in)    :: delta_icetable ! Mass of H2O that condensed/sublimated at the ice table
    424 real(dp), dimension(:,:), intent(in)    :: h2o_ice, d_h2oice
    425 real(qp),                 intent(out)   :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm
    426 type(stopFlags),          intent(inout) :: stopcrit
    427 
    428 ! LOCAL VARIABLES
    429 ! ---------------
    430 integer(di) :: i, islope
    431 
    432 ! CODE
    433 ! ----
    434 ! Check to escape the subroutine (not relevant in 1D or if the PEM should stop already)
    435 if (ngrid == 1 .or. stopcrit%stop_code() > 0) then
    436     S_atm_2_h2o = 1._qp
    437     S_h2o_2_atm = S_atm_2_h2o
    438     S_atm_2_h2oice = 1._qp
    439     S_h2oice_2_atm = S_atm_2_h2oice
    440     return
    441 end if
    442 
    443 ! We compute the amount of water going out from the atmosphere (S_atm_2_h2o) and going into the atmophere (S_h2o_2_atm)
    444 S_atm_2_h2o = 0._qp
    445 S_h2o_2_atm = 0._qp
    446 S_atm_2_h2oice = 0._qp
    447 S_h2oice_2_atm = 0._qp
    448 do i = 1,ngrid
    449     if (delta_h2o_ads(i) > 0._dp) then
    450         S_atm_2_h2o = S_atm_2_h2o + delta_h2o_ads(i)*cell_area(i)
    451     else
    452         S_h2o_2_atm = S_h2o_2_atm - delta_h2o_ads(i)*cell_area(i)
    453     end if
    454     if (delta_icetable(i) > 0._dp) then
    455         S_atm_2_h2o = S_atm_2_h2o + delta_icetable(i)*cell_area(i)
    456     else
    457         S_h2o_2_atm = S_h2o_2_atm - delta_icetable(i)*cell_area(i)
    458     end if
    459     do islope = 1,nslope
    460         if (d_h2oice(i,islope) > 0._dp) then
    461             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)
    462             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)
    463         else if (d_h2oice(i,islope) < 0._dp .and. h2o_ice(i,islope) > 0._dp) then
    464             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)
    465             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)
    466         end if
    467     end do
    468 end do
    469 
    470 ! Since relative atmospheric water is kept constant, we need to equate condensing reservoirs to the sublimating ones.
    471 ! It is not possible if one of them is missing.
    472 if ((abs(S_atm_2_h2oice) < minieps_qp .or. abs(S_h2oice_2_atm) < minieps_qp)) then
    473     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)
    474     call print_msg("Amount of condensing ice  [kg] = "//real2str(S_atm_2_h2oice),LVL_WRN)
    475     call print_msg("Amount of sublimating ice [kg] = "//real2str(S_h2oice_2_atm),LVL_WRN)
    476     stopcrit%zero_dh2oice = .true.
    477 end if
    478 
    479 END SUBROUTINE stopping_crit_h2o
    480 !=======================================================================
    481 
    482392END MODULE stopping_crit

  • trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90

    r4170 r4174  
    1616! DEPENDENCIES
    1717! ------------
    18 use numerics, only: dp, qp, di, k4
     18use numerics, only: dp, qp, di, k4, minieps
    1919
    2020! DECLARATION
    … …  
    2424! PARAMETERS
    2525! ----------
     26integer(di),                              parameter :: WEIGHT_UNIFORM   = 1_di ! Equal weight for all adjustable points
     27integer(di),                              parameter :: WEIGHT_ABSFLUX   = 2_di ! Weight by absolute local H2O surface flux
     28integer(di),                              parameter :: WEIGHT_CAPACITY  = 3_di ! Weight by available local adjustment range
     29integer(di),                              parameter :: WEIGHT_COMBINED  = 4_di ! Capacity weighted by current absolute flux
     30integer(di),                              parameter :: WEIGHT_RESERVOIR = 5_di ! Weight by locally available H2O ice mass
     31integer(di),                              parameter :: WEIGHT_DIRECTION = 6_di ! Weight by directional capacity (up/down)
     32integer(di),                              private   :: weight_method = WEIGHT_DIRECTION ! Default method for balancing
    2633real(dp),                                 parameter :: rho_co2ice = 1650._dp ! Density of CO2 ice [kg.m-3]
    2734real(dp),                                 parameter :: rho_h2oice = 920._dp  ! Density of H2O ice [kg.m-3]
    … …  
    338345! DEPENDENCIES
    339346! ------------
    340 use geometry,      only: ngrid, nslope
     347use geometry,      only: ngrid
    341348use evolution,     only: dt
    342 use stopping_crit, only: stopping_crit_h2o, stopFlags
     349use stopping_crit, only: stopFlags
    343350use display,       only: print_msg, LVL_NFO
    344351
    … …  
    354361real(dp), dimension(:,:), intent(out)   :: zshift_surf
    355362
    356 ! LOCAL VARIABLES
    357 ! ---------------
    358 real(qp)                          :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! H2O balance variables
    359 real(dp), dimension(ngrid,nslope) :: d_h2oice_new                                             ! Tendencies computed to keep balance
    360 
    361363! CODE
    362364! ----
    … …  
    366368    where (d_h2oice(:,:) < 0._dp .and. abs(d_h2oice(:,:)*dt) > h2o_ice(:,:)) ! Not enough ice to sublimate everything
    367369        ! We sublimate what we can
    368         d_h2oice_new(:,:) = h2o_ice(:,:)/dt
    369         ! It means the tendency is zero next time
    370         d_h2oice(:,:) = 0._dp
    371     else where
    372         d_h2oice_new(:,:) = d_h2oice(:,:)
     370        d_h2oice(:,:) = -h2o_ice(:,:)/dt
    373371    end where
    374372else ! In 3D
    375     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)
     373    call balance_h2o_fluxes(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,stopcrit)
    376374    if (stopcrit%stop_code() > 0) return
    377     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)
    378375end if
    379376
    380 h2o_ice(:,:) = h2o_ice(:,:) + d_h2oice_new(:,:)*dt
    381 zshift_surf(:,:) = zshift_surf(:,:) + d_h2oice_new(:,:)*dt/rho_h2oice
     377h2o_ice(:,:) = h2o_ice(:,:) + d_h2oice(:,:)*dt
     378zshift_surf(:,:) = zshift_surf(:,:) + d_h2oice(:,:)*dt/rho_h2oice
    382379
    383380END SUBROUTINE evolve_h2oice
    … …  
    385382
    386383!=======================================================================
    387 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)
    388 !-----------------------------------------------------------------------
    389 ! NAME
    390 !     balance_h2oice_reservoirs
    391 !
    392 ! DESCRIPTION
    393 !     Balance H2O flux from and into atmosphere across reservoirs.
    394 !
    395 ! AUTHORS & DATE
    396 !     JB Clement, 2025
    397 !
    398 ! NOTES
    399 !
    400 !-----------------------------------------------------------------------
    401 
    402 ! DEPENDENCIES
    403 ! ------------
    404 use evolution, only: dt
    405 use geometry,  only: ngrid, nslope
    406 
    407 ! DECLARATION
    408 ! -----------
    409 implicit none
    410 
    411 ! ARGUMENTS
    412 ! ---------
    413 real(qp),                 intent(in)    :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm
    414 real(dp), dimension(:,:), intent(in)    :: h2o_ice
    415 real(dp), dimension(:,:), intent(inout) :: d_h2oice
    416 real(dp), dimension(:,:), intent(out)   :: d_balanced
     384SUBROUTINE balance_h2o_fluxes(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,stopcrit,weight_type)
     385!-----------------------------------------------------------------------
     386! NAME
     387!     balance_h2o_fluxes
     388!
     389! DESCRIPTION
     390!     Balance H2O fluxes from and into atmosphere across reservoirs.
     391!
     392! AUTHORS & DATE
     393!     JB Clement, 04/2025
     394!
     395! NOTES
     396!
     397!-----------------------------------------------------------------------
     398
     399! DEPENDENCIES
     400! ------------
     401use evolution,     only: dt
     402use geometry,      only: ngrid, nslope, cell_area
     403use stopping_crit, only: stopFlags
     404use utility,       only: real2str, int2str
     405use display,       only: print_msg, LVL_WRN
     406
     407! DECLARATION
     408! -----------
     409implicit none
     410
     411! ARGUMENTS
     412! ---------
     413real(dp), dimension(:),   intent(in)           :: delta_h2o_ads, delta_icetable
     414real(dp), dimension(:,:), intent(in)           :: h2o_ice
     415integer(di),              intent(in), optional :: weight_type
     416real(dp), dimension(:,:), intent(inout)        :: d_h2oice
     417type(stopFlags),          intent(inout)        :: stopcrit
    417418
    418419! LOCAL VARIABLES
    419420! ---------------
    420 integer(di) :: i, islope
    421 real(qp)    :: S_corr_subl, S_corr_cond, S_target_subl, S_target_cond, S_ghostice ! Balance variables
    422 real(dp)    :: d_target
    423 
    424 ! CODE
    425 ! ----
    426 ! Compute global targets
    427 S_corr_cond = (S_h2o_2_atm - S_atm_2_h2o)/2._qp ! Correction = deviation to the mean
    428 S_corr_subl = (S_atm_2_h2o - S_h2o_2_atm)/2._qp ! Correction = deviation to the mean
    429 S_target_cond = abs(S_atm_2_h2oice + S_corr_cond)
    430 S_target_subl = abs(S_h2oice_2_atm + S_corr_subl)
    431 
    432 ! Compute balanced tendencies with positivity limiter
    433 d_balanced(:,:) = 0._dp
    434 S_ghostice = 0._qp
     421integer(di), parameter :: max_iter = 50_di ! Maximum number of iterations for the balancing procedure
     422real(dp),    parameter :: eps = 1.e-12_dp ! Small number to prevent division by zero in weights normalization
     423real(dp),    parameter :: tiny_corr = minieps ! Minimum correction to consider that progress is made in the balancing procedure
     424integer(di)                              :: i, islope, iter
     425integer(di)                              :: method_used
     426real(dp)                                 :: flux_scale, tol_flux
     427real(qp)                                 :: B_flux, B_flux_surfice, B_flux_nonsurfice, Delta, Delta_used, wsum
     428real(dp),    dimension(:,:), allocatable :: w, d_corr, flux_new
     429real(dp),    dimension(:,:), allocatable :: flux_min, flux_max
     430real(dp),    dimension(:,:), allocatable :: capacity, capacity_up, capacity_down
     431logical(k4), dimension(:,:), allocatable :: active, adjustable
     432logical(k4), dimension(:,:), allocatable :: clipped_high, clipped_low
     433
     434! CODE
     435! ----
     436! Initialization
     437allocate(w(ngrid,nslope),d_corr(ngrid,nslope),flux_new(ngrid,nslope))
     438allocate(flux_min(ngrid,nslope),flux_max(ngrid,nslope))
     439allocate(capacity(ngrid,nslope),capacity_up(ngrid,nslope),capacity_down(ngrid,nslope))
     440allocate(active(ngrid,nslope),adjustable(ngrid,nslope))
     441allocate(clipped_high(ngrid,nslope),clipped_low(ngrid,nslope))
     442
     443! Build physical bounds and geometry weights
    435444do i = 1,ngrid
    436445    do islope = 1,nslope
    437         if (d_h2oice(i,islope) > 0._dp) then ! Condensation
    438             d_balanced(i,islope) = d_h2oice(i,islope)*real(S_target_cond/S_atm_2_h2oice,dp)
    439         else if (d_h2oice(i,islope) < 0._dp) then ! Sublimation
    440             d_target = d_h2oice(i,islope)*real(S_target_subl/S_h2oice_2_atm,dp)
    441             if (abs(d_target*dt) <= h2o_ice(i,islope)) then ! Enough ice to sublimate
    442                 d_balanced(i,islope) = d_target
    443             else ! Not enough ice to sublimate
    444                 ! Sublimate all the available ice
    445                 d_balanced(i,islope) = -h2o_ice(i,islope)/dt
    446                 ! If fully depleted, zero tendency for next time
    447                 d_h2oice(i,islope) = 0._dp
    448                 ! Compute the amount of H2O ice unable to sublimate
    449                 S_ghostice = S_ghostice + real(abs(d_target*dt),qp) - real(h2o_ice(i,islope),qp)
    450             end if
     446        flux_min(i,islope) = -h2o_ice(i,islope)/dt
     447        flux_max(i,islope) = huge(1._dp)
     448        if (d_h2oice(i,islope) > 0._dp) then
     449            flux_min(i,islope) = 0._dp
     450        else if (d_h2oice(i,islope) < 0._dp) then
     451            flux_max(i,islope) = 0._dp
    451452        end if
     453        capacity(i,islope) = flux_max(i,islope) - flux_min(i,islope)
    452454    end do
    453455end do
    454 
    455 ! Enforce conservation removing the ghost ice from places where ice condensed
    456 where (d_balanced(:,:) > 0._dp) d_balanced(:,:) = d_balanced(:,:)*real((S_target_cond - S_ghostice)/S_target_cond,dp)
    457 
    458 END SUBROUTINE balance_h2oice_reservoirs
     456active(:,:) = .true.
     457
     458! Choose weighting method
     459method_used = weight_method
     460if (present(weight_type)) then
     461    if (weight_type >= WEIGHT_UNIFORM .and. weight_type <= WEIGHT_DIRECTION) then
     462        weight_method = weight_type
     463        method_used = weight_type
     464    else
     465        call print_msg('Unknown H2O flux weighting method. Falling back to WEIGHT_DIRECTION.',LVL_WRN)
     466        method_used = WEIGHT_DIRECTION
     467    end if
     468end if
     469
     470! Compute initial imbalance
     471B_flux_nonsurfice = 0._qp
     472B_flux_surfice = 0._qp
     473do i = 1,ngrid
     474    B_flux_nonsurfice = B_flux_nonsurfice + real(delta_h2o_ads(i) + delta_icetable(i),qp)/real(cell_area(i),qp)/real(dt,qp)
     475    do islope = 1,nslope
     476        B_flux_surfice = B_flux_surfice + real(d_h2oice(i,islope),qp)
     477    end do
     478end do
     479B_flux = B_flux_nonsurfice + B_flux_surfice
     480Delta = -B_flux
     481flux_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)
     482tol_flux = max(1.e-12_dp,1.e-10_dp*flux_scale) ! Tolerance for closing the flux balance, scaled to the initial imbalance
     483
     484! Iterative correction
     485do iter = 1,max_iter
     486    ! Stopping criterion
     487    if (abs(Delta) <= real(tol_flux,qp)) exit
     488
     489    ! Prepare correction
     490    w(:,:) = 0._dp
     491    ! Compute directional capacities
     492    capacity_up(:,:) = flux_max(:,:) - d_h2oice(:,:)
     493    capacity_down(:,:) = d_h2oice(:,:) - flux_min(:,:)
     494
     495    ! Compute adjustable mask
     496    adjustable(:,:) = active(:,:)
     497    if (Delta > 0._qp) then
     498        where (capacity_up <= eps) adjustable = .false.
     499    else
     500        where (capacity_down <= eps) adjustable = .false.
     501    end if
     502    if (.not. any(adjustable)) exit
     503
     504    ! Compute weights
     505    select case (method_used)
     506        case (WEIGHT_UNIFORM) ! Flat redistribution across adjustable cells
     507            where (adjustable) w = 1._dp
     508        case (WEIGHT_ABSFLUX) ! Prioritize cells already carrying strong fluxes
     509            where (adjustable) w = abs(d_h2oice)
     510        case (WEIGHT_CAPACITY) ! Prioritize cells with large admissible adjustment room
     511            where (adjustable) w = capacity
     512        case (WEIGHT_COMBINED) ! Mix current activity (abs flux) and admissible room
     513            where (adjustable) w = capacity*max(abs(d_h2oice),eps)
     514        case (WEIGHT_RESERVOIR) ! Favor points with larger local ice reservoirs
     515            where (adjustable) w = max(h2o_ice,eps)
     516        case (WEIGHT_DIRECTION) ! Follow directional room needed by the sign of Delta
     517            if (Delta > 0._qp) then
     518                where (adjustable) w = max(capacity_up,eps)
     519            else
     520                where (adjustable) w = max(capacity_down,eps)
     521            end if
     522        case default
     523            where (adjustable) w = 1._dp
     524    end select
     525
     526    wsum = sum(real(w,qp),mask = adjustable)
     527    if (wsum <= 0._qp) exit ! Fallback all weights are zero
     528
     529    ! Compute correction
     530    d_corr(:,:) = 0._dp
     531    where (adjustable) d_corr = real(Delta*real(w,qp)/wsum,dp)
     532    flux_new(:,:) = d_h2oice(:,:) + d_corr(:,:)
     533
     534    ! Apply bounds (clipping)
     535    clipped_high(:,:) = .false.
     536    clipped_low(:,:) = .false.
     537    where (adjustable .and. flux_new > flux_max)
     538        d_corr = flux_max - d_h2oice
     539        d_h2oice = flux_max
     540        clipped_high = .true.
     541    end where
     542    where (adjustable .and. flux_new < flux_min)
     543        d_corr = flux_min - d_h2oice
     544        d_h2oice = flux_min
     545        clipped_low = .true.
     546    end where
     547    where (adjustable .and. .not. (clipped_high .or. clipped_low))
     548        d_h2oice = flux_new
     549    end where
     550    where (clipped_high .or. clipped_low)
     551        active = .false.
     552    end where
     553
     554    ! Update imbalance
     555    Delta_used = sum(real(d_corr,qp))
     556    Delta = Delta - Delta_used
     557
     558    if (abs(Delta_used) < tiny_corr) exit
     559end do
     560
     561! Diagnostics
     562if (abs(Delta) > real(tol_flux,qp)) then
     563    call print_msg('Unable to close H2O flux balance on surface-ice reservoirs.',LVL_WRN)
     564    call print_msg('Weighting method used = '//int2str(method_used),LVL_WRN)
     565    call print_msg('Residual imbalance [kg/m2/y] = '//real2str(Delta),LVL_WRN)
     566    call print_msg('Initial imbalance  [kg/m2/y] = '//real2str(B_flux),LVL_WRN)
     567    call print_msg('Initial surface ice flux     [kg/m2/y] = '//real2str(B_flux_surfice),LVL_WRN)
     568    call print_msg('Initial non-surface ice flux [kg/m2/y] = '//real2str(B_flux_nonsurfice),LVL_WRN)
     569    stopcrit%h2o_flux_balance_unclosed = .true.
     570end if
     571
     572! Cleanup
     573deallocate(w,d_corr,flux_new,flux_min,flux_max)
     574deallocate(capacity,capacity_up,capacity_down)
     575deallocate(active,adjustable,clipped_high,clipped_low)
     576
     577END SUBROUTINE balance_h2o_fluxes
    459578!=======================================================================
    460579
Note: See TracChangeset for help on using the changeset viewer.