compute_ice_totals_efp Subroutine

private subroutine compute_ice_totals_efp(wet_T, areaT, part_size, m_ice, ncat, nghost, wet_area_efp, ci_area_efp, hi_area_efp)

EFP twin of compute_ice_totals. 2D-only (no k-slab blocking needed – one “slab” per accumulator), so a single efp_summands_guard call suffices. THREE SEPARATE single-pass reductions (one per accumulator, 7 reduction scalars each) rather than one kernel combining all 21 – ice diagnostics are a status-cadence cold path (SS11.11: “the EFP path costs ~nz times more kernel launches… unmeasurable at status cadence”), so the extra category-sum pass for ci/hi is free.

GPU bug this works around (reproduced on a V100, NVHPC 25.5, validation_examples/ocean/sea_ice_pack/sea_ice_pack.nml, the default &ocean_diag_nml reproducing_sums = .true. path): the original implementation ran ONE !$acc parallel loop reduction(...) combining all 21 scalar accumulators (ew*/ec*/eh*) with ~17 more private scratch scalars AND an un-annotated inner do c = 1, ncat category-gather loop ahead of the efp_decompose_impl call (an !$acc routine seq helper with SIX intent(out) arguments). On device this silently corrupted the reductions: wet_area_efp%poison came back a nonzero, RUN-INVARIANT garbage value on every status line (not a data-dependent one) – efp_to_real then returns NaN for a poisoned total, the console’s g_wet_area > tiny(0.0_wp) guard is FALSE for NaN (IEEE comparisons with NaN are always false), and mean_ci/mean_hi silently kept their pre-set 0.0_wp default – the observed “Ice: conc 0.0000 thick 0.0000” on every GPU status line, even step 0 off a 100%-covered IC. Splitting into three 7-accumulator passes (matching compute_total_h_efp’s per-k-slab kernel, which never showed this defect) fixed the wet-area pass outright (no inner category loop there), but the ci/hi passes – which DO have the inner do c = 1, ncat gather – stayed poisoned until that loop got an explicit !$acc loop seq (see below): an un-annotated serial loop nested in a !$acc parallel loop reduction(...) region, followed by a multi-out !$acc routine seq call that feeds the reduction, is what NVHPC mis-schedules. Both the split AND the explicit loop seq are required; see test_ocean_console_stats_efp’s test_ice_totals_efp_many_accumulators for the regression gate (CPU-portable: the defect is GPU-codegen-specific, but the test pins the VALUES, which must match compute_ice_totals on every backend). Gather logic + extent clamping copied verbatim from compute_ice_totals.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: wet_T(:,:)
real(kind=wp), intent(in) :: areaT(:,:)
real(kind=wp), intent(in) :: part_size(:,:,0:)
real(kind=wp), intent(in) :: m_ice(:,:,:)
integer, intent(in) :: ncat
integer, intent(in) :: nghost
type(efp_t), intent(out) :: wet_area_efp
type(efp_t), intent(out) :: ci_area_efp
type(efp_t), intent(out) :: hi_area_efp

Calls

proc~~compute_ice_totals_efp~~CallsGraph proc~compute_ice_totals_efp compute_ice_totals_efp proc~efp_carry efp_carry proc~compute_ice_totals_efp->proc~efp_carry proc~efp_decompose_impl efp_decompose_impl proc~compute_ice_totals_efp->proc~efp_decompose_impl proc~efp_summands_guard efp_summands_guard proc~compute_ice_totals_efp->proc~efp_summands_guard error error proc~efp_summands_guard->error

Called by

proc~~compute_ice_totals_efp~~CalledByGraph proc~compute_ice_totals_efp compute_ice_totals_efp proc~ocean_console_stats_report ocean_console_stats_report proc~ocean_console_stats_report->proc~compute_ice_totals_efp proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~ocean_console_stats_report proc~driver_run driver_run proc~driver_run->proc~driver_run_ocean

Variables

Type Visibility Attributes Name Initial
integer, private :: c
real(kind=wp), private :: ci
integer(kind=int64), private :: dw1
integer(kind=int64), private :: dw2
integer(kind=int64), private :: dw3
integer(kind=int64), private :: dw4
integer(kind=int64), private :: dw5
integer(kind=int64), private :: dw6
integer(kind=int64), private :: dwp
integer(kind=int64), private :: ec1
integer(kind=int64), private :: ec2
integer(kind=int64), private :: ec3
integer(kind=int64), private :: ec4
integer(kind=int64), private :: ec5
integer(kind=int64), private :: ec6
integer(kind=int64), private :: ecp
integer(kind=int64), private :: eh1
integer(kind=int64), private :: eh2
integer(kind=int64), private :: eh3
integer(kind=int64), private :: eh4
integer(kind=int64), private :: eh5
integer(kind=int64), private :: eh6
integer(kind=int64), private :: ehp
integer(kind=int64), private :: ew1
integer(kind=int64), private :: ew2
integer(kind=int64), private :: ew3
integer(kind=int64), private :: ew4
integer(kind=int64), private :: ew5
integer(kind=int64), private :: ew6
integer(kind=int64), private :: ewp
integer, private :: i
integer, private :: i_hi
integer, private :: i_lo
integer, private :: j
integer, private :: j_hi
integer, private :: j_lo
real(kind=wp), private :: mice
integer, private :: nx
integer, private :: ny
real(kind=real64), private :: val_w

Source Code

   subroutine compute_ice_totals_efp(wet_T, areaT, part_size, m_ice, ncat, nghost, &
                                     wet_area_efp, ci_area_efp, hi_area_efp)
      !! EFP twin of `compute_ice_totals`.  2D-only (no k-slab blocking
      !! needed -- one "slab" per accumulator), so a single
      !! `efp_summands_guard` call suffices.  THREE SEPARATE single-pass
      !! reductions (one per accumulator, 7 reduction scalars each)
      !! rather than one kernel combining all 21 -- ice diagnostics are a
      !! status-cadence cold path (SS11.11: "the EFP path costs ~nz times
      !! more kernel launches... unmeasurable at status cadence"), so the
      !! extra category-sum pass for `ci`/`hi` is free.
      !!
      !! **GPU bug this works around** (reproduced on a V100, NVHPC 25.5,
      !! `validation_examples/ocean/sea_ice_pack/sea_ice_pack.nml`, the
      !! default `&ocean_diag_nml reproducing_sums = .true.` path): the
      !! original implementation ran ONE `!$acc parallel loop
      !! reduction(...)` combining all 21 scalar accumulators
      !! (`ew*`/`ec*`/`eh*`) with ~17 more private scratch scalars AND an
      !! un-annotated inner `do c = 1, ncat` category-gather loop ahead of
      !! the `efp_decompose_impl` call (an `!$acc routine seq` helper with
      !! SIX `intent(out)` arguments). On device this silently corrupted
      !! the reductions: `wet_area_efp%poison` came back a nonzero,
      !! RUN-INVARIANT garbage value on every status line (not a
      !! data-dependent one) -- `efp_to_real` then returns NaN for a
      !! poisoned total, the console's `g_wet_area > tiny(0.0_wp)` guard
      !! is FALSE for NaN (IEEE comparisons with NaN are always false),
      !! and `mean_ci`/`mean_hi` silently kept their pre-set `0.0_wp`
      !! default -- the observed "Ice: conc 0.0000 thick 0.0000" on every
      !! GPU status line, even step 0 off a 100%-covered IC.  Splitting
      !! into three 7-accumulator passes (matching `compute_total_h_efp`'s
      !! per-k-slab kernel, which never showed this defect) fixed the
      !! wet-area pass outright (no inner category loop there), but the
      !! `ci`/`hi` passes -- which DO have the inner `do c = 1, ncat`
      !! gather -- stayed poisoned until that loop got an explicit
      !! `!$acc loop seq` (see below): an un-annotated serial loop nested
      !! in a `!$acc parallel loop reduction(...)` region, followed by a
      !! multi-out `!$acc routine seq` call that feeds the reduction, is
      !! what NVHPC mis-schedules. Both the split AND the explicit `loop
      !! seq` are required; see `test_ocean_console_stats_efp`'s
      !! `test_ice_totals_efp_many_accumulators` for the regression gate
      !! (CPU-portable: the defect is GPU-codegen-specific, but the test
      !! pins the VALUES, which must match `compute_ice_totals` on every
      !! backend).  Gather logic + extent clamping copied verbatim from
      !! `compute_ice_totals`.
      real(wp), intent(in) :: wet_T(:, :), areaT(:, :)
      real(wp), intent(in) :: part_size(:, :, 0:)
      real(wp), intent(in) :: m_ice(:, :, :)
      integer, intent(in) :: ncat, nghost
      type(efp_t), intent(out) :: wet_area_efp, ci_area_efp, hi_area_efp
      integer :: i, j, c, nx, ny, i_lo, i_hi, j_lo, j_hi
      integer(int64) :: ew1, ew2, ew3, ew4, ew5, ew6, ewp
      integer(int64) :: ec1, ec2, ec3, ec4, ec5, ec6, ecp
      integer(int64) :: eh1, eh2, eh3, eh4, eh5, eh6, ehp
      integer(int64) :: dw1, dw2, dw3, dw4, dw5, dw6, dwp
      real(real64) :: val_w
      real(wp) :: ci, mice

      nx = min(size(wet_T, 1), size(areaT, 1), size(m_ice, 1))
      ny = min(size(wet_T, 2), size(areaT, 2), size(m_ice, 2))
      i_lo = nghost + 1
      i_hi = nx - nghost
      j_lo = nghost + 1
      j_hi = ny - nghost
      call efp_summands_guard(i_hi - i_lo + 1, j_hi - j_lo + 1, "compute_ice_totals_efp")

      ! ---- Pass 1: wet area -- Sigma wet_T*areaT. Identical in both
      ! ncat modes (does not touch part_size/m_ice).
      ew1 = 0_int64
      ew2 = 0_int64
      ew3 = 0_int64
      ew4 = 0_int64
      ew5 = 0_int64
      ew6 = 0_int64
      ewp = 0_int64
      !$acc parallel loop collapse(2) reduction(+:ew1,ew2,ew3,ew4,ew5,ew6,ewp) &
      !$acc&    private(val_w, dw1, dw2, dw3, dw4, dw5, dw6, dwp) present(wet_T, areaT)
      do j = j_lo, j_hi
         do i = i_lo, i_hi
            val_w = real(wet_T(i, j)*areaT(i, j), real64)
            call efp_decompose_impl(val_w, dw1, dw2, dw3, dw4, dw5, dw6, dwp)
            ew1 = ew1 + dw1
            ew2 = ew2 + dw2
            ew3 = ew3 + dw3
            ew4 = ew4 + dw4
            ew5 = ew5 + dw5
            ew6 = ew6 + dw6
            ewp = ewp + dwp
         end do
      end do
      wet_area_efp%v = [ew1, ew2, ew3, ew4, ew5, ew6]
      wet_area_efp%poison = ewp
      call efp_carry(wet_area_efp%v)

      ! ---- Pass 2: ice-covered area -- Sigma ci*areaT. `dw1..dwp` reused
      ! as the per-cell decompose scratch (renamed `dc*` would only add
      ! more private-list entries for the same purpose).
      ec1 = 0_int64
      ec2 = 0_int64
      ec3 = 0_int64
      ec4 = 0_int64
      ec5 = 0_int64
      ec6 = 0_int64
      ecp = 0_int64
      if (ncat == 1) then
         !$acc parallel loop collapse(2) reduction(+:ec1,ec2,ec3,ec4,ec5,ec6,ecp) &
         !$acc&    private(ci, val_w, dw1, dw2, dw3, dw4, dw5, dw6, dwp) present(wet_T, areaT, m_ice)
         do j = j_lo, j_hi
            do i = i_lo, i_hi
               ci = 0.0_wp
               if (wet_T(i, j) > 0.5_wp .and. m_ice(i, j, 1) > 0.0_wp) ci = 1.0_wp
               val_w = real(ci*areaT(i, j), real64)
               call efp_decompose_impl(val_w, dw1, dw2, dw3, dw4, dw5, dw6, dwp)
               ec1 = ec1 + dw1
               ec2 = ec2 + dw2
               ec3 = ec3 + dw3
               ec4 = ec4 + dw4
               ec5 = ec5 + dw5
               ec6 = ec6 + dw6
               ecp = ecp + dwp
            end do
         end do
      else
         !$acc parallel loop collapse(2) reduction(+:ec1,ec2,ec3,ec4,ec5,ec6,ecp) &
         !$acc&    private(c, ci, val_w, dw1, dw2, dw3, dw4, dw5, dw6, dwp) &
         !$acc&    present(wet_T, areaT, part_size)
         do j = j_lo, j_hi
            do i = i_lo, i_hi
               ci = 0.0_wp
               if (wet_T(i, j) > 0.5_wp) then
                  !$acc loop seq
                  do c = 1, ncat
                     ci = ci + part_size(i, j, c)
                  end do
                  ci = min(1.0_wp, ci)
               end if
               val_w = real(ci*areaT(i, j), real64)
               call efp_decompose_impl(val_w, dw1, dw2, dw3, dw4, dw5, dw6, dwp)
               ec1 = ec1 + dw1
               ec2 = ec2 + dw2
               ec3 = ec3 + dw3
               ec4 = ec4 + dw4
               ec5 = ec5 + dw5
               ec6 = ec6 + dw6
               ecp = ecp + dwp
            end do
         end do
      end if
      ci_area_efp%v = [ec1, ec2, ec3, ec4, ec5, ec6]
      ci_area_efp%poison = ecp
      call efp_carry(ci_area_efp%v)

      ! ---- Pass 3: ice volume (as area) -- Sigma (mice/ICE_RHO_ICE)*areaT.
      eh1 = 0_int64
      eh2 = 0_int64
      eh3 = 0_int64
      eh4 = 0_int64
      eh5 = 0_int64
      eh6 = 0_int64
      ehp = 0_int64
      if (ncat == 1) then
         !$acc parallel loop collapse(2) reduction(+:eh1,eh2,eh3,eh4,eh5,eh6,ehp) &
         !$acc&    private(mice, val_w, dw1, dw2, dw3, dw4, dw5, dw6, dwp) present(wet_T, areaT, m_ice)
         do j = j_lo, j_hi
            do i = i_lo, i_hi
               mice = 0.0_wp
               if (wet_T(i, j) > 0.5_wp .and. m_ice(i, j, 1) > 0.0_wp) mice = m_ice(i, j, 1)
               val_w = real((mice/ICE_RHO_ICE)*areaT(i, j), real64)
               call efp_decompose_impl(val_w, dw1, dw2, dw3, dw4, dw5, dw6, dwp)
               eh1 = eh1 + dw1
               eh2 = eh2 + dw2
               eh3 = eh3 + dw3
               eh4 = eh4 + dw4
               eh5 = eh5 + dw5
               eh6 = eh6 + dw6
               ehp = ehp + dwp
            end do
         end do
      else
         !$acc parallel loop collapse(2) reduction(+:eh1,eh2,eh3,eh4,eh5,eh6,ehp) &
         !$acc&    private(c, mice, val_w, dw1, dw2, dw3, dw4, dw5, dw6, dwp) &
         !$acc&    present(wet_T, areaT, part_size, m_ice)
         do j = j_lo, j_hi
            do i = i_lo, i_hi
               mice = 0.0_wp
               if (wet_T(i, j) > 0.5_wp) then
                  !$acc loop seq
                  do c = 1, ncat
                     mice = mice + part_size(i, j, c)*m_ice(i, j, c)
                  end do
               end if
               val_w = real((mice/ICE_RHO_ICE)*areaT(i, j), real64)
               call efp_decompose_impl(val_w, dw1, dw2, dw3, dw4, dw5, dw6, dwp)
               eh1 = eh1 + dw1
               eh2 = eh2 + dw2
               eh3 = eh3 + dw3
               eh4 = eh4 + dw4
               eh5 = eh5 + dw5
               eh6 = eh6 + dw6
               ehp = ehp + dwp
            end do
         end do
      end if
      hi_area_efp%v = [eh1, eh2, eh3, eh4, eh5, eh6]
      hi_area_efp%poison = ehp
      call efp_carry(hi_area_efp%v)
   end subroutine compute_ice_totals_efp