tracer_advect_vertical_one_impl Subroutine

private pure subroutine tracer_advect_vertical_one_impl(nx, ny, nz, dt, h, w_interface, hTr, F_face, budget)

Flat-impl first-order upwind-in-z vertical advection for one tracer. Two passes:

Pass 1: fill F_face(:, :, k) for k = 1..nz+1 with the upwind tracer flux at interface k. Bed (k=1) and surface (k=nz+1) faces get zero. Reads hTr only.

Pass 2: hTr(k) += dt * (F_face(k) - F_face(k+1)). Reads F_face only.

Two-pass avoids the do-concurrent race a single-pass version would have (each layer reads its neighbour’s hTr while another iteration writes it). The shared F_face scratch lives on ocean_vertical_advection_t and is reused across all tracers in one call.

Vanishing-layer donor side: if h_donor <= 0, T_donor falls back to zero (no flux).

budget (optional): Phase D v2 contributor slot. When present, the per-cell dt · (F(k) - F(k+1)) increment is also accumulated into the slot (alongside the hTr update) for budget closure.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in) :: h(nx,ny,nz)
real(kind=wp), intent(in) :: w_interface(nx,ny,nz+1)
real(kind=wp), intent(inout) :: hTr(nx,ny,nz)
real(kind=wp), intent(inout) :: F_face(nx,ny,nz+1)
real(kind=wp), intent(inout), optional :: budget(nx,ny,nz)

Calls

proc~~tracer_advect_vertical_one_impl~~CallsGraph proc~tracer_advect_vertical_one_impl tracer_advect_vertical_one_impl local local proc~tracer_advect_vertical_one_impl->local

Called by

proc~~tracer_advect_vertical_one_impl~~CalledByGraph proc~tracer_advect_vertical_one_impl tracer_advect_vertical_one_impl proc~tracer_advect_vertical tracer_advect_vertical proc~tracer_advect_vertical->proc~tracer_advect_vertical_one_impl proc~run_continuity_chain run_continuity_chain proc~run_continuity_chain->proc~tracer_advect_vertical proc~run_gm_step run_gm_step proc~run_gm_step->proc~tracer_advect_vertical proc~run_stage run_stage proc~run_stage->proc~tracer_advect_vertical 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_gm_step proc~run_stage_split run_stage_split proc~ocean_dyn_step_split->proc~run_stage_split proc~run_stage_split->proc~run_continuity_chain 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 :: T_donor
real(kind=wp), private :: h_donor
integer, private :: i
integer, private :: j
integer, private :: k
real(kind=wp), private :: w_at_k

Source Code

   pure subroutine tracer_advect_vertical_one_impl(nx, ny, nz, dt, &
                                                   h, w_interface, hTr, F_face, budget)
      !! Flat-impl first-order upwind-in-z vertical advection for
      !! one tracer.  Two passes:
      !!
      !!   Pass 1: fill `F_face(:, :, k)` for k = 1..nz+1 with the
      !!     upwind tracer flux at interface k.  Bed (k=1) and
      !!     surface (k=nz+1) faces get zero.  Reads `hTr` only.
      !!
      !!   Pass 2: `hTr(k) += dt * (F_face(k) - F_face(k+1))`.
      !!     Reads `F_face` only.
      !!
      !! Two-pass avoids the do-concurrent race a single-pass
      !! version would have (each layer reads its neighbour's hTr
      !! while another iteration writes it).  The shared `F_face`
      !! scratch lives on `ocean_vertical_advection_t` and is
      !! reused across all tracers in one call.
      !!
      !! Vanishing-layer donor side: if h_donor <= 0, `T_donor`
      !! falls back to zero (no flux).
      !!
      !! `budget` (optional): Phase D v2 contributor slot.  When
      !! present, the per-cell `dt · (F(k) - F(k+1))` increment is
      !! also accumulated into the slot (alongside the hTr update)
      !! for budget closure.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: dt
      real(wp), intent(in) :: h(nx, ny, nz)
      real(wp), intent(in) :: w_interface(nx, ny, nz + 1)
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(inout) :: F_face(nx, ny, nz + 1)
      real(wp), intent(inout), optional :: budget(nx, ny, nz)

      integer :: i, j, k
      real(wp) :: w_at_k, T_donor, h_donor

      ! ---- Pass 1: per-interface upwind tracer flux ----
      do concurrent(k=2:nz, j=1:ny, i=1:nx) &
         local(w_at_k, T_donor, h_donor)
         w_at_k = w_interface(i, j, k)
         if (w_at_k >= 0.0_wp) then
            h_donor = h(i, j, k - 1)
            if (h_donor > 0.0_wp) then
               T_donor = hTr(i, j, k - 1)/h_donor
            else
               T_donor = 0.0_wp
            end if
         else
            h_donor = h(i, j, k)
            if (h_donor > 0.0_wp) then
               T_donor = hTr(i, j, k)/h_donor
            else
               T_donor = 0.0_wp
            end if
         end if
         F_face(i, j, k) = w_at_k*T_donor
      end do
      ! Bed and surface fluxes: zero by closed BC.
      do concurrent(j=1:ny, i=1:nx)
         F_face(i, j, 1) = 0.0_wp
         F_face(i, j, nz + 1) = 0.0_wp
      end do

      ! ---- Pass 2: apply flux divergence ----
      if (present(budget)) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            hTr(i, j, k) = hTr(i, j, k) + &
                           dt*(F_face(i, j, k) - F_face(i, j, k + 1))
            budget(i, j, k) = budget(i, j, k) + &
                              dt*(F_face(i, j, k) - F_face(i, j, k + 1))
         end do
      else
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            hTr(i, j, k) = hTr(i, j, k) + &
                           dt*(F_face(i, j, k) - F_face(i, j, k + 1))
         end do
      end if
   end subroutine tracer_advect_vertical_one_impl