compute_distributed_drag Subroutine

private pure subroutine compute_distributed_drag(du_drag, dv_drag, u_face, v_face, h_layer, wet_mask, k_bot_u, k_bot_v, variant, r, c_d, hbbl, bg_vel, bbl_min, bed_factor, dt_imp, nx_u, ny_u, nx_v, ny_v, nx, ny, nz)

HBBL-distributed bottom drag. Mirrors MOM6’s LINEAR_DRAG and quadratic-with-HBBL formulations: the drag stress is spread across the bottom hbbl metres rather than dumped into the bed-most layer.

Linear branch (variant == BDRAG_LINEAR, r > 0): du_k/dt = -r · u_k · (h_in_bbl_k / h_face_k) · f_k Each layer that overlaps the BBL band gets a damping rate proportional to its fractional BBL coverage. Integrating over k recovers the bulk MOM6 result r · U_bbl.

Quadratic branch (variant == BDRAG_QUADRATIC, c_d > 0): 1. First pass per face: compute BBL-mean velocity U_bbl = Σ_k u_k · h_in_bbl_k / max(Σ_k h_in_bbl_k, bbl_min). 2. Effective speed |U_eff| = max(bg_vel, |U_bbl|). 3. Per-layer apply: du_k/dt = -c_d · |U_eff| · u_k · (h_in_bbl_k / h_face_k) · f_k / max(h_in_bbl_total, bbl_min)

The band walk starts at the face’s first LIVE layer counting up from the bed, kb = k_bot_u/v (≡ 1 off z_fixed ⇒ the historical k = 1 walk), so under z_fixed the inert bed fillers below it neither take a share of the drag nor count toward the band thickness.

f_k = bed_factor for k = kb and f_k = 1 otherwise — lets the bed layer carry stronger drag than the rest of the BBL while preserving HBBL-distribution shape for layers k>=2. Default bed_factor = 1.0 ⇒ f_k ≡ 1 ⇒ bit-identical to the pre-knob path.

Flat-impl: explicit-shape dummies, no derived-type derefs in the device kernels.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout) :: du_drag(nx_u,ny_u,nz)
real(kind=wp), intent(inout) :: dv_drag(nx_v,ny_v,nz)
real(kind=wp), intent(in) :: u_face(nx_u,ny_u,nz)
real(kind=wp), intent(in) :: v_face(nx_v,ny_v,nz)
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(in) :: wet_mask(nx,ny)
integer, intent(in) :: k_bot_u(nx_u,ny_u)

ms%k_bot_u/v — first live layer of the face, counting up.
integer, intent(in) :: k_bot_v(nx_v,ny_v)

ms%k_bot_u/v — first live layer of the face, counting up.
integer, intent(in) :: variant
real(kind=wp), intent(in) :: r
real(kind=wp), intent(in) :: c_d
real(kind=wp), intent(in) :: hbbl
real(kind=wp), intent(in) :: bg_vel
real(kind=wp), intent(in) :: bbl_min
real(kind=wp), intent(in) :: bed_factor
real(kind=wp), intent(in) :: dt_imp

Implicit timestep: dt for backward-Euler drag, 0 for explicit.
integer, intent(in) :: nx_u
integer, intent(in) :: ny_u
integer, intent(in) :: nx_v
integer, intent(in) :: ny_v
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz

Calls

proc~~compute_distributed_drag~~CallsGraph proc~compute_distributed_drag compute_distributed_drag local local proc~compute_distributed_drag->local

Called by

proc~~compute_distributed_drag~~CalledByGraph proc~compute_distributed_drag compute_distributed_drag proc~ocean_bottom_drag_compute_tendencies ocean_bottom_drag_compute_tendencies proc~ocean_bottom_drag_compute_tendencies->proc~compute_distributed_drag proc~run_stage run_stage proc~run_stage->proc~ocean_bottom_drag_compute_tendencies proc~run_stage_split run_stage_split proc~run_stage_split->proc~ocean_bottom_drag_compute_tendencies proc~ocean_dyn_step ocean_dyn_step proc~ocean_dyn_step->proc~run_stage proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_stage_split proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step proc~engine_step->proc~ocean_dyn_step_split proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_step proc~rdb_ocean_step rdb_ocean_step proc~rdb_ocean_step->proc~engine_step

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: abs_U_eff
real(kind=wp), private :: cumul_h
real(kind=wp), private :: f_k
real(kind=wp), private :: h_eff_denom
real(kind=wp), private :: h_face_k
real(kind=wp), private :: h_in_bbl
real(kind=wp), private :: h_in_bbl_total
integer, private :: i
integer, private :: j
integer, private :: k
integer, private :: kb
real(kind=wp), private :: mask_face
real(kind=wp), private :: u_at_v
real(kind=wp), private :: u_bbl
real(kind=wp), private :: u_bbl_int
real(kind=wp), private :: v_at_u
real(kind=wp), private :: v_bbl
real(kind=wp), private :: v_bbl_int

Source Code

   pure subroutine compute_distributed_drag(du_drag, dv_drag, u_face, v_face, &
                                            h_layer, wet_mask, k_bot_u, k_bot_v, &
                                            variant, r, c_d, hbbl, bg_vel, bbl_min, &
                                            bed_factor, dt_imp, &
                                            nx_u, ny_u, nx_v, ny_v, nx, ny, nz)
      !! HBBL-distributed bottom drag.  Mirrors MOM6's LINEAR_DRAG
      !! and quadratic-with-HBBL formulations: the drag stress is
      !! spread across the bottom `hbbl` metres rather than dumped
      !! into the bed-most layer.
      !!
      !! Linear branch (`variant == BDRAG_LINEAR, r > 0`):
      !!   du_k/dt = -r · u_k · (h_in_bbl_k / h_face_k) · f_k
      !! Each layer that overlaps the BBL band gets a damping rate
      !! proportional to its fractional BBL coverage.  Integrating
      !! over k recovers the bulk MOM6 result `r · U_bbl`.
      !!
      !! Quadratic branch (`variant == BDRAG_QUADRATIC, c_d > 0`):
      !!   1. First pass per face: compute BBL-mean velocity
      !!      `U_bbl = Σ_k u_k · h_in_bbl_k / max(Σ_k h_in_bbl_k, bbl_min)`.
      !!   2. Effective speed `|U_eff| = max(bg_vel, |U_bbl|)`.
      !!   3. Per-layer apply:
      !!      du_k/dt = -c_d · |U_eff| · u_k · (h_in_bbl_k / h_face_k) · f_k /
      !!                 max(h_in_bbl_total, bbl_min)
      !!
      !! The band walk starts at the face's first LIVE layer counting up
      !! from the bed, `kb = k_bot_u/v` (`≡ 1` off `z_fixed` ⇒ the
      !! historical `k = 1` walk), so under `z_fixed` the inert bed
      !! fillers below it neither take a share of the drag nor count
      !! toward the band thickness.
      !!
      !! `f_k = bed_factor` for `k = kb` and `f_k = 1` otherwise — lets
      !! the bed layer carry stronger drag than the rest of the BBL
      !! while preserving HBBL-distribution shape for layers k>=2.
      !! Default `bed_factor = 1.0` ⇒ `f_k ≡ 1` ⇒ bit-identical to
      !! the pre-knob path.
      !!
      !! Flat-impl: explicit-shape dummies, no derived-type derefs in
      !! the device kernels.
      integer, intent(in) :: nx_u, ny_u, nx_v, ny_v, nx, ny, nz, variant
      real(wp), intent(in) :: r, c_d, hbbl, bg_vel, bbl_min, bed_factor
      real(wp), intent(in) :: dt_imp
         !! Implicit timestep: dt for backward-Euler drag, 0 for explicit.
      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)
      integer, intent(in)     :: k_bot_u(nx_u, ny_u), k_bot_v(nx_v, ny_v)
         !! `ms%k_bot_u/v` — first live layer of the face, counting up.
      real(wp), intent(inout) :: du_drag(nx_u, ny_u, nz), dv_drag(nx_v, ny_v, nz)

      integer :: i, j, k, kb
      real(wp) :: cumul_h, h_face_k, h_in_bbl, mask_face, f_k
      real(wp) :: h_in_bbl_total, u_bbl_int, v_bbl_int, u_bbl, v_bbl
      real(wp) :: u_at_v, v_at_u, abs_U_eff, h_eff_denom

      ! Linear: per-layer independent — single pass walks k=kb upward,
      ! accumulating BBL thickness and applying the per-layer drag.
      ! Bed layer (k=kb) scaled by `bed_factor`; layers above unchanged.
      if (variant == BDRAG_LINEAR .and. r > 0.0_wp) then
         do concurrent(j=1:ny, i=2:nx) &
            local(k, kb, cumul_h, h_face_k, h_in_bbl, mask_face, f_k)
            mask_face = min(wet_mask(i - 1, j), wet_mask(i, j))
            kb = k_bot_u(i, j)
            cumul_h = 0.0_wp
            do k = kb, nz
               if (cumul_h >= hbbl) exit
               h_face_k = 0.5_wp*(h_layer(i - 1, j, k) + h_layer(i, j, k))
               if (h_face_k > 0.0_wp) then
                  h_in_bbl = max(0.0_wp, min(h_face_k, hbbl - cumul_h))
                  if (k == kb) then
                     f_k = bed_factor
                  else
                     f_k = 1.0_wp
                  end if
                  du_drag(i, j, k) = mask_face*(-r*u_face(i, j, k)*h_in_bbl/h_face_k)*f_k &
                                     /(1.0_wp + dt_imp*r*(h_in_bbl/h_face_k)*f_k)
                  cumul_h = cumul_h + h_face_k
               else
                  exit
               end if
            end do
         end do
         do concurrent(j=2:ny, i=1:nx) &
            local(k, kb, cumul_h, h_face_k, h_in_bbl, mask_face, f_k)
            mask_face = min(wet_mask(i, j - 1), wet_mask(i, j))
            kb = k_bot_v(i, j)
            cumul_h = 0.0_wp
            do k = kb, nz
               if (cumul_h >= hbbl) exit
               h_face_k = 0.5_wp*(h_layer(i, j - 1, k) + h_layer(i, j, k))
               if (h_face_k > 0.0_wp) then
                  h_in_bbl = max(0.0_wp, min(h_face_k, hbbl - cumul_h))
                  if (k == kb) then
                     f_k = bed_factor
                  else
                     f_k = 1.0_wp
                  end if
                  dv_drag(i, j, k) = mask_face*(-r*v_face(i, j, k)*h_in_bbl/h_face_k)*f_k &
                                     /(1.0_wp + dt_imp*r*(h_in_bbl/h_face_k)*f_k)
                  cumul_h = cumul_h + h_face_k
               else
                  exit
               end if
            end do
         end do
         return
      end if

      ! Quadratic with HBBL: need U_bbl, then apply.  Two-pass per
      ! face encoded as one DC with a sequential inner k-loop that
      ! does both passes (sum then write).
      if (variant == BDRAG_QUADRATIC .and. c_d > 0.0_wp) then
         do concurrent(j=2:ny, i=2:nx) &
            local(k, kb, cumul_h, h_face_k, h_in_bbl, mask_face, f_k, &
                  h_in_bbl_total, u_bbl_int, v_bbl_int, u_bbl, v_bbl, &
                  v_at_u, abs_U_eff, h_eff_denom)
            mask_face = min(wet_mask(i - 1, j), wet_mask(i, j))
            kb = k_bot_u(i, j)
            ! Pass 1: integrate u_face * h_in_bbl over the BBL band.
            cumul_h = 0.0_wp
            u_bbl_int = 0.0_wp
            v_bbl_int = 0.0_wp
            h_in_bbl_total = 0.0_wp
            do k = kb, nz
               if (cumul_h >= hbbl) exit
               h_face_k = 0.5_wp*(h_layer(i - 1, j, k) + h_layer(i, j, k))
               if (h_face_k <= 0.0_wp) exit
               h_in_bbl = max(0.0_wp, min(h_face_k, hbbl - cumul_h))
               u_bbl_int = u_bbl_int + u_face(i, j, k)*h_in_bbl
               v_at_u = 0.25_wp*( &
                        v_face(i - 1, j, k) + v_face(i, j, k) + &
                        v_face(i - 1, j + 1, k) + v_face(i, j + 1, k))
               v_bbl_int = v_bbl_int + v_at_u*h_in_bbl
               h_in_bbl_total = h_in_bbl_total + h_in_bbl
               cumul_h = cumul_h + h_face_k
            end do
            h_eff_denom = max(h_in_bbl_total, bbl_min)
            u_bbl = u_bbl_int/h_eff_denom
            v_bbl = v_bbl_int/h_eff_denom
            abs_U_eff = max(bg_vel, sqrt(u_bbl*u_bbl + v_bbl*v_bbl))
            ! Pass 2: apply stress per-layer.  Bed layer (k=kb) scaled
            ! by `bed_factor`; layers above unchanged.
            cumul_h = 0.0_wp
            do k = kb, nz
               if (cumul_h >= hbbl) exit
               h_face_k = 0.5_wp*(h_layer(i - 1, j, k) + h_layer(i, j, k))
               if (h_face_k <= 0.0_wp) exit
               h_in_bbl = max(0.0_wp, min(h_face_k, hbbl - cumul_h))
               if (k == kb) then
                  f_k = bed_factor
               else
                  f_k = 1.0_wp
               end if
               du_drag(i, j, k) = mask_face*( &
                                  -c_d*abs_U_eff*u_face(i, j, k)* &
                                  (h_in_bbl/h_face_k)/h_eff_denom)*f_k &
                                  /(1.0_wp + dt_imp*c_d*abs_U_eff*(h_in_bbl/h_face_k)/h_eff_denom*f_k)
               cumul_h = cumul_h + h_face_k
            end do
         end do
         do concurrent(j=2:ny, i=2:nx) &
            local(k, kb, cumul_h, h_face_k, h_in_bbl, mask_face, f_k, &
                  h_in_bbl_total, u_bbl_int, v_bbl_int, u_bbl, v_bbl, &
                  u_at_v, abs_U_eff, h_eff_denom)
            mask_face = min(wet_mask(i, j - 1), wet_mask(i, j))
            kb = k_bot_v(i, j)
            cumul_h = 0.0_wp
            u_bbl_int = 0.0_wp
            v_bbl_int = 0.0_wp
            h_in_bbl_total = 0.0_wp
            do k = kb, nz
               if (cumul_h >= hbbl) exit
               h_face_k = 0.5_wp*(h_layer(i, j - 1, k) + h_layer(i, j, k))
               if (h_face_k <= 0.0_wp) exit
               h_in_bbl = max(0.0_wp, min(h_face_k, hbbl - cumul_h))
               v_bbl_int = v_bbl_int + v_face(i, j, k)*h_in_bbl
               u_at_v = 0.25_wp*( &
                        u_face(i, j - 1, k) + u_face(i + 1, j - 1, k) + &
                        u_face(i, j, k) + u_face(i + 1, j, k))
               u_bbl_int = u_bbl_int + u_at_v*h_in_bbl
               h_in_bbl_total = h_in_bbl_total + h_in_bbl
               cumul_h = cumul_h + h_face_k
            end do
            h_eff_denom = max(h_in_bbl_total, bbl_min)
            u_bbl = u_bbl_int/h_eff_denom
            v_bbl = v_bbl_int/h_eff_denom
            abs_U_eff = max(bg_vel, sqrt(u_bbl*u_bbl + v_bbl*v_bbl))
            cumul_h = 0.0_wp
            do k = kb, nz
               if (cumul_h >= hbbl) exit
               h_face_k = 0.5_wp*(h_layer(i, j - 1, k) + h_layer(i, j, k))
               if (h_face_k <= 0.0_wp) exit
               h_in_bbl = max(0.0_wp, min(h_face_k, hbbl - cumul_h))
               if (k == kb) then
                  f_k = bed_factor
               else
                  f_k = 1.0_wp
               end if
               dv_drag(i, j, k) = mask_face*( &
                                  -c_d*abs_U_eff*v_face(i, j, k)* &
                                  (h_in_bbl/h_face_k)/h_eff_denom)*f_k &
                                  /(1.0_wp + dt_imp*c_d*abs_U_eff*(h_in_bbl/h_face_k)/h_eff_denom*f_k)
               cumul_h = cumul_h + h_face_k
            end do
         end do
      end if
   end subroutine compute_distributed_drag