ocean_accumulate_mass_out Subroutine

   subroutine ocean_accumulate_mass_out(ms, flux_h_layer, areaT, nghost, dt, weight)
      !! Accumulate the net mass (kg) that left the domain this RK stage into
      !! `ms%mass_out`.  `flux_h_layer` is the total horizontal divergence
      !! (`h_layer -= dt·flux_h_layer`), so `-Σ_interior(−flux_h_layer)·areaT`
      !! is the boundary outflux (interior faces telescope — divergence
      !! theorem), and it is the SAME field the thickness update consumes, so
      !! with the RK2 stage weight this closes the mass budget to round-off.
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: flux_h_layer(:, :, :)
      real(wp), intent(in) :: areaT(:, :)
      integer, intent(in) :: nghost
      real(wp), intent(in) :: dt, weight
      integer :: i, j, k, nz, nx, ny, i_lo, i_hi, j_lo, j_hi
      real(wp) :: acc
      integer(int64) :: e1, e2, e3, e4, e5, e6, epoison, d1, d2, d3, d4, d5, d6, dpoison
      integer(int64) :: slab_e(EFP_DIGITS)
      real(real64) :: val, scale, colsum

      nx = min(size(flux_h_layer, 1), size(areaT, 1))
      ny = min(size(flux_h_layer, 2), size(areaT, 2))
      nz = size(flux_h_layer, 3)
      i_lo = nghost + 1
      i_hi = nx - nghost
      j_lo = nghost + 1
      j_hi = ny - nghost
      acc = 0.0_wp
      !$acc parallel loop collapse(3) reduction(+:acc) present(flux_h_layer, areaT)
      do k = 1, nz
         do j = j_lo, j_hi
            do i = i_lo, i_hi
               acc = acc + flux_h_layer(i, j, k)*areaT(i, j)
            end do
         end do
      end do
      ! h_layer -= dt·flux_h_layer ⇒ mass change = −dt·ρ·Σ; outflux (positive
      ! = leaving) is its negative.
      ms%mass_out = ms%mass_out + weight*dt*RHO_WATER*acc
      ms%mass_out_tracked = .true.
      if (ms%mass_out_efp_on) then
         ! Order-invariant twin.  Each COLUMN's contribution is formed in a
         ! fixed k order (the same on every decomposition), decomposed into
         ! fixed-point bins, and the bins added exactly -- one 2-D reduction
         ! per call, so the running sum depends on neither the decomposition
         ! nor the reduction order.
         scale = real(weight*dt*RHO_WATER, real64)
         e1 = 0_int64
         e2 = 0_int64
         e3 = 0_int64
         e4 = 0_int64
         e5 = 0_int64
         e6 = 0_int64
         epoison = 0_int64
         !$acc parallel loop collapse(2) reduction(+:e1,e2,e3,e4,e5,e6,epoison) &
         !$acc&  private(val, colsum, d1, d2, d3, d4, d5, d6, dpoison) present(flux_h_layer, areaT)
         do j = j_lo, j_hi
            do i = i_lo, i_hi
               colsum = 0.0_real64
               !$acc loop seq
               do k = 1, nz
                  colsum = colsum + real(flux_h_layer(i, j, k), real64)
               end do
               val = scale*colsum*real(areaT(i, j), real64)
               call efp_decompose_impl(val, d1, d2, d3, d4, d5, d6, dpoison)
               e1 = e1 + d1
               e2 = e2 + d2
               e3 = e3 + d3
               e4 = e4 + d4
               e5 = e5 + d5
               e6 = e6 + d6
               epoison = epoison + dpoison
            end do
         end do
         slab_e = [e1, e2, e3, e4, e5, e6]
         call efp_carry(slab_e)
         ms%mass_out_efp = ms%mass_out_efp + slab_e
         ms%mass_out_efp_poison = ms%mass_out_efp_poison + epoison
         call efp_carry(ms%mass_out_efp)
      end if
   end subroutine ocean_accumulate_mass_out