vmix_split_ddiff_impl Subroutine

   pure subroutine vmix_split_ddiff_impl(nx, ny, nzp1, kt, ks, temp_h, salt_h, &
                                         h_layer, alpha_T, beta_S, strat_param_max, &
                                         kappa_s, exp1, exp2, param1, param2, param3, &
                                         mol_diff, use_k90)
      !! Double-diffusion split.  Replaces `ks := kt` with the asymmetric
      !!   ks = kt_pre + kd_extra_s ;  kt = kt_pre + kd_extra_t
      !! both from the SAME pre-split kt (read into `kt_pre` before either
      !! write).  CVMix `cvmix_coeffs_ddiff` algebra (Large et al. 1994
      !! fingering; Marmorino-Caldwell 1976 / Kelley 1990 convection).
      !!
      !! Interior interfaces k=2..nz only (boundary interfaces k=1/nzp1
      !! carry no fingering and stay at kt's closed-BC zero, reproducing
      !! `ks`'s boundary zeros -- same all-k structure as the plain split).
      !! Interface k sits at the bottom of layer k, between layer k-1
      !! (lower) and layer k (upper, toward-surface): bottom-up convention.
      !!
      !! Branch on the SIGNED alpha*dT / beta*dS (never a pre-divided
      !! R_rho) so R_rho<=0, R_rho->inf and 0/0 never arise; R_rho is
      !! formed only inside a branch whose denominator is provably nonzero.
      !! No floor/ceiling here -- `vmix_assemble` is the single gate.
      !! Explicit-shape dummies (per the "vmix incident" descriptor-walk
      !! rule); linear-EOS constant alpha_T/beta_S passed by value.
      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)
         !! Temperature thickness-integral hTr = T*h_layer (degC*m).
      real(wp), intent(in) :: salt_h(nx, ny, nzp1 - 1)
         !! Salinity thickness-integral hTr = S*h_layer (PSU*m).
      real(wp), intent(in) :: h_layer(nx, ny, nzp1 - 1)
      real(wp), intent(in) :: alpha_T, beta_S
         !! |d rho/dT| and d rho/dS (kg/m^3 per degC / per PSU).
      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

      nz = nzp1 - 1

      do concurrent(k=1:nzp1, j=1:ny, i=1:nx) &
         local(kt_pre, adT, bdS, rrho, ddiff, kd_t, kd_s, hu, hl)
         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
            if (hu > H_VANISHED .and. hl > H_VANISHED) then
               ! alpha*dT and beta*dS across interface k (upper - lower)
               adT = alpha_T*(temp_h(i, j, k)/hu - temp_h(i, j, k - 1)/hl)
               bdS = beta_S*(salt_h(i, j, k)/hu - salt_h(i, j, k - 1)/hl)

               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 subroutine vmix_split_ddiff_impl