epbl_column_kernel Subroutine

   pure subroutine epbl_column_kernel(grid, this, ms, hT, hS, ss, dt, &
                                      Q_heat_field, inv_rho0_cp, &
                                      Q_salt_field, inv_rho0, &
                                      sw_src_field, sw_ctke_active, sw_pen_frac, &
                                      sw_R, sw_zeta1, sw_zeta2, wet_mask_field, &
                                      p_top_in_eos, nx_arg, ny_arg)
      !! Per-column EPBL solve.  One `do concurrent (j, i)` with the
      !! serial work in k inside (j -> i -> k ordering); ALL sweep
      !! state is carried in scalars (design doc D6) — the only
      !! column arrays are the six iteration-invariant workspaces
      !! filled by the prep sweep.
      !!
      !! Kinematic surface fluxes are read per-column inside the DC:
      !!   q_t_kin = Q_heat_field(i,j) / (rho0 · cp)  [degC·m/s]
      !!   q_s_kin = Q_salt_field(i,j) / rho0           [PSU·m/s]
      !! For the constant-fill default both fields are uniform, giving
      !! arithmetic identical to the former scalar-broadcast path.
      !!
      !! Algorithm:
      !!   prep:  T0, S0 and the pressure-weighted PE / steric
      !!          sensitivities per layer (downward pressure sum).
      !!   outer: MLD root-find (false position / bisection) —
      !!          mstar and the mixing-length shape depend on MLD.
      !!   sweep: interfaces Ki = nz..2 downward.  Decay mech TKE
      !!          across the layer above; rotation-reduce the
      !!          convective reservoir; closed-form energy solve for
      !!          the largest affordable Kd (the gravity-wave
      !!          column-height correction folds into PEc_core);
      !!          advance the embedded tridiagonal forward
      !!          elimination (hp_a / dX_to_dPE_a / Th_a recursions).
      type(hgrid_t), intent(in) :: grid
      type(ocean_epbl_t), intent(inout) :: this
      type(multilayer_state_t), intent(in) :: ms
      ! assumed-shape-ok: tracer registry outer-shim — caller host-dereferences
      ! ms%tracers(idx)%hTr before passing; size varies per tracer slot
      ! (see CLAUDE.md "outer-shim + flat-impl" pattern); thermo cadence.
      real(wp), intent(in) :: hT(:, :, :)
         !! Temperature tracer hTr (degC*m), host-dereferenced.
      real(wp), intent(in) :: hS(:, :, :)  ! assumed-shape-ok: tracer registry outer-shim; thermo cadence
         !! Salinity tracer hTr (PSU*m), host-dereferenced.
      type(ocean_surface_stress_t), intent(in) :: ss
      real(wp), intent(in) :: dt
      real(wp), intent(in) :: Q_heat_field(:, :)   ! assumed-shape-ok: outer-shim pass; thermo cadence
         !! 2D heat-flux field (W/m²), host-dereferenced from sf%Q_heat.
      real(wp), intent(in) :: inv_rho0_cp
         !! Precomputed 1/(rho0·cp) multiplier.
      real(wp), intent(in) :: Q_salt_field(:, :)   ! assumed-shape-ok: outer-shim pass; thermo cadence
         !! 2D salt-flux field (kg/m²/s), host-dereferenced from sf%Q_salt.
      real(wp), intent(in) :: inv_rho0
         !! Precomputed 1/rho0 multiplier.
      real(wp), intent(in) :: sw_src_field(:, :)   ! assumed-shape-ok: outer-shim pass; thermo cadence
         !! (PR-21) 2D irradiance source (W/m², >= 0 on the q_sw path),
         !! host-selected from sf%Q_heat or sf%q_sw.  Read only when
         !! `sw_ctke_active`; otherwise the legacy sf%Q_heat is passed and
         !! ignored.
      logical, intent(in) :: sw_ctke_active
         !! (PR-21) Charge the TKE ledger for penetrating SW (host-side
         !! `epbl_sw_ctke .and. sf%has_sw`).  False ⇒ the surface
         !! energetics + sweep run the unmodified legacy lines.
      real(wp), intent(in) :: sw_pen_frac, sw_R, sw_zeta1, sw_zeta2
         !! (PR-21) Two-band SW parameters (from sf), by value.
      real(wp), intent(in) :: wet_mask_field(:, :)   ! assumed-shape-ok: outer-shim pass; thermo cadence
         !! (PR-21) Wet mask — mirrors the deposition kernel's I0 gate so
         !! the ledger and the tracer field account the same heat.
      logical, intent(in) :: p_top_in_eos
         !! (E3) Seed the column pressure stack at `ms%p_top(i,j)` rather
         !! than at 0 Pa — `&ocean_psurf_nml in_eos`, by value.  `.false.`
         !! ⇒ the pre-E3 arithmetic, character for character.
      integer, intent(in) :: nx_arg, ny_arg
         !! Grid extents — used to bounds-check Q_* indexing.

      integer :: i, j, nx, ny, nz
      ! prep locals
      integer :: k
      real(wp) :: q_t_kin, q_s_kin
      real(wp) :: hk, hk_eff, inv_h, t0k, s0k, dmass, dpres, p_mid, pres
      real(wp) :: dsv_dt_k, dsv_ds_k, dsv_dt_sfc, dsv_ds_sfc, h_sum
      ! forcing / environment locals
      real(wp) :: ustar, absf, idecay, mech_in
      real(wp) :: b0, ctke_sfc
      ! MLD iteration locals
      integer :: obl_it, n_its
      real(wp) :: min_mld, max_mld, mld_guess, mld_found
      real(wp) :: dmld_min, dmld_max
      logical :: have_min, have_max
      real(wp) :: mstar_val, mech_tke, conv_perel, forcing_clip
      ! sweep carried scalars
      integer :: ki, ka, kb
      real(wp) :: htot, z_int, pres_int, mld_output
      logical :: sfc_connected, sfc_disconnect
      real(wp) :: hp_a, dpe_t_a, dpe_s_a, dch_t_a, dch_s_a
      real(wp) :: th_a, sh_a, te_lag, se_lag, kddt_prev, kddt_cur
      real(wp) :: exp_kh, nstar_fc, tot_tke, dt_h, tke_here, vstar
      real(wp) :: h_ka, h_kb, hb_hs, shape_fn, hbs, mixlen, kd_g0
      real(wp) :: hp_b, th_b, sh_b, hps, bdt1, dt_c, ds_c
      real(wp) :: pec_core, colht_core, dkddt, pe_max
      real(wp) :: kd_val, tke_used, frac_bl, pe_g0, dpe_conv
      real(wp) :: b1, c1, r_reduc, te_new, se_new, surf_scale
      ! Langmuir (LF17) wave state + Langmuir number
      real(wp) :: lt_u10, lt_ustokes, lt_kphil, la_val
      ! TKE budget ledger (final iteration's values survive)
      real(wp) :: d_wind, d_conv, d_forcing, d_mixing, d_mdecay, d_cdecay
      ! (PR-21) penetrating-SW TKE ledger locals
      real(wp) :: i0_col, d_top_sw, d_bot_sw, heat_sw1, heat_sw2, q_nonpen_kin
      real(wp) :: ctke_sw_kb, r_sw

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

      do concurrent(j=1:ny, i=1:nx) &
         local(k, hk, hk_eff, inv_h, t0k, s0k, dmass, dpres, p_mid, pres, &
               dsv_dt_k, dsv_ds_k, dsv_dt_sfc, dsv_ds_sfc, h_sum, &
               ustar, absf, idecay, mech_in, b0, ctke_sfc, &
               obl_it, n_its, min_mld, max_mld, mld_guess, mld_found, &
               dmld_min, dmld_max, have_min, have_max, &
               mstar_val, mech_tke, conv_perel, forcing_clip, &
               ki, ka, kb, htot, z_int, pres_int, mld_output, &
               sfc_connected, sfc_disconnect, &
               hp_a, dpe_t_a, dpe_s_a, dch_t_a, dch_s_a, &
               th_a, sh_a, te_lag, se_lag, kddt_prev, kddt_cur, &
               exp_kh, nstar_fc, tot_tke, dt_h, tke_here, vstar, &
               h_ka, h_kb, hb_hs, shape_fn, hbs, mixlen, kd_g0, &
               hp_b, th_b, sh_b, hps, bdt1, dt_c, ds_c, &
               pec_core, colht_core, dkddt, pe_max, &
               kd_val, tke_used, frac_bl, pe_g0, dpe_conv, &
               b1, c1, r_reduc, te_new, se_new, surf_scale, &
               lt_u10, lt_ustokes, lt_kphil, la_val, &
               d_wind, d_conv, d_forcing, d_mixing, d_mdecay, d_cdecay, &
               q_t_kin, q_s_kin, &
               i0_col, d_top_sw, d_bot_sw, heat_sw1, heat_sw2, q_nonpen_kin, &
               ctke_sw_kb, r_sw)

         ! ---- Column prep: T0/S0 + PE/steric weights (downward) ----
         ! (E3) The top of this column.  0 Pa at a free surface; the
         ! ice-shelf load + surface pressure `ms%p_top(i,j)` under a lid,
         ! when `&ocean_psurf_nml in_eos` is set.  Everything below
         ! accumulates `g*rho_0*h` downward from here, so `p_mid` is a
         ! true per-layer hydrostatic pressure either way — which is the
         ! test the `p_top` seam contract applies to a joining builder.
         ! It reaches BOTH consumers of the stack (the in-situ EOS
         ! argument and the PE weight `dmass*p_mid*dsv`), because they
         ! are the same pressure: the load a layer's centre of mass has
         ! to lift, ice included.
         pres = 0.0_wp
         if (p_top_in_eos) pres = ms%p_top(i, j)
         h_sum = 0.0_wp
         dsv_dt_sfc = 0.0_wp
         dsv_ds_sfc = 0.0_wp
         ! (PR-21) penetrating-SW column irradiance + running depth of the
         ! layer TOP below the free surface (0 at k=nz, grows downward to
         ! the opaque bed at k=1).  `i0_col` mirrors the deposition
         ! kernel's `I0 = sw_pen_frac·sw_src·wet_mask` exactly.
         if (sw_ctke_active) then
            i0_col = sw_pen_frac*sw_src_field(i, j)*wet_mask_field(i, j)
         else
            i0_col = 0.0_wp
         end if
         d_top_sw = 0.0_wp
         do k = nz, 1, -1
            hk = ms%h_layer(i, j, k)
            if (hk > 0.0_wp) then
               inv_h = 1.0_wp/(hk + H_NEGLECT)
               t0k = hT(i, j, k)*inv_h
               s0k = hS(i, j, k)*inv_h
            else
               t0k = 0.0_wp
               s0k = 0.0_wp
            end if
            dmass = this%rho0*hk
            dpres = GRAVITY*dmass
            p_mid = pres + 0.5_wp*dpres
            call eos_specvol_derivs(this%eos, t0k, s0k, p_mid, dsv_dt_k, dsv_ds_k)
            this%t0%data(i, j, k) = t0k
            this%s0%data(i, j, k) = s0k
            this%dpe_t%data(i, j, k) = dmass*p_mid*dsv_dt_k
            this%dpe_s%data(i, j, k) = dmass*p_mid*dsv_ds_k
            this%dcolht_t%data(i, j, k) = dmass*dsv_dt_k
            this%dcolht_s%data(i, j, k) = dmass*dsv_ds_k
            ! (PR-21) Per-layer penetrating-SW TKE cost.  Absorbed heat
            ! per band (degC·m) is the two-band difference form with an
            ! opaque bed at k = k_bot — the first LIVE layer counting up,
            ! `1` off z_fixed; the inert bed fillers below it absorb (and
            ! cost) nothing — (identical to the deposition kernel, so
            ! Σ_k Σ_n heat = I0·dt/(rho0·cp) exactly).  The in-layer PE
            ! cost of homogenising that exponentially-distributed heating
            ! is `Phi(h/zeta) <= 1` of the skin cost; the skin identity
            ! `rho0²·h·dsv_dt ≡ rho0·dcolht_t` makes this division-free.
            if (sw_ctke_active) then
               d_bot_sw = d_top_sw + hk
               if (k > ms%k_bot(i, j)) then
                  heat_sw1 = i0_col*sw_R*inv_rho0_cp*dt* &
                             (exp(-d_top_sw/sw_zeta1) - exp(-d_bot_sw/sw_zeta1))
                  heat_sw2 = i0_col*(1.0_wp - sw_R)*inv_rho0_cp*dt* &
                             (exp(-d_top_sw/sw_zeta2) - exp(-d_bot_sw/sw_zeta2))
               else
                  ! opaque bed: absorb the whole remaining irradiance.
                  heat_sw1 = i0_col*sw_R*inv_rho0_cp*dt*exp(-d_top_sw/sw_zeta1)
                  heat_sw2 = i0_col*(1.0_wp - sw_R)*inv_rho0_cp*dt*exp(-d_top_sw/sw_zeta2)
               end if
               this%ctke_sw%data(i, j, k) = &
                  -0.5_wp*GRAVITY*this%rho0*this%dcolht_t%data(i, j, k)* &
                  (heat_sw1*sw_pe_cost_shape(hk/sw_zeta1) + &
                   heat_sw2*sw_pe_cost_shape(hk/sw_zeta2))
               ! Below the opaque bed row: nothing arrives.  `k_bot ≡ 1`
               ! off z_fixed, so this never fires there (bit-identical).
               if (k < ms%k_bot(i, j)) this%ctke_sw%data(i, j, k) = 0.0_wp
               d_top_sw = d_bot_sw
            end if
            if (k == nz) then
               dsv_dt_sfc = dsv_dt_k
               dsv_ds_sfc = dsv_ds_k
            end if
            pres = pres + dpres
            h_sum = h_sum + hk
         end do
         h_sum = h_sum + H_NEGLECT

         ! ---- Surface forcing energetics ----
         ! Kinematic fluxes at this column:
         !   q_T_kin = Q_heat(i,j) / (rho0·cp)  [degC·m/s]
         !   q_S_kin = Q_salt(i,j) / rho0         [PSU·m/s]
         q_t_kin = inv_rho0_cp*Q_heat_field(i, j)
         q_s_kin = inv_rho0*Q_salt_field(i, j)
         ! B0 = g rho0 (dSV_dT q_T + dSV_dS q_S); > 0 stabilizing.
         ! B0 drives mstar / Langmuir (surface-buoyancy scaling); it uses
         ! the FULL net heat flux, unchanged by this PR (MOM6 does the
         ! same — penetrating SW enters ONLY through the per-layer TKE
         ! ledger, not the mstar buoyancy scale).
         b0 = GRAVITY*this%rho0*(dsv_dt_sfc*q_t_kin + dsv_ds_sfc*q_s_kin)
         this%b0(i, j) = b0   ! persist for the Bodner MLE convective velocity scale
         ! cTKE_sfc: PE released (> 0) / required (< 0) to homogenize
         ! the freshly applied SKIN fluxes through the surface layer.
         ! (PR-21) When penetrating SW is charged to the ledger, the skin
         ! carries only the NON-penetrating heat `q_nonpen = q_heat - I0`;
         ! the SW absorbed within the surface layer is the separately
         ! computed `ctke_sw(nz)` (Phi-weighted, <= the all-skin charge).
         ! I0 = 0 ⇒ q_nonpen = q_heat and ctke_sw(nz) = 0 ⇒ the else
         ! branch is character-for-character the legacy line.
         if (sw_ctke_active) then
            q_nonpen_kin = inv_rho0_cp*(Q_heat_field(i, j) - i0_col)
            ctke_sfc = -0.5_wp*GRAVITY*this%rho0**2*ms%h_layer(i, j, nz)* &
                       (q_nonpen_kin*dt*dsv_dt_sfc + q_s_kin*dt*dsv_ds_sfc) &
                       + this%ctke_sw%data(i, j, nz)
         else
            ctke_sfc = -0.5_wp*GRAVITY*this%rho0**2*ms%h_layer(i, j, nz)* &
                       (q_t_kin*dt*dsv_dt_sfc + q_s_kin*dt*dsv_ds_sfc)
         end if

         ! PR-12 dedup: |tau| at cell centres is a shared field
         ! (ocean_surface_stress_set_derived) — the inner sqrt(tau_xc^2 +
         ! tau_yc^2) below IS stress_mag, computed with the identical FP
         ! op order, so this substitution is bit-identical (§7.5).
         ! Phase 4b: `stress_shelf` adds the ICE-SHELF base stress, which
         ! is not in `tau` (the cover mask zeroes the wind there).  The
         ! supports are disjoint, the field is always allocated, and it
         ! is the zero array without a cavity — `x + 0.0` is `x`.  See
         ! the contract in `rdb_ocean_surface_stress`.
         ustar = max(sqrt((ss%stress_mag(i, j) + ss%stress_shelf(i, j))/ &
                          this%rho0), this%ustar_min)
         absf = sqrt((1.0_wp - this%omega_frac)*this%f_centre(i, j)**2 + &
                     this%omega_frac*4.0_wp*this%omega**2)
         idecay = this%tke_decay*absf/ustar
         mech_in = dt*this%rho0*ustar**3

         ! LF17 wave state is BLD-independent: compute once per
         ! column; only the surface-layer average inside the MLD
         ! iteration depends on the guess.
         la_val = 0.0_wp
         lt_u10 = 0.0_wp
         lt_ustokes = 0.0_wp
         lt_kphil = 0.0_wp
         if (this%use_lt) then
            call epbl_lf17_wave_state(ustar, this%rho0, lt_u10, lt_ustokes, lt_kphil)
         end if

         d_wind = 0.0_wp
         d_conv = 0.0_wp
         d_forcing = 0.0_wp
         d_mixing = 0.0_wp
         d_mdecay = 0.0_wp
         d_cdecay = 0.0_wp
         mld_found = 0.0_wp

         if (ms%wet_mask(i, j) <= 0.0_wp .or. h_sum <= 2.0_wp*H_NEGLECT) then
            ! Dry / land column: no mixing.
            do k = 1, nz + 1
               this%kd_int(i, j, k) = 0.0_wp
            end do
         else

            ! ---- Outer MLD iteration ----
            min_mld = 0.0_wp
            max_mld = h_sum
            dmld_min = 0.0_wp
            dmld_max = 0.0_wp
            have_min = .false.
            have_max = .false.
            mld_guess = 0.5_wp*(min_mld + max_mld)
            if (this%mld_use_prev_guess .and. this%mld(i, j) > 0.0_wp) then
               mld_guess = min(this%mld(i, j), max_mld)
            end if
            n_its = 1
            if (this%mld_iteration) n_its = this%mld_max_its

            do obl_it = 1, n_its
               ! Budget ledger restarts each iteration; the last
               ! iteration's values are the ones reported.
               d_wind = 0.0_wp
               d_conv = 0.0_wp
               d_forcing = 0.0_wp
               d_mixing = 0.0_wp
               d_mdecay = 0.0_wp
               d_cdecay = 0.0_wp

               ! (A) mstar at the current MLD guess.
               call epbl_find_mstar(this%mstar_scheme, this%mstar_const, &
                                    this%mstar_cap, this%mstar_coef1, &
                                    this%c_ek, this%mstar_conv_adj, &
                                    this%rh18_cn1, this%rh18_cn2, this%rh18_cn3, &
                                    this%rh18_cs1, this%rh18_cs2, &
                                    b0, ustar, mld_guess, absf, mstar_val)
               if (this%use_lt) then
                  la_val = epbl_lf17_la(ustar, this%la_frac_hbl*mld_guess, &
                                        lt_ustokes, lt_kphil)
                  call epbl_lt_enhance(this%lt_scheme, this%lt_enhance_coef, &
                                       this%lt_enhance_exp, this%lt_max_enhance, &
                                       this%von_karman, &
                                       this%lt_lac1, this%lt_lac2, this%lt_lac3, &
                                       this%lt_lac4, this%lt_lac5, &
                                       la_val, b0, ustar, mld_guess, absf, mstar_val)
               end if
               mech_tke = mstar_val*mech_in
               d_wind = mech_tke

               ! (B) seed the reservoirs from the surface forcing.
               if (ctke_sfc <= 0.0_wp) then
                  forcing_clip = max(ctke_sfc, -mech_tke)
                  mech_tke = mech_tke + forcing_clip
                  conv_perel = 0.0_wp
                  d_forcing = forcing_clip
               else
                  conv_perel = ctke_sfc
                  d_conv = this%nstar*ctke_sfc
               end if

               ! (C) sweep initialization at the surface layer.
               h_ka = ms%h_layer(i, j, nz) + H_NEGLECT
               hp_a = h_ka
               dpe_t_a = this%dpe_t%data(i, j, nz)
               dpe_s_a = this%dpe_s%data(i, j, nz)
               dch_t_a = this%dcolht_t%data(i, j, nz)
               dch_s_a = this%dcolht_s%data(i, j, nz)
               th_a = h_ka*this%t0%data(i, j, nz)
               sh_a = h_ka*this%s0%data(i, j, nz)
               te_lag = 0.0_wp
               se_lag = 0.0_wp
               kddt_prev = 0.0_wp
               htot = ms%h_layer(i, j, nz)
               z_int = ms%h_layer(i, j, nz)
               pres_int = GRAVITY*this%rho0*ms%h_layer(i, j, nz)
               mld_output = ms%h_layer(i, j, nz)
               sfc_connected = .true.
               this%kd_int(i, j, nz + 1) = 0.0_wp

               ! (D) downward sweep over interfaces Ki = nz .. 2.
               ! Layer above the interface: ka = Ki; below: kb = Ki-1.
               do ki = nz, 2, -1
                  ka = ki
                  kb = ki - 1
                  h_ka = ms%h_layer(i, j, ka) + H_NEGLECT
                  h_kb = ms%h_layer(i, j, kb) + H_NEGLECT
                  sfc_disconnect = .false.

                  ! (1) mechanical TKE decays across the layer above
                  !     (Ekman-scale e-folding; no decay at f = 0).
                  exp_kh = exp(-h_ka*idecay)
                  d_mdecay = d_mdecay + (1.0_wp - exp_kh)*mech_tke
                  mech_tke = mech_tke*exp_kh

                  ! (2) per-layer convective forcing accrual (PR-21).
                  !     The surface layer's cTKE seeds the reservoirs in
                  !     (B); here the sub-surface layer `kb` accrues its
                  !     penetrating-SW cost.  A POSITIVE ctke_sw(kb) (SW
                  !     cooling: only reachable on the legacy net-heat
                  !     path at night) releases convective PE — into
                  !     conv_perel, posted to the convective ledger like
                  !     the surface seed.  A NEGATIVE ctke_sw(kb) (solar
                  !     heating: the physical case) is a TKE SINK and
                  !     drains the reservoirs in (3b), after tot_tke is
                  !     formed.
                  if (sw_ctke_active) then
                     ctke_sw_kb = this%ctke_sw%data(i, j, kb)
                     if (ctke_sw_kb > 0.0_wp) then
                        conv_perel = conv_perel + ctke_sw_kb
                        d_conv = d_conv + this%nstar*ctke_sw_kb
                     end if
                  else
                     ctke_sw_kb = 0.0_wp
                  end if

                  ! (3) rotation-reduced convective efficiency.
                  nstar_fc = this%nstar
                  if (conv_perel > 0.0_wp .and. absf > 0.0_wp) then
                     nstar_fc = this%nstar*conv_perel/ &
                                (conv_perel + 0.2_wp* &
                                 sqrt(0.5_wp*dt*this%rho0*(absf*htot)**3*conv_perel))
                  end if
                  tot_tke = mech_tke + nstar_fc*conv_perel

                  ! (3b) penetrating-SW TKE drain (PR-21).  A negative
                  !      ctke_sw(kb) (solar heating stratifies below the
                  !      interface) must be paid before any mixing here —
                  !      discard that TKE to homogenise the SW through the
                  !      next denser cell.  Mechanical + convective
                  !      reservoirs drain PROPORTIONATELY (Reichl &
                  !      Hallberg 2018).  The consumed energy is booked to
                  !      d_mixing (the PE-raising work the SW stratification
                  !      demands), and the nstar-vs-nstar_fc slop of the
                  !      drained convective reservoir to d_cdecay — mirroring
                  !      the proven (8b) closed-form drain below, so the
                  !      column TKE budget still closes to round-off.  (Note:
                  !      posting to d_mixing OR d_forcing balances the ledger;
                  !      posting to BOTH double-counts.)
                  if (sw_ctke_active) then
                     if (ctke_sw_kb < 0.0_wp) then
                        if (ctke_sw_kb + tot_tke < 0.0_wp) then
                           ! SW cost exhausts all the TKE at this interface.
                           d_mixing = d_mixing + tot_tke
                           d_cdecay = d_cdecay + (this%nstar - nstar_fc)*conv_perel
                           tot_tke = 0.0_wp
                           mech_tke = 0.0_wp
                           conv_perel = 0.0_wp
                        else
                           r_sw = (tot_tke + ctke_sw_kb)/tot_tke
                           d_mixing = d_mixing - ctke_sw_kb
                           d_cdecay = d_cdecay + &
                                      (1.0_wp - r_sw)*(this%nstar - nstar_fc)*conv_perel
                           tot_tke = r_sw*tot_tke
                           mech_tke = r_sw*mech_tke
                           conv_perel = r_sw*conv_perel
                        end if
                     end if
                  end if

                  ! (5) static-stability short-circuit: no energy and
                  !     a stable interface => no mixing here.
                  if (tot_tke <= 0.0_wp .and. &
                      0.0_wp <= (this%dcolht_t%data(i, j, kb) + this%dcolht_t%data(i, j, ka))* &
                      (this%t0%data(i, j, ka) - this%t0%data(i, j, kb)) + &
                      (this%dcolht_s%data(i, j, kb) + this%dcolht_s%data(i, j, ka))* &
                      (this%s0%data(i, j, ka) - this%s0%data(i, j, kb))) then
                     kd_val = 0.0_wp
                     sfc_disconnect = .true.
                  else
                     ! (6) velocity scale, mixing length, first-guess Kd.
                     dt_h = dt/max(0.5_wp*(h_ka + h_kb), 1.0e-15_wp*h_sum)
                     tke_here = mech_tke + this%wstar_ustar_coef*conv_perel
                     hb_hs = (h_sum - z_int)/h_sum
                     shape_fn = epbl_mixlen_shape(z_int, mld_guess, &
                                                  this%translay_scale, &
                                                  this%mixlen_exponent, &
                                                  this%mld_iteration)
                     hbs = min(hb_hs, shape_fn)
                     if (tke_here > 0.0_wp) then
                        if (this%vstar_scheme == EPBL_VSTAR_RH18) then
                           surf_scale = max(0.05_wp, 1.0_wp - htot/mld_guess)
                           vstar = this%vstar_scale_fac*surf_scale* &
                                   (this%vstar_surf_fac*ustar + &
                                    (this%wstar_ustar_coef*conv_perel/ &
                                     (dt*this%rho0))**(1.0_wp/3.0_wp))
                        else
                           vstar = this%vstar_scale_fac* &
                                   (tke_here/(dt*this%rho0))**(1.0_wp/3.0_wp)
                        end if
                        if (this%mld_iteration) then
                           mixlen = max(this%min_mix_len, &
                                        (htot*hbs*vstar)/ &
                                        (this%ekman_scale_coef*absf*htot*hbs + vstar))
                           kd_g0 = vstar*this%von_karman*mixlen
                        else
                           kd_g0 = vstar*this%von_karman*(htot*hbs*vstar)/ &
                                   (this%ekman_scale_coef*absf*htot*hbs + vstar)
                        end if
                     else
                        vstar = 0.0_wp
                        kd_g0 = 0.0_wp
                     end if

                     ! (7) pivot quantities for the layer below.
                     hp_b = h_kb
                     th_b = h_kb*this%t0%data(i, j, kb)
                     sh_b = h_kb*this%s0%data(i, j, kb)

                     ! (8) closed-form energy solve (direct path).
                     hps = hp_a + hp_b
                     bdt1 = hp_a*hp_b
                     dt_c = hp_a*th_b - hp_b*th_a
                     ds_c = hp_a*sh_b - hp_b*sh_a
                     pec_core = hp_b*(dpe_t_a*dt_c + dpe_s_a*ds_c) - &
                                hp_a*(this%dpe_t%data(i, j, kb)*dt_c + &
                                      this%dpe_s%data(i, j, kb)*ds_c)
                     colht_core = hp_b*(dch_t_a*dt_c + dch_s_a*ds_c) - &
                                  hp_a*(this%dcolht_t%data(i, j, kb)*dt_c + &
                                        this%dcolht_s%data(i, j, kb)*ds_c)
                     ! Gravity-wave radiation correction: a shrinking
                     ! column radiates energy that cannot drive mixing.
                     if (colht_core < 0.0_wp) then
                        pec_core = pec_core - pres_int*colht_core
                     end if
                     dkddt = kd_g0*dt_h
                     pe_max = pec_core/(bdt1*hps)

                     if (pe_max < 0.0_wp) then
                        ! (8a) convectively unstable: mixing RELEASES
                        ! PE.  Recompute vstar with the released
                        ! energy included; Kd from the mixing length
                        ! (not energy-limited); bank the release.
                        tke_here = mech_tke + &
                                   this%wstar_ustar_coef*(conv_perel - pe_max)
                        if (tke_here > 0.0_wp) then
                           if (this%vstar_scheme == EPBL_VSTAR_RH18) then
                              surf_scale = max(0.05_wp, 1.0_wp - htot/mld_guess)
                              vstar = this%vstar_scale_fac*surf_scale* &
                                      (this%vstar_surf_fac*ustar + &
                                       (this%wstar_ustar_coef*conv_perel/ &
                                        (dt*this%rho0))**(1.0_wp/3.0_wp))
                           else
                              vstar = this%vstar_scale_fac* &
                                      (tke_here/(dt*this%rho0))**(1.0_wp/3.0_wp)
                           end if
                           if (this%mld_iteration) then
                              mixlen = max(this%min_mix_len, &
                                           (htot*hbs*vstar)/ &
                                           (this%ekman_scale_coef*absf*htot*hbs + vstar))
                              kd_val = vstar*this%von_karman*mixlen
                           else
                              kd_val = vstar*this%von_karman*(htot*hbs*vstar)/ &
                                       (this%ekman_scale_coef*absf*htot*hbs + vstar)
                           end if
                        else
                           vstar = 0.0_wp
                           kd_val = 0.0_wp
                        end if
                        pe_g0 = pec_core*dkddt/(bdt1*(bdt1 + dkddt*hps))
                        dpe_conv = pec_core*(kd_val*dt_h)/ &
                                   (bdt1*(bdt1 + (kd_val*dt_h)*hps))
                        if (dpe_conv > 0.0_wp) then
                           kd_val = kd_g0
                           dpe_conv = pe_g0
                        end if
                        ! dpe_conv < 0 => the reservoir grows; the
                        ! reservoirs are NOT proportionally drained on
                        ! this branch.
                        conv_perel = conv_perel - dpe_conv
                        d_conv = d_conv - this%nstar*dpe_conv
                        if (sfc_connected) then
                           mld_output = mld_output + ms%h_layer(i, j, kb)
                        end if
                     else
                        ! (8b) stable: direct closed-form Kd from the
                        ! energy budget; drain the reservoirs.
                        if ((pec_core*dkddt <= &
                             tot_tke*(bdt1*(bdt1 + dkddt*hps))) .or. &
                            (pec_core <= 0.0_wp)) then
                           kd_val = kd_g0
                           tke_used = pec_core*dkddt/(bdt1*(bdt1 + dkddt*hps))
                           frac_bl = 1.0_wp
                        else
                           kd_val = (bdt1**2*tot_tke)/ &
                                    (dt_h*(pec_core - bdt1*hps*tot_tke))
                           tke_used = tot_tke
                           frac_bl = tot_tke*(bdt1*(bdt1 + dkddt*hps))/ &
                                     (pec_core*dkddt)
                        end if
                        if (sfc_connected) then
                           mld_output = mld_output + frac_bl*ms%h_layer(i, j, kb)
                        end if
                        if (frac_bl < 1.0_wp) sfc_disconnect = .true.
                        r_reduc = 0.0_wp
                        if (tot_tke > 0.0_wp .and. tot_tke > tke_used) then
                           r_reduc = (tot_tke - tke_used)/tot_tke
                        end if
                        d_mixing = d_mixing + tke_used
                        d_cdecay = d_cdecay + &
                                   (1.0_wp - r_reduc)*(this%nstar - nstar_fc)*conv_perel
                        mech_tke = r_reduc*mech_tke
                        conv_perel = r_reduc*conv_perel
                     end if
                  end if

                  this%kd_int(i, j, ki) = kd_val
                  kddt_cur = kd_val*dt_h

                  ! (9) advance the embedded forward elimination.
                  b1 = 1.0_wp/(hp_a + kddt_cur)
                  c1 = kddt_cur*b1
                  if (ki == nz) then
                     te_lag = b1*(h_ka*this%t0%data(i, j, ka))
                     se_lag = b1*(h_ka*this%s0%data(i, j, ka))
                  else
                     te_new = b1*(h_ka*this%t0%data(i, j, ka) + kddt_prev*te_lag)
                     se_new = b1*(h_ka*this%s0%data(i, j, ka) + kddt_prev*se_lag)
                     te_lag = te_new
                     se_lag = se_new
                  end if
                  hp_a = h_kb + (hp_a*b1)*kddt_cur
                  dpe_t_a = this%dpe_t%data(i, j, kb) + c1*dpe_t_a
                  dpe_s_a = this%dpe_s%data(i, j, kb) + c1*dpe_s_a
                  dch_t_a = this%dcolht_t%data(i, j, kb) + c1*dch_t_a
                  dch_s_a = this%dcolht_s%data(i, j, kb) + c1*dch_s_a
                  th_a = h_kb*this%t0%data(i, j, kb) + kddt_cur*te_lag
                  sh_a = h_kb*this%s0%data(i, j, kb) + kddt_cur*se_lag
                  kddt_prev = kddt_cur
                  if (sfc_disconnect) then
                     htot = ms%h_layer(i, j, kb)
                     sfc_connected = .false.
                  else
                     htot = htot + ms%h_layer(i, j, kb)
                  end if
                  z_int = z_int + ms%h_layer(i, j, kb)
                  pres_int = pres_int + GRAVITY*this%rho0*ms%h_layer(i, j, kb)
               end do
               this%kd_int(i, j, 1) = 0.0_wp

               ! Leftover stocks at the bed count as dissipated.
               d_mdecay = d_mdecay + mech_tke
               d_cdecay = d_cdecay + this%nstar*conv_perel

               mld_found = mld_output
               if (.not. this%mld_iteration) exit
               if (abs(mld_found - mld_guess) < this%mld_tol) exit
               if (obl_it == n_its) exit

               ! Bracket update + next guess (false position with a
               ! bisection fallback; bisection-only when configured).
               if (mld_found > mld_guess) then
                  min_mld = mld_guess
                  dmld_min = mld_found - mld_guess
                  have_min = .true.
               else
                  max_mld = mld_guess
                  dmld_max = mld_found - mld_guess
                  have_max = .true.
               end if
               if (this%mld_bisection) then
                  mld_guess = 0.5_wp*(min_mld + max_mld)
               else if (have_min .and. have_max .and. obl_it > 2 .and. &
                        mod(obl_it - 1, 4) > 0) then
                  mld_guess = min_mld + dmld_min*(max_mld - min_mld)/ &
                              (dmld_min - dmld_max)
               else if (mld_found > min_mld .and. mld_found < max_mld) then
                  mld_guess = mld_found
               else
                  mld_guess = 0.5_wp*(min_mld + max_mld)
               end if
            end do

         end if

         this%mld(i, j) = mld_found
         this%la(i, j) = la_val
         if (this%tke_diags) then
            this%tke_wind(i, j) = d_wind/dt
            this%tke_conv(i, j) = d_conv/dt
            this%tke_forcing(i, j) = d_forcing/dt
            this%tke_mixing(i, j) = d_mixing/dt
            this%tke_mech_decay(i, j) = d_mdecay/dt
            this%tke_conv_decay(i, j) = d_cdecay/dt
         end if
      end do
   end subroutine epbl_column_kernel