vmix_split_ddiff_eos_impl Subroutine

   pure subroutine vmix_split_ddiff_eos_impl(nx, ny, nzp1, kt, ks, temp_h, salt_h, &
                                             h_layer, p_top, eos, rho0, p_top_in_eos, &
                                             strat_param_max, kappa_s, exp1, exp2, &
                                             param1, param2, param3, mol_diff, use_k90)
      !! `buoyancy_coeffs = "eos"` twin of `vmix_split_ddiff_impl` — the
      !! SAME CVMix closed forms, the same branch structure, the same
      !! outputs; the only change is where `α` and `β` come from.
      !!
      !! The constant version forms `adT = α·ΔT` and `bdS = β·ΔS` with one
      !! scalar pair for the whole domain.  Here both are evaluated at the
      !! INTERFACE's own state — `eos_buoyancy_coeffs` at the mean of the
      !! two abutting layer (T, S) and at the in-situ hydrostatic pressure
      !! there.  The stratification parameter the closure actually
      !! branches on is the density ratio `R_ρ = α·ΔT / β·ΔS`, and under a
      !! nonlinear EOS α varies by a factor of several between a 25 degC
      !! surface and a −1.9 degC cavity, and again with depth, so a single
      !! α can put an interface in the wrong REGIME (fingering vs
      !! diffusive convection), not merely off by a coefficient.
      !!
      !! **Pressure.** Seeded from `p_top` (the E3 surface-load seam) and
      !! accumulated DOWNWARD as `g·ρ₀·h`, exactly the way EPBL seeds and
      !! walks its `p_mid` stack — a true per-interface hydrostatic
      !! pressure, which is the test the `p_top` seam contract applies to
      !! any joining builder.  `eos%p_ref` deliberately does NOT enter: an
      !! interior interface has a real depth of its own, and adding a
      !! potential-density reference on top of it would double-count.
      !! (KPP's `B_0` is the other way round — it is a SURFACE flux with
      !! no depth of its own, so it falls back to `p_ref`.)
      !!
      !! **Why the loop nest differs from the constant twin.** The
      !! pressure is a running sum down the column, so `k` cannot be part
      !! of the concurrent index set.  The constant path keeps its fully
      !! collapsed `(i, j, k)` launch untouched; this one is
      !! `do concurrent(j, i)` with a serial surface→bed `k` sweep, the
      !! same shape as `epbl_column_kernel` and `ks_solve_column`.
      !!
      !! Interface `k` sits at the BOTTOM of layer `k`, between layer `k`
      !! (upper, surfaceward — `hu`) and layer `k-1` (lower, bedward —
      !! `hl`); `k = nzp1` is the free surface and `k = 1` the bed, and
      !! both are left at `kd = 0` exactly as the constant twin leaves
      !! them, which is what preserves the closed-BC invariant.
      !!
      !! Explicit-shape dummies throughout (no assumed-shape in a
      !! `do concurrent`); `eos` is a flat POD passed by value, so the
      !! device copy is register-resident and `eos_buoyancy_coeffs` is
      !! reachable as an `!$acc routine seq` from the same shared library.
      integer, intent(in) :: nx, ny, nzp1
      real(wp), intent(inout) :: kt(nx, ny, nzp1)
      real(wp), intent(inout) :: ks(nx, ny, nzp1)
      real(wp), intent(in) :: temp_h(nx, ny, nzp1 - 1)
         !! `hTr` of temperature (degC·m).
      real(wp), intent(in) :: salt_h(nx, ny, nzp1 - 1)
         !! `hTr` of salinity (PSU·m).
      real(wp), intent(in) :: h_layer(nx, ny, nzp1 - 1)
      real(wp), intent(in) :: p_top(nx, ny)
         !! Surface load (Pa) — `multilayer_state_t%p_top`, the zero array
         !! unless `&ocean_psurf_nml in_eos`.
      type(eos_t), intent(in) :: eos
         !! Active EOS handle, by value.
      real(wp), intent(in) :: rho0
         !! Boussinesq reference density for the hydrostatic accumulation
         !! (the one configured ρ₀ of record, via `vmix%rho0`).
      logical, intent(in) :: p_top_in_eos
         !! Whether to seed the stack from `p_top` — mirrors EPBL's gate.
      real(wp), intent(in) :: strat_param_max, kappa_s, exp1, exp2
      real(wp), intent(in) :: param1, param2, param3, mol_diff
      logical, intent(in) :: use_k90

      integer :: i, j, k, nz
      real(wp) :: kt_pre, adT, bdS, rrho, ddiff, kd_t, kd_s, hu, hl
      real(wp) :: p_int, t_u, t_l, s_u, s_l, alpha_i, beta_i

      nz = nzp1 - 1

      do concurrent(j=1:ny, i=1:nx) &
         local(k, kt_pre, adT, bdS, rrho, ddiff, kd_t, kd_s, hu, hl, &
               p_int, t_u, t_l, s_u, s_l, alpha_i, beta_i)
         p_int = 0.0_wp
         if (p_top_in_eos) p_int = p_top(i, j)

         do k = nzp1, 1, -1
            ! Walk down from the free surface: stepping from interface
            ! k+1 to interface k crosses layer k, so charge its weight
            ! first.  At k = nzp1 nothing has been crossed yet and p_int
            ! is still the surface load.
            if (k <= nz) p_int = p_int + GRAVITY*rho0*h_layer(i, j, k)

            kt_pre = kt(i, j, k)
            kd_t = 0.0_wp
            kd_s = 0.0_wp

            if (k >= 2 .and. k <= nz) then
               hu = h_layer(i, j, k)      ! upper layer (toward surface)
               hl = h_layer(i, j, k - 1)  ! lower layer
               ! vanished-ok: a PAIRWISE gate — the double-diffusive term needs BOTH layers
               ! live and is skipped otherwise; a per-layer `rdb_vl_conc` would
               ! silently make the difference finite.
               if (hu > H_VANISHED .and. hl > H_VANISHED) then
                  t_u = temp_h(i, j, k)/hu
                  t_l = temp_h(i, j, k - 1)/hl
                  s_u = salt_h(i, j, k)/hu
                  s_l = salt_h(i, j, k - 1)/hl
                  ! α, β at the interface state: the mean of the two
                  ! abutting layer centres (the same two-point average
                  ! the ΔT / ΔS below are differences of), at the in-situ
                  ! interface pressure.
                  call eos_buoyancy_coeffs(eos, 0.5_wp*(t_u + t_l), &
                                           0.5_wp*(s_u + s_l), p_int, &
                                           alpha_i, beta_i)
                  ! alpha*dT and beta*dS across interface k (upper - lower)
                  adT = alpha_i*(t_u - t_l)
                  bdS = beta_i*(s_u - s_l)

                  if (adT >= bdS .and. bdS > 0.0_wp) then
                     ! ---- salt fingering (R_rho >= 1) ----
                     rrho = adT/bdS
                     if (rrho < strat_param_max) then
                        ddiff = (1.0_wp - ((rrho - 1.0_wp)/ &
                                           (strat_param_max - 1.0_wp))**exp1)**exp2
                        kd_s = kappa_s*ddiff
                     end if
                     kd_t = 0.7_wp*kd_s
                  else if (adT >= bdS .and. adT < 0.0_wp) then
                     ! ---- diffusive convection (0 < R_rho < 1) ----
                     rrho = adT/bdS
                     if (use_k90) then
                        ddiff = mol_diff*8.7_wp*rrho**1.1_wp
                     else
                        ddiff = mol_diff*param1* &
                                exp(param2*exp(param3*(1.0_wp/rrho - 1.0_wp)))
                     end if
                     kd_t = ddiff
                     if (rrho < 0.5_wp) then
                        kd_s = 0.15_wp*rrho*ddiff
                     else
                        kd_s = (1.85_wp*rrho - 0.85_wp)*ddiff
                     end if
                  end if
               end if
            end if

            ks(i, j, k) = kt_pre + kd_s
            kt(i, j, k) = kt_pre + kd_t
         end do
      end do
   end subroutine vmix_split_ddiff_eos_impl