h_face_upstream_open_impl Subroutine

   pure subroutine h_face_upstream_open_impl(nx, ny, nz, use_por, h_layer, &
                                             u_face, v_face, open_u, open_v, &
                                             por_u, por_v, dy_cu, dx_cv, &
                                             dy_cu_bt, dx_cv_bt, h_up_x, h_up_y)
      !! OPEN-column upstream face thickness, for `upstream_h_face` under
      !! `&vcoord_nml zfixed_closed_faces`.
      !!
      !! ### What the barotropic transport has to equal
      !!
      !! Everything else on the barotropic path under closed faces already
      !! means "the OPEN column": `derive_bt_from_layers` builds `ubt` as
      !! `Σ_k u_k·h_up,k·open_k / Σ_k h_up,k·open_k` (upstream pick under
      !! this knob), and the renormaliser hands the fast loop's `uhbt` to
      !! the layers with weight `dy_cu·por·open` on the upstream PPM
      !! thickness.  For a depth-uniform open-layer velocity `u` the layers
      !! therefore carry `dy_cu·Σ_k h_up,k·por_k·open_k·u`, and the fast
      !! loop must transport EXACTLY that, or the renormaliser's `du` is
      !! not zero: `u_av = r·u` with `r` = (barotropic face depth)/(open
      !! upstream depth), and the slow tendencies of the next stage are
      !! evaluated on a velocity the barotropic mode never had.
      !!
      !! ### Why the full-column sum is wrong by O(1), not by round-off
      !!
      !! The fast loop transports `h_face_up·ubt·dy_cu_bt`, and
      !! `dy_cu_bt = dy_cu·φ_c` carries the CENTRED open fraction
      !! `φ_c = Σ h_c·por·open / Σ h_c` (`closed_faces_update_bt_widths`).
      !! The full-column upstream sum `H_up` times `φ_c` is the open depth
      !! only when `H_up = Σ h_c` — a flat face.  At a staircase face
      !! between a deep column `H_D` and a shallow one `H_S` (`h_c` of a
      !! layer live on one side only is `≈ h/2`):
      !!
      !! ```
      !! φ_c        = H_S / ((H_D + H_S)/2)
      !! from deep:    r = H_D·φ_c / H_S = 2·H_D/(H_D + H_S)    -> 2
      !! from shallow: r = H_S·φ_c / H_S = 2·H_S/(H_D + H_S)    -> 0
      !! ```
      !!
      !! (prototype: `python_prototypes/bt_upstream_zfixed/`).  Measured on
      !! the 1-degree Southern Ocean (`zfixed_audit/bt_upstream_h_face`):
      !! barotropic velocity 2.7 m/s by step 13 against 0.84 m/s knob-off,
      !! then the `maxvel` clamp, then NaN at step 309; with the open
      !! column it runs the 10 days at `En 5.709E-04` against the knob-off
      !! `5.493E-04`.
      !!
      !! ### What is stored
      !!
      !! The substep multiplies `h_face_up` by the NARROWED width
      !! `dy_cu_bt`, which already carries `φ_c`.  Storing the open sum
      !! `s = Σ_k h_up,k·por_k·open_k` itself would count the open fraction
      !! twice, so the stored value is the full-column EQUIVALENT
      !!
      !! ```
      !! h_face_up = s · dy_cu / dy_cu_bt     ⇒   h_face_up·dy_cu_bt = s·dy_cu
      !! ```
      !!
      !! which makes the fast-loop transport the open-column upstream
      !! transport to round-off, on BOTH substep kernels, without touching
      !! either (the hot nonlinear kernel receives `dy_cu_bt` as its only
      !! width).  Read against the WIDTH the substep will actually use, so
      !! it is exact whichever step `dy_cu_bt` was refreshed at.  A face
      !! with `dy_cu_bt = 0` (every layer closed, or land) transports
      !! nothing whatever is stored, and stores `s` (`= 0` when every layer
      !! is closed).
      !!
      !! Porous barriers (`use_por`) enter `s` the way they enter the
      !! renormaliser's weight; `por_u`/`por_v` are never indexed when
      !! `.false.` (the caller hands over the mask as an inert,
      !! device-present stand-in — see `compute_h_face_upstream`).
      integer, intent(in) :: nx, ny, nz
      logical, intent(in) :: use_por
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: u_face(nx + 1, ny, nz), v_face(nx, ny + 1, nz)
      real(wp), intent(in) :: open_u(nx + 1, ny, nz), open_v(nx, ny + 1, nz)
      real(wp), intent(in) :: por_u(nx + 1, ny, nz), por_v(nx, ny + 1, nz)
      real(wp), intent(in) :: dy_cu(nx + 1, ny), dx_cv(nx, ny + 1)
      real(wp), intent(in) :: dy_cu_bt(nx + 1, ny), dx_cv_bt(nx, ny + 1)
      real(wp), intent(inout) :: h_up_x(nx + 1, ny), h_up_y(nx, ny + 1)

      integer :: i, j, k
      real(wp) :: s, h_k, w

      do concurrent(j=1:ny, i=1:nx + 1) local(k, s, h_k, w)
         s = 0.0_wp
         do k = 1, nz
            if (i == 1) then
               h_k = h_layer(1, j, k)
            else if (i == nx + 1) then
               h_k = h_layer(nx, j, k)
            else if (u_face(i, j, k) >= 0.0_wp) then
               h_k = h_layer(i - 1, j, k)
            else
               h_k = h_layer(i, j, k)
            end if
            w = open_u(i, j, k)
            if (use_por) w = w*por_u(i, j, k)
            s = s + h_k*w
         end do
         if (dy_cu_bt(i, j) > 0.0_wp) then
            h_up_x(i, j) = s*(dy_cu(i, j)/dy_cu_bt(i, j))
         else
            h_up_x(i, j) = s
         end if
      end do

      do concurrent(j=1:ny + 1, i=1:nx) local(k, s, h_k, w)
         s = 0.0_wp
         do k = 1, nz
            if (j == 1) then
               h_k = h_layer(i, 1, k)
            else if (j == ny + 1) then
               h_k = h_layer(i, ny, k)
            else if (v_face(i, j, k) >= 0.0_wp) then
               h_k = h_layer(i, j - 1, k)
            else
               h_k = h_layer(i, j, k)
            end if
            w = open_v(i, j, k)
            if (use_por) w = w*por_v(i, j, k)
            s = s + h_k*w
         end do
         if (dx_cv_bt(i, j) > 0.0_wp) then
            h_up_y(i, j) = s*(dx_cv(i, j)/dx_cv_bt(i, j))
         else
            h_up_y(i, j) = s
         end if
      end do
   end subroutine h_face_upstream_open_impl