ice_snow_ice_flood Subroutine

   pure subroutine ice_snow_ice_flood(nk, m_lay, enthalpy, salin, rho_ratio, snow_to_ice)
      !! Archimedes freeboard snow-ice flooding — the FINAL substantive
      !! block of `ice_resize_SIS2` (SIS2_ice_thm.F90:1303-1320; PR 27).
      !! Standard closure: Leppäranta (1983), *A growth model for black
      !! ice, snow ice and snow thickness in subarctic basins*, Nordic
      !! Hydrology 14, 59-70; Fichefet & Morales Maqueda (1997), JGR 102,
      !! 12609-12646 §2.3.
      !!
      !! The column floats, displacing its own mass of seawater. If the
      !! ice alone cannot support the snow load (`m_submerged =
      !! (m_i+m_s)*rho_ratio > m_i`, `rho_ratio = ICE_RHO_ICE/
      !! ICE_RHO_OCEAN`), the snow-ice interface is below the waterline
      !! ("flooded"): convert `snow_to_ice = min(m_submerged - m_i,
      !! m_lay(0))` kg/m^2 of snow into the TOP ice layer (local index 1
      !! — TOP-DOWN column, TRAP #2 in `rdb_ice_column`'s module
      !! docstring; the caller flips this to the state's bottom-up
      !! `enth_ice(...,nk)` at the scatter boundary). One non-iterative
      !! step suffices: the conversion is 1:1 in mass, so `m_i+m_lay(0)`
      !! is invariant and `m_submerged` does not move, so the interface
      !! lands EXACTLY on the waterline (see PLAN_PR27_snow_ice_
      !! flooding.md §3 for the algebraic proof + a worked golden).
      !!
      !! Roundabout has no ponds (`m_pond` is a hardwired dead local in
      !! `rdb_ice_column`), so the `min(...)` clamp is PROVABLY DEAD
      !! CODE here: it binds iff `m_lay(0) < -m_i`, impossible for
      !! non-negative masses (see the plan §9.2 sweep test). Kept for
      !! SIS2 parity — ponds would resurrect it.
      !!
      !! Mass-, enthalpy- and salt-conserving WITHIN THE COLUMN, exactly
      !! (`(m_0-s)*E_0 + (m_1+s)*[(m_1*E_1+s*E_0)/(m_1+s)] = m_0*E_0 +
      !! m_1*E_1`, identity; the salinity dilution `S_1*m_1/(m_1+s)` IS
      !! the mass-weighted mix since `S_snow == 0`, `salin(0)`).
      !! `enthalpy(1)` and `salin(1)` both DECREASE — snow's cold, fresh
      !! mass lands in the ice.
      !!
      !! **Known physical incompleteness, inherited from SIS2 (not fixed
      !! here — true seawater flooding is a materially larger, separate
      !! closure, see the plan's §12/§14 Q1).** Real flooding draws
      !! SEAWATER into the pore space, which then refreezes, releasing
      !! latent heat and rejecting brine. SIS2 instead moves snow mass
      !! with the SNOW's enthalpy and ZERO salinity — snow enthalpy
      !! (<= -L_f ~ -3.34e5 J/kg) is far more negative than near-freezing
      !! seawater, so the new ice is too COLD; `S_snow == 0` means it is
      !! too FRESH. `snow_to_ice` is diagnostic-only (SIS2's SN2IC):
      !! flooding exchanges NOTHING with the ocean (SIS2_ice_thm.F90:1273,
      !! "There are no further heat or mass losses or gains by the
      !! ice+snow").
      !$acc routine seq
      integer, intent(in) :: nk
         !! Number of ice layers (declared first — decl-order).
      real(wp), intent(inout) :: m_lay(0:nk)
         !! Layer masses (kg/m²), 0 = snow, 1..nk = ice.
      real(wp), intent(inout) :: enthalpy(0:nk)
         !! Layer specific enthalpies (J/kg).
      real(wp), intent(inout) :: salin(0:nk)
         !! Layer bulk salinities (PSU).
      real(wp), intent(in) :: rho_ratio
         !! `ICE_RHO_ICE/ICE_RHO_OCEAN` (nondim, < 1) — computed by the
         !! caller (`rdb_ice_column`), which already has both density
         !! parameters in scope; keeps this module free of a circular
         !! `use rdb_ice_column` (that module `use`s this one).
      real(wp), intent(out) :: snow_to_ice
         !! Mass converted from snow to the top ice layer this call
         !! (kg/m²), >= 0 — SIS2 `SN2IC`.

      real(wp) :: m_i, m_submerged

      m_i = sum(m_lay(1:nk))
      m_submerged = (m_i + m_lay(0))*rho_ratio

      snow_to_ice = 0.0_wp
      if (m_submerged > m_i) then
         snow_to_ice = min(m_submerged - m_i, m_lay(0))
      end if

      ! Denominator `m_lay(1) + snow_to_ice` is provably > 0 whenever
      ! this branch is taken: entering it forces `m_lay(0) > 0` (since
      ! rho_ratio < 1, m_submerged > m_i needs m_lay(0) > 0), hence
      ! snow_to_ice = min(m_lay(0), ...) > 0 — even if m_lay(1) == 0
      ! (top layer melted out this step). The inner guard below is free
      ! and makes the no-snow no-op provable by inspection; it is NOT
      ! `H_DIV_EPS` armour (CLAUDE.md taxonomy — the branch guard above
      ! already rules out the zero case).
      if (snow_to_ice > 0.0_wp) then
         m_lay(0) = m_lay(0) - snow_to_ice
         enthalpy(1) = (m_lay(1)*enthalpy(1) + snow_to_ice*enthalpy(0))/ &
                       (m_lay(1) + snow_to_ice)
         salin(1) = salin(1)*m_lay(1)/(m_lay(1) + snow_to_ice)
         m_lay(1) = m_lay(1) + snow_to_ice
      end if
   end subroutine ice_snow_ice_flood