vmix_assemble Subroutine

   subroutine vmix_assemble(grid, this, ms, geolat, status)
      !! The single downstream gate of the vmix diffusivity assembly —
      !! the `set_diffusivity`-style stage that every interior and
      !! overlay closure feeds into before vdiff consumes kv/kt/ks.
      !!
      !! CONTRACT.  Every interior / overlay closure CONTRIBUTES into
      !! `kv` / `kt` / `ks` upstream of this call: PP81 writes the
      !! interior background + Richardson term, KPP / EPBL overlay (max
      !! or add), kappa-shear merges additively, KV_ML_INVZ2 augments
      !! the surface band.  `vmix_assemble` is the SINGLE place the
      !! assembled fields are gated.  Future Area-C contributors (tidal
      !! mixing, double-diffusion, geothermal-adjacent floors, depth-
      !! varying background profiles) plug in as additional upstream
      !! contributors — they do NOT add their own floor/clip; they rely
      !! on this stage.  Applied here, in order, on interior interfaces
      !! k = 2..nz (the boundary interfaces k=1/bed and k=nz+1/surface
      !! stay at the closed-BC zero, untouched):
      !!   (a) optional debug-gated negative/NaN guard (`vmix_guard`,
      !!       default off) — runs FIRST, on the raw closure output before
      !!       any floor or clip.  The guard must see raw values because the
      !!       subsequent clip would launder NaN/negatives: NVHPC -O2
      !!       evaluates min(max(NaN, bg), huge) to bg, converting NaN
      !!       diffusivities into silent plausible mixing.  With `status`
      !!       present the routine returns a non-zero code (testable path);
      !!       without it a tripped guard `error stop`s.
      !!   (b) background floors — `kv ≥ kv_bg`, `kt ≥ kt_bg`,
      !!       `ks ≥ ks_bg`.  Defaults match `pp81_nu_bg` / `pp81_kappa_bg`
      !!       (structural invariant set in `init`) so the floor is a no-op
      !!       for the shipped closure path (bit-identical).  This is the
      !!       single place the constant background is enforced going forward.
      !!       Two MUTUALLY EXCLUSIVE opt-in variants replace the scalar
      !!       floor: `bkgnd_profile` (Bryan-Lewis per-interface `kd_bg`
      !!       field, kv floor = `bkgnd_prandtl·kd_bg`) and `bkgnd_henyey`
      !!       (Henyey latitude factor on the scalar tracer floors,
      !!       `max(bkgnd_kd_min, kt_bg·L(phi))`).  Configure refuses both
      !!       at once, so this dispatch is a three-way `if/else if/else`.
      !!   (c) ceilings — `kv ≤ kv_max`, `kt,ks ≤ kd_max`.  Defaults
      !!       `huge(1.0)` ⇒ no clip ⇒ bit-identical.
      !!   (d) optional 1-2-1 horizontal smoothing of kv/kt/ks at interfaces
      !!       (`kd_smooth_iterations`, default 0 = off).  Wet-mask aware:
      !!       contributions from dry neighbours are excluded and the stencil
      !!       weight is renormalised over wet cells only; a dry centre column
      !!       is left unchanged.  NOTE: under future MPI, multi-pass
      !!       smoothing requires nghost >= kd_smooth_iterations and a halo
      !!       refresh between passes; with nghost = 2 only 1-2 passes are
      !!       safe without the refresh.
      !!
      !! `ks` is derived from `kt` by `vmix_split_kd_heat_salt`, called
      !! upstream of this gate (after the last kv/kt contributor), and
      !! consumed by `vdiff_apply_tracers` (salinity + every passive
      !! tracer).  It equals `kt` until a double-diffusion contributor
      !! lands (PR-33).
      type(hgrid_t), intent(in) :: grid
      type(ocean_vmix_t), intent(inout) :: this
      type(multilayer_state_t), intent(in) :: ms
      real(wp), intent(in), optional :: geolat(grid%nx_total, grid%ny_total)
         !! T-point geographic latitude (degN, `ocean_metrics_t%geolatT` —
         !! the T stagger specifically; `geolatBu` is corner-shaped
         !! `(nx+1, ny+1)` and is NOT interchangeable here).  Explicit-shape
         !! so the value flows into `vmix_assemble_clip_henyey_impl`'s
         !! `do concurrent` without a descriptor walk.  Only read when
         !! `bkgnd_henyey`.  Required in that case —
         !! absent then `error stop`s (a caller forgot to thread `metrics`
         !! through).  Every production call site (`vmix_apply_in_stage`)
         !! has `metrics` in scope and always passes it; test harnesses that
         !! never enable `bkgnd_henyey` may omit it.
      integer, intent(out), optional :: status
         !! 0 = ok; 1 = guard tripped (negative or NaN K).  Only
         !! written when `vmix_guard` is on.  When absent and the guard
         !! trips, the routine `error stop`s instead.

      integer :: nx, ny, nz, it
      integer :: bad_count

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

      if (present(status)) status = 0

      ! (a) optional negative/NaN guard — BEFORE clip so laundering
      ! (NaN→bg via min/max) cannot suppress a real closure error.
      if (this%vmix_guard) then
         call vmix_guard_impl(nx, ny, nz + 1, this%kv, this%kt, this%ks, bad_count)
         if (bad_count > 0) then
            if (present(status)) then
               status = 1
            else
               error stop "vmix_assemble: negative or NaN diffusivity detected"
            end if
         end if
      end if

      ! (b)+(c) floors + ceilings on interior interfaces.  One DC over
      ! (i,j,k) with explicit-shape args via the _impl helper (assumed-
      ! shape dummies in a do concurrent make NVHPC walk descriptors
      ! per launch).  Two floor paths:
      !   * default (bkgnd_profile off): the SCALAR kv_bg/kt_bg/ks_bg
      !     floor exactly as before — bit-identical.
      !   * C7 Bryan-Lewis on: fill the per-interface kd_bg depth profile
      !     from the current column thicknesses (correct under any vcoord)
      !     and floor with that field; kv floor = bkgnd_prandtl * kd_bg.
      if (this%bkgnd_profile) then
         call vmix_bkgnd_fill_impl(nx, ny, nz + 1, this%kd_bg, ms%h_layer, &
                                   this%bkgnd_kd_sfc, this%bkgnd_kd_deep, &
                                   this%bkgnd_z0, this%bkgnd_delta)
         call vmix_assemble_clip_profile_impl(nx, ny, nz + 1, this%kv, this%kt, this%ks, &
                                              this%kd_bg, this%bkgnd_prandtl, &
                                              this%kv_max, this%kd_max)
      else if (this%bkgnd_henyey) then
         ! C7 Henyey: the latitude factor scales the SCALAR kt_bg/ks_bg
         ! floors, floored at `bkgnd_kd_min` — see the `bkgnd_henyey`
         ! docstring + `henyey_lat_factor_impl` for the exact form.  This
         ! branch is unreachable together with `bkgnd_profile`: the two are
         ! mutually exclusive (configure refuses both), so the factor never
         ! touches the Bryan-Lewis deep asymptote.
         !
         ! `geolat` is required here — a caller that turned this knob on but
         ! never threaded `metrics%geolatT` through is a wiring bug, not a
         ! silent no-op.  Deliberately a SEPARATE `_impl` from the plain
         ! scalar clip (rather than one routine with optional Henyey args):
         ! with no `present()` inside the kernel there is no way for a
         ! partially-supplied argument list to degrade into a silent no-op,
         ! and the default path's routine body is literally unchanged from
         ! before Henyey existed.
         if (.not. present(geolat)) then
            error stop "vmix_assemble: bkgnd_henyey requires the geolat "// &
               "argument (thread ocean_metrics_t%geolatT through the caller)"
         end if
         call vmix_assemble_clip_henyey_impl(nx, ny, nz + 1, this%kv, this%kt, this%ks, &
                                             this%kv_bg, this%kt_bg, this%ks_bg, &
                                             this%kv_max, this%kd_max, &
                                             this%bkgnd_henyey_n0_2omega, &
                                             this%bkgnd_henyey_max_lat, &
                                             vmix_resolve_kd_min(this%bkgnd_kd_min, &
                                                                 this%kt_bg), &
                                             geolat)
      else
         call vmix_assemble_clip_impl(nx, ny, nz + 1, this%kv, this%kt, this%ks, &
                                      this%kv_bg, this%kt_bg, this%ks_bg, &
                                      this%kv_max, this%kd_max)
      end if

      ! (d) optional 1-2-1 horizontal smoothing of kv / kt.  scratch is a
      ! persistent slot buffer mapped in enter_data — no lazy attach.
      ! Wet-mask aware: dry neighbours are excluded, stencil renormalised.
      if (this%kd_smooth_iterations > 0) then
         do it = 1, this%kd_smooth_iterations
            call vmix_smooth_121_impl(nx, ny, nz + 1, this%kv, this%smooth_scratch, ms%wet_mask)
            call vmix_smooth_121_impl(nx, ny, nz + 1, this%kt, this%smooth_scratch, ms%wet_mask)
            call vmix_smooth_121_impl(nx, ny, nz + 1, this%ks, this%smooth_scratch, ms%wet_mask)
         end do
      end if
   end subroutine vmix_assemble