ocean_diag_step Subroutine

   subroutine ocean_diag_step(this, state_handle, dt, t)
      !! Advance every registered variable.  Behaviour by `time_op`:
      !!
      !!   * `INSTANT`: skip until `dt_accum >= dt_out`, then fill once
      !!     (snapshot) and emit.
      !!   * `MEAN` / `INTEGRAL` / `MAX` / `MIN`: fill EVERY step, fold
      !!     the sample into the accumulator.  When `dt_accum >= dt_out`,
      !!     finalise the accumulator into `output_buffer` (see
      !!     `finalise_accumulator`), emit, reset.
      !!
      !! Fill writes into `output_buffer` for LAYER vgrid or into
      !! `layer_buffer` (then remap to `output_buffer`) for Z_FIXED.
      !! The post-remap value is the unit of accumulation — so the
      !! mean of T at fixed z-levels is computed in z-space, not
      !! layer-space (matters when h_layer drifts).
      !!
      !! Time stamping (CF-1.8 convention):
      !!   * INSTANT, MAX, MIN -> emit at `t` (end of window).
      !!   * MEAN, INTEGRAL    -> emit at the centre of the window,
      !!     `t - dt_accum / 2`.  The centre lines up with the
      !!     temporal centroid of the time-weighted average so
      !!     downstream tools (CF-aware analysis) plot it correctly.
      !!
      !! The emitted `[diag]` line's min / max / mean are taken over the
      !! FINITE cells only (`diag_field_stats`) — land columns and
      !! vanished layers carry the NaN missing-data sentinel — and the
      !! line gains a `missing=<excluded>/<total>` suffix whenever any
      !! cell was excluded.
      class(ocean_diag_t), intent(inout) :: this
      class(*), intent(in) :: state_handle
      real(wp), intent(in) :: dt
      real(wp), intent(in) :: t
      integer :: i, n_valid, n_total
      real(wp) :: vmin, vmax, vmean, t_emit
      logical :: fire
      character(len=256) :: line
      character(len=64) :: missing_tag

      if (.not. this%is_init) return
      if (.not. this%enabled) return

      do i = 1, this%nvars
         associate (v => this%vars(i))
            v%dt_accum = v%dt_accum + dt
            fire = (v%dt_accum + 1.0e-9_wp >= v%dt_out)

            if (.not. associated(v%fill)) cycle
            if (.not. v%enabled) cycle

            select case (v%time_op)
            case (DIAG_OP_INSTANT)
               if (.not. fire) cycle
               call fill_and_remap(v, state_handle, this%z_out, this%rho_out, &
                                   this%sigma_out, this%zstar_out)
            case (DIAG_OP_MEAN, DIAG_OP_INTEGRAL, DIAG_OP_MAX, DIAG_OP_MIN)
               call fill_and_remap(v, state_handle, this%z_out, this%rho_out, &
                                   this%sigma_out, this%zstar_out)
               call fold_sample(v, dt)
               if (.not. fire) cycle
               call finalise_accumulator(v)
            case default
               error stop "ocean_diag%step: unknown v%time_op"
            end select

            select case (v%time_op)
            case (DIAG_OP_MEAN, DIAG_OP_INTEGRAL)
               t_emit = t - 0.5_wp*v%dt_accum
            case default
               t_emit = t
            end select

            ! Single device→host pull at cadence fire — the fills, fold,
            ! and finalise above all ran on device, so `output_buffer`
            ! lives on the device.  Pull it once here so both the log-line
            ! reductions and the NetCDF write below see fresh host data.
            !
            ! `if_present` keeps lightweight unit tests (which skip
            ! `ocean_state_enter_data` and run on host) working — when
            ! the buffer hasn't been attached, the directive is a no-op
            ! and the host buffer is already current from the host-side
            ! fills.
            ! Statistics first, where the data already is.  On a GPU build
            ! the buffer is device-resident, so one fused device reduction
            ! returns a handful of scalars instead of dragging the whole
            ! array through three host passes.  Without device residency
            ! (unit tests that skip `enter_data`) the host path is the only
            ! correct one -- a device read of an unmapped buffer under
            ! `-gpu=mem:separate` returns garbage without failing.  Both
            ! paths skip the NaN missing-data sentinel; see
            ! `diag_field_stats`.
            call diag_field_stats(v%output_buffer, this%on_device, &
                                  vmin, vmax, vmean, n_valid, n_total)

            ! The NetCDF write below needs the array itself on the host.
            !$acc update self(v%output_buffer) if_present
            v%fire_count = v%fire_count + 1
            write (line, "(A,F12.2,A,A,A,A,A,ES13.5,A,ES13.5,A,ES13.5)") &
               "[diag] t=", t_emit, " ", trim(v%name), " [", trim(v%units), &
               "]  min=", vmin, "  max=", vmax, "  mean=", vmean
            ! Only appended when cells were actually excluded, so a run with
            ! no masked cells emits the byte-identical line it always did.
            if (n_valid < n_total) then
               write (missing_tag, "(A,I0,A,I0)") "  missing=", n_total - n_valid, &
                  "/", n_total
               line = trim(line)//trim(missing_tag)
            end if
            call logger%info(trim(line))
            if (associated(this%emit_post_fire)) then
               call this%emit_post_fire(this, i, t_emit)
            end if
            v%dt_accum = 0.0_wp
            if (v%time_op /= DIAG_OP_INSTANT) then
               call reset_accumulator(v)
               v%n_accum = 0
            end if
         end associate
      end do
      this%dt_last_eval = t
   end subroutine ocean_diag_step