apply_factored_tracer Subroutine

private pure subroutine apply_factored_tracer(nx, ny, nz, h_layer, hTr, a_diag, b_diag, c_diag, rhs, budget)

Apply the pre-factored tridiagonal (from build_factorize_tracer_matrix) to one tracer: form T = hTr/h, run the Thomas RHS sweep + back-substitution against the stored pivots, and reconstitute hTr = T_new * h_layer. Operates on concentration so it conserves cell-centred T; h_layer and the factored a/b/c are read-only (shared across every tracer).

When budget is present, accumulate the per-cell increment (T_new·h − hTr_old) into the contributor slot before overwriting.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(inout) :: hTr(nx,ny,nz)
real(kind=wp), intent(in) :: a_diag(nx,ny,nz)
real(kind=wp), intent(in) :: b_diag(nx,ny,nz)
real(kind=wp), intent(in) :: c_diag(nx,ny,nz)
real(kind=wp), intent(inout) :: rhs(nx,ny,nz)
real(kind=wp), intent(inout), optional :: budget(nx,ny,nz)

Calls

proc~~apply_factored_tracer~~CallsGraph proc~apply_factored_tracer apply_factored_tracer local local proc~apply_factored_tracer->local

Called by

proc~~apply_factored_tracer~~CalledByGraph proc~apply_factored_tracer apply_factored_tracer proc~vdiff_apply_tracers vdiff_apply_tracers proc~vdiff_apply_tracers->proc~apply_factored_tracer proc~vmix_apply_in_stage vmix_apply_in_stage proc~vmix_apply_in_stage->proc~vdiff_apply_tracers proc~run_stage run_stage proc~run_stage->proc~vmix_apply_in_stage proc~run_stage_split run_stage_split proc~run_stage_split->proc~vmix_apply_in_stage 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

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: hTr_new
integer, private :: i
integer, private :: j
integer, private :: k

Source Code

   pure subroutine apply_factored_tracer(nx, ny, nz, h_layer, hTr, &
                                         a_diag, b_diag, c_diag, rhs, budget)
      !! Apply the pre-factored tridiagonal (from
      !! `build_factorize_tracer_matrix`) to one tracer: form `T = hTr/h`,
      !! run the Thomas RHS sweep + back-substitution against the stored
      !! pivots, and reconstitute `hTr = T_new * h_layer`.  Operates on
      !! concentration so it conserves cell-centred T; `h_layer` and the
      !! factored `a/b/c` are read-only (shared across every tracer).
      !!
      !! When `budget` is present, accumulate the per-cell increment
      !! (T_new·h − hTr_old) into the contributor slot before overwriting.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(in) :: a_diag(nx, ny, nz)
      real(wp), intent(in) :: b_diag(nx, ny, nz)
      real(wp), intent(in) :: c_diag(nx, ny, nz)
      real(wp), intent(inout) :: rhs(nx, ny, nz)
      real(wp), intent(inout), optional :: budget(nx, ny, nz)

      integer :: i, j, k
      real(wp) :: hTr_new

      do concurrent(j=1:ny, i=1:nx) local(k, hTr_new)
         ! RHS = T = hTr/h̃, where h̃ = max(h, H_VANISHED) is the SAME floored
         ! thickness used by build_factorize_tracer_matrix.  A vanishing /
         ! collapsed layer (h ≤ H_VANISHED, e.g. a thin z* surface layer
         ! driven ≤ 0 by combined Fox-Kemper + resolved transport on an
         ! intermediate RK2 stage) is decoupled (identity row in the matrix)
         ! and reconstituted with the same h̃ ⇒ hTr is preserved EXACTLY
         ! rather than zeroed (the previous `h ≤ 0 → T = 0 → hTr = 0` branch
         ! silently dropped its frozen tracer mass → ~5%/day leak under
         ! FK × windowed advect).  For h ≫ H_VANISHED (sigma / double_gyre)
         ! h̃ = h exactly ⇒ bit-identical.
         do k = 1, nz
            rhs(i, j, k) = hTr(i, j, k)/max(h_layer(i, j, k), H_VANISHED)
         end do

         ! ---- Thomas RHS forward sweep against the stored pivots ----
         rhs(i, j, 1) = rhs(i, j, 1)/b_diag(i, j, 1)
         do k = 2, nz
            rhs(i, j, k) = (rhs(i, j, k) - a_diag(i, j, k)*rhs(i, j, k - 1))/b_diag(i, j, k)
         end do

         ! ---- Back-substitution (rhs now holds T_new) ----
         do k = nz - 1, 1, -1
            rhs(i, j, k) = rhs(i, j, k) - c_diag(i, j, k)*rhs(i, j, k + 1)
         end do

         ! ---- Reconstitute hTr = T_new * h ----
         ! Budget write gated INSIDE the single do concurrent (splitting
         ! present() into two loops makes NVHPC emit a far slower kernel
         ! for one branch — see the remap fix).
         if (present(budget)) then
            do k = 1, nz
               hTr_new = rhs(i, j, k)*max(h_layer(i, j, k), H_VANISHED)
               budget(i, j, k) = budget(i, j, k) + (hTr_new - hTr(i, j, k))
               hTr(i, j, k) = hTr_new
            end do
         else
            do k = 1, nz
               hTr(i, j, k) = rhs(i, j, k)*max(h_layer(i, j, k), H_VANISHED)
            end do
         end if
      end do
   end subroutine apply_factored_tracer