compute_total_ke Function

   function compute_total_ke(h_layer, u_face, v_face, areaT, nghost) result(total)
      !! Σ 0.5·h·(u_c²+v_c²)·areaT over PHYSICAL cells (ghosts excluded),
      !! using cell-centred face averages. Direct OpenACC reduction.
      real(wp), intent(in) :: h_layer(:, :, :)
      real(wp), intent(in) :: u_face(:, :, :), v_face(:, :, :)
      real(wp), intent(in) :: areaT(:, :)
      integer, intent(in) :: nghost
      real(wp) :: total
      real(wp) :: acc, uc, vc
      integer :: i, j, k, nx, ny, nz, i_lo, i_hi, j_lo, j_hi
      nx = min(size(h_layer, 1), size(u_face, 1) - 1, size(v_face, 1), size(areaT, 1))
      ny = min(size(h_layer, 2), size(u_face, 2), size(v_face, 2) - 1, size(areaT, 2))
      nz = min(size(h_layer, 3), size(u_face, 3), size(v_face, 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) &
      !$acc&         private(uc, vc) present(h_layer, u_face, v_face, areaT)
      do k = 1, nz
         do j = j_lo, j_hi
            do i = i_lo, i_hi
               uc = 0.5_wp*(u_face(i, j, k) + u_face(i + 1, j, k))
               vc = 0.5_wp*(v_face(i, j, k) + v_face(i, j + 1, k))
               acc = acc + 0.5_wp*h_layer(i, j, k)*(uc*uc + vc*vc)*areaT(i, j)
            end do
         end do
      end do
      total = acc
   end function compute_total_ke