min_thickness_target_column Subroutine

   pure subroutine min_thickness_target_column(nz, h_old, h_floor, h_new, grounded)
      !$acc routine seq
      !$omp declare target
      !! Build the floor-only conservative target thickness column.
      !!
      !! `grounded` (out): `.true.` iff any layer is strictly below `floor`.
      !! When `.false.`, `h_new` is an exact copy of `h_old` (no-op) and the
      !! caller must leave the column untouched.
      !!
      !! When grounded, each sub-floor layer is inflated to `floor` and the
      !! required volume is removed from the surplus layers (`h_old > floor`)
      !! in proportion to their surplus, so `Sigma h_new == Sigma h_old` and
      !! every returned `h_new(k) >= floor` (feasible whenever the column can
      !! hold `nz*floor`).  The largest-surplus layer absorbs the summation
      !! round-off so the column total is exact to a single ULP.
      !!
      !! Degenerate fallback: if the whole column cannot hold `nz*floor`
      !! (never reached in a deep isopycnal ocean), the mass is spread
      !! uniformly — still conservative, floor not guaranteed.
      integer, intent(in) :: nz
      real(wp), intent(in) :: h_old(nz)
      real(wp), intent(in) :: h_floor
      real(wp), intent(out) :: h_new(nz)
      logical, intent(out) :: grounded

      integer :: k, k_donor
      real(wp) :: total, surplus_sum, deficit_sum, surplus_k
      real(wp) :: give, sum_new, best_surplus

      grounded = .false.
      total = 0.0_wp
      do k = 1, nz
         h_new(k) = h_old(k)
         total = total + h_old(k)
         if (h_old(k) < h_floor) grounded = .true.
      end do
      if (.not. grounded) return

      ! Degenerate: not enough water to floor every layer -> spread uniformly.
      if (total <= real(nz, wp)*h_floor) then
         do k = 1, nz
            h_new(k) = total/real(nz, wp)
         end do
         return
      end if

      ! Surplus available above the floor and total deficit below it.
      surplus_sum = 0.0_wp
      deficit_sum = 0.0_wp
      do k = 1, nz
         surplus_k = h_old(k) - h_floor
         if (surplus_k > 0.0_wp) then
            surplus_sum = surplus_sum + surplus_k
         else
            deficit_sum = deficit_sum - surplus_k
         end if
      end do
      ! total > nz*floor guarantees surplus_sum >= deficit_sum (feasible).

      ! Inflate sub-floor layers to the floor; draw the deficit from surplus
      ! layers proportionally to their surplus.
      do k = 1, nz
         surplus_k = h_old(k) - h_floor
         if (surplus_k > 0.0_wp) then
            give = (surplus_k/surplus_sum)*deficit_sum
            h_new(k) = h_old(k) - give
         else
            h_new(k) = h_floor
         end if
      end do

      ! Absorb summation round-off into the largest-surplus layer so the
      ! column total is preserved to a single ULP (exact conservation of h).
      k_donor = 1
      best_surplus = -1.0_wp
      sum_new = 0.0_wp
      do k = 1, nz
         sum_new = sum_new + h_new(k)
         surplus_k = h_old(k) - h_floor
         if (surplus_k > best_surplus) then
            best_surplus = surplus_k
            k_donor = k
         end if
      end do
      h_new(k_donor) = h_new(k_donor) + (total - sum_new)
   end subroutine min_thickness_target_column