ice_transport_compress_cell Subroutine

   pure subroutine ice_transport_compress_cell(part_size, m_ice, m_snow, enth_ice, enth_snow, &
                                               sal_ice, mh_lim, ncat, nk, ok)
      !! KEEP IN SYNC with `ice_compress_cell_inline` (the DEVICE production
      !! twin, directly above). Same excess/ratio/compaction algorithm; this
      !! HOST twin exists only as the directly unit-testable seam (small
      !! standalone per-cell arrays, no `!$acc routine seq`), while the
      !! device twin takes full device-present arrays + a scalar `(i,j)`
      !! index. Any change to the compaction logic here MUST be mirrored
      !! there (the fused device path is what production runs).
      !!
      !! SIS2 `compress_ice` (`SIS_transport.F90:898-1100`), no-ridge
      !! default, ponds dropped — single-CELL kernel (public: directly
      !! unit-testable, SPEC §7 test 3).  `part_size(0:ncat)`,
      !! `m_ice`/`m_snow`/`enth_snow` `(ncat)`, `enth_ice`/`sal_ice`
      !! `(ncat, nk)`.  `part_size(0)` may be negative on entry (the
      !! open-water deficit from Phase 3); no-op (returns `ok=.true.`)
      !! when it is already `>= 0`.
      !!
      !! HOST-ONLY test entry point (no `!$acc routine seq` — this is
      !! deliberately NOT called from device kernels; a dummy-`ncat`-sized
      !! automatic array (`mca_c`/`msnow_c` below) is fine on the host but
      !! would violate the fixed-size-locals rule for a device routine).
      !! The PRODUCTION per-cell body is `ice_compress_cell_inline`
      !! (device-safe: full arrays + scalar `(i,j)` index, no per-category
      !! `mca_c`/`msnow_c` snapshot).  The two are CONSERVING +
      !! SIS2-roundoff-equivalent, but NOT bitwise identical across
      !! CHAINED cascades: this host body carries a pre-materialized
      !! `mca_c(ncat)` array (SIS2's own running array), whereas the
      !! device body recomputes `mca_this = part*m_ice` fresh per
      !! iteration — after a c-1 -> c transfer they differ at round-off,
      !! though both conserve `Σ_c part*m` exactly.  On a SINGLE transfer
      !! (the unit-test hand-check inputs) they ARE bit-identical and
      !! reproduce the Python prototype exactly.  Duplicated rather than
      !! shared, per CLAUDE.md "duplicate explicitly".
      !!
      !! Documented divergence (D8): the tracer merge is done INSIDE the
      !! thinnest-first cascade at each transfer site, whereas SIS2 defers
      !! it to `advect_tracers_thicker` AFTER the category k-loop
      !! (`SIS_tracer_advect.F90`) — equivalent to round-off since each
      !! boundary is visited once in a fixed order with running masses
      !! (same argument as the PR-4a `ice_adjust_categories` merge).
      !!
      !! Algorithm (thinnest-first, `c = 1..ncat-1`):
      !!   `excess = -part_size(0)`; `part_size(0) = 0`.
      !!   `mca_c(c) = part_size(c)*m_ice(c)`, `msnow_c(c) =
      !!   part_size(c)*m_snow(c)` (the CAS-absent branch, SIS2 `:986-992`
      !!   — recomputed here, not reused from Phase 2's `mca_ice`; SIS2
      !!   documents the same roundoff freedom, `:986`).
      !!   Per category `c` (while `excess > 0 .and. mca_c(c) > 0`):
      !!     `ratio = m_ice(c)/mh_lim(c+1)`.
      !!     If `part_size(c)*(1-ratio) >= excess`: IN-PLACE compaction —
      !!       `f = part_size(c)/(part_size(c)-excess)`;
      !!       `m_ice(c) *= f`; `m_snow(c) *= f`; `part_size(c) -= excess`;
      !!       `excess = 0`.
      !!     Else: `excess -= part_size(c)*(1-ratio)`; if
      !!       `mca_c(c) > MASS_NEGLECT_ICE_TRANSPORT`: transfer c -> c+1
      !!       (`part_size(c+1) += part_size(c)*ratio`; mass-weighted
      !!       merge of `mca_c`, `msnow_c`, mass-weighted tracer merge of
      !!       `enth_ice`/`sal_ice`/`enth_snow`; three-branch thickness
      !!       underflow guard on `m_ice(c+1)`/`part_size(c+1)`); zero
      !!       category `c`.
      !!   Top category (`c = ncat`): if `excess > 0`, consistency check
      !!   (`part_size(ncat) <= 1 .and. excess > 2*ncat*epsilon` => FATAL
      !!   in SIS2; here `ok = .false.`), else in-place compaction (same
      !!   formula, no transfer branch).
      !!   Cell epilogue: `part_size(0) = max(1 - Sum_c part_size(c), 0)`.
      integer, intent(in) :: ncat, nk
      real(wp), intent(inout) :: part_size(0:ncat)
      real(wp), intent(inout) :: m_ice(ncat)
      real(wp), intent(inout) :: m_snow(ncat)
      real(wp), intent(inout) :: enth_ice(ncat, nk)
      real(wp), intent(inout) :: enth_snow(ncat, 1)
      real(wp), intent(inout) :: sal_ice(ncat, nk)
      real(wp), intent(in) :: mh_lim(ncat + 1)
      logical, intent(out) :: ok

      integer :: c, l
      real(wp) :: excess, ratio, f, part_sum
      real(wp) :: mca_c(ncat), msnow_c(ncat)
      real(wp) :: mca_old, trans, mnew, msnew

      ok = .true.
      if (part_size(0) >= 0.0_wp) return

      excess = -part_size(0)
      part_size(0) = 0.0_wp
      do c = 1, ncat
         mca_c(c) = part_size(c)*m_ice(c)
         msnow_c(c) = part_size(c)*m_snow(c)
      end do

      do c = 1, ncat - 1
         if (excess > 0.0_wp .and. mca_c(c) > 0.0_wp) then
            ratio = m_ice(c)/mh_lim(c + 1)
            if (part_size(c)*(1.0_wp - ratio) >= excess) then
               f = part_size(c)/(part_size(c) - excess)
               m_ice(c) = m_ice(c)*f
               m_snow(c) = m_snow(c)*f
               part_size(c) = part_size(c) - excess
               excess = 0.0_wp
            else
               excess = excess - part_size(c)*(1.0_wp - ratio)
               if (mca_c(c) > MASS_NEGLECT_ICE_TRANSPORT) then
                  part_size(c + 1) = part_size(c + 1) + part_size(c)*ratio
                  mca_old = mca_c(c + 1)
                  trans = mca_c(c)
                  mca_c(c + 1) = mca_c(c + 1) + mca_c(c)
                  if (part_size(c + 1) > MASS_NEGLECT_ICE_TRANSPORT) then
                     m_ice(c + 1) = mca_c(c + 1)/part_size(c + 1)
                  else if (trans > mca_old) then
                     part_size(c + 1) = mca_c(c + 1)/mh_lim(c + 1)
                     m_ice(c + 1) = mh_lim(c + 1)
                  else
                     part_size(c + 1) = mca_c(c + 1)/m_ice(c + 1)
                  end if
                  msnow_c(c + 1) = msnow_c(c + 1) + msnow_c(c)
                  if (part_size(c + 1) > 0.0_wp) then
                     m_snow(c + 1) = msnow_c(c + 1)/part_size(c + 1)
                  else
                     m_snow(c + 1) = 0.0_wp
                  end if
                  mnew = trans + mca_old
                  if (trans > 0.0_wp .and. mnew > 0.0_wp) then
                     do l = 1, nk
                        enth_ice(c + 1, l) = (trans*enth_ice(c, l) + mca_old*enth_ice(c + 1, l))/mnew
                        sal_ice(c + 1, l) = (trans*sal_ice(c, l) + mca_old*sal_ice(c + 1, l))/mnew
                     end do
                  end if
                  msnew = msnow_c(c + 1)
                  if (msnow_c(c) > 0.0_wp .and. msnew > 0.0_wp) then
                     enth_snow(c + 1, 1) = ((msnew - msnow_c(c))*enth_snow(c + 1, 1) &
                                            + msnow_c(c)*enth_snow(c, 1))/msnew
                  end if
               end if
               mca_c(c) = 0.0_wp
               msnow_c(c) = 0.0_wp
               m_ice(c) = 0.0_wp
               m_snow(c) = 0.0_wp
               part_size(c) = 0.0_wp
            end if
         end if
      end do

      if (excess > 0.0_wp) then
         c = ncat
         ! Same top-category fail-loud as `ice_compress_cell_inline`,
         ! incl. the `part(ncat) <= excess` negative-or-zero-denominator guard
         ! (the `==` tie is the exact zero-denominator divide).
         if ((part_size(c) <= 1.0_wp .and. excess > 2.0_wp*real(ncat, wp)*epsilon(1.0_wp)) &
             .or. part_size(c) <= excess) then
            ok = .false.
            return
         end if
         f = part_size(c)/(part_size(c) - excess)
         m_ice(c) = m_ice(c)*f
         m_snow(c) = m_snow(c)*f
         part_size(c) = part_size(c) - excess
      end if

      part_sum = 0.0_wp
      do c = 1, ncat
         part_sum = part_sum + part_size(c)
      end do
      part_size(0) = max(1.0_wp - part_sum, 0.0_wp)
   end subroutine ice_transport_compress_cell