ice_compress_cell_inline Subroutine

   pure subroutine ice_compress_cell_inline(part_size, m_ice, m_snow, enth_ice, enth_snow, &
                                            sal_ice, mh_lim, i, j, ncat, nk, nx, ny, ok)
      !$acc routine seq
      !! KEEP IN SYNC with `ice_transport_compress_cell` (the HOST tested
      !! seam twin, directly below). Same excess/ratio/compaction algorithm;
      !! this twin exists only for a different argument shape — full
      !! device-present arrays + scalar `(i,j)` indices for the fused device
      !! path (fixed-size device locals), vs the twin's small standalone
      !! per-cell arrays for the unit test. Any change to the compaction
      !! logic here MUST be mirrored there (the test only drives the twin).
      !!
      !! Device-safe per-cell compress body: takes the FULL state arrays
      !! plus scalar `(i,j)` indices (no array-section slicing — every
      !! access below is a direct `(i,j,c[,l])` element read/write, so
      !! this is safe to call from `!$acc parallel loop`/`do concurrent`
      !! with the full arrays already device-present).  Same algorithm as
      !! `ice_transport_compress_cell` (see that routine's docstring for
      !! the algorithm narrative); duplicated rather than shared because
      !! the two need different argument shapes (device-safe
      !! full-array-plus-index here; small standalone per-cell arrays
      !! there for the unit test) — CLAUDE.md "duplicate explicitly" (no
      !! include-style body sharing).  Conservation is exact and the two
      !! are SIS2-roundoff-equivalent, but NOT bitwise identical across
      !! CHAINED cascades: this device body recomputes `mca_this =
      !! part_size(c)*m_ice(c)` fresh per iteration, so after a c-1 -> c
      !! transfer (which just re-derived `m_ice(c)` by division) its
      !! `mca_this` differs at round-off from the standalone's running
      !! `mca_c(c)` array (and from SIS2's own running array) — the
      !! transferred/removed masses still telescope to the same total, so
      !! `Σ_c part*m` is conserved regardless.
      integer, intent(in) :: i, j, ncat, nk, 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) :: enth_snow(nx, ny, ncat, 1)
      real(wp), intent(inout) :: sal_ice(nx, ny, 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_this, msnow_this, mca_next, msnow_next
      real(wp) :: mca_old, trans, mnew, msnew

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

      excess = -part_size(i, j, 0)
      part_size(i, j, 0) = 0.0_wp

      ! No pre-materialized per-category mca/msnow arrays (ncat is not a
      ! compile-time bound — CLAUDE.md fixed-size-locals rule for
      ! `!$acc routine seq`).  Each iteration recomputes `mca_this`
      ! (category c) and `mca_next` (category c+1) FRESH from
      ! `part_size*m_ice`, BEFORE this iteration's own mutations — since
      ! the cascade only ever reads/writes categories c and c+1 at step
      ! c, and category c+1 has not yet been touched by any earlier
      ! iteration, this reproduces the array-based algorithm exactly
      ! (verified against the Python prototype, both branches, to full
      ! precision).
      do c = 1, ncat - 1
         mca_this = part_size(i, j, c)*m_ice(i, j, c)
         if (excess > 0.0_wp .and. mca_this > 0.0_wp) then
            msnow_this = part_size(i, j, c)*m_snow(i, j, c)
            ratio = m_ice(i, j, c)/mh_lim(c + 1)
            if (part_size(i, j, c)*(1.0_wp - ratio) >= excess) then
               f = part_size(i, j, c)/(part_size(i, j, c) - excess)
               m_ice(i, j, c) = m_ice(i, j, c)*f
               m_snow(i, j, c) = m_snow(i, j, c)*f
               part_size(i, j, c) = part_size(i, j, c) - excess
               excess = 0.0_wp
            else
               excess = excess - part_size(i, j, c)*(1.0_wp - ratio)
               if (mca_this > MASS_NEGLECT_ICE_TRANSPORT) then
                  mca_next = part_size(i, j, c + 1)*m_ice(i, j, c + 1)
                  msnow_next = part_size(i, j, c + 1)*m_snow(i, j, c + 1)
                  part_size(i, j, c + 1) = part_size(i, j, c + 1) + part_size(i, j, c)*ratio
                  mca_old = mca_next
                  trans = mca_this
                  mnew = mca_next + mca_this
                  if (part_size(i, j, c + 1) > MASS_NEGLECT_ICE_TRANSPORT) then
                     m_ice(i, j, c + 1) = mnew/part_size(i, j, c + 1)
                  else if (trans > mca_old) then
                     part_size(i, j, c + 1) = mnew/mh_lim(c + 1)
                     m_ice(i, j, c + 1) = mh_lim(c + 1)
                  else
                     part_size(i, j, c + 1) = mnew/m_ice(i, j, c + 1)
                  end if
                  msnew = msnow_next + msnow_this
                  if (part_size(i, j, c + 1) > 0.0_wp) then
                     m_snow(i, j, c + 1) = msnew/part_size(i, j, c + 1)
                  else
                     m_snow(i, j, c + 1) = 0.0_wp
                  end if
                  if (trans > 0.0_wp .and. mnew > 0.0_wp) then
                     do l = 1, nk
                        enth_ice(i, j, c + 1, l) = (trans*enth_ice(i, j, c, l) &
                                                    + mca_old*enth_ice(i, j, c + 1, l))/mnew
                        sal_ice(i, j, c + 1, l) = (trans*sal_ice(i, j, c, l) &
                                                   + mca_old*sal_ice(i, j, c + 1, l))/mnew
                     end do
                  end if
                  if (msnow_this > 0.0_wp .and. msnew > 0.0_wp) then
                     enth_snow(i, j, c + 1, 1) = ((msnew - msnow_this)*enth_snow(i, j, c + 1, 1) &
                                                  + msnow_this*enth_snow(i, j, c, 1))/msnew
                  end if
               end if
               m_ice(i, j, c) = 0.0_wp
               m_snow(i, j, c) = 0.0_wp
               part_size(i, j, c) = 0.0_wp
            end if
         end if
      end do

      if (excess > 0.0_wp) then
         c = ncat
         ! Fail-loud (D7) when the top category cannot absorb the leftover
         ! excess: (a) SIS2's own consistency check (`part(ncat) <= 1 .and.
         ! excess > 2*ncat*eps`), PLUS (b) `part(ncat) <= excess` — the
         ! in-place compaction `f = part/(part-excess)` has a NEGATIVE (or,
         ! at the exact `part(ncat) == excess` tie, ZERO) denominator there
         ! and would flip `m_ice(ncat)` large-negative or divide by zero
         ! (SIS2 has the same hole; guarding it is strictly safer and the
         ! Phase-2 mca-only reduction cannot catch a compress-produced
         ! negative).
         if ((part_size(i, j, c) <= 1.0_wp .and. excess > 2.0_wp*real(ncat, wp)*epsilon(1.0_wp)) &
             .or. part_size(i, j, c) <= excess) then
            ok = .false.
            return
         end if
         f = part_size(i, j, c)/(part_size(i, j, c) - excess)
         m_ice(i, j, c) = m_ice(i, j, c)*f
         m_snow(i, j, c) = m_snow(i, j, c)*f
         part_size(i, j, c) = part_size(i, j, c) - excess
      end if

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