ice_frazil_uptake_multicat_impl Subroutine

   pure subroutine ice_frazil_uptake_multicat_impl(hTr_T, hTr_S, h_layer, wet_mask, eos, &
                                                   frazil_heat, part_size, m_ice, m_snow, &
                                                   enth_ice, sal_ice, &
                                                   m_frozen_diag, salt_flux_diag, &
                                                   dt_therm, nghost, ncat, nk, nz, nx, ny)
      !! ncat>1 SIS2 ITD-mode frazil spend. Port of SIS2
      !! `SIS_slow_thermo.F90:1121-1146 + 1181-1186` (non-filling mode —
      !! SIS2's default `SIS2_FILLING_FRAZIL=.true.` thin-category fill is
      !! DEFERRED, see module docstring). Per banked cell (same
      !! wet/non-vanished/bank>0 gate as `ice_frazil_uptake_impl`):
      !!   1. k_merge scan (SIS2:1124-1129): first category `c` with
      !!      `part(0) + part(c) > 0.01`; falls back to `k_merge = 1` if
      !!      no category qualifies (SIS2's `k_merge` default).
      !!   2. Open-water annexation (SIS2:1131-1145): if `part(0) > 0`,
      !!      dilute category `k_merge`'s thickness at CONSTANT MASS —
      !!      `m_ice(k_merge)`/`m_snow(k_merge)` scale by
      !!      `part(k_merge)/(part(k_merge)+part(0))`, `part(k_merge)`
      !!      absorbs all of `part(0)`, `part(0)` -> 0. `enth`/`sal` are
      !!      per-MASS intensive — untouched by an area-only dilution.
      !!   3. Per-ice-area spend (SIS2:1181-1186):
      !!      `frazil_col = frazil_heat/part(k_merge)` (J per m² of
      !!      category area; the denominator is > 0 by construction —
      !!      either an occupied category was found, or step 2 just grew
      !!      `part(k_merge)` from the Σpart=1 invariant), then the SAME
      !!      UNCHANGED `ice_frazil_uptake_column` spends it on category
      !!      `k_merge`'s column.
      !!   4. Diags (per cell, part-weighted back to CELL-area units to
      !!      match the ncat==1 diag convention that the couplers and
      !!      `rdb_ice_thermo_driver` consume): `m_frozen_diag =
      !!      part(k_merge)*m_frozen_pt`; `salt_flux_diag =
      !!      part(k_merge)*(m_frozen_pt*s_surf - salt_to_ice_pt)/dt_therm`;
      !!      `frazil_heat` reset to 0 (fully spent). Both diags are
      !!      zeroed UNCONDITIONALLY at loop top, same contract as
      !!      `ice_frazil_uptake_impl`.
      !!
      !! Decl-order: all integer dims declared before the explicit-shape
      !! arrays that use them. Inner `if` gate only, serial `do cat` scan
      !! for k_merge — never a masked `do concurrent` header.
      integer, intent(in) :: nghost, ncat, nk, nz, nx, ny
      real(wp), intent(in) :: hTr_T(nx, ny, nz)
      real(wp), intent(in) :: hTr_S(nx, ny, nz)
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: wet_mask(nx, ny)
      type(eos_t), intent(in) :: eos
      real(wp), intent(inout) :: frazil_heat(nx, ny)
      real(wp), intent(inout) :: part_size(nx, ny, 0:ncat)
      real(wp), intent(inout) :: m_ice(nx, ny, ncat)
      real(wp), intent(inout) :: m_snow(nx, ny, ncat)
      real(wp), intent(inout) :: enth_ice(nx, ny, ncat, nk)
      real(wp), intent(inout) :: sal_ice(nx, ny, ncat, nk)
      real(wp), intent(inout) :: m_frozen_diag(nx, ny)
      real(wp), intent(inout) :: salt_flux_diag(nx, ny)
      real(wp), intent(in) :: dt_therm

      integer :: i, j, i_lo, i_hi, j_lo, j_hi
      integer :: cat, k_merge
      real(wp) :: h, sst, s_surf, tfw, i_part
      real(wp) :: m_snow_pt, m_ice_pt, m_frozen_pt, salt_to_ice_pt, frazil_col
      real(wp) :: enth_ice_col(ICE_NK_MAX), sal_ice_col(ICE_NK_MAX)

      i_lo = nghost + 1
      i_hi = nx - nghost
      j_lo = nghost + 1
      j_hi = ny - nghost

      do concurrent(j=j_lo:j_hi, i=i_lo:i_hi) local(h, sst, s_surf, tfw, i_part, &
                                                    cat, k_merge, &
                                                    m_snow_pt, m_ice_pt, m_frozen_pt, &
                                                    salt_to_ice_pt, frazil_col, &
                                                    enth_ice_col, sal_ice_col)
         m_frozen_diag(i, j) = 0.0_wp
         salt_flux_diag(i, j) = 0.0_wp
         h = h_layer(i, j, nz)
         if (wet_mask(i, j) > 0.5_wp .and. h > H_VANISHED .and. frazil_heat(i, j) > 0.0_wp) then
            sst = hTr_T(i, j, nz)/h
            s_surf = hTr_S(i, j, nz)/h
            tfw = eos_freezing_point(eos, s_surf, 0.0_wp)

            ! ---- 1. k_merge scan (SIS2:1124-1129) ----
            k_merge = 1
            do cat = 1, ncat
               if (part_size(i, j, 0) + part_size(i, j, cat) > 0.01_wp) then
                  k_merge = cat
                  exit
               end if
            end do

            ! ---- 2. Open-water annexation (SIS2:1131-1145) ----
            if (part_size(i, j, 0) > 0.0_wp) then
               i_part = 1.0_wp/(part_size(i, j, k_merge) + part_size(i, j, 0))
               m_snow(i, j, k_merge) = (m_snow(i, j, k_merge)*part_size(i, j, k_merge))*i_part
               m_ice(i, j, k_merge) = (m_ice(i, j, k_merge)*part_size(i, j, k_merge))*i_part
               part_size(i, j, k_merge) = part_size(i, j, k_merge) + part_size(i, j, 0)
               part_size(i, j, 0) = 0.0_wp
            end if

            ! ---- 3. Per-ice-area spend (SIS2:1181-1186) ----
            frazil_col = frazil_heat(i, j)/part_size(i, j, k_merge)

            m_snow_pt = m_snow(i, j, k_merge)
            m_ice_pt = m_ice(i, j, k_merge)
            enth_ice_col(1:nk) = enth_ice(i, j, k_merge, 1:nk)
            sal_ice_col(1:nk) = sal_ice(i, j, k_merge, 1:nk)

            call ice_frazil_uptake_column(nk, frazil_col, tfw, sst, s_surf, &
                                          m_snow_pt, m_ice_pt, enth_ice_col(1:nk), &
                                          sal_ice_col(1:nk), m_frozen_pt, salt_to_ice_pt)

            m_ice(i, j, k_merge) = m_ice_pt
            enth_ice(i, j, k_merge, 1:nk) = enth_ice_col(1:nk)
            sal_ice(i, j, k_merge, 1:nk) = sal_ice_col(1:nk)

            ! ---- 4. Diags (per cell) ----
            m_frozen_diag(i, j) = part_size(i, j, k_merge)*m_frozen_pt
            salt_flux_diag(i, j) = part_size(i, j, k_merge) &
                                   *(m_frozen_pt*s_surf - salt_to_ice_pt)/dt_therm
            frazil_heat(i, j) = 0.0_wp
         end if
      end do
   end subroutine ice_frazil_uptake_multicat_impl