cavity_two_equation Subroutine

   pure subroutine cavity_two_equation(T_w, S_w, p_b, gamma_t, ice, eos, const, &
                                       T_b, S_b, m_mass, ierr)
      !! Two-equation variant: the interface salinity is the FAR-FIELD
      !! salinity, `S_b = S_w` exactly, i.e. the `gamma_s -> infinity`
      !! limit of the three-equation form (NOT the `gamma_s -> 0` limit,
      !! which sends `T_b -> T_w` and the melt rate to zero).  Holland &
      !! Jenkins (1999) section 2b(2) p. 1791; Jenkins, Nicholls & Corr
      !! (2010) eq. (6) p. 2302, where the single transfer coefficient is
      !! `Gamma_TS ~ 0.006` and the freezing point is evaluated at the
      !! far-field salinity.
      !!
      !! **IDENTITY WITH THE SHIPPED SEA-ICE BASAL FLUX.**  This form, at
      !! `gamma_t*dt = h`, IS `rdb_ice_basal_flux`'s
      !! "relax the top layer to the freezing point in one thermo step"
      !! law: that kernel forms
      !! `fb = rho*c_p*max(0, SST - T_f)*h/dt_therm`, and setting
      !! `gamma_t = h/dt` in `q_ocean = rho_w*c_w*gamma_t*(T_w - T_b)`
      !! with `T_b = T_f(S_w, p)` reproduces it term for term.  It is an
      !! exact identity, not an asymptote — the two differences are that
      !! the sea-ice kernel clamps the supercooled branch away with a
      !! `max(0, .)` (frazil owns that side) and evaluates the liquidus at
      !! `p = 0`.  The practical reading: the shipped sea-ice law already
      !! has an implied exchange velocity, and it is resolution- and
      !! timestep-dependent by construction, which is exactly the
      !! far-field-sampling problem the cavity literature is about.
      !$acc routine seq
      real(wp), intent(in) :: T_w
         !! Far-field temperature (degC).
      real(wp), intent(in) :: S_w
         !! Far-field salinity (g/kg).
      real(wp), intent(in) :: p_b
         !! Interface pressure (Pa).
      real(wp), intent(in) :: gamma_t
         !! Heat exchange velocity (m/s), >= 0.
      type(ocean_cavity_ice_t), intent(in) :: ice
         !! Ice-conduction bundle.
      type(eos_t), intent(in) :: eos
         !! Shared EOS handle — the liquidus.
      type(ocean_cavity_const_t), intent(in) :: const
         !! Constants bundle.
      real(wp), intent(out) :: T_b
         !! Interface temperature (degC) = `T_f(S_w, p_b)`.
      real(wp), intent(out) :: S_b
         !! Interface salinity (g/kg) = `S_w` exactly.
      real(wp), intent(out) :: m_mass
         !! Melt mass flux (kg/m^2/s), > 0 melting.
      integer, intent(out) :: ierr
         !! `CAVITY_MELT_*` status.
      real(wp) :: T_ice, c_i_eff, kh, l_eff, q_oc, q_ic
      logical :: melting

      call cavity_safe_state(eos, S_w, p_b, T_b, S_b, m_mass)
      T_ice = cavity_t_ice(ice)

      ierr = CAVITY_MELT_OK
      if (.not. (ieee_is_finite(T_w) .and. ieee_is_finite(S_w) .and. &
                 ieee_is_finite(p_b) .and. ieee_is_finite(gamma_t) .and. &
                 ieee_is_finite(T_ice))) then
         ierr = CAVITY_MELT_NONFINITE_INPUT
         return
      end if
      if (gamma_t < 0.0_wp) then
         ierr = CAVITY_MELT_BAD_INPUT
         return
      end if

      melting = (T_w - T_b) > 0.0_wp
      call cavity_ice_terms(ice, melting, const, c_i_eff, kh, ierr)
      if (ierr /= CAVITY_MELT_OK) return

      l_eff = const%L_f + c_i_eff*(T_b - T_ice)
      if (l_eff <= 0.0_wp) then
         ierr = CAVITY_MELT_NO_PHYSICAL_ROOT
         return
      end if
      q_oc = const%rho_w*const%c_w*gamma_t*(T_w - T_b)
      q_ic = kh*(T_b - T_ice)
      m_mass = (q_oc - q_ic)/l_eff
      if (.not. ieee_is_finite(m_mass)) then
         call cavity_safe_state(eos, S_w, p_b, T_b, S_b, m_mass)
         ierr = CAVITY_MELT_NONFINITE_STATE
      end if
   end subroutine cavity_two_equation