vdiff_apply_tracers Subroutine

   subroutine vdiff_apply_tracers(grid, this, ms, dt, kt_source, ks_source)
      !! Backward-Euler vertical diffusivity on every registered
      !! tracer.  Each tracer is converted to `T = hTr/h`, the
      !! tridiagonal solve runs, then `hTr = T*h` is reconstituted.
      !! `h_layer` is untouched.  Per-tracer
      !! `do_vertical_diffusion` flag gates participation.
      !!
      !! Diffusivity dispatch — three cases:
      !!   1. `ks_source` absent -> single-source legacy path: one
      !!      factorize from `kt_source` (3D, interface-located) if
      !!      present, else the scalar `K_v_tracer` fallback; ALL
      !!      tracers use it.  Every existing caller lands here ⇒
      !!      bit-identical to the pre-PR-20 behaviour.
      !!   2. `ks_source` present, `kt_source` absent -> programming
      !!      error, fail loud (a caller that supplies salt but not
      !!      heat has a bug — MOM6 guards the same pairing).
      !!   3. Both present -> two passes, same buffers reused: pass 1
      !!      factorizes `kt_source` and applies it to temperature
      !!      only; pass 2 factorizes `ks_source` (overwriting the
      !!      same `a_diag_t`/`b_diag_t`/`c_diag_t` buffers — no new
      !!      scratch) and applies it to salinity AND every other
      !!      registered passive tracer.  This is MOM6's `Kd_salt`
      !!      convention — "the diapycnal diffusivity of salt AND
      !!      PASSIVE TRACERS" — so passives follow salt,
      !!      not heat.
      type(hgrid_t), intent(in) :: grid
      type(ocean_vdiff_t), intent(inout) :: this
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      real(wp), intent(in), optional :: kt_source(:, :, :)
         !! Temperature diffusivity, (nx, ny, nz+1), interface-located.
      real(wp), intent(in), optional :: ks_source(:, :, :)
         !! Salinity + passive-tracer diffusivity, same shape.  MUST
         !! be accompanied by `kt_source` (case 2 fails loud).

      integer :: it, nx, ny, nz
      logical :: use_source, two_pass

      if (present(ks_source) .and. .not. present(kt_source)) then
         call logger%error("vdiff_apply_tracers: ks_source requires kt_source")
         error stop "vdiff_apply_tracers: ks_source requires kt_source"
      end if

      use_source = present(kt_source)
      two_pass = present(ks_source)
      if (.not. use_source .and. this%K_v_tracer <= 0.0_wp) return
      if (.not. allocated(ms%tracers)) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml

      if (.not. use_source) then
         call fill_kv_scalar_buf(this%kv_scalar_buf%data, &
                                 this%K_v_tracer, nx, ny, nz)
      end if

      ! The tridiagonal is tracer-independent (depends only on kv/h/dt),
      ! so build + Thomas-factorize it ONCE per diffusivity field; every
      ! tracer sharing that field then reuses the factored coefficients.
      ! Saves rebuilding the matrix (the face_thick-heavy part) per
      ! tracer — the win grows with the tracer count.  Case 3 pays two
      ! factorizations (heat, then salt+passives) instead of one; case 1
      ! is the literal pre-PR-20 single-factorize path, unchanged.
      if (use_source) then
         call build_factorize_tracer_matrix(nx, ny, nz, this%use_harmonic, dt, &
                                            kt_source, ms%h_layer, &
                                            this%a_diag_t%data, this%b_diag_t%data, &
                                            this%c_diag_t%data)
      else
         call build_factorize_tracer_matrix(nx, ny, nz, this%use_harmonic, dt, &
                                            this%kv_scalar_buf%data, ms%h_layer, &
                                            this%a_diag_t%data, this%b_diag_t%data, &
                                            this%c_diag_t%data)
      end if

      do it = 1, size(ms%tracers)
         if (.not. ms%tracers(it)%do_vertical_diffusion) cycle
         if (two_pass .and. it /= ms%idx_temperature) cycle
         select case (ms%tracers(it)%budget_id)
         case (TRACER_BUDGET_HEAT)
            call apply_factored_tracer(nx, ny, nz, ms%h_layer, ms%tracers(it)%hTr, &
                                       this%a_diag_t%data, this%b_diag_t%data, &
                                       this%c_diag_t%data, this%rhs_t%data, &
                                       budget=ms%heat_budget_vdiff)
         case (TRACER_BUDGET_SALT)
            call apply_factored_tracer(nx, ny, nz, ms%h_layer, ms%tracers(it)%hTr, &
                                       this%a_diag_t%data, this%b_diag_t%data, &
                                       this%c_diag_t%data, this%rhs_t%data, &
                                       budget=ms%salt_budget_vdiff)
         case default
            call apply_factored_tracer(nx, ny, nz, ms%h_layer, ms%tracers(it)%hTr, &
                                       this%a_diag_t%data, this%b_diag_t%data, &
                                       this%c_diag_t%data, this%rhs_t%data)
         end select
      end do

      if (.not. two_pass) return

      ! ---- Pass 2: ks_source -> salinity + every passive tracer ----
      call build_factorize_tracer_matrix(nx, ny, nz, this%use_harmonic, dt, &
                                         ks_source, ms%h_layer, &
                                         this%a_diag_t%data, this%b_diag_t%data, &
                                         this%c_diag_t%data)

      do it = 1, size(ms%tracers)
         if (.not. ms%tracers(it)%do_vertical_diffusion) cycle
         if (it == ms%idx_temperature) cycle
         if (it == ms%idx_salinity) then
            call apply_factored_tracer(nx, ny, nz, ms%h_layer, ms%tracers(it)%hTr, &
                                       this%a_diag_t%data, this%b_diag_t%data, &
                                       this%c_diag_t%data, this%rhs_t%data, &
                                       budget=ms%salt_budget_vdiff)
         else
            call apply_factored_tracer(nx, ny, nz, ms%h_layer, ms%tracers(it)%hTr, &
                                       this%a_diag_t%data, this%b_diag_t%data, &
                                       this%c_diag_t%data, this%rhs_t%data)
         end if
      end do
   end subroutine vdiff_apply_tracers