MODULE surf_ice !----------------------------------------------------------------------- ! NAME ! surf_ice ! ! DESCRIPTION ! Surface ice management. ! ! AUTHORS & DATE ! JB Clement, 12/2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use numerics, only: dp, qp, di, k4, minieps ! DECLARATION ! ----------- implicit none ! 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] real(dp), protected :: threshold_h2oice_cap ! Threshold to consider H2O ice as infinite reservoir [kg/m2] real(dp), protected :: h2oice_huge_ini ! Initial value for huge reservoirs of H2O ice [kg/m^2] logical(k4), dimension(:), allocatable, protected :: is_h2o_perice_PCM ! Flag for the presence of perennial H2O ice in the PCM at the beginning [kg/m2] real(dp), dimension(:,:), allocatable, protected :: co2_perice_PCM ! Perennial CO2 ice in the PCM at the beginning [kg/m2] contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE ini_surf_ice() !----------------------------------------------------------------------- ! NAME ! ini_surf_ice ! ! DESCRIPTION ! Initialize surface ice arrays. ! ! AUTHORS & DATE ! JB Clement 12/2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope ! DECLARATION ! ----------- implicit none ! CODE ! ---- allocate(is_h2o_perice_PCM(ngrid)) allocate(co2_perice_PCM(ngrid,nslope)) is_h2o_perice_PCM(:) = .false. co2_perice_PCM(:,:) = 0._dp END SUBROUTINE ini_surf_ice !======================================================================= !======================================================================= SUBROUTINE end_surf_ice() !----------------------------------------------------------------------- ! NAME ! end_surf_ice ! ! DESCRIPTION ! Deallocate surface ice arrays. ! ! AUTHORS & DATE ! JB Clement, 12/2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! CODE ! ---- if (allocated(is_h2o_perice_PCM)) deallocate(is_h2o_perice_PCM) if (allocated(co2_perice_PCM)) deallocate(co2_perice_PCM) END SUBROUTINE end_surf_ice !======================================================================= !======================================================================= SUBROUTINE set_surf_ice_config(threshold_h2oice_cap_in,h2oice_huge_ini_in) !----------------------------------------------------------------------- ! NAME ! set_surf_ice_config ! ! DESCRIPTION ! Setter for "surf_ice" configuration parameters. ! ! AUTHORS & DATE ! JB Clement, 02/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use utility, only: real2str use display, only: print_msg, LVL_NFO use stoppage, only: stop_clean ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), intent(in) :: threshold_h2oice_cap_in, h2oice_huge_ini_in ! CODE ! ---- threshold_h2oice_cap = threshold_h2oice_cap_in h2oice_huge_ini = h2oice_huge_ini_in call print_msg('threshold_h2oice_cap = '//real2str(threshold_h2oice_cap),LVL_NFO) call print_msg('h2oice_huge_ini = '//real2str(h2oice_huge_ini),LVL_NFO) if (threshold_h2oice_cap < 0._dp) call stop_clean(__FILE__,__LINE__,'''threshold_h2oice_cap'' must be positive or zero!',1) if (h2oice_huge_ini < 0._dp) call stop_clean(__FILE__,__LINE__,'''h2oice_huge_ini'' must be positive or zero!',1) END SUBROUTINE set_surf_ice_config !======================================================================= !======================================================================= SUBROUTINE set_co2_perice_PCM(co2_perice_PCM_in) !----------------------------------------------------------------------- ! NAME ! set_co2_perice_PCM ! ! DESCRIPTION ! Setter for 'co2_perice_PCM'. ! ! AUTHORS & DATE ! JB Clement, 12/2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:), intent(in) :: co2_perice_PCM_in ! CODE ! ---- co2_perice_PCM(:,:) = co2_perice_PCM_in(:,:) END SUBROUTINE set_co2_perice_PCM !======================================================================= !======================================================================= SUBROUTINE set_is_h2o_perice_PCM(is_h2o_perice_PCM_in) !----------------------------------------------------------------------- ! NAME ! set_is_h2o_perice_PCM ! ! DESCRIPTION ! Setter for 'is_h2o_perice_PCM'. ! ! AUTHORS & DATE ! JB Clement, 12/2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:), intent(in) :: is_h2o_perice_PCM_in ! CODE ! ---- where(is_h2o_perice_PCM_in > 0.5_dp) is_h2o_perice_PCM = .true. else where is_h2o_perice_PCM = .false. end where END SUBROUTINE set_is_h2o_perice_PCM !======================================================================= !======================================================================= SUBROUTINE build4PCM_perice(h2o_ice,co2_ice,is_h2o_perice,h2o_ice4PCM,co2_ice4PCM) !----------------------------------------------------------------------- ! NAME ! build4PCM_perice ! ! DESCRIPTION ! Reconstructs perennial ice and frost for the PCM. ! ! AUTHORS & DATE ! JB Clement, 12/2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid use frost, only: h2o_frost4PCM use slopes, only: subslope_dist, def_slope_mean use maths, only: pi use display, only: print_msg, LVL_NFO ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:), intent(inout) :: h2o_ice, co2_ice logical(k4), dimension(:), intent(out) :: is_h2o_perice ! H2O perennial ice flag real(dp), dimension(:,:), intent(out) :: h2o_ice4PCM, co2_ice4PCM ! Ice for PCM ! LOCAL VARIABLES ! --------------- integer(di) :: i ! CODE ! ---- call print_msg('> Building surface ice/frost for the PCM',LVL_NFO) co2_ice4PCM(:,:) = co2_ice(:,:) h2o_ice4PCM(:,:) = 0._dp ! Because in the Mars PCM, only the variation of perennial H2O ice is monitored, not the absolute quantity do i = 1,ngrid ! Is H2O ice still considered as an infinite reservoir for the PCM? if (sum(h2o_ice(i,:)*subslope_dist(i,:)/cos(pi*def_slope_mean(:)/180._dp)) > threshold_h2oice_cap) then ! There is enough ice to be considered as an infinite reservoir is_h2o_perice(i) = .true. else ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost is_h2o_perice(i) = .false. h2o_frost4PCM(i,:) = h2o_frost4PCM(i,:) + h2o_ice(i,:) h2o_ice(i,:) = 0._dp end if end do END SUBROUTINE build4PCM_perice !======================================================================= !======================================================================= SUBROUTINE evolve_co2ice(co2_ice,d_co2ice,zshift_surf) !----------------------------------------------------------------------- ! NAME ! evolve_co2ice ! ! DESCRIPTION ! Compute the evolution of CO2 ice over one year. ! ! AUTHORS & DATE ! R. Vandemeulebrouck ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope use evolution, only: dt use display, only: print_msg, LVL_NFO ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:), intent(inout) :: co2_ice real(dp), dimension(:,:), intent(inout) :: d_co2ice real(dp), dimension(:,:), intent(out) :: zshift_surf ! LOCAL VARIABLES ! --------------- real(dp), dimension(ngrid,nslope) :: co2_ice_old ! Old density of CO2 ice ! CODE ! ---- ! Evolution of CO2 ice for each physical point call print_msg('> Evolution of CO2 ice',LVL_NFO) co2_ice_old(:,:) = co2_ice(:,:) co2_ice(:,:) = co2_ice(:,:) + d_co2ice(:,:)*dt where (co2_ice < 0._dp) co2_ice(:,:) = 0._dp d_co2ice(:,:) = -co2_ice_old(:,:)/dt end where zshift_surf(:,:) = zshift_surf(:,:) + d_co2ice(:,:)*dt/rho_co2ice END SUBROUTINE evolve_co2ice !======================================================================= !======================================================================= SUBROUTINE evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit) !----------------------------------------------------------------------- ! NAME ! evolve_h2oice ! ! DESCRIPTION ! Compute the evolution of H2O ice over one year. ! ! AUTHORS & DATE ! R. Vandemeulebrouck ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid use evolution, only: dt use stopping_crit, only: stopFlags use display, only: print_msg, LVL_NFO ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:), intent(in) :: delta_h2o_ads, delta_icetable real(dp), dimension(:,:), intent(inout) :: h2o_ice, d_h2oice type(stopFlags), intent(inout) :: stopcrit real(dp), dimension(:,:), intent(out) :: zshift_surf ! CODE ! ---- call print_msg('> Evolution of H2O ice',LVL_NFO) if (ngrid == 1) then ! In 1D 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(:,:) = -h2o_ice(:,:)/dt end where else ! In 3D call balance_h2o_fluxes(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,stopcrit) if (stopcrit%stop_code() > 0) return end if h2o_ice(:,:) = h2o_ice(:,:) + d_h2oice(:,:)*dt zshift_surf(:,:) = zshift_surf(:,:) + d_h2oice(:,:)*dt/rho_h2oice END SUBROUTINE evolve_h2oice !======================================================================= !======================================================================= 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), 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 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 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 !======================================================================= END MODULE surf_ice