top_drag_stress_mag_impl Subroutine

   pure subroutine top_drag_stress_mag_impl(stress_top, u_face, v_face, &
                                            h_layer, wet_mask, cover_frac, &
                                            k_top, h_vanished, &
                                            variant, r, c_d, h_floor, &
                                            htbl, bg_vel, tbl_min, rho0, &
                                            nx_u, ny_u, nx_v, ny_v, nx, ny, nz)
      !! Cell-centred magnitude of the top stress (N/m^2), for the
      !! later `ustar_shelf` consumer:
      !!
      !!   quadratic  `|tau| = rho_0 * C_d * |U_tbl_eff|^2`
      !!   linear     `|tau| = rho_0 * r * |U_tbl| * h_tbl`
      !!
      !! Both are `rho_0 * (drag acceleration) * (band thickness)`, so the
      !! two forms are one definition, and `sqrt(|tau|/rho_0)` is the
      !! friction velocity either way.  Cell-centred velocities come from
      !! the ordinary 2-point face averages; no cover interpolation is
      !! needed because `cover_frac` IS cell-centred here.
      integer, intent(in) :: nx_u, ny_u, nx_v, ny_v, nx, ny, nz, variant
      real(wp), intent(in) :: r, c_d, h_floor, htbl, bg_vel, tbl_min, rho0
      real(wp), intent(in) :: u_face(nx_u, ny_u, nz), v_face(nx_v, ny_v, nz)
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: wet_mask(nx, ny), cover_frac(nx, ny)
      integer, intent(in) :: k_top(nx, ny)
         !! `ms%k_top` -- the first LIVE layer, `nz` off a rigid top.
         !! `stress_top` is the one number the boundary-layer schemes
         !! turn into `u_*` under the shelf (through
         !! `ocean_surface_stress_t%stress_shelf`), so a band mean built
         !! from a filler is a wrong `u_*` in BOTH KPP and EPBL.
      real(wp), intent(in) :: h_vanished
         !! `H_VANISHED`; see `top_drag_tendencies_impl`.
      real(wp), intent(out) :: stress_top(nx, ny)

      integer :: i, j, k, kt
      logical :: layer_only, quad
      real(wp) :: cumul_h, h_k, h_in, u_int, v_int, h_in_total
      real(wp) :: u_c, v_c, u_tbl, v_tbl, h_eff, abs_u_eff, spd

      layer_only = (htbl <= 0.0_wp)
      quad = (variant == TDRAG_QUADRATIC)

      do concurrent(j=1:ny, i=1:nx) &
         local(k, kt, cumul_h, h_k, h_in, u_int, v_int, h_in_total, &
               u_c, v_c, u_tbl, v_tbl, h_eff, abs_u_eff, spd)
         kt = k_top(i, j)
         cumul_h = 0.0_wp
         u_int = 0.0_wp
         v_int = 0.0_wp
         h_in_total = 0.0_wp
         do k = kt, 1, -1
            if (layer_only .and. k < kt) exit
            if ((.not. layer_only) .and. cumul_h >= htbl) exit
            h_k = h_layer(i, j, k)
            if (h_k <= h_vanished) exit
            if (layer_only) then
               h_in = h_k
            else
               h_in = max(0.0_wp, min(h_k, htbl - cumul_h))
            end if
            u_c = 0.5_wp*(u_face(i, j, k) + u_face(i + 1, j, k))
            v_c = 0.5_wp*(v_face(i, j, k) + v_face(i, j + 1, k))
            u_int = u_int + u_c*h_in
            v_int = v_int + v_c*h_in
            h_in_total = h_in_total + h_in
            cumul_h = cumul_h + h_k
         end do
         h_eff = max(max(h_in_total, tbl_min), h_floor)
         u_tbl = u_int/h_eff
         v_tbl = v_int/h_eff
         spd = sqrt(u_tbl*u_tbl + v_tbl*v_tbl)
         abs_u_eff = max(bg_vel, spd)
         if (quad) then
            stress_top(i, j) = rho0*c_d*abs_u_eff*abs_u_eff* &
                               wet_mask(i, j)*cover_frac(i, j)
         else
            stress_top(i, j) = rho0*r*spd*h_eff* &
                               wet_mask(i, j)*cover_frac(i, j)
         end if
      end do
   end subroutine top_drag_stress_mag_impl