mle_compute_transports Subroutine

   subroutine mle_compute_transports(grid, metrics, this, ms, epbl, ss, dt_limit, bc)
      !! Fill `uhml`/`vhml` (m^3/s) with the FK MLE overturning transport.
      !! Run once per outer step at thermo cadence, before the continuity
      !! divergence.  Steps: (1) b_ml + htot_ml at cell centres (surface→bed
      !! band to mld, partial-weight the straddling layer); (2) uDml/vDml at
      !! faces from grad b_bar, timescale, H_vel²; (2b) optional per-layer
      !! availability cap (a scalar shrink keeping sum_k a(k)=0); (3) fold
      !! the mu profile a(k) → uhml/vhml.  No-op when `enable=.false.` or
      !! the slot / state arrays are absent.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(ocean_mle_t), intent(inout) :: this
      type(multilayer_state_t), intent(in) :: ms
      type(ocean_epbl_t), intent(in) :: epbl
      type(ocean_surface_stress_t), intent(in), optional :: ss
      real(wp), intent(in), optional :: dt_limit
         !! Window (s) the FK transport is integrated over (= `dt_therm`).
         !! When present and > 0, the per-layer availability cap bounds the
         !! transport so the windowed drain's `hprev` reconstruction stays
         !! non-negative in thin surface layers.  Absent / ≤ 0 ⇒ no cap.
         !! PR-8: this cap is UNCONDITIONAL in production — both production
         !! call sites (`rdb_ocean_dyn.F90`) always pass
         !! `dt_limit=dyn%therm_dt(dt)` > 0, so `do_limit` is always
         !! `.true.` on the live path.  The former `&ocean_foxkemper_nml
         !! apply_cfl_limit` knob was deleted as dead/incoherent — it could
         !! only ever be set to a thing the code already always does; do
         !! not re-add a knob that toggles this cap.
      type(ocean_bc_state_t), intent(in), optional :: bc
         !! Per-edge OBC tags.  Masks the FK transport on closed
         !! (non-periodic) physical wall faces so no MLE overturning crosses
         !! a land boundary.  Absent ⇒ array-edge zeroing only.

      integer :: i, j, k, nx, ny, nz, nghost, nx_phys, ny_phys
      real(wp) :: ce_l, f_floor_l, rho0_l, g_over_rho0
      real(wp) :: cr_l, mstar_l, nstar_l, minw2_l
      logical :: use_mr_l, use_bodner_l, has_ustar, do_limit, do_filter, n_seam, s_seam, w_seam, e_seam
      logical :: use_open
      real(wp) :: h_remain, w, htot, rho_int
      real(wp) :: db, h_vel, f_abs, ustar, ts, uDml, vDml, i4dt, h_av, h_open
      real(wp) :: a_stack(NZ_STACK_MAX), hf_stack(NZ_STACK_MAX)
      real(wp) :: a_fac, b_fac, mld_inst, mld_use

      if (.not. this%is_init) return
      if (.not. this%enable) return
      if (.not. allocated(ms%rho_layer)) return
      if (.not. allocated(ms%h_layer)) return
      if (.not. allocated(epbl%mld)) return
      if (.not. allocated(epbl%f_centre)) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      nghost = grid%nghost
      nx_phys = grid%nx_phys
      ny_phys = grid%ny_phys

      ce_l = this%ce
      f_floor_l = this%f_floor
      use_mr_l = this%use_mom_mixrate
      use_bodner_l = this%use_bodner
      cr_l = this%cr
      mstar_l = this%bodner_mstar
      nstar_l = this%bodner_nstar
      minw2_l = this%min_wstar2
      ! z-level closed faces (module header): open-column ML walk + face
      ! overturning.  Host scalar; `.false.` ⇒ the legacy arithmetic.
      use_open = metrics%use_closed_faces
      ! Bodner needs the surface buoyancy flux (epbl%b0); if EPBL didn't
      ! persist it, there is nothing to restratify with -> no-op.
      if (use_bodner_l .and. .not. allocated(epbl%b0)) return
      rho0_l = epbl%rho0
      g_over_rho0 = GRAVITY/rho0_l
      ! Availability-cap setup (FK transport CFL limiter).  I4dt = 1/(4*dt)
      ! so a single donor face can evacuate at most 1/4 of the donor layer
      ! volume over the window — guarantees the windowed-drain hprev
      ! reconstruction stays positive in thin layers.
      do_limit = .false.
      i4dt = 0.0_wp
      if (present(dt_limit)) then
         if (dt_limit > 0.0_wp) then
            do_limit = .true.
            i4dt = 1.0_wp/(4.0_wp*dt_limit)
         end if
      end if
      ! u* only needed for the FK11 mixrate form; surface stress is a
      ! Pa wind stress tau, u* = sqrt(|tau|/rho0).  Bare form ignores it.
      has_ustar = .false.
      if ((use_mr_l .or. use_bodner_l) .and. present(ss)) has_ustar = .true.

      ! ---- Running-mean MLD filter ----------
      ! Psi ~ MLD^2, so a thinning/oscillating EPBL MLD injects spiky
      ! transport.  Damp with a running mean that resets instantly to a
      ! deeper MLD but decays over `mld_decay_time` when it retreats:
      !     aFac = T/(dt+T),  bFac = dt/(dt+T)
      !     MLD_filt = max( MLD, bFac*MLD + aFac*MLD_filt )
      ! `dt` = the FK call cadence (= `dt_limit`).  Filter off (default)
      ! reads the instantaneous MLD ⇒ bit-identical.  Fox-Kemper et al. (2011).
      do_filter = .false.
      if (this%mld_decay_time > 0.0_wp .and. present(dt_limit)) then
         if (dt_limit > 0.0_wp) do_filter = .true.
      end if
      if (do_filter) then
         a_fac = this%mld_decay_time/(dt_limit + this%mld_decay_time)
         b_fac = dt_limit/(dt_limit + this%mld_decay_time)
         do concurrent(j=1:ny, i=1:nx) local(mld_inst)
            mld_inst = epbl%mld(i, j)
            if (this%mld_filtered(i, j) < 0.0_wp) then
               ! Unseeded: start the running mean from the true MLD.
               this%mld_filtered(i, j) = mld_inst
            else
               this%mld_filtered(i, j) = max(mld_inst, &
                                             b_fac*mld_inst + a_fac*this%mld_filtered(i, j))
            end if
         end do
      end if

      ! ---- 1. ML-averaged buoyancy + clamped MLD at cell centres ----
      ! Walk k=nz (surface) -> bed; partial-weight the layer straddling
      ! the MLD base so htot reaches mld exactly.  b = -(g/rho0)*rho_bar.
      ! `mld_use` is the filtered MLD when the decay-time filter is on,
      ! else the instantaneous EPBL MLD (bit-identical legacy path).
      ! Closed faces (`use_open`): a filler (in the bed or the ice draft)
      ! is not water — it is skipped, so the walk starts at the first live
      ! layer from the top and its EOS-reference rho never enters b_bar.
      do concurrent(j=1:ny, i=1:nx) local(k, h_remain, w, htot, rho_int, mld_use)
         if (do_filter) then
            mld_use = this%mld_filtered(i, j)
         else
            mld_use = epbl%mld(i, j)
         end if
         htot = 0.0_wp
         rho_int = 0.0_wp
         do k = nz, 1, -1
            h_remain = mld_use - htot
            if (h_remain <= 0.0_wp) exit
            if (use_open) then
               if (.not. rdb_vl_is_live(ms%h_layer(i, j, k))) cycle
            end if
            w = min(ms%h_layer(i, j, k), h_remain)   ! partial weight
            htot = htot + w
            rho_int = rho_int + ms%rho_layer(i, j, k)*w
         end do
         this%htot_ml(i, j) = htot
         this%b_ml(i, j) = -g_over_rho0*(rho_int/(htot + H_NEGLECT))
      end do

      ! ---- 2+3. u-face transports.  Interior faces i=2..nx; wall faces
      ! (i=1, i=nx+1) get zero transport so FK injects no flux through
      ! walls (mirrors the continuity wall convention).
      ! uDml = timescale * dy_cu * (b_E - b_W) * H_vel^2.
      ! Positive (b_E - b_W) (light to the east) => positive uDml; the
      ! surface a(k)<0 then drives SURFACE transport WESTWARD toward the
      ! dense column => the front slumps / restratifies (see module head).
      do concurrent(k=1:nz, j=1:ny)
         this%uhml(1, j, k) = 0.0_wp
         this%uhml(nx + 1, j, k) = 0.0_wp
      end do
      ! Closed faces: stage the OPEN face thickness in `uhml` (the kernel
      ! below reads it, then overwrites it with the transport, column by
      ! column — no cross-iteration dependence).  Inline and host-gated:
      ! with the knob off this kernel is never launched, so `open_u` (the
      ! `(1,1,1)` placeholder) is never present-checked over the face range.
      if (use_open) then
         do concurrent(k=1:nz, j=1:ny, i=2:nx)
            if (metrics%open_u(i, j, k) > 0.5_wp .and. &
                rdb_vl_is_live(ms%h_layer(i - 1, j, k)) .and. &
                rdb_vl_is_live(ms%h_layer(i, j, k))) then
               this%uhml(i, j, k) = 0.5_wp*(ms%h_layer(i - 1, j, k) + ms%h_layer(i, j, k))
            else
               this%uhml(i, j, k) = 0.0_wp
            end if
         end do
      end if
      do concurrent(j=1:ny, i=2:nx) &
         local(k, db, h_vel, f_abs, ustar, ts, uDml, h_av, h_open, a_stack, hf_stack)
         db = this%b_ml(i, j) - this%b_ml(i - 1, j)        ! b_E - b_W
         h_vel = 0.5_wp*(this%htot_ml(i - 1, j) + this%htot_ml(i, j))
         if (use_open) then
            ! Open face column (module header): the masked thickness, and a
            ! cell no deeper than the open column so sigma reaches -1.
            h_open = 0.0_wp
            do k = 1, nz
               hf_stack(k) = this%uhml(i, j, k)
               h_open = h_open + hf_stack(k)
            end do
            h_vel = min(h_vel, h_open)
         else
            do k = 1, nz
               hf_stack(k) = 0.5_wp*(ms%h_layer(i - 1, j, k) + ms%h_layer(i, j, k))
            end do
         end if
         f_abs = abs(0.5_wp*(epbl%f_centre(i - 1, j) + epbl%f_centre(i, j)))
         ustar = 0.0_wp
         if (has_ustar) ustar = mle_face_ustar_x(ss, rho0_l, i, j)
         if (use_bodner_l) then
            ts = mle_bodner_timescale(f_abs, ustar, h_vel, &
                                      0.5_wp*(epbl%b0(i - 1, j) + epbl%b0(i, j)), &
                                      sqrt(0.5_wp*(metrics%dxCu(i, j)**2 + metrics%dyCu(i, j)**2)), &
                                      cr_l, mstar_l, nstar_l, minw2_l)
         else
            ts = mle_timescale(f_abs, ustar, h_vel, ce_l, f_floor_l, use_mr_l)
         end if
         ! db is the raw buoyancy DIFFERENCE b_E - b_W; idxCu = 1/dxCu turns
         ! it into the gradient db/dx, so dy_cu*idxCu = dyCu/dxCu is the
         ! transport aspect ratio (FK streamfunction * face width). Omitting
         ! idxCu made uDml ~dxCu too large (the FK over-amplification bug).
         uDml = ts*metrics%dy_cu(i, j)*metrics%idxCu(i, j)*db*h_vel*h_vel
         call mle_layer_weights(hf_stack, nz, h_vel, a_stack)
         ! Per-layer availability cap (MOM6).  Donor is the WEST cell (i-1)
         ! for positive transport a(k)*uDml > 0, the EAST cell (i) for
         ! negative.  Shrinking the scalar uDml preserves sum_k a(k)=0.
         if (do_limit) then
            do k = 1, nz
               if (a_stack(k)*uDml > 0.0_wp) then
                  h_av = max(i4dt*metrics%areaT(i - 1, j)* &
                             (ms%h_layer(i - 1, j, k) - MLE_H_AVAIL_MIN), 0.0_wp)
                  if (a_stack(k)*uDml > h_av) uDml = h_av/a_stack(k)
               else if (a_stack(k)*uDml < 0.0_wp) then
                  h_av = max(i4dt*metrics%areaT(i, j)* &
                             (ms%h_layer(i, j, k) - MLE_H_AVAIL_MIN), 0.0_wp)
                  if (-a_stack(k)*uDml > h_av) uDml = -h_av/a_stack(k)
               end if
            end do
         end if
         do k = 1, nz
            this%uhml(i, j, k) = a_stack(k)*uDml
         end do
         if (use_open) then
            ! Exactly zero off the open set (a(k) = 0 there already; this
            ! also holds should uDml ever be non-finite).
            do k = 1, nz
               if (.not. (hf_stack(k) > 0.0_wp)) this%uhml(i, j, k) = 0.0_wp
            end do
         end if
      end do

      ! ---- 2+3. v-face transports.  Mirror; wall faces j=1, j=ny+1 = 0.
      do concurrent(k=1:nz, i=1:nx)
         this%vhml(i, 1, k) = 0.0_wp
         this%vhml(i, ny + 1, k) = 0.0_wp
      end do
      ! Closed faces: stage the open face thickness in `vhml` (as above).
      if (use_open) then
         do concurrent(k=1:nz, j=2:ny, i=1:nx)
            if (metrics%open_v(i, j, k) > 0.5_wp .and. &
                rdb_vl_is_live(ms%h_layer(i, j - 1, k)) .and. &
                rdb_vl_is_live(ms%h_layer(i, j, k))) then
               this%vhml(i, j, k) = 0.5_wp*(ms%h_layer(i, j - 1, k) + ms%h_layer(i, j, k))
            else
               this%vhml(i, j, k) = 0.0_wp
            end if
         end do
      end if
      do concurrent(j=2:ny, i=1:nx) &
         local(k, db, h_vel, f_abs, ustar, ts, vDml, h_av, h_open, a_stack, hf_stack)
         db = this%b_ml(i, j) - this%b_ml(i, j - 1)        ! b_N - b_S
         h_vel = 0.5_wp*(this%htot_ml(i, j - 1) + this%htot_ml(i, j))
         if (use_open) then
            h_open = 0.0_wp
            do k = 1, nz
               hf_stack(k) = this%vhml(i, j, k)
               h_open = h_open + hf_stack(k)
            end do
            h_vel = min(h_vel, h_open)
         else
            do k = 1, nz
               hf_stack(k) = 0.5_wp*(ms%h_layer(i, j - 1, k) + ms%h_layer(i, j, k))
            end do
         end if
         f_abs = abs(0.5_wp*(epbl%f_centre(i, j - 1) + epbl%f_centre(i, j)))
         ustar = 0.0_wp
         if (has_ustar) ustar = mle_face_ustar_y(ss, rho0_l, i, j)
         if (use_bodner_l) then
            ts = mle_bodner_timescale(f_abs, ustar, h_vel, &
                                      0.5_wp*(epbl%b0(i, j - 1) + epbl%b0(i, j)), &
                                      sqrt(0.5_wp*(metrics%dxCv(i, j)**2 + metrics%dyCv(i, j)**2)), &
                                      cr_l, mstar_l, nstar_l, minw2_l)
         else
            ts = mle_timescale(f_abs, ustar, h_vel, ce_l, f_floor_l, use_mr_l)
         end if
         ! idyCv = 1/dyCv turns the raw db = b_N - b_S into db/dy; the
         ! dx_cv*idyCv = dxCv/dyCv aspect ratio mirrors the u-face above.
         vDml = ts*metrics%dx_cv(i, j)*metrics%idyCv(i, j)*db*h_vel*h_vel
         call mle_layer_weights(hf_stack, nz, h_vel, a_stack)
         ! Per-layer availability cap.  Donor is the SOUTH cell (i,j-1) for
         ! positive transport, the NORTH cell (i,j) for negative.
         if (do_limit) then
            do k = 1, nz
               if (a_stack(k)*vDml > 0.0_wp) then
                  h_av = max(i4dt*metrics%areaT(i, j - 1)* &
                             (ms%h_layer(i, j - 1, k) - MLE_H_AVAIL_MIN), 0.0_wp)
                  if (a_stack(k)*vDml > h_av) vDml = h_av/a_stack(k)
               else if (a_stack(k)*vDml < 0.0_wp) then
                  h_av = max(i4dt*metrics%areaT(i, j)* &
                             (ms%h_layer(i, j, k) - MLE_H_AVAIL_MIN), 0.0_wp)
                  if (-a_stack(k)*vDml > h_av) vDml = -h_av/a_stack(k)
               end if
            end do
         end if
         do k = 1, nz
            this%vhml(i, j, k) = a_stack(k)*vDml
         end do
         if (use_open) then
            do k = 1, nz
               if (.not. (hf_stack(k) > 0.0_wp)) this%vhml(i, j, k) = 0.0_wp
            end do
         end if
      end do

      ! ---- Physical closed-wall face mask -------
      ! No MLE transport may cross a closed (land) boundary.  Array-edge
      ! zeroing above misses the PHYSICAL wall faces (at i=nghost+1 /
      ! i=nghost+nx_phys+1, interior to the array); a nonzero uhml/vhml
      ! there leaks tracer across the wall (worst under the windowed drain).
      ! Zero the FK transport on every non-periodic physical edge, mirroring
      ! the continuity wall convention.  Periodic edges are seams, not walls,
      ! and so is an MPI subdomain edge (`has_* = .false.`): the face there is
      ! an interior face the neighbour rank computes identically.
      if (present(bc)) then
         if (.not. bc%periodic_x) then
            ! host scalars: never deref bc on device
            w_seam = .not. bc%has_west
            e_seam = .not. bc%has_east
            do concurrent(k=1:nz, j=1:ny)
               if (.not. w_seam) this%uhml(nghost + 1, j, k) = 0.0_wp
               if (.not. e_seam) this%uhml(nghost + nx_phys + 1, j, k) = 0.0_wp
            end do
         end if
         if (.not. bc%periodic_y) then
            ! A tripolar north fold is a seam too: its fold-line face keeps
            ! the FK transport (projected antisymmetric with the resolved
            ! mass flux in `continuity_tracer_step_split`).  `north_fold` is
            ! rank-local (the north-edge rank only), hence the `has_north`.
            s_seam = .not. bc%has_south
            n_seam = bc%north_fold .or. .not. bc%has_north
            do concurrent(k=1:nz, i=1:nx)
               if (.not. s_seam) this%vhml(i, nghost + 1, k) = 0.0_wp
               if (.not. n_seam) this%vhml(i, nghost + ny_phys + 1, k) = 0.0_wp
            end do
         end if
      end if
   end subroutine mle_compute_transports