check_h_positive_or_die Subroutine

   subroutine check_h_positive_or_die(grid, ms, label, stage, outer_step, check_layers)
      !! `&ocean_isopycnal_nml check_h_positive` guard.  Abort on the FIRST
      !! negative layer thickness, naming the pipeline stage that produced it,
      !! the offending `(i,j,k)`, and that column's full thickness profile.
      !!
      !! Why this exists: a negative `h` surfaces much later and somewhere else
      !! — as a console-stats NaN, or as an hourly diagnostic minimum — by
      !! which point the producing kernel cannot be identified.  Under
      !! `conservative_floor` it should not be reachable at all: continuity is
      !! handed `h_min = 0`, so its `max(h - dt*div, h_min)` clamps at zero,
      !! and the conservative borrow only redistributes WITHIN a column.  A
      !! trip here therefore localises a real defect rather than a tuning
      !! problem.
      !!
      !! Cheap on the healthy path: one device-side min-reduction, no H<-D
      !! copy.  The host walk + column dump run only when already aborting.
      type(hgrid_t), intent(in) :: grid
      type(multilayer_state_t), intent(inout) :: ms
      character(len=*), intent(in) :: label
         !! Pipeline stage that last wrote `h_layer` (e.g. "after continuity").
      integer, intent(in) :: stage, outer_step
      logical, intent(in) :: check_layers
         !! `.true.` => also abort on a negative single LAYER.  Pass `.false.`
         !! between the raw continuity update and the conservative borrow,
         !! where a transiently negative layer is legal; the column-total
         !! check runs unconditionally either way.
      integer :: i, j, k, ig, i0, i1, j0, j1, ib, jb, kb
      real(wp) :: h_min, col_sum, col_min, worst, acc

      ig = grid%nghost
      i0 = ig + 1
      i1 = grid%nx_total - ig
      j0 = ig + 1
      j1 = grid%ny_total - ig

      ! Column TOTAL is the invariant that always holds.  A negative single
      ! LAYER is legal between the raw continuity update and the conservative
      ! borrow (see the call-site comment), but no column may ever reach a
      ! non-positive total — and if one does, `min_thickness_target_column`'s
      ! degenerate branch (`total <= nz*h_floor`) spreads that negative total
      ! uniformly over every layer, poisoning the whole column silently.
      ! Catching the total here names the stage that evacuated the column.
      col_min = huge(1.0_wp)
      !$acc parallel loop collapse(2) reduction(min:col_min) &
      !$acc   private(k, col_sum) present(ms%h_layer)
      do j = j0, j1
         do i = i0, i1
            col_sum = 0.0_wp
            do k = 1, ms%nz_ml
               col_sum = col_sum + ms%h_layer(i, j, k)
            end do
            col_min = min(col_min, col_sum)
         end do
      end do

      h_min = huge(1.0_wp)
      if (check_layers) then
         !$acc parallel loop collapse(3) reduction(min:h_min) present(ms%h_layer)
         do k = 1, ms%nz_ml
            do j = j0, j1
               do i = i0, i1
                  h_min = min(h_min, ms%h_layer(i, j, k))
               end do
            end do
         end do
      end if
      if (h_min >= 0.0_wp .and. col_min > 0.0_wp) return

      ! Abort path.  Pull the COMPONENT array, never the aggregate `ms` — an
      ! aggregate D->H copy overwrites the host allocatable descriptors with
      ! DEVICE addresses and the next host read segfaults.
      !$acc update self(ms%h_layer)

      ! Locate the worst column: if a column total went non-positive that is
      ! the primary failure, so report the most-negative TOTAL.  Otherwise
      ! report the column holding the most-negative single layer.
      ib = i0
      jb = j0
      worst = huge(1.0_wp)
      if (col_min <= 0.0_wp) then
         do j = j0, j1
            do i = i0, i1
               acc = 0.0_wp
               do k = 1, ms%nz_ml
                  acc = acc + ms%h_layer(i, j, k)
               end do
               if (acc < worst) then
                  worst = acc
                  ib = i
                  jb = j
               end if
            end do
         end do
      else
         do k = 1, ms%nz_ml
            do j = j0, j1
               do i = i0, i1
                  if (ms%h_layer(i, j, k) < worst) then
                     worst = ms%h_layer(i, j, k)
                     ib = i
                     jb = j
                  end if
               end do
            end do
         end do
      end if
      ! Thinnest layer within the reported column.
      kb = 1
      do k = 1, ms%nz_ml
         if (ms%h_layer(ib, jb, k) < ms%h_layer(ib, jb, kb)) kb = k
      end do

      write (output_unit, "(a)") repeat("=", 68)
      if (col_min <= 0.0_wp) then
         write (output_unit, '("[h-guard] NON-POSITIVE COLUMN TOTAL (column evacuated)")')
      else
         write (output_unit, '("[h-guard] NEGATIVE LAYER THICKNESS")')
      end if
      write (output_unit, '("[h-guard]   stage      : ", a)') trim(label)
      write (output_unit, '("[h-guard]   outer step : ", i0, "   rk2 stage : ", i0)') &
         outer_step, stage
      write (output_unit, '("[h-guard]   at (i,j,k) : ", i0, ", ", i0, ", ", i0, &
                          &"   (physical i0/j0 = ", i0, "/", i0, ")")') &
         ib, jb, kb, i0, j0
      write (output_unit, '("[h-guard]   thinnest k : ", i0, "   h = ", es13.5)') &
         kb, ms%h_layer(ib, jb, kb)
      write (output_unit, '("[h-guard]   min layer h over interior  : ", es13.5)') h_min
      write (output_unit, '("[h-guard]   min column total over interior : ", es13.5)') col_min
      col_sum = 0.0_wp
      do k = 1, ms%nz_ml
         col_sum = col_sum + ms%h_layer(ib, jb, k)
      end do
      write (output_unit, '("[h-guard]   column sum : ", es13.5)') col_sum
      write (output_unit, '("[h-guard]   column profile (k, h):")')
      do k = 1, ms%nz_ml
         write (output_unit, '("[h-guard]     ", i4, "  ", es14.6)') k, ms%h_layer(ib, jb, k)
      end do
      write (output_unit, "(a)") repeat("=", 68)
      flush (output_unit)
      error stop "h-guard: negative layer thickness (see [h-guard] block above)"
   end subroutine check_h_positive_or_die