top_drag_tendencies_impl Subroutine

public pure subroutine top_drag_tendencies_impl(du_drag, dv_drag, lambda_u, lambda_v, u_face, v_face, h_layer, wet_mask, cover_u, cover_v, k_top_u, k_top_v, h_vanished, variant, r, c_d, h_floor, htbl, bg_vel, tbl_min, dt_imp, fold, nx_u, ny_u, nx_v, ny_v, nx, ny, nz)

Flat device kernel: explicit-shape dummies, no derived-type dereference inside the do concurrent.

One code path covers both modes. htbl <= 0 is the LAYER-ONLY mode: the band is layer k_top alone and h_in/h_face == 1, which reduces the formulae below to the exact algebraic mirror of ocean_bottom_drag_compute_tendencies’ bed-only branch (at the default bg_vel = 0). htbl > 0 spreads the stress over the top htbl metres, the mirror of compute_distributed_drag.

Per face, two sequential passes over k = nz downward: 1. band-mean velocity U_tbl = sum_k u_k*h_in_k / max(sum_k h_in_k, tbl_min) and the band thickness; 2. per-layer rate and tendency.

Linear: rate_k = r * (h_in_k/h_face_k) Quadratic: rate_k = C_d * |U_eff| * (h_in_k/h_face_k) / h_tbl with |U_eff| = max(bg_vel, |U_tbl|), and the tendency -rate_k*u_k/(1 + dt_imp*rate_k) — dt_imp = 0 gives the explicit form bit-identically.

h_in_k/h_face_k is bounded by 1 by construction (h_in_k = min(h_face_k, ...)) and the loop exits on h_face_k <= 0, so the ratio needs no epsilon.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(out) :: du_drag(nx_u,ny_u,nz)
real(kind=wp), intent(out) :: dv_drag(nx_v,ny_v,nz)
real(kind=wp), intent(out) :: lambda_u(nx+1,ny)
real(kind=wp), intent(out) :: lambda_v(nx,ny+1)
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)
real(kind=wp), intent(in) :: cover_u(nx+1,ny)
real(kind=wp), intent(in) :: cover_v(nx,ny+1)
integer, intent(in) :: k_top_u(nx+1,ny)

ms%k_top_u / k_top_v – the first layer LIVE on BOTH sides of the face (min of the two columns’ own k_top), nz wherever nothing vanishes against the top, so the walks below start exactly where they do today on every coordinate but z_fixed under a rigid top.
integer, intent(in) :: k_top_v(nx,ny+1)

ms%k_top_u / k_top_v – the first layer LIVE on BOTH sides of the face (min of the two columns’ own k_top), nz wherever nothing vanishes against the top, so the walks below start exactly where they do today on every coordinate but z_fixed under a rigid top.
real(kind=wp), intent(in) :: h_vanished

H_VANISHED. The band walks exit on a face thickness at or below this instead of at or below ZERO: a filler has h = zstar_h_min > 0, so the old <= 0 gate let it into the band with h_in/h_face = 1 – FULL drag rate on a massless layer – while contributing nothing to cumul_h.
integer, intent(in) :: variant
real(kind=wp), intent(in) :: r
real(kind=wp), intent(in) :: c_d
real(kind=wp), intent(in) :: h_floor
real(kind=wp), intent(in) :: htbl
real(kind=wp), intent(in) :: bg_vel
real(kind=wp), intent(in) :: tbl_min
real(kind=wp), intent(in) :: dt_imp
logical, intent(in) :: fold
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~~top_drag_tendencies_impl~~CallsGraph proc~top_drag_tendencies_impl top_drag_tendencies_impl local local proc~top_drag_tendencies_impl->local

Called by

proc~~top_drag_tendencies_impl~~CalledByGraph proc~top_drag_tendencies_impl top_drag_tendencies_impl proc~ocean_top_drag_compute_tendencies ocean_top_drag_compute_tendencies proc~ocean_top_drag_compute_tendencies->proc~top_drag_tendencies_impl proc~run_stage run_stage proc~run_stage->proc~ocean_top_drag_compute_tendencies proc~run_stage_split run_stage_split proc~run_stage_split->proc~ocean_top_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 :: frac
real(kind=wp), private :: h_eff
real(kind=wp), private :: h_face_k
real(kind=wp), private :: h_in
real(kind=wp), private :: h_in_total
integer, private :: i
integer, private :: j
integer, private :: k
integer, private :: kt
logical, private :: layer_only
real(kind=wp), private :: mask_face
logical, private :: quad
real(kind=wp), private :: rate
real(kind=wp), private :: u_at_v
real(kind=wp), private :: u_int
real(kind=wp), private :: u_tbl
real(kind=wp), private :: v_at_u
real(kind=wp), private :: v_int
real(kind=wp), private :: v_tbl

Source Code

   pure subroutine top_drag_tendencies_impl(du_drag, dv_drag, lambda_u, lambda_v, &
                                            u_face, v_face, h_layer, wet_mask, &
                                            cover_u, cover_v, k_top_u, k_top_v, &
                                            h_vanished, &
                                            variant, r, c_d, h_floor, &
                                            htbl, bg_vel, tbl_min, dt_imp, fold, &
                                            nx_u, ny_u, nx_v, ny_v, &
                                            nx, ny, nz)
      !! Flat device kernel: explicit-shape dummies, no derived-type
      !! dereference inside the `do concurrent`.
      !!
      !! One code path covers both modes.  `htbl <= 0` is the LAYER-ONLY
      !! mode: the band is layer `k_top` alone and `h_in/h_face == 1`, which
      !! reduces the formulae below to the exact algebraic mirror of
      !! `ocean_bottom_drag_compute_tendencies`' bed-only branch (at the
      !! default `bg_vel = 0`).  `htbl > 0` spreads the stress over the
      !! top `htbl` metres, the mirror of `compute_distributed_drag`.
      !!
      !! Per face, two sequential passes over `k = nz` downward:
      !!   1. band-mean velocity `U_tbl = sum_k u_k*h_in_k / max(sum_k
      !!      h_in_k, tbl_min)` and the band thickness;
      !!   2. per-layer rate and tendency.
      !!
      !! Linear:     `rate_k = r * (h_in_k/h_face_k)`
      !! Quadratic:  `rate_k = C_d * |U_eff| * (h_in_k/h_face_k) / h_tbl`
      !! with `|U_eff| = max(bg_vel, |U_tbl|)`, and the tendency
      !! `-rate_k*u_k/(1 + dt_imp*rate_k)` — `dt_imp = 0` gives the
      !! explicit form bit-identically.
      !!
      !! `h_in_k/h_face_k` is bounded by 1 by construction (`h_in_k =
      !! min(h_face_k, ...)`) and the loop exits on `h_face_k <= 0`, so
      !! the ratio needs no epsilon.
      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, dt_imp
      logical, intent(in) :: fold
      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)
      real(wp), intent(in) :: cover_u(nx + 1, ny), cover_v(nx, ny + 1)
      integer, intent(in) :: k_top_u(nx + 1, ny), k_top_v(nx, ny + 1)
         !! `ms%k_top_u` / `k_top_v` -- the first layer LIVE on BOTH
         !! sides of the face (`min` of the two columns' own `k_top`),
         !! `nz` wherever nothing vanishes against the top, so the walks
         !! below start exactly where they do today on every coordinate
         !! but `z_fixed` under a rigid top.
      real(wp), intent(in) :: h_vanished
         !! `H_VANISHED`.  The band walks `exit` on a face thickness at
         !! or below this instead of at or below ZERO: a filler has
         !! `h = zstar_h_min > 0`, so the old `<= 0` gate let it into the
         !! band with `h_in/h_face = 1` -- FULL drag rate on a massless
         !! layer -- while contributing nothing to `cumul_h`.
      real(wp), intent(out) :: du_drag(nx_u, ny_u, nz), dv_drag(nx_v, ny_v, nz)
      real(wp), intent(out) :: lambda_u(nx + 1, ny), lambda_v(nx, ny + 1)

      integer :: i, j, k, kt
      logical :: layer_only, quad
      real(wp) :: cumul_h, h_face_k, h_in, mask_face, frac
      real(wp) :: h_in_total, u_int, v_int, u_tbl, v_tbl, u_at_v, v_at_u
      real(wp) :: abs_u_eff, h_eff, rate

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

      ! ---- Zero every level + both rate fields first ----
      do concurrent(k=1:nz, j=1:ny_u, i=1:nx_u)
         du_drag(i, j, k) = 0.0_wp
      end do
      do concurrent(k=1:nz, j=1:ny_v, i=1:nx_v)
         dv_drag(i, j, k) = 0.0_wp
      end do
      do concurrent(j=1:ny, i=1:nx + 1)
         lambda_u(i, j) = 0.0_wp
      end do
      do concurrent(j=1:ny + 1, i=1:nx)
         lambda_v(i, j) = 0.0_wp
      end do

      if (variant == TDRAG_LINEAR .and. r <= 0.0_wp) return
      if (quad .and. c_d <= 0.0_wp) return
      if (.not. (quad .or. variant == TDRAG_LINEAR)) return

      ! ---- East (u) faces ----
      do concurrent(j=1:ny, i=2:nx) &
         local(k, kt, cumul_h, h_face_k, h_in, mask_face, frac, &
               h_in_total, u_int, v_int, u_tbl, v_tbl, v_at_u, &
               abs_u_eff, h_eff, rate)
         mask_face = min(wet_mask(i - 1, j), wet_mask(i, j))*cover_u(i, j)
         kt = k_top_u(i, j)
         ! Pass 1: band mean.
         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_face_k = 0.5_wp*(h_layer(i - 1, j, k) + h_layer(i, j, k))
            if (h_face_k <= h_vanished) exit
            if (layer_only) then
               h_in = h_face_k
            else
               h_in = max(0.0_wp, min(h_face_k, htbl - cumul_h))
            end if
            u_int = u_int + u_face(i, j, k)*h_in
            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_int = v_int + v_at_u*h_in
            h_in_total = h_in_total + h_in
            cumul_h = cumul_h + h_face_k
         end do
         h_eff = max(h_in_total, tbl_min)
         h_eff = max(h_eff, h_floor)
         u_tbl = u_int/h_eff
         v_tbl = v_int/h_eff
         abs_u_eff = max(bg_vel, sqrt(u_tbl*u_tbl + v_tbl*v_tbl))
         ! Pass 2: per-layer rate + tendency.
         cumul_h = 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_face_k = 0.5_wp*(h_layer(i - 1, j, k) + h_layer(i, j, k))
            if (h_face_k <= h_vanished) exit
            if (layer_only) then
               h_in = h_face_k
            else
               h_in = max(0.0_wp, min(h_face_k, htbl - cumul_h))
            end if
            frac = h_in/h_face_k
            if (quad) then
               rate = mask_face*c_d*abs_u_eff*frac/h_eff
            else
               rate = mask_face*r*frac
            end if
            du_drag(i, j, k) = -rate*u_face(i, j, k)/(1.0_wp + dt_imp*rate)
            if (fold .and. k == kt) lambda_u(i, j) = rate
            cumul_h = cumul_h + h_face_k
         end do
      end do

      ! ---- North (v) faces ----
      do concurrent(j=2:ny, i=1:nx) &
         local(k, kt, cumul_h, h_face_k, h_in, mask_face, frac, &
               h_in_total, u_int, v_int, u_tbl, v_tbl, u_at_v, &
               abs_u_eff, h_eff, rate)
         mask_face = min(wet_mask(i, j - 1), wet_mask(i, j))*cover_v(i, j)
         kt = k_top_v(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_face_k = 0.5_wp*(h_layer(i, j - 1, k) + h_layer(i, j, k))
            if (h_face_k <= h_vanished) exit
            if (layer_only) then
               h_in = h_face_k
            else
               h_in = max(0.0_wp, min(h_face_k, htbl - cumul_h))
            end if
            v_int = v_int + v_face(i, j, k)*h_in
            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_int = u_int + u_at_v*h_in
            h_in_total = h_in_total + h_in
            cumul_h = cumul_h + h_face_k
         end do
         h_eff = max(h_in_total, tbl_min)
         h_eff = max(h_eff, h_floor)
         u_tbl = u_int/h_eff
         v_tbl = v_int/h_eff
         abs_u_eff = max(bg_vel, sqrt(u_tbl*u_tbl + v_tbl*v_tbl))
         cumul_h = 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_face_k = 0.5_wp*(h_layer(i, j - 1, k) + h_layer(i, j, k))
            if (h_face_k <= h_vanished) exit
            if (layer_only) then
               h_in = h_face_k
            else
               h_in = max(0.0_wp, min(h_face_k, htbl - cumul_h))
            end if
            frac = h_in/h_face_k
            if (quad) then
               rate = mask_face*c_d*abs_u_eff*frac/h_eff
            else
               rate = mask_face*r*frac
            end if
            dv_drag(i, j, k) = -rate*v_face(i, j, k)/(1.0_wp + dt_imp*rate)
            if (fold .and. k == kt) lambda_v(i, j) = rate
            cumul_h = cumul_h + h_face_k
         end do
      end do
   end subroutine top_drag_tendencies_impl