ice_frazil_uptake_column Subroutine

   pure subroutine ice_frazil_uptake_column(nk, frazil, tfw, sst, s_surf, m_snow, &
                                            m_ice_tot, enth_ice_bu, sal_ice_bu, &
                                            m_frozen, salt_to_ice)
      !! Per-(cell,category-1) column worker: gather + flip
      !! (bottom-up -> top-down, TRAP #2 discipline), spend the whole
      !! `frazil` bank evenly over the `nk` layers (SIS2:1402), rebalance,
      !! scatter + flip back. Port of `add_frazil_SIS2`
      !! (SIS2_ice_thm.F90:1342-1444), bulk-salinity mode
      !! (`salin_freeze = ICE_BULK_SALINITY`, SIS_slow_thermo.F90:
      !! 1207-1209) — the `ice_rel_salin` mode is NOT ported.
      !!
      !! Snow slot (`m_snow`, index 0 in the local top-down column) is
      !! carried through untouched — the frazil spend loop only ever
      !! touches `k = 1..nk` (ice layers); `enthalpy(0)`/`salin(0)` are
      !! placeholder zeros that `ice_rebalance_layers` neither reads nor
      !! writes.
      !!
      !! Empty-column degeneracy: `m_lay(k) = 0` on entry makes the
      !! mass-weighted mix reduce exactly to
      !! `enthalpy(k) = enth_frazil`, `salin(k) = ICE_BULK_SALINITY` —
      !! SIS2's post-rebalance `mH_ice == 0 => enthalpy_liquid_freeze`
      !! reset (SIS2_ice_thm.F90:1420-1424) is unreachable here (a spend
      !! only runs when `frazil > 0`, which always deposits mass) and is
      !! deliberately NOT ported.
      !$acc routine seq
      integer, intent(in) :: nk
         !! Number of ice layers (declared first — decl-order).
      real(wp), intent(in) :: frazil
         !! Whole per-cell frazil bank to spend this call (J/m² of cell).
      real(wp), intent(in) :: tfw
         !! Seawater freezing temperature at the surface (degC).
      real(wp), intent(in) :: sst
         !! Sea-surface temperature (degC) — feeds the TRAP-#1 liquid
         !! ocean enthalpy.
      real(wp), intent(in) :: s_surf
         !! Sea-surface salinity (PSU) — feeds the TRAP-#1 liquid ocean
         !! enthalpy (unused by the linear formula, kept for call-site
         !! parity with `ice_enthalpy_liquid`).
      real(wp), intent(in) :: m_snow
         !! Snow mass per unit CELL area (kg/m²) — untouched, carried
         !! through only to seed the local column's slot 0.
      real(wp), intent(inout) :: m_ice_tot
         !! Total category-1 ice mass per unit CELL area (kg/m²); in =
         !! prior step, out = post-freeze.
      real(wp), intent(inout) :: enth_ice_bu(nk)
         !! BOTTOM-UP ice specific enthalpies (J/kg) — state order,
         !! `enth_ice_bu(1)` = ice bottom.
      real(wp), intent(inout) :: sal_ice_bu(nk)
         !! BOTTOM-UP ice bulk salinities (PSU) — state order.
      real(wp), intent(out) :: m_frozen
         !! Total new-ice mass formed this call (kg/m² of cell).
      real(wp), intent(out) :: salt_to_ice
         !! Salt content retained by the new ice (kg/m² of cell) —
         !! `m_frozen * ICE_BULK_SALINITY`.

      real(wp) :: m_lay(0:ICE_NK_MAX), enthalpy(0:ICE_NK_MAX), salin(0:ICE_NK_MAX)
      real(wp) :: frazil_per_layer, enth_ocean, min_denth, t_frazil, enth_frazil, m_frazil
      real(wp) :: mtot_ice
      integer :: k

      ! ---- Gather + flip (TRAP #2): bottom-up state -> top-down local ----
      do k = 1, nk
         enthalpy(k) = enth_ice_bu(nk + 1 - k)
         salin(k) = sal_ice_bu(nk + 1 - k)
      end do
      m_lay(0) = m_snow
      enthalpy(0) = 0.0_wp
      salin(0) = 0.0_wp
      do k = 1, nk
         m_lay(k) = m_ice_tot/real(nk, wp)
      end do

      ! ---- Spend the whole bank, evenly split, top-down (SIS2:1402) ----
      frazil_per_layer = frazil/real(nk, wp)
      enth_ocean = ice_enthalpy_liquid(sst, s_surf)   ! TRAP #1
      min_denth = ICE_LAT_FUS*(1.0_wp - ICE_LIQ_LIM)
      m_frozen = 0.0_wp
      salt_to_ice = 0.0_wp
      do k = 1, nk
         t_frazil = min(tfw, ice_t_freeze(salin(k)) - ICE_FRAZIL_T_OFFSET)
         enth_frazil = min(ice_enth_from_ts(t_frazil, salin(k)), enth_ocean - min_denth)
         m_frazil = frazil_per_layer/(enth_ocean - enth_frazil)
         enthalpy(k) = (m_lay(k)*enthalpy(k) + m_frazil*enth_frazil)/(m_lay(k) + m_frazil)
         salin(k) = (m_lay(k)*salin(k) + m_frazil*ICE_BULK_SALINITY)/(m_lay(k) + m_frazil)
         m_lay(k) = m_lay(k) + m_frazil
         m_frozen = m_frozen + m_frazil
         salt_to_ice = salt_to_ice + m_frazil*ICE_BULK_SALINITY
      end do

      call ice_rebalance_layers(nk, m_lay(0:nk), enthalpy(0:nk), salin(0:nk), mtot_ice)

      ! ---- Scatter + flip back: top-down local -> bottom-up state ----
      m_ice_tot = mtot_ice
      do k = 1, nk
         enth_ice_bu(nk + 1 - k) = enthalpy(k)
         sal_ice_bu(nk + 1 - k) = salin(k)
      end do
   end subroutine ice_frazil_uptake_column