MODULE surf_ice !----------------------------------------------------------------------- ! NAME ! surf_ice ! ! DESCRIPTION ! Surface ice management. ! ! AUTHORS & DATE ! JB Clement, 12/2025 & 04/2026 ! C. Metz, 04/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use numerics, only: dp, qp, di, k4, eps, tol ! 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_COMBINED ! 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 :: threshold_co2ice_cap ! Threshold to consider CO2 ice as infinite reservoir [kg/m2] real(dp), protected :: h2oice_huge_ini ! Initial value for huge reservoirs of H2O ice [kg/m^2] real(dp), protected :: co2_ice_huge_ini ! Initial value for huge reservoirs of CO2 ice [kg/m^2] contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE set_surf_ice_config(threshold_h2oice_cap_in,h2oice_huge_ini_in,threshold_co2ice_cap_in,co2_ice_huge_ini_in) !----------------------------------------------------------------------- ! NAME ! set_surf_ice_config ! ! DESCRIPTION ! Setter for "surf_ice" configuration parameters. ! ! AUTHORS & DATE ! JB Clement, 02/2026 & 04/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, threshold_co2ice_cap_in, co2_ice_huge_ini_in ! CODE ! ---- threshold_h2oice_cap = threshold_h2oice_cap_in h2oice_huge_ini = h2oice_huge_ini_in threshold_co2ice_cap = threshold_co2ice_cap_in co2_ice_huge_ini = co2_ice_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) call print_msg('threshold_co2ice_cap = '//real2str(threshold_co2ice_cap),LVL_NFO) call print_msg('co2_ice_huge_ini = '//real2str(co2_ice_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) if (threshold_co2ice_cap < 0._dp) call stop_clean(__FILE__,__LINE__,'''threshold_co2ice_cap'' must be positive or zero!',1) if (co2_ice_huge_ini < 0._dp) call stop_clean(__FILE__,__LINE__,'''co2_ice_huge_ini'' must be positive or zero!',1) END SUBROUTINE set_surf_ice_config !======================================================================= !======================================================================= SUBROUTINE build4PCM_perice(h2o_ice,co2_ice,h2o_ice4PCM,co2_ice4PCM) !----------------------------------------------------------------------- ! NAME ! build4PCM_perice ! ! DESCRIPTION ! Reconstructs perennial ice and frost for the PCM. ! ! AUTHORS & DATE ! JB Clement, 12/2025 & 04/2026 ! C. Metz, 04/2026 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use geometry, only: ngrid, nslope use frost, only: h2o_frost4PCM, co2_frost4PCM use display, only: print_msg, LVL_NFO ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:), intent(inout) :: h2o_ice, co2_ice real(dp), dimension(:,:), intent(out) :: h2o_ice4PCM, co2_ice4PCM ! Ice for PCM ! LOCAL VARIABLES ! --------------- integer(di) :: i, islope ! CODE ! ---- call print_msg('> Building surface ice/frost for the PCM',LVL_NFO) do i = 1,ngrid ! cmetz we track the real mass of the ice in the PEM, not the PCM ! So we give the PCM an infinite value of ice wherever we have perennial ice, and just give it the seasonal frost everywhere else do islope = 1,nslope if (h2o_ice(i,islope) > threshold_h2oice_cap) then h2o_ice4PCM(i,islope) = h2oice_huge_ini else h2o_ice4PCM(i,islope) = h2o_frost4PCM(i,islope) + h2o_ice(i,islope) endif if (co2_ice(i,islope) > threshold_co2ice_cap) then co2_ice4PCM(i,islope) = co2_ice_huge_ini else co2_ice4PCM(i,islope) = co2_frost4PCM(i,islope) + co2_ice(i,islope) endif enddo 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(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(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(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(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 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) :: 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 integer(di) :: i, islope, iter integer(di) :: method_used real(dp) :: flux_scale, tol_flux real(qp) :: B_flux, 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 = 0._qp do i = 1,ngrid do islope = 1,nslope B_flux = B_flux + real(d_h2oice(i,islope),qp) end do end do 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 ! Test if the balancing procedure progresses if (abs(Delta_used) < tol) 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) 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