ocean_bottom_drag_compute_tendencies Subroutine

   pure subroutine ocean_bottom_drag_compute_tendencies(grid, this, ms, dt)
      !! Fill `du_drag` / `dv_drag` with the bottom-layer drag
      !! acceleration on each face's first LIVE layer `k_bot_u/v`
      !! (`≡ 1` off `z_fixed`); every other layer gets zero in the
      !! bed-only mode (default).  When `hbbl > 0` the stress is distributed
      !! across the bottom-most `hbbl` metres — every layer with
      !! `cumulative_depth_from_bed_top ≤ hbbl` gets a proportional
      !! share of the drag tendency.
      !!
      !! `implicit=.true.` (knob, default .false.) makes the drag
      !! BACKWARD-EULER in the dragged velocity: a per-face/-layer rate
      !! `λ` (= `c_d·|U|/h` etc.) gives `u^{n+1} = u/(1+dt·λ)`, formed by
      !! the tendency `-λ·u/(1+dt·λ)` so the standalone `u += dt·du_drag`
      !! apply reproduces it exactly.  Unconditionally stable for ANY h
      !! (matches MOM6's implicit bottom-BC drag); the explicit form
      !! (default) is conditionally unstable on thin bottom layers
      !! (`λ·dt > 1`).  `dt` is unused in the explicit branch.
      !!
      !! Wall faces (i=1, i=nx+1 for u; j=1, j=ny+1 for v) get a
      !! zero tendency — they don't move under drag because they
      !! don't move at all under any kernel here.
      type(hgrid_t), intent(in) :: grid
      type(ocean_bottom_drag_t), intent(inout) :: this
      type(multilayer_state_t), intent(in) :: ms
      real(wp), intent(in) :: dt
         !! Outer-step length (s); only read when `this%implicit`.

      integer :: i, j, k, nx, ny, nz, kb
      real(wp) :: r, c_d, h_floor, u_bot, v_bot, h_face, u_at_v, v_at_u
      real(wp) :: speed_at_u, speed_at_v, lam, dt_imp
      real(wp) :: hbbl, bg_vel, bbl_min
      logical :: implicit_drag, fold

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      r = this%r_linear
      c_d = this%c_drag
      h_floor = this%h_min
      hbbl = this%hbbl
      bg_vel = this%drag_bg_vel
      bbl_min = this%bbl_thick_min
      if (bbl_min <= 0.0_wp) bbl_min = h_floor
      implicit_drag = this%implicit
      dt_imp = merge(dt, 0.0_wp, implicit_drag)  ! 0 ⇒ explicit (bit-identical)
      fold = this%implicit_fold

      ! ---- Zero every level first; only k = k_bot will get filled ----
      do concurrent(k=1:nz, j=1:ny, i=1:nx + 1)
         this%du_drag%data(i, j, k) = 0.0_wp
      end do
      do concurrent(k=1:nz, j=1:ny + 1, i=1:nx)
         this%dv_drag%data(i, j, k) = 0.0_wp
      end do
      ! Always re-zero the implicit-fold bed-rate fields so a face that
      ! left the wet stencil between steps doesn't carry a stale rate.
      do concurrent(j=1:ny, i=1:nx + 1)
         this%lambda_bot_u(i, j) = 0.0_wp
      end do
      do concurrent(j=1:ny + 1, i=1:nx)
         this%lambda_bot_v(i, j) = 0.0_wp
      end do

      ! ---- Bed-only implicit-fold rate fields (k = k_bot) ----
      ! λ is the Rayleigh RATE the vdiff diagonal consumes (`+dt·λ`): r
      ! (linear) or c_d·|U_bbl|/h_kb (quadratic, |U| frozen at uⁿ), all
      ! read on the face's first LIVE layer `kb = k_bot_u/v` — the SAME
      ! row `diffuse_velocity_columns_impl` adds it to.  Filled here from
      ! the SAME stencil that forms du_drag so there is one drag algebra.
      ! Configure forbids hbbl>0 + fold, so the bed-only 2D field is
      ! sufficient.  Zero where the face touches land (mask).
      if (fold) then
         if (this%variant == BDRAG_LINEAR .and. r > 0.0_wp) then
            do concurrent(j=1:ny, i=2:nx)
               this%lambda_bot_u(i, j) = min(ms%wet_mask(i - 1, j), ms%wet_mask(i, j))*r
            end do
            do concurrent(j=2:ny, i=1:nx)
               this%lambda_bot_v(i, j) = min(ms%wet_mask(i, j - 1), ms%wet_mask(i, j))*r
            end do
         else if (this%variant == BDRAG_QUADRATIC .and. c_d > 0.0_wp) then
            do concurrent(j=1:ny, i=2:nx) local(kb, u_bot, v_at_u, h_face, speed_at_u)
               kb = ms%k_bot_u(i, j)
               u_bot = ms%u_face_x_layer(i, j, kb)
               v_at_u = 0.25_wp*( &
                        ms%v_face_y_layer(i - 1, j, kb) + ms%v_face_y_layer(i, j, kb) + &
                        ms%v_face_y_layer(i - 1, j + 1, kb) + ms%v_face_y_layer(i, j + 1, kb))
               h_face = max(0.5_wp*(ms%h_layer(i - 1, j, kb) + ms%h_layer(i, j, kb)), h_floor)
               speed_at_u = sqrt(u_bot*u_bot + v_at_u*v_at_u)
               this%lambda_bot_u(i, j) = min(ms%wet_mask(i - 1, j), ms%wet_mask(i, j))* &
                                         c_d*speed_at_u/h_face
            end do
            do concurrent(j=2:ny, i=1:nx) local(kb, v_bot, u_at_v, h_face, speed_at_v)
               kb = ms%k_bot_v(i, j)
               v_bot = ms%v_face_y_layer(i, j, kb)
               u_at_v = 0.25_wp*( &
                        ms%u_face_x_layer(i, j - 1, kb) + ms%u_face_x_layer(i + 1, j - 1, kb) + &
                        ms%u_face_x_layer(i, j, kb) + ms%u_face_x_layer(i + 1, j, kb))
               h_face = max(0.5_wp*(ms%h_layer(i, j - 1, kb) + ms%h_layer(i, j, kb)), h_floor)
               speed_at_v = sqrt(v_bot*v_bot + u_at_v*u_at_v)
               this%lambda_bot_v(i, j) = min(ms%wet_mask(i, j - 1), ms%wet_mask(i, j))* &
                                         c_d*speed_at_v/h_face
            end do
         end if
      end if

      ! ---- HBBL-distributed branch (MOM6 LINEAR_DRAG / BBL_THICK_MIN) ----
      if (hbbl > 0.0_wp) then
         call compute_distributed_drag(this%du_drag%data, this%dv_drag%data, &
                                       ms%u_face_x_layer, ms%v_face_y_layer, &
                                       ms%h_layer, ms%wet_mask, ms%k_bot_u, ms%k_bot_v, &
                                       this%variant, r, c_d, hbbl, bg_vel, bbl_min, &
                                       this%bed_factor, dt_imp, &
                                       size(ms%u_face_x_layer, 1), size(ms%u_face_x_layer, 2), &
                                       size(ms%v_face_y_layer, 1), size(ms%v_face_y_layer, 2), &
                                       nx, ny, nz)
         return
      end if

      ! Face wet-mask: drag only fires at faces between two ocean cells.
      ! `min(wet_left, wet_right)` zeros the tendency at any face that
      ! touches land.  All-1.0 mask (analytical tests) is a no-op.

      ! ---- Linear branch ----
      ! Implicit (dt_imp>0): u^{n+1}=u/(1+dt·r) via tendency -r·u/(1+dt·r).
      if (this%variant == BDRAG_LINEAR .and. r > 0.0_wp) then
         do concurrent(j=1:ny, i=2:nx) local(kb, u_bot)
            kb = ms%k_bot_u(i, j)
            u_bot = ms%u_face_x_layer(i, j, kb)
            this%du_drag%data(i, j, kb) = &
               min(ms%wet_mask(i - 1, j), ms%wet_mask(i, j))*(-r*u_bot/(1.0_wp + dt_imp*r))
         end do
         do concurrent(j=2:ny, i=1:nx) local(kb, v_bot)
            kb = ms%k_bot_v(i, j)
            v_bot = ms%v_face_y_layer(i, j, kb)
            this%dv_drag%data(i, j, kb) = &
               min(ms%wet_mask(i, j - 1), ms%wet_mask(i, j))*(-r*v_bot/(1.0_wp + dt_imp*r))
         end do
         return
      end if

      ! ---- Quadratic branch ----
      if (this%variant == BDRAG_QUADRATIC .and. c_d > 0.0_wp) then
         ! du/dt = -C_d * |U| * u / h_bot, where |U| = sqrt(u^2 + v^2)
         ! evaluated at the same face.  For u-faces we average v from
         ! the four surrounding v-faces (standard C-grid stencil); for
         ! v-faces we average u from the four surrounding u-faces.
         do concurrent(j=1:ny, i=2:nx) &
            local(kb, u_bot, v_at_u, h_face, speed_at_u)
            kb = ms%k_bot_u(i, j)
            u_bot = ms%u_face_x_layer(i, j, kb)
            v_at_u = 0.25_wp*( &
                     ms%v_face_y_layer(i - 1, j, kb) + ms%v_face_y_layer(i, j, kb) + &
                     ms%v_face_y_layer(i - 1, j + 1, kb) + ms%v_face_y_layer(i, j + 1, kb))
            h_face = 0.5_wp*(ms%h_layer(i - 1, j, kb) + ms%h_layer(i, j, kb))
            h_face = max(h_face, h_floor)
            speed_at_u = sqrt(u_bot*u_bot + v_at_u*v_at_u)
            ! Implicit: denom h_face → h_face + dt·c_d·|U| ⇒ u/(1+dt·c_d·|U|/h).
            this%du_drag%data(i, j, kb) = &
               min(ms%wet_mask(i - 1, j), ms%wet_mask(i, j))* &
               (-c_d*speed_at_u*u_bot/(h_face + dt_imp*c_d*speed_at_u))
         end do
         do concurrent(j=2:ny, i=1:nx) &
            local(kb, v_bot, u_at_v, h_face, speed_at_v)
            kb = ms%k_bot_v(i, j)
            v_bot = ms%v_face_y_layer(i, j, kb)
            u_at_v = 0.25_wp*( &
                     ms%u_face_x_layer(i, j - 1, kb) + ms%u_face_x_layer(i + 1, j - 1, kb) + &
                     ms%u_face_x_layer(i, j, kb) + ms%u_face_x_layer(i + 1, j, kb))
            h_face = 0.5_wp*(ms%h_layer(i, j - 1, kb) + ms%h_layer(i, j, kb))
            h_face = max(h_face, h_floor)
            speed_at_v = sqrt(v_bot*v_bot + u_at_v*u_at_v)
            this%dv_drag%data(i, j, kb) = &
               min(ms%wet_mask(i, j - 1), ms%wet_mask(i, j))* &
               (-c_d*speed_at_v*v_bot/(h_face + dt_imp*c_d*speed_at_v))
         end do
      end if
   end subroutine ocean_bottom_drag_compute_tendencies