meke_step Subroutine

   subroutine meke_step(grid, metrics, this, gm, varmix, wavespeed, ms, dt, ke_diss_ext)
      !! Advance the MEKE field one thermo step (Strang split), update the
      !! derived diffusivity `kh_diff`, and feed the geometric-mean kh into
      !! VarMix's per-face KhTh/KhTr (the GM↔MEKE feedback).  Run once per
      !! outer step at thermo cadence, after `varmix_compute` and before
      !! `gm_compute_transports` (MEKE reads `gm%gm_src` from the previous
      !! thermo step — a one-step lag).  No-op when `enable=.false.`,
      !! uninitialised, or the GM slot is absent.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(ocean_meke_t), intent(inout) :: this
      type(ocean_gm_t), intent(in) :: gm
      type(ocean_varmix_t), intent(inout), optional :: varmix
      type(ocean_wave_speed_t), intent(in), optional :: wavespeed
      type(multilayer_state_t), intent(in) :: ms
      real(wp), intent(in) :: dt
      real(wp), intent(in), optional :: ke_diss_ext(:, :)
         !! hvisc KE-dissipation rate `(nx,ny)` for the frictional source;
         !! absent ⇒ the source is inert (staged to 0).

      integer :: nx, ny, nz
      real(wp) :: sdt, sdt_damp, damp_step
      logical :: have_rho, kh_flux_enabled, have_ws, have_vm

      if (.not. this%is_init) return
      if (.not. this%enable) return
      if (.not. gm%is_init) return
      if (.not. allocated(ms%h_layer)) return
      if (.not. allocated(gm%gm_src)) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      if (this%nz_ml /= nz) return

      sdt = dt*this%dtscale
      damp_step = 1.0_wp
      if (this%kh >= 0.0_wp .or. this%k4 >= 0.0_wp) damp_step = 0.5_wp
      sdt_damp = sdt*damp_step
      kh_flux_enabled = (this%kh >= 0.0_wp)
      have_rho = allocated(ms%rho_layer)
      have_ws = .false.
      if (present(wavespeed)) have_ws = wavespeed%is_init .and. allocated(wavespeed%rd_over_dx)
      have_vm = .false.
      if (present(varmix)) have_vm = varmix%is_init .and. allocated(varmix%sn_u)

      ! ---- 0. stage the optional inputs into device-resident workspaces.
      ! Copy ON-DEVICE from the slot arrays (which are device-resident) so the
      ! length-scale kernel reads explicit-shape slot args, never an
      ! optional-or-host-staged actual.  When the source slot is absent the
      ! workspace is zeroed (Ldeform->0 / Eady inert).  `f_centre` is the
      ! slot's own |f| (filled at setup via set_f_centre); when never filled
      ! it is zero ⇒ beta=0 ⇒ Rhines inert (alpha_rhines default 0).
      if (have_ws) then
         call meke_stage_rd(nx, ny, wavespeed%rd_over_dx, this%rd_ws)
      else
         call meke_zero_2d(nx, ny, this%rd_ws)
      end if
      if (have_vm) then
         call meke_stage_sn(nx, ny, varmix%sn_u, varmix%sn_v, &
                            this%sn_u_ws, this%sn_v_ws)
      else
         call meke_zero_2d(nx + 1, ny, this%sn_u_ws)
         call meke_zero_2d(nx, ny + 1, this%sn_v_ws)
      end if

      ! ---- 1. column mass + depth (I_mass, depth_tot, mass_ws). ----
      ! Pass rho_layer only when allocated; otherwise feed h_layer as a
      ! harmless placeholder (gated off by have_rho ⇒ mass uses gm%rho0).
      if (have_rho) then
         call meke_mass(nx, ny, nz, gm%rho0, .true., ms%h_layer, ms%rho_layer, &
                        this%i_mass, this%depth_tot, this%mass_ws)
      else
         call meke_mass(nx, ny, nz, gm%rho0, .false., ms%h_layer, ms%h_layer, &
                        this%i_mass, this%depth_tot, this%mass_ws)
      end if

      ! ---- 2. structure factors (bottomFac2, barotrFac2) + Lmix into le. ----
      call meke_length_scales(nx, ny, this%cd_scale, this%cb, this%ct, &
                              this%min_gamma2, this%cdrag, &
                              this%alpha_deform, this%alpha_rhines, &
                              this%alpha_eady, this%alpha_frict, this%alpha_grid, &
                              metrics%areaT, metrics%idxT, metrics%idyT, &
                              this%f_centre, &
                              this%depth_tot, this%meke, &
                              this%rd_ws, this%sn_u_ws, this%sn_v_ws, &
                              this%bottom_fac2, this%barotr_fac2, this%le)

      ! ---- 3. explicit source bump: E += sdt*src. ----
      ! Stage the hvisc KE-dissipation rate for the frictional source
      ! (-frcoeff·i_mass·ke_diss); 0 when the seam is not wired ⇒ inert.
      if (present(ke_diss_ext)) then
         call meke_stage_rd(nx, ny, ke_diss_ext, this%ke_diss_ws)
      else
         call meke_zero_2d(nx, ny, this%ke_diss_ws)
      end if
      call meke_source(nx, ny, this%bgsrc, this%gmcoeff, this%frcoeff, sdt, &
                       this%i_mass, gm%gm_src, this%ke_diss_ws, this%src, this%meke)

      ! ---- 4. implicit drag half: E <- E/(1+sdt_damp*damp_rate). ----
      ! Resolved bed-layer eddy velocity (MOM6 drag_rate_visc); 0 when off.
      if (this%use_bbl_drag) then
         call meke_bbl_speed2(nx, ny, nz, ms%u_face_x_layer, ms%v_face_y_layer, &
                              ms%k_bot_u, ms%k_bot_v, this%u_bbl2)
      end if
      call meke_drag(nx, ny, sdt_damp, this%damping, this%cdrag, this%uscale, gm%rho0, &
                     this%i_mass, this%bottom_fac2, this%u_bbl2, this%meke)

      ! ---- 5. lateral diffusion (+ biharmonic). ----
      if (kh_flux_enabled .or. this%k4 >= 0.0_wp) then
         call meke_lateral(nx, ny, sdt, this%kh, this%k4, this%khmeke_fac, &
                           kh_flux_enabled, &
                           metrics%dy_cu, metrics%dx_cv, metrics%idxCu, &
                           metrics%idyCv, metrics%iareaT, &
                           this%i_mass, this%mass_ws, &
                           this%kh_diff, this%uflux, this%vflux, this%del2, &
                           this%meke)
      end if

      ! ---- 5b. upwind advection of E by the barotropic transport. ----
      ! Self-contained transport stage: baroHu = Sum_k (mass-weighted face
      ! transport), upwind flux E*baroHu, single conservative divergence.
      ! `advection_factor = 0` (default) ⇒ no-op ⇒ bit-identical.
      if (this%advection_factor > 0.0_wp) then
         if (have_rho) then
            call meke_baro_transport(nx, ny, nz, gm%rho0, .true., &
                                     ms%mass_flux_x_layer, ms%mass_flux_y_layer, &
                                     ms%rho_layer, this%baro_hu, this%baro_hv)
         else
            call meke_baro_transport(nx, ny, nz, gm%rho0, .false., &
                                     ms%mass_flux_x_layer, ms%mass_flux_y_layer, &
                                     ms%rho_layer, this%baro_hu, this%baro_hv)
         end if
         call meke_advect(nx, ny, sdt, this%advection_factor, metrics%iareaT, &
                          this%i_mass, this%baro_hu, this%baro_hv, &
                          this%uflux, this%vflux, this%meke)
      end if

      ! ---- 6. implicit drag half (only when Strang-split, damp_step=0.5). ----
      ! `u_bbl2` already filled in step 4 (or held at 0 when off).
      if (this%kh >= 0.0_wp .or. this%k4 >= 0.0_wp) then
         call meke_drag(nx, ny, sdt_damp, this%damping, this%cdrag, this%uscale, gm%rho0, &
                        this%i_mass, this%bottom_fac2, this%u_bbl2, this%meke)
      end if

      ! ---- 7. MEKE -> KhTh closure: kh = khcoeff*sqrt(2*gamma_t2*E)*Lmix. ----
      call meke_kh_closure(nx, ny, this%khcoeff, &
                           this%barotr_fac2, this%meke, this%le, this%kh_diff)

      ! ---- 7b. MEKE -> Ku backscatter coefficient (v1, harmonic only). ----
      ! ku = visc_coeff_ku*sqrt(2*gamma_t2*E)*Lmix; same eddy-velocity ×
      ! mixing-length form as kh.  Filled only when `backscatter`; the
      ! driver subtracts a face-average of `ku` from the resolved harmonic
      ! viscosity (stability-floored) ⇒ a negative-viscosity energy return.
      ! Default off ⇒ ku stays 0 ⇒ bit-identical.
      call meke_ku_closure(nx, ny, this%backscatter, this%visc_coeff_ku, &
                           this%meke, this%le, this%ku)

      ! ---- 8. feedback seam: add geom-mean kh into VarMix face KhTh/KhTr. ----
      if (present(varmix)) then
         if (varmix%is_init .and. (this%khth_fac /= 0.0_wp .or. &
                                   this%khtr_fac /= 0.0_wp)) then
            call meke_feed_khth(nx, ny, this%khth_fac, this%khtr_fac, &
                                this%kh_diff, varmix%khth_u, varmix%khth_v, &
                                varmix%khtr_u, varmix%khtr_v)
         end if
      end if
   end subroutine meke_step