cavity_solve_melt_f Subroutine

private pure subroutine cavity_solve_melt_f(T_w, S_w, p_b, u_star, S_i, par, f_cor, ice, eos, const, sol, ierr)

Solve the three-equation system with any implemented exchange law, and return the full interface state plus the fluxes a coupling seam will consume.

For an EXPLICIT law (const_gamma) this is one call to the closed-form quadratic.

For a STRATIFICATION-DEPENDENT law (hj99, yung25) the system is implicit — the exchange coefficients depend on the interfacial buoyancy flux, which depends on the melt rate they produce — and the outer iteration is BISECTION on x = ln(L+) over [CAVITY_LP_X_LO, CAVITY_LP_X_HI]. It is guaranteed convergent because the residual G(x) = ln(L+_new(x)) - x provably changes sign across that bracket:

  • x -> x_lo (maximal suppression): gamma -> 0, so |B_b| -> 0 and L+_new -> +infinity, giving G > 0;
  • x -> x_hi (neutral): gamma sits at its cap, |B_b| is maximal and L+_new is finite, giving G < 0.

UNLESS the interface is DESTABILISING (B_b >= 0, i.e. freezing or a strongly cooling interface), in which case L+_new is +infinity everywhere and the neutral limit IS the fixed point — short-circuited, not iterated. That is the branch KINK: H&J99 sets eta* := 1 there (p. 1792) and Yung et al. (2025) falls back to the ConstCoeff values (pp. 5832-5833). The melt rate is continuous across it; its derivative is not, which is exactly why this is bisection and not Newton.

After convergence the whole column is re-evaluated ONCE at the converged scalar, so the returned exchange velocities and the returned (T_b, S_b, m_mass) belong to ONE value of L+.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: T_w

Far-field temperature (degC).
real(kind=wp), intent(in) :: S_w

Far-field salinity (g/kg). Must exceed S_i.
real(kind=wp), intent(in) :: p_b

Interface pressure (Pa).
real(kind=wp), intent(in) :: u_star

Friction velocity (m/s), strictly positive — from cavity_ustar.
real(kind=wp), intent(in) :: S_i

Ice salinity (g/kg), >= 0.
type(ocean_cavity_exchange_t), intent(in) :: par

Exchange-law bundle (its f_cor member is NOT read).
real(kind=wp), intent(in) :: f_cor

Coriolis parameter (1/s) for this column; hj99 only. See cavity_exchange_velocities_f for why it is not par%f_cor.
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.
type(ocean_cavity_solution_t), intent(out) :: sol

Interface state, fluxes and solver diagnostics. On any non-OK status this carries the safe state (m_mass exactly zero).
integer, intent(out) :: ierr

CAVITY_MELT_* status.

Calls

proc~~cavity_solve_melt_f~~CallsGraph proc~cavity_solve_melt_f cavity_solve_melt_f proc~cavity_heat_fluxes cavity_heat_fluxes proc~cavity_solve_melt_f->proc~cavity_heat_fluxes proc~cavity_law_is_implicit cavity_law_is_implicit proc~cavity_solve_melt_f->proc~cavity_law_is_implicit proc~cavity_obukhov_length cavity_obukhov_length proc~cavity_solve_melt_f->proc~cavity_obukhov_length proc~cavity_outer_residual cavity_outer_residual proc~cavity_solve_melt_f->proc~cavity_outer_residual proc~cavity_safe_state cavity_safe_state proc~cavity_solve_melt_f->proc~cavity_safe_state proc~cavity_solution_reset cavity_solution_reset proc~cavity_solve_melt_f->proc~cavity_solution_reset proc~cavity_state_at_x cavity_state_at_x proc~cavity_solve_melt_f->proc~cavity_state_at_x proc~eos_freezing_point eos_freezing_point proc~cavity_solve_melt_f->proc~eos_freezing_point proc~cavity_ice_terms cavity_ice_terms proc~cavity_heat_fluxes->proc~cavity_ice_terms proc~cavity_t_ice cavity_t_ice proc~cavity_heat_fluxes->proc~cavity_t_ice proc~cavity_l_plus_is_neutral cavity_l_plus_is_neutral proc~cavity_outer_residual->proc~cavity_l_plus_is_neutral proc~cavity_safe_state->proc~eos_freezing_point proc~cavity_state_at_x->proc~cavity_safe_state proc~cavity_buoyancy_flux cavity_buoyancy_flux proc~cavity_state_at_x->proc~cavity_buoyancy_flux proc~cavity_exchange_velocities_f cavity_exchange_velocities_f proc~cavity_state_at_x->proc~cavity_exchange_velocities_f proc~cavity_l_plus_from_state cavity_l_plus_from_state proc~cavity_state_at_x->proc~cavity_l_plus_from_state proc~cavity_three_equation cavity_three_equation proc~cavity_state_at_x->proc~cavity_three_equation proc~cavity_gamma_hj99 cavity_gamma_hj99 proc~cavity_exchange_velocities_f->proc~cavity_gamma_hj99 proc~cavity_gamma_yung25 cavity_gamma_yung25 proc~cavity_exchange_velocities_f->proc~cavity_gamma_yung25 proc~cavity_three_equation->proc~cavity_safe_state proc~cavity_three_equation->proc~eos_freezing_point proc~cavity_three_equation->proc~cavity_ice_terms proc~cavity_three_equation->proc~cavity_t_ice proc~cavity_gamma_hj99->proc~cavity_l_plus_is_neutral proc~cavity_gamma_yung25->proc~cavity_l_plus_is_neutral

Called by

proc~~cavity_solve_melt_f~~CalledByGraph proc~cavity_solve_melt_f cavity_solve_melt_f proc~cavity_melt_point_gamma_f cavity_melt_point_gamma_f proc~cavity_melt_point_gamma_f->proc~cavity_solve_melt_f proc~cavity_solve_melt cavity_solve_melt proc~cavity_solve_melt->proc~cavity_solve_melt_f proc~cavity_melt_columns_2d cavity_melt_columns_2d proc~cavity_melt_columns_2d->proc~cavity_melt_point_gamma_f proc~cavity_melt_point cavity_melt_point proc~cavity_melt_point->proc~cavity_solve_melt proc~cavity_melt_point_gamma cavity_melt_point_gamma proc~cavity_melt_point_gamma->proc~cavity_melt_point_gamma_f proc~cavity_melt_columns cavity_melt_columns proc~cavity_melt_columns->proc~cavity_melt_point proc~ocean_cavity_flux_step ocean_cavity_flux_step proc~ocean_cavity_flux_step->proc~cavity_melt_columns_2d proc~engine_step_finalize engine_step_finalize proc~engine_step_finalize->proc~ocean_cavity_flux_step proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_step_finalize proc~rdb_ocean_step rdb_ocean_step proc~rdb_ocean_step->proc~engine_step_finalize

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: S_b
real(kind=wp), private :: T_b
real(kind=wp), private :: b_flux
logical, private :: converged
real(kind=wp), private :: g_hi
real(kind=wp), private :: g_lo
real(kind=wp), private :: g_mid
real(kind=wp), private :: gamma_s
real(kind=wp), private :: gamma_t
real(kind=wp), private :: hi
integer, private :: it
real(kind=wp), private :: l_plus
real(kind=wp), private :: lo
real(kind=wp), private :: lp_new
real(kind=wp), private :: m_mass
real(kind=wp), private :: mid
logical, private :: settled
real(kind=wp), private :: xstar

Source Code

   pure subroutine cavity_solve_melt_f(T_w, S_w, p_b, u_star, S_i, par, f_cor, ice, eos, &
                                       const, sol, ierr)
      !! Solve the three-equation system with any implemented exchange
      !! law, and return the full interface state plus the fluxes a
      !! coupling seam will consume.
      !!
      !! For an EXPLICIT law (`const_gamma`) this is one call to the
      !! closed-form quadratic.
      !!
      !! For a STRATIFICATION-DEPENDENT law (`hj99`, `yung25`) the system
      !! is implicit — the exchange coefficients depend on the interfacial
      !! buoyancy flux, which depends on the melt rate they produce — and
      !! the outer iteration is BISECTION on `x = ln(L+)` over
      !! `[CAVITY_LP_X_LO, CAVITY_LP_X_HI]`.  It is guaranteed convergent
      !! because the residual `G(x) = ln(L+_new(x)) - x` provably changes
      !! sign across that bracket:
      !!
      !!   * `x -> x_lo` (maximal suppression): `gamma -> 0`, so
      !!     `|B_b| -> 0` and `L+_new -> +infinity`, giving `G > 0`;
      !!   * `x -> x_hi` (neutral): `gamma` sits at its cap, `|B_b|` is
      !!     maximal and `L+_new` is finite, giving `G < 0`.
      !!
      !! UNLESS the interface is DESTABILISING (`B_b >= 0`, i.e. freezing
      !! or a strongly cooling interface), in which case `L+_new` is
      !! `+infinity` everywhere and the neutral limit IS the fixed point —
      !! short-circuited, not iterated.  That is the branch KINK: H&J99
      !! sets `eta* := 1` there (p. 1792) and Yung et al. (2025) falls
      !! back to the ConstCoeff values (pp. 5832-5833).  The melt rate is
      !! continuous across it; its derivative is not, which is exactly why
      !! this is bisection and not Newton.
      !!
      !! After convergence the whole column is re-evaluated ONCE at the
      !! converged scalar, so the returned exchange velocities and the
      !! returned `(T_b, S_b, m_mass)` belong to ONE value of `L+`.
      !$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).
      real(wp), intent(in) :: u_star
         !! Friction velocity (m/s), strictly positive — from
         !! `cavity_ustar`.
      real(wp), intent(in) :: S_i
         !! Ice salinity (g/kg), >= 0.
      type(ocean_cavity_exchange_t), intent(in) :: par
         !! Exchange-law bundle (its `f_cor` member is NOT read).
      real(wp), intent(in) :: f_cor
         !! Coriolis parameter (1/s) for this column; `hj99` only.  See
         !! `cavity_exchange_velocities_f` for why it is not `par%f_cor`.
      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.
      type(ocean_cavity_solution_t), intent(out) :: sol
         !! Interface state, fluxes and solver diagnostics.  On any non-OK
         !! status this carries the safe state (`m_mass` exactly zero).
      integer, intent(out) :: ierr
         !! `CAVITY_MELT_*` status.
      real(wp) :: T_b, S_b, m_mass, gamma_t, gamma_s, b_flux, lp_new, l_plus
      real(wp) :: lo, hi, mid, xstar, g_lo, g_hi, g_mid
      integer :: it
      logical :: converged, settled

      ierr = CAVITY_MELT_OK
      call cavity_solution_reset(sol)
      sol%u_star = u_star
      sol%l_plus = CAVITY_L_PLUS_NEUTRAL
      sol%L_obukhov = CAVITY_L_PLUS_NEUTRAL
      call cavity_safe_state(eos, S_w, p_b, sol%T_b, sol%S_b, sol%m_mass)

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

      converged = .true.
      it = 0
      l_plus = CAVITY_L_PLUS_NEUTRAL

      ! The neutral evaluation is both the answer for an explicit law and
      ! the upper bracket for an implicit one.
      call cavity_state_at_x(CAVITY_LP_X_HI, par, f_cor, ice, eos, const, u_star, T_w, S_w, &
                             p_b, S_i, T_b, S_b, m_mass, gamma_t, gamma_s, b_flux, &
                             lp_new, ierr)
      if (ierr /= CAVITY_MELT_OK) then
         call cavity_safe_state(eos, S_w, p_b, sol%T_b, sol%S_b, sol%m_mass)
         return
      end if
      l_plus = lp_new

      if (cavity_law_is_implicit(par%law)) then
         g_hi = cavity_outer_residual(lp_new, CAVITY_LP_X_HI)
         ! `cavity_l_plus_is_neutral(lp_new)` is the destabilising
         ! short-circuit and makes `g_hi = +huge >= 0`; the explicit
         ! `g_hi >= 0` test below therefore covers both accept branches.
         if (g_hi < 0.0_wp) then
            call cavity_state_at_x(CAVITY_LP_X_LO, par, f_cor, ice, eos, const, u_star, T_w, &
                                   S_w, p_b, S_i, T_b, S_b, m_mass, gamma_t, gamma_s, &
                                   b_flux, lp_new, ierr)
            if (ierr /= CAVITY_MELT_OK) then
               call cavity_safe_state(eos, S_w, p_b, sol%T_b, sol%S_b, sol%m_mass)
               return
            end if
            g_lo = cavity_outer_residual(lp_new, CAVITY_LP_X_LO)
            if (g_lo <= 0.0_wp) then
               ! The bracket argument failed — refuse rather than return
               ! whichever end happens to look plausible.
               call cavity_safe_state(eos, S_w, p_b, sol%T_b, sol%S_b, sol%m_mass)
               ierr = CAVITY_MELT_NOT_CONVERGED
               return
            end if

            lo = CAVITY_LP_X_LO
            hi = CAVITY_LP_X_HI
            converged = .false.
            do it = 1, CAVITY_MAX_ITER
               mid = 0.5_wp*(lo + hi)
               settled = (mid <= lo) .or. (mid >= hi)
               if (.not. settled) then
                  call cavity_state_at_x(mid, par, f_cor, ice, eos, const, u_star, T_w, S_w, &
                                         p_b, S_i, T_b, S_b, m_mass, gamma_t, gamma_s, &
                                         b_flux, lp_new, ierr)
                  if (ierr /= CAVITY_MELT_OK) then
                     call cavity_safe_state(eos, S_w, p_b, sol%T_b, sol%S_b, sol%m_mass)
                     return
                  end if
                  g_mid = cavity_outer_residual(lp_new, mid)
                  if (g_mid == 0.0_wp) then
                     settled = .true.
                  else if (g_mid > 0.0_wp) then
                     lo = mid
                  else
                     hi = mid
                  end if
                  if ((.not. settled) .and. (hi - lo < CAVITY_LP_X_TOL)) settled = .true.
               end if
               if (settled) then
                  converged = .true.
                  exit
               end if
            end do
            if (.not. converged) then
               call cavity_safe_state(eos, S_w, p_b, sol%T_b, sol%S_b, sol%m_mass)
               ierr = CAVITY_MELT_NOT_CONVERGED
               return
            end if

            ! One final evaluation at the converged scalar, so the
            ! returned gammas and the returned interface state belong to
            ! the SAME `L+`.
            xstar = 0.5_wp*(lo + hi)
            call cavity_state_at_x(xstar, par, f_cor, ice, eos, const, u_star, T_w, S_w, p_b, &
                                   S_i, T_b, S_b, m_mass, gamma_t, gamma_s, b_flux, &
                                   lp_new, ierr)
            if (ierr /= CAVITY_MELT_OK) then
               call cavity_safe_state(eos, S_w, p_b, sol%T_b, sol%S_b, sol%m_mass)
               return
            end if
            l_plus = exp(xstar)
         end if
      end if

      sol%T_b = T_b
      sol%S_b = S_b
      sol%m_mass = m_mass
      sol%gamma_t = gamma_t
      sol%gamma_s = gamma_s
      sol%b_flux = b_flux
      sol%l_plus = l_plus
      sol%L_obukhov = cavity_obukhov_length(const, u_star, b_flux)
      sol%T_star = T_w - eos_freezing_point(eos, S_w, p_b)
      sol%S_star = S_w - S_b
      call cavity_heat_fluxes(T_w, T_b, m_mass, gamma_t, ice, const, &
                              sol%q_ocean, sol%q_ice, sol%q_latent)
      sol%n_iter = it
      sol%converged = converged
   end subroutine cavity_solve_melt_f