bbl_faces_impl Subroutine

   pure subroutine bbl_faces_impl(nu, nv, nx, ny, nz, x_face, u_x, v_y, h, wet, &
                                  conc_t, conc_s, ncx, ncy, ncz, f_corner, use_eos, eos, &
                                  form, cd, hbbl, bg_vel, thick_min, rino_cap, rho0, &
                                  kv_bbl, bbl_thick)
      !! Per-face body of `vdiff_set_viscous_bbl` (see there).  `x_face`
      !! selects the u-faces `(nx+1, ny)` (cells `i-1`, `i`) or the v-faces
      !! `(nx, ny+1)` (cells `j-1`, `j`).
      integer, intent(in) :: nu, nv, nx, ny, nz, ncx, ncy, ncz
      logical, intent(in) :: x_face
      real(wp), intent(in) :: u_x(nx + 1, ny, nz)
      real(wp), intent(in) :: v_y(nx, ny + 1, nz)
      real(wp), intent(in) :: h(nx, ny, nz)
      real(wp), intent(in) :: wet(nx, ny)
      real(wp), intent(in) :: conc_t(ncx, ncy, ncz)
      real(wp), intent(in) :: conc_s(ncx, ncy, ncz)
         !! Cell T/S concentrations; read only when `use_eos` (then
         !! `(nx, ny, nz)`; a `(1,1,1)` placeholder otherwise).
      real(wp), intent(in) :: f_corner(nx + 1, ny + 1)
      logical, intent(in) :: use_eos
      type(eos_t), intent(in) :: eos
      integer, intent(in) :: form
      real(wp), intent(in) :: cd, hbbl, bg_vel, thick_min, rho0
      logical, intent(in) :: rino_cap
      real(wp), intent(inout) :: kv_bbl(nu, nv)
      real(wp), intent(inout) :: bbl_thick(nu, nv)

      real(wp), parameter :: ANGSTROM = 1.0e-10_wp
         !! MOM6 `ANGSTROM` default (m): the thinnest layer counted in the
         !! HBBL average (`hweight < 1.5·Angstrom` is skipped, line 666).
      real(wp), parameter :: EPS_NEG = 1.0e-30_wp
         !! MOM6 `h_neglect` / `dz_neglect` (Boussinesq, m).
      integer :: i, j, k, il, ir, jl, jr, im, ip, jm, jp
      real(wp) :: hav(NZ_STACK_MAX), tv(NZ_STACK_MAX), sv(NZ_STACK_MAX)
      real(wp) :: hl, hr, un, press, cd_sqrt, u2_bg, ustar, ustarsq
      real(wp) :: htot_vel, hwtot, hutot, thtot, shtot, hweight, vt
      real(wp) :: wa, wb, wc, wd, wsum, t_eos, s_eos
      real(wp) :: drdt, drds, htot, oldfn, dfn, dh, c2f, root, thick

      cd_sqrt = sqrt(cd)
      u2_bg = bg_vel*bg_vel
      do concurrent(j=1:nv, i=1:nu) &
         local(il, ir, jl, jr, im, ip, jm, jp, k, hav, tv, sv, hl, hr, un, press, &
               ustar, ustarsq, htot_vel, hwtot, hutot, thtot, shtot, hweight, vt, &
               wa, wb, wc, wd, wsum, t_eos, s_eos, drdt, drds, &
               htot, oldfn, dfn, dh, c2f, root, thick)
         if (x_face) then
            il = max(1, i - 1)
            ir = min(nx, i)
            jl = j
            jr = j
         else
            il = i
            ir = i
            jl = max(1, j - 1)
            jr = min(ny, j)
         end if
         kv_bbl(i, j) = 0.0_wp
         bbl_thick(i, j) = hbbl
         if (wet(il, jl)*wet(ir, jr) > 0.0_wp) then
            ! ---- 1. face thickness, T, S, bottom pressure (lines 491-516, 791-803)
            press = 0.0_wp
            do k = 1, nz
               hl = h(il, jl, k)
               hr = h(ir, jr, k)
               if (x_face) then
                  un = u_x(i, j, k)
               else
                  un = v_y(i, j, k)
               end if
               if (un*(hr - hl) >= 0.0_wp) then
                  hav(k) = 2.0_wp*hl*hr/(hl + hr + EPS_NEG)
               else
                  hav(k) = 0.5_wp*(hl + hr)
               end if
               press = press + rho0*GRAVITY*0.5_wp*(hl + hr)
               tv(k) = 0.0_wp
               sv(k) = 0.0_wp
               if (use_eos) then
                  tv(k) = 0.5_wp*(conc_t(il, jl, k) + conc_t(ir, jr, k))
                  sv(k) = 0.5_wp*(conc_s(il, jl, k) + conc_s(ir, jr, k))
               end if
            end do

            ! ---- 2. u_bbl and u* over the bottom HBBL ----
            htot_vel = 0.0_wp
            hwtot = 0.0_wp
            hutot = 0.0_wp
            thtot = 0.0_wp
            shtot = 0.0_wp
            do k = 1, nz
               if (htot_vel >= hbbl) exit
               hweight = min(hbbl - htot_vel, hav(k))
               if (hweight < 1.5_wp*ANGSTROM + EPS_NEG) cycle
               htot_vel = htot_vel + hav(k)
               hwtot = hwtot + hweight
               if (form == BBL_FORM_QUADRATIC) then
                  ! `set_v_at_u` / `set_u_at_v`: masked, thickness-weighted
                  ! mean of the four transverse faces at this layer.
                  if (x_face) then
                     un = u_x(i, j, k)
                     jm = max(1, j - 1)
                     jp = min(ny, j + 1)
                     wa = max(0.0_wp, h(il, jm, k) + h(il, j, k))*wet(il, jm)*wet(il, j)
                     wb = max(0.0_wp, h(ir, jm, k) + h(ir, j, k))*wet(ir, jm)*wet(ir, j)
                     wc = max(0.0_wp, h(il, j, k) + h(il, jp, k))*wet(il, j)*wet(il, jp)
                     wd = max(0.0_wp, h(ir, j, k) + h(ir, jp, k))*wet(ir, j)*wet(ir, jp)
                     wsum = (wa + wd) + (wb + wc)
                     vt = 0.0_wp
                     if (wsum > 0.0_wp) then
                        vt = ((wa*v_y(il, j, k) + wd*v_y(ir, j + 1, k)) + &
                              (wb*v_y(ir, j, k) + wc*v_y(il, j + 1, k)))/wsum
                     end if
                  else
                     un = v_y(i, j, k)
                     im = max(1, i - 1)
                     ip = min(nx, i + 1)
                     wa = max(0.0_wp, h(im, jl, k) + h(i, jl, k))*wet(im, jl)*wet(i, jl)
                     wb = max(0.0_wp, h(im, jr, k) + h(i, jr, k))*wet(im, jr)*wet(i, jr)
                     wc = max(0.0_wp, h(i, jl, k) + h(ip, jl, k))*wet(i, jl)*wet(ip, jl)
                     wd = max(0.0_wp, h(i, jr, k) + h(ip, jr, k))*wet(i, jr)*wet(ip, jr)
                     wsum = (wa + wd) + (wb + wc)
                     vt = 0.0_wp
                     if (wsum > 0.0_wp) then
                        vt = ((wa*u_x(i, jl, k) + wd*u_x(i + 1, jr, k)) + &
                              (wb*u_x(i, jr, k) + wc*u_x(i + 1, jl, k)))/wsum
                     end if
                  end if
                  hutot = hutot + hweight*sqrt(un*un + vt*vt + u2_bg)
               end if
               thtot = thtot + hweight*tv(k)
               shtot = shtot + hweight*sv(k)
            end do
            if (hwtot <= 0.0_wp .or. form == BBL_FORM_LINEAR) then
               ustar = cd_sqrt*bg_vel
            else
               ustar = cd_sqrt*hutot/hwtot
            end if

            ! ---- 3. stratification-limited thickness h_N ----
            drdt = 0.0_wp
            drds = 0.0_wp
            if (use_eos) then
               t_eos = 0.0_wp
               s_eos = 0.0_wp
               if (hwtot > 0.0_wp) then
                  t_eos = thtot/hwtot
                  s_eos = shtot/hwtot
               end if
               ! MOM6 `calculate_density_derivs` (closed-form ∂ρ/∂T, ∂ρ/∂S).
               call eos_density_derivs(eos, t_eos, s_eos, press, drdt, drds)
            end if
            ! The 400 is Ci² (KW99 eq. 2.22); Boussinesq `Rho0x400_G`.
            ustarsq = 400.0_wp*rho0/GRAVITY*ustar*ustar
            htot = 0.0_wp
            thtot = 0.0_wp
            shtot = 0.0_wp
            oldfn = 0.0_wp
            do k = 1, nz - 1
               if (hav(k) <= 0.0_wp) cycle
               ! Δρ·h of the BBL with everything below homogenised.
               oldfn = drdt*(thtot - tv(k)*htot) + drds*(shtot - sv(k)*htot)
               if (oldfn >= ustarsq) exit
               dfn = (drdt*(tv(k) - tv(k + 1)) + drds*(sv(k) - sv(k + 1)))*(hav(k) + htot)
               if (oldfn + dfn <= ustarsq) then
                  dh = hav(k)
               else
                  dh = hav(k)*sqrt((ustarsq - oldfn)/dfn)
               end if
               htot = htot + dh
               thtot = thtot + tv(k)*dh
               shtot = shtot + sv(k)*dh
            end do
            ! The top layer might be part of the BBL.
            if (oldfn < ustarsq .and. hav(nz) > 0.0_wp) then
               if (drdt*(thtot - tv(nz)*htot) + drds*(shtot - sv(nz)*htot) < ustarsq) then
                  htot = htot + hav(nz)
               end if
            end if

            ! ---- 4. rotation (Killworth & Edwards 1999 eq. 2.20) + caps ----
            if (x_face) then
               c2f = f_corner(i, j) + f_corner(i, j + 1)
            else
               c2f = f_corner(i, j) + f_corner(i + 1, j)
            end if
            if (cd*u2_bg <= 0.0_wp) then
               root = sqrt(0.25_wp*ustar*ustar + (htot*c2f)**2)
               if (htot*ustar <= (thick_min + EPS_NEG)*(0.5_wp*ustar + root)) then
                  thick = thick_min
               else
                  thick = (htot*ustar)/(0.5_wp*ustar + root)
               end if
            else
               thick = htot/(0.5_wp + sqrt(0.25_wp + htot*htot*c2f*c2f/(ustar*ustar)))
               if (thick < thick_min) thick = thick_min
            end if
            if (rino_cap .and. thick > 0.5_wp*hbbl) thick = 0.5_wp*hbbl

            ! ---- 5. the viscosity that carries CDRAG·u_bbl² ----
            kv_bbl(i, j) = cd_sqrt*ustar*thick
            bbl_thick(i, j) = thick
         end if
      end do
   end subroutine bbl_faces_impl