cavity_three_equation Subroutine

   pure subroutine cavity_three_equation(T_w, S_w, p_b, gamma_t, gamma_s, S_i, &
                                         ice, eos, const, T_b, S_b, m_mass, ierr)
      !! Closed-form solve of (E1)-(E3) on the linear liquidus carried by
      !! `eos`.  Returns the interface state and the canonical melt mass
      !! flux (kg/m^2/s, > 0 melting).  See the derivation block above for
      !! the root selection, the cancellation-safe quadratic and the
      !! pre-solve melt/freeze branch.
      !!
      !! `gamma_s > 0` is REQUIRED: the `gamma_s -> infinity` limit is the
      !! two-equation form and has its own entry point
      !! (`cavity_two_equation`), while `gamma_s = 0` is not a limit of
      !! this system at all.
      !$acc routine seq
      real(wp), intent(in) :: T_w
         !! Far-field temperature (degC).
      real(wp), intent(in) :: S_w
         !! Far-field salinity (g/kg).  Must exceed `S_i`.
      real(wp), intent(in) :: p_b
         !! Interface pressure (Pa) — `g*rho_i*h_ice`, the undiluted ice
         !! load.
      real(wp), intent(in) :: gamma_t
         !! Heat exchange velocity (m/s), >= 0.
      real(wp), intent(in) :: gamma_s
         !! Salt exchange velocity (m/s), > 0.
      real(wp), intent(in) :: S_i
         !! Ice salinity (g/kg), >= 0.  Holland & Jenkins (1999) p. 1789
         !! treats marine ice as fresh ("we can treat S_I as zero
         !! always"); the argument is kept because the root-bracketing
         !! proof is stated for general `S_i < S_w`.
      type(ocean_cavity_ice_t), intent(in) :: ice
         !! Ice-conduction bundle.
      type(eos_t), intent(in) :: eos
         !! Shared EOS handle — the liquidus.  By value, flat POD.
      type(ocean_cavity_const_t), intent(in) :: const
         !! Constants bundle.
      real(wp), intent(out) :: T_b
         !! Interface temperature (degC).
      real(wp), intent(out) :: S_b
         !! Interface salinity (g/kg).
      real(wp), intent(out) :: m_mass
         !! Melt mass flux (kg/m^2/s), > 0 melting.
      integer, intent(out) :: ierr
         !! `CAVITY_MELT_*` status.  On anything but OK the three outputs
         !! are the safe state (see `cavity_safe_state`).
      real(wp) :: T_ice, T_star, c_i_eff, kh
      real(wp) :: a, b, d_w, e0, e1, rc, h0, h1, p1, qa, qb, qc, disc, sq, q
      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(gamma_s) .and. ieee_is_finite(S_i) .and. &
                 ieee_is_finite(T_ice))) then
         ierr = CAVITY_MELT_NONFINITE_INPUT
         return
      end if
      if (gamma_t < 0.0_wp .or. gamma_s <= 0.0_wp .or. S_i < 0.0_wp) then
         ierr = CAVITY_MELT_BAD_INPUT
         return
      end if
      if (.not. (S_w > S_i)) then
         ! The whole root-bracketing argument rests on `f(S_i) > 0`, which
         ! needs `S_w > S_i`.
         ierr = CAVITY_MELT_BAD_INPUT
         return
      end if

      T_star = T_w - eos_freezing_point(eos, S_w, p_b)
      melting = T_star > 0.0_wp
      call cavity_ice_terms(ice, melting, const, c_i_eff, kh, ierr)
      if (ierr /= CAVITY_MELT_OK) return

      a = eos%tfr_s
      b = eos%tfr_0 + eos%tfr_p*p_b
      d_w = T_w - b
      e0 = const%L_f + c_i_eff*(b - T_ice)
      e1 = c_i_eff*a
      rc = const%rho_w*const%c_w*gamma_t
      h0 = rc*d_w - kh*(b - T_ice)
      h1 = -a*(rc + kh)
      p1 = const%rho_w*gamma_s

      qa = -p1*e1 - h1
      qb = p1*(S_w*e1 - e0) - h0 + h1*S_i
      qc = p1*S_w*e0 + h0*S_i

      if (qa == 0.0_wp) then
         ! Degenerate: a salinity-independent liquidus (`lambda1 = 0`), or
         ! the pathological `c_w*gamma_t + kh/rho_w == c_i*gamma_s`.  The
         ! system is then linear in `S_b`.
         if (qb == 0.0_wp) then
            ierr = CAVITY_MELT_NO_PHYSICAL_ROOT
            return
         end if
         S_b = -qc/qb
      else
         disc = qb*qb - 4.0_wp*qa*qc
         if (.not. ieee_is_finite(disc)) then
            ierr = CAVITY_MELT_NONFINITE_STATE
            return
         end if
         if (disc < 0.0_wp) then
            ! Cannot happen for a well-posed set (`S_i` lies between the
            ! roots, so `disc > 0`).  Treat it as corruption and refuse —
            ! NEVER clamp it to zero, which would manufacture a double
            ! root and a plausible melt rate out of a broken column.
            ierr = CAVITY_MELT_NO_PHYSICAL_ROOT
            return
         end if
         sq = sqrt(disc)
         if (qb >= 0.0_wp) then
            q = -0.5_wp*(qb + sq)
            S_b = q/qa
         else
            q = -0.5_wp*(qb - sq)
            ! `q == 0` requires `B == 0` and `disc == 0`, excluded above.
            S_b = qc/q
         end if
      end if

      if (.not. ieee_is_finite(S_b)) then
         call cavity_safe_state(eos, S_w, p_b, T_b, S_b, m_mass)
         ierr = CAVITY_MELT_NONFINITE_STATE
         return
      end if
      if (S_b <= S_i) then
         call cavity_safe_state(eos, S_w, p_b, T_b, S_b, m_mass)
         ierr = CAVITY_MELT_NO_PHYSICAL_ROOT
         return
      end if

      T_b = eos_freezing_point(eos, S_b, p_b)
      m_mass = const%rho_w*gamma_s*(S_w - S_b)/(S_b - S_i)
      if (.not. (ieee_is_finite(m_mass) .and. ieee_is_finite(T_b))) 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_three_equation