rdb_ocean_cavity_melt Module


Uses

  • module~~rdb_ocean_cavity_melt~~UsesGraph module~rdb_ocean_cavity_melt rdb_ocean_cavity_melt ieee_arithmetic ieee_arithmetic module~rdb_ocean_cavity_melt->ieee_arithmetic module~rdb_constants rdb_constants module~rdb_ocean_cavity_melt->module~rdb_constants module~rdb_eos rdb_eos module~rdb_ocean_cavity_melt->module~rdb_eos pic_types pic_types module~rdb_constants->pic_types module~rdb_eos->module~rdb_constants module~rdb_grid rdb_grid module~rdb_eos->module~rdb_grid module~rdb_grid->module~rdb_constants

Used by

  • module~~rdb_ocean_cavity_melt~~UsedByGraph module~rdb_ocean_cavity_melt rdb_ocean_cavity_melt module~rdb_ocean_cavity_flux rdb_ocean_cavity_flux module~rdb_ocean_cavity_flux->module~rdb_ocean_cavity_melt module~rdb_ocean_setup rdb_ocean_setup module~rdb_ocean_setup->module~rdb_ocean_cavity_melt module~rdb_ocean_dyn rdb_ocean_dyn module~rdb_ocean_setup->module~rdb_ocean_dyn module~rdb_ocean_state rdb_ocean_state module~rdb_ocean_setup->module~rdb_ocean_state proc~configure_ocean_cavity_melt configure_ocean_cavity_melt proc~configure_ocean_cavity_melt->module~rdb_ocean_cavity_melt proc~validate_config validate_config proc~validate_config->module~rdb_ocean_cavity_melt module~rdb_ocean_dyn->module~rdb_ocean_cavity_flux module~rdb_ocean_engine rdb_ocean_engine module~rdb_ocean_engine->module~rdb_ocean_cavity_flux module~rdb_ocean_engine->module~rdb_ocean_setup module~rdb_ocean_engine->module~rdb_ocean_dyn module~rdb_ocean_engine->module~rdb_ocean_state module~rdb_ocean_diag_derived rdb_ocean_diag_derived module~rdb_ocean_engine->module~rdb_ocean_diag_derived module~rdb_ocean_diag_fills rdb_ocean_diag_fills module~rdb_ocean_engine->module~rdb_ocean_diag_fills module~rdb_ocean_state->module~rdb_ocean_cavity_flux module~rdb_ocean_state->module~rdb_ocean_dyn module~rdb_driver rdb_driver module~rdb_driver->module~rdb_ocean_dyn module~rdb_driver->module~rdb_ocean_engine module~rdb_driver->module~rdb_ocean_state module~rdb_handle rdb_handle module~rdb_handle->module~rdb_ocean_engine module~rdb_handle->module~rdb_ocean_state module~rdb_ocean_api rdb_ocean_api module~rdb_ocean_api->module~rdb_ocean_dyn module~rdb_ocean_api->module~rdb_ocean_engine module~rdb_ocean_api->module~rdb_handle module~rdb_ocean_api->module~rdb_ocean_diag_derived module~rdb_ocean_api->module~rdb_ocean_diag_fills module~rdb_ocean_diag_derived->module~rdb_ocean_state module~rdb_ocean_diag_derived->module~rdb_ocean_diag_fills module~rdb_ocean_diag_fills->module~rdb_ocean_state

Variables

Type Visibility Attributes Name Initial
real(kind=wp), public, parameter :: CAVITY_CD_ISOMIP = 2.5e-3_wp

ISOMIP+ top drag coefficient C_D,top — Asay-Davis et al. (2016) Table 4 p. 2483. The least constrained number in the whole subject: the literature spans 1.5e-3 (Holland & Jenkins 1999 Table 1 p. 1790) to 9.7e-3 (Jenkins, Nicholls & Corr 2010 Table 2 p. 2309) — a factor of 6.5 — and Yung et al. (2025) p. 5831 reports order-of-magnitude variation within a single crevasse. Must be a namelist knob.
integer, public, parameter :: CAVITY_FW_INVALID = 0

Unrecognised freshwater string.
integer, public, parameter :: CAVITY_FW_MASS = 2

Real Boussinesq volume source on the top layer, dh = m*dt/rho_0.
integer, public, parameter :: CAVITY_FW_VIRTUAL = 1

Fixed column mass; the dilution is emulated by the exact fixed-mass equivalent salt flux -m*(S_far - s_ice). Default.
real(kind=wp), public, parameter :: CAVITY_GAMMA_RATIO_ISOMIP = 35.0_wp

Gamma_T/Gamma_S, the ratio the ISOMIP+ protocol adopts — Asay-Davis et al. (2016) Table 4 p. 2483. The 35 is Jenkins, Nicholls & Corr (2010) p. 2309: the ratio “should lie somewhere in the range 35-70. Adopting a value at the lower end of this range…”. Named because the coupling layer resolves the &ocean_cavity_melt_nml gamma_s “unset” sentinel through it, and a second literal 35 in a second file is how the two drift.
real(kind=wp), public, parameter :: CAVITY_GAMMA_S_ISOMIP = CAVITY_GAMMA_T_ISOMIP/CAVITY_GAMMA_RATIO_ISOMIP

ISOMIP+ salt-transfer coefficient, Gamma_S = Gamma_T/35 — Asay-Davis et al. (2016) Table 4 p. 2483. The 35 is Jenkins, Nicholls & Corr (2010) p. 2309: the ratio “should lie somewhere in the range 35-70. Adopting a value at the lower end of this range…”.
real(kind=wp), public, parameter :: CAVITY_GAMMA_T_ISOMIP = 2.2e-2_wp

ISOMIP+ heat-transfer coefficient Gamma_T, dimensionless — Asay-Davis et al. (2016) section 3.2.1 p. 2487, which derives it from the Stanton number sqrt(C_D,top)*Gamma_T = 1.1e-3 “suggesting that Gamma_T = 2.2 x 10^-2 might be a good initial guess”. A STARTING GUESS, NOT A CONSTANT OF NATURE: the protocol has participants TUNE it until the Ocean0 mean melt lands in the prescribed band, and Yung et al. (2026) Table 2 p. 2058 shows the twelve ISOMIP+ submissions landing anywhere from 0.011 to 0.2. A namelist knob in the coupling PR, never a hard-wired number.
integer, public, parameter :: CAVITY_ICE_ADV_DIFF = 2

SHIPS. Holland & Jenkins (1999) constant-vertical-advection + diffusion in the linearised eq. (31) p. 1794 form, where the amplification factor is replaced by its asymptote Pi = w_I H_I/kappa_I for MELTING and Pi = 0 for FREEZING. Substituting (31) into (26)+(6) cancels H_I and kappa_I identically and collapses the whole ice column to

q_ice = m_mass*c_i*(T_b - T_ice) (melting), 0 (freezing)

i.e. the melt flux must additionally warm the incorporated ice from T_ice to T_b. It is LINEAR in m_mass, so the closed-form quadratic survives — it merely replaces L_f by an effective latent heat L_eff = L_f + c_i*(T_b - T_ice). The ice thickness and diffusivity are NOT needed in this mode. The zeroing on freezing IS what the paper prescribes (eq. 31), not an approximation added here: p. 1796, “The model with constant vertical heat advection in the ice shelf has no effect unless the mixed layer is warmer than the freezing point.”
integer, public, parameter :: CAVITY_ICE_DIFFUSIVE = 3

RESERVED. Holland & Jenkins (1999) section 2d(2) eq. (21) p. 1793, q_ice = k_ice*(T_b - T_ice)/h_ice: steady, no advection, linear ice profile. Not shipped because it also changes the BRANCH LOGIC — q_ice is then independent of m_mass, so sign(m) = sign(T*) no longer holds and a strictly POSITIVE thermal driving is required for zero melt (their p. 1796: “causes a net shift toward freezing”). Wiring it up means revisiting the pre-solve melt/freeze branch AND adding k_ice/h_ice to the parameter bundles.
integer, public, parameter :: CAVITY_ICE_INSULATING = 1

SHIPS. Perfect insulator, q_ice = 0 — Holland & Jenkins (1999) section 2d(1) p. 1793. PRESCRIBED by the ISOMIP+ protocol (Asay-Davis et al. 2016 Table 4 p. 2483 sets kappa_i = 0, and p. 2485 explicitly instructs participants NOT to use the H&J99 advection-diffusion scheme). T_ice is then unread.
integer, public, parameter :: CAVITY_ICE_INVALID = 0

Unrecognised ice-conduction string.
integer, public, parameter :: CAVITY_LAW_BURCHARD22 = 8

RESERVED. Burchard et al. (2022) resolution-robust log-layer law — the only law in the set that treats the far-field sampling depth as a physical input, and the natural control for a vertical-coordinate study.
integer, public, parameter :: CAVITY_LAW_CONST_GAMMA = 1

SHIPS. gamma = Gamma*u* — Jenkins, Nicholls & Corr (2010) eqs. (1), (2), (5) p. 2300; the ISOMIP+ form, Asay-Davis et al. (2016) eqs. (24), (26) p. 2485. The recommended default.
integer, public, parameter :: CAVITY_LAW_HJ99 = 2

SHIPS. Holland & Jenkins (1999) eqs. (14)-(18) p. 1792 — turbulent + molecular sublayer with the McPhee (1981) stability parameter eta*. Stratification-dependent, hence implicit in the melt rate (outer bisection).
integer, public, parameter :: CAVITY_LAW_INVALID = 0

Unrecognised law string. parse_cavity_exchange_law returns this rather than falling back to a default — a mistyped exchange law is a silent physics change, exactly like a mistyped liquidus.
integer, public, parameter :: CAVITY_LAW_JENKINS21 = 9

RESERVED, and probably permanently: Jenkins (2021) is a 1-D boundary-CURRENT model, not a per-column transfer law. Its only separable piece is the two-equation constant-Stanton law already available as cavity_two_equation.
integer, public, parameter :: CAVITY_LAW_JENKINS91 = 3

RESERVED. Kader & Yaglom smooth-wall form, Holland & Jenkins (1999) eqs. (11)-(12) p. 1792. Explicit in the melt rate, but it carries a free boundary-layer thickness h for which the paper gives no value.
integer, public, parameter :: CAVITY_LAW_MK18 = 7

RESERVED. McConnochie & Kerr convective floor, reachable only through Yung et al. (2025) eqs. (10)-(13) p. 5836. It is the ONLY law that reads the interface state (T_b, S_b) — which is why those two arguments are in cavity_exchange_velocities’ argument list already.
integer, public, parameter :: CAVITY_LAW_ROSEVEAR22 = 5

RESERVED. Rosevear, Gayen & Galton-Fenzi (2022) eqs. (27)-(28) p. 2601. Needs two documented choices the paper does not make (natural vs base-10 logarithm; behaviour below L+ = 2500).
integer, public, parameter :: CAVITY_LAW_VT19 = 6

RESERVED. Vreugdenhil & Taylor (2019) Monin-Obukhov recipe. Carries a free reference height z_inf with no defensible default, and its authors disown their own salt branch.
integer, public, parameter :: CAVITY_LAW_YUNG25 = 4

SHIPS. Yung et al. (2025) “StratFeedback”, eqs. (7)-(8) p. 5832 — two power laws in the viscous Obukhov scale L+, capped at the Vreugdenhil & Taylor (2019) passive-scalar maxima, so it reduces EXACTLY to a constant-Gamma law in the shear-dominated limit. Stratification-dependent (implicit).
real(kind=wp), public, parameter :: CAVITY_L_PLUS_NEUTRAL = huge(1.0_wp)

Sentinel for “neutral / unsuppressed”, i.e. the L+ -> +infinity limit, carried in the same scalar as a real L+ so that the exchange-law seam is ONE argument list. Any l_plus that is non-finite, non-positive or >= CAVITY_L_PLUS_NEUTRAL is treated as neutral (see cavity_l_plus_is_neutral), which is exactly what every stratification-dependent law prescribes for a destabilising or vanishing buoyancy flux.
integer, public, parameter :: CAVITY_MELT_BAD_INPUT = 3

Finite but outside the solver’s domain: a negative drag or friction-velocity floor, a negative exchange velocity, a negative ice salinity, u_star <= 0, gamma_s <= 0 in the three-equation form, or S_w <= S_i (which is what the root-bracketing argument rests on).
integer, public, parameter :: CAVITY_MELT_LAW_DOMAIN = 9

An exchange law was evaluated outside its own domain: a non-positive logarithm argument, a non-positive eta* argument, or a non-positive Gamma_Turb + Gamma_Mole denominator (H&J99’s Gamma_Turb can go strongly negative at very small u*; the molecular terms normally dominate, but a negative exchange velocity is refused rather than returned).
integer, public, parameter :: CAVITY_MELT_LAW_INVALID = 7

The law / ice-mode code is not one of the enum values at all (e.g. CAVITY_LAW_INVALID straight out of a mistyped namelist string). Distinct from NOT_IMPLEMENTED: this one is a typo.
integer, public, parameter :: CAVITY_MELT_NONFINITE_INPUT = 1

An INPUT was NaN or +/-Inf. Every guard fires BEFORE any min/max/clamp, so a non-finite input can never come back out as a plausible clamp bound (the nvfortran relaxed-FP hazard).
integer, public, parameter :: CAVITY_MELT_NONFINITE_STATE = 2

An INTERMEDIATE went non-finite (overflow in the discriminant, a non-finite root or melt rate) from finite inputs.
integer, public, parameter :: CAVITY_MELT_NOT_CONVERGED = 5

The outer stratification bisection did not bracket or did not converge within CAVITY_MAX_ITER.
integer, public, parameter :: CAVITY_MELT_NOT_IMPLEMENTED = 6

A RESERVED exchange law or ice-conduction mode — the enum value exists so the dispatch seam is stable, the physics does not ship yet. See the per-enum !! notes for what the prototype has.
integer, public, parameter :: CAVITY_MELT_NO_CORIOLIS = 8

The Holland & Jenkins (1999) law was asked for at f = 0. Their eq. (15) p. 1792 takes ln(u* xi_N eta*^2 / (|f| h_nu)) and eq. (18) divides by f L_O, so the law simply does not exist on the equator — the same fail-loud stance this repository already takes for the Henyey background-mixing latitude factor on a cartesian grid.
integer, public, parameter :: CAVITY_MELT_NO_PHYSICAL_ROOT = 4

The quadratic has no admissible root: a negative discriminant (impossible for a well-posed set, so it is treated as corruption and never clamped to zero), a fully degenerate A = B = 0, a root that failed S_b > S_i, or a non-positive effective latent heat in the two-equation form.
integer, public, parameter :: CAVITY_MELT_OK = 0

Solved. The outputs are physics.
real(kind=wp), public, parameter :: CAVITY_USTAR_MIN_YUNG25 = 1.0e-4_wp

Friction-velocity floor (m/s) — Yung et al. (2025) eq. (14) p. 5836 with the value from their Table 2 p. 5838. It exists because “a friction velocity of zero (perhaps created by initialising the model at rest) will result in identically zero melt … which would be inconsistent with the presence of heat available for melting” (p. 5836).
real(kind=wp), public, parameter :: CAVITY_U_TIDE_ISOMIP = 1.0e-2_wp

ISOMIP+ RMS tidal velocity u_tidal (m/s) entering the MELT friction velocity — Asay-Davis et al. (2016) Table 4 p. 2483 and eq. (27) p. 2485, after Jenkins, Nicholls & Corr (2010) eq. (10) p. 2309 u*^2 = C_d (U^2 + <U_T^2>). TRAP: the protocol applies it to the melt u* ONLY, not to the momentum drag (p. 2486, “The computation of top and bottom drag do not incorporate utidal”).
integer, public, parameter :: CAVITY_VC_INVALID = 0

Unrecognised volume_compensation string.
integer, public, parameter :: CAVITY_VC_NONE = 1

No compensation; a closed domain gains the melt volume.
integer, public, parameter :: CAVITY_VC_UNIFORM_OPEN = 2

Remove the domain-integrated melt volume again, uniformly per unit area over the wet cells the ice does NOT cover.
real(kind=wp), private, parameter :: CAVITY_LP_X_HI = 46.051701859880914_wp

Upper bracket, ln(1e20), and the NEUTRAL evaluation point: at or above it the trial L+ is taken to be +infinity, every law sits at its own neutral limit, |B_b| is maximal and L+_new is finite, so G < 0.
real(kind=wp), private, parameter :: CAVITY_LP_X_LO = -18.420680743952367_wp

Lower bracket, ln(1e-8). As x -> x_lo the exchange velocities go to zero, so |B_b| -> 0 and the re-diagnosed L+ -> +infinity: the residual G = ln(L+_new) - x is > 0. (Spelled as a literal, not log(1.0e-8_wp), because a transcendental is not a Fortran constant expression.)
real(kind=wp), private, parameter :: CAVITY_LP_X_TOL = 1.0e-13_wp

Bracket width in ln(L+) at which the bisection is declared converged. At x ~ 46 the double-precision spacing is ~7e-15, so this is ~15 representable steps above the floor — tight enough that the returned L+ is good to ~1e-13 relative, loose enough to terminate on every toolchain.
integer, private, parameter :: CAVITY_MAX_ITER = 200

Hard cap on the outer bisection. A 64.5-wide bracket halved 200 times is far beyond exhausting double precision, so the loop normally exits on CAVITY_LP_X_TOL or on mid <= lo; the cap exists so a corrupted residual cannot spin forever on device.
real(kind=wp), private, parameter :: HJ99_H_NU_COEFF = 5.0_wp

Viscous sublayer thickness h_nu = 5*nu/u*, Holland & Jenkins (1999) eq. (17) p. 1792 — their hydraulically-smooth replacement for McPhee, Maykut & Morison (1987) eq. (10) p. 7029’s roughness length z_0.
real(kind=wp), private, parameter :: HJ99_MOLE_OFFSET = 6.0_wp

The - 6 of the same equation. (Malyarenko et al. (2020) Table B.1 records this additive constant appearing as -10.1, -8.68, -9 and -6 across the lineage; we use H&J99’s own rendering, which is what the oracle was generated with.)
real(kind=wp), private, parameter :: HJ99_MOLE_SLOPE = 12.5_wp

Gamma_Mole = 12.5*(Pr,Sc)^(2/3) - 6, attributed by Holland & Jenkins (1999) eq. (16) p. 1792 to Kader & Yaglom (1972).
real(kind=wp), private, parameter :: Y25_A_S = -4.30_wp

Gamma_S = 10^A_S * (L+)^n_S prefactor exponent, eq. (8) p. 5832.
real(kind=wp), private, parameter :: Y25_A_T = -3.21_wp

Gamma_T = 10^A_T * (L+)^n_T prefactor exponent, eq. (7) p. 5832.
real(kind=wp), private, parameter :: Y25_GAMMA_S_CC = 3.9e-4_wp

Constant-coefficient cap on Gamma_S, Yung et al. (2025) Table 1 p. 5832. The forced crossover to the caps is APPROXIMATE: with the published exponents the power laws reach 0.011967 and 3.876e-4 at L+ = 1e4, so the min() actually engages at L+ = 1.04e4 (heat) and 1.13e4 (salt). Immaterial physically; material to any test that asserts equality AT 1e4.
real(kind=wp), private, parameter :: Y25_GAMMA_T_CC = 0.012_wp

Constant-coefficient cap on Gamma_T (Yung et al. 2025 Table 1 p. 5832) — the Vreugdenhil & Taylor (2019) passive-scalar maximum. NOTE this law’s neutral limit is its OWN cap, not the configured Gamma_T; par%gamma_t_coeff is ignored by design.
real(kind=wp), private, parameter :: Y25_N_S = 0.223_wp

Gamma_S power-law slope in L+, eq. (8) p. 5832.
real(kind=wp), private, parameter :: Y25_N_T = 0.322_wp

Gamma_T power-law slope in L+, eq. (7) p. 5832.

Derived Types

type, public ::  ocean_cavity_const_t

Thermodynamic + turbulence constants. All SI.

Read more…

Components

Type Visibility Attributes Name Initial
real(kind=wp), public :: L_f = 3.34e5_wp

Latent heat of fusion (J/kg) — Holland & Jenkins (1999) Table 1 p. 1790; Asay-Davis et al. (2016) Table 4 p. 2483; Yung et al. (2025) Table 1 p. 5832. (Burchard et al. (2022) Table 1 p. 8 uses 3.335e5.)
real(kind=wp), public :: Pr = 13.8_wp

Molecular Prandtl number — Holland & Jenkins (1999) Table 1 p. 1790, confirmed by McPhee, Maykut & Morison (1987) p. 7029.
real(kind=wp), public :: R_c = 0.20_wp

Critical flux Richardson number — Holland & Jenkins (1999) Table 1 p. 1790.
real(kind=wp), public :: Sc = 2432.0_wp

Molecular Schmidt number — same two sources.
real(kind=wp), public :: alpha_T = 3.733e-5_wp

FRACTIONAL thermal expansion coefficient (1/degC) of the ISOMIP+ linear EOS — Asay-Davis et al. (2016) Table 4 p. 2483. Feeds the interfacial buoyancy flux only. NOTE eos_t%alpha_T in this repository is DIMENSIONAL (kg/m^3 per degC) = rho0 times this; convert at the coupling seam. OPEN QUESTION, recorded not resolved: near -2 degC the true thermal expansion under a nonlinear EOS is near zero or NEGATIVE, which flips the sign of the (small) temperature term in the buoyancy flux and therefore of the stratification feedback near neutrality. No paper in the set addresses it.
real(kind=wp), public :: beta_S = 7.843e-4_wp

FRACTIONAL haline contraction coefficient (1/(g/kg)) of the ISOMIP+ linear EOS — Asay-Davis et al. (2016) Table 4 p. 2483.
real(kind=wp), public :: c_i = 2009.0_wp

Specific heat capacity of ice (J/kg/K) — Holland & Jenkins (1999) Table 1 p. 1790. Read only by CAVITY_ICE_ADV_DIFF.
real(kind=wp), public :: c_w = 3974.0_wp

Specific heat capacity of seawater (J/kg/K) — unanimous across Holland & Jenkins (1999), Jenkins et al. (2010), Asay-Davis et al. (2016) and Yung et al. (2025).
real(kind=wp), public :: g = 9.81_wp

Gravitational acceleration (m/s^2), as used by the buoyancy flux — Asay-Davis et al. (2016) Table 4 p. 2483.
real(kind=wp), public :: kappa_vk = 0.40_wp

Von Karman constant — Holland & Jenkins (1999) Table 1 p. 1790. Also the kappa of the Obukhov scale.
real(kind=wp), public :: nu = 1.95e-6_wp

Kinematic viscosity of seawater (m^2/s) — Holland & Jenkins (1999) Table 1 p. 1790.
real(kind=wp), public :: rho_fw = 1000.0_wp

Freshwater density (kg/m^3), REPORTING ONLY — the density ISOMIP+ reports its m_w melt rate with (Asay-Davis et al. (2016) eq. (24) p. 2485).
real(kind=wp), public :: rho_i = 918.0_wp

Ice density (kg/m^3), REPORTING ONLY — Asay-Davis et al. (2016) p. 2479 and Yung et al. (2025).
real(kind=wp), public :: rho_w = 1028.0_wp

Seawater density multiplying the turbulent fluxes (kg/m^3) — Asay-Davis et al. (2016) p. 2479. The papers span 1025 (H&J99) to 1030 (Jenkins et al. 2010), a 0.5% spread that lands directly on the melt rate.
real(kind=wp), public :: xi_N = 0.052_wp

McPhee stability constant xi_N — Holland & Jenkins (1999) Table 1 p. 1790; McPhee, Maykut & Morison (1987) p. 7029; Yung et al. (2025) Table A1 p. 5849. (NOT 0.13: no paper in the set contains a zeta_N or that value.)

type, public ::  ocean_cavity_exchange_t

Exchange-law selector + its parameters. ONE bundle for every law, so the dispatch is a single select case and a later do concurrent kernel can call it without reshaping. Members a given law does not read are ignored.

Components

Type Visibility Attributes Name Initial
real(kind=wp), public :: f_cor = -1.4e-4_wp

Coriolis parameter (1/s) — read by CAVITY_LAW_HJ99 only, and only as |f|. Holland & Jenkins (1999) Table 1 p. 1790 prints f = -1.0e-4 (Southern Hemisphere), but their eq. (15) takes ln(.../f h_nu) and eq. (18) needs f L_O > 0: both only make sense with the magnitude. Neither that paper nor McPhee, Maykut & Morison (1987) says |f| — it is inferred here, and recorded as inferred.
real(kind=wp), public :: gamma_s_coeff = CAVITY_GAMMA_S_ISOMIP

Dimensionless Gamma_S.
real(kind=wp), public :: gamma_t_coeff = CAVITY_GAMMA_T_ISOMIP

Dimensionless Gamma_T of gamma_t = Gamma_T*u*. Named *_coeff to keep it distinct from the exchange VELOCITY gamma_t (m/s) the laws return.
integer, public :: law = CAVITY_LAW_CONST_GAMMA

CAVITY_LAW_*. Default = the recommended const_gamma.

type, public ::  ocean_cavity_ice_t

Ice-side conduction selector + its parameter.

Components

Type Visibility Attributes Name Initial
real(kind=wp), public :: T_ice = -25.0_wp

Ice interior / surface temperature (degC), read by CAVITY_ICE_ADV_DIFF only — Holland & Jenkins (1999) Table 1 p. 1790 uses T_S ~ -25.0. Under CAVITY_ICE_INSULATING this member is IGNORED and the kernel substitutes exactly zero, so a stale or absent ice temperature cannot leak into an insulating run.
integer, public :: mode = CAVITY_ICE_INSULATING

CAVITY_ICE_*. Default = insulating, which is what the ISOMIP+ protocol prescribes.

type, public ::  ocean_cavity_solution_t

Everything cavity_solve_melt returns: the interface state, the fluxes a coupling seam needs, and the solver diagnostics.

Read more…

Components

Type Visibility Attributes Name Initial
real(kind=wp), public :: L_obukhov

Dimensional Obukhov length (m), > 0 stabilising.
real(kind=wp), public :: S_b

Interface salinity (g/kg).
real(kind=wp), public :: S_star

Haline driving S_w - S_b (g/kg).
real(kind=wp), public :: T_b

Interface temperature (degC), on the liquidus by construction.
real(kind=wp), public :: T_star

Thermal driving T_w - T_f(S_w, p_b) (degC).
real(kind=wp), public :: b_flux

Interfacial buoyancy flux (m^2/s^3), < 0 stabilising.
logical, public :: converged

Outer iteration converged (always .true. for an explicit law).
real(kind=wp), public :: gamma_s

Salt exchange velocity actually used (m/s).
real(kind=wp), public :: gamma_t

Heat exchange velocity actually used (m/s).
real(kind=wp), public :: l_plus

Viscous Obukhov scale (dimensionless), > 0 stabilising; CAVITY_L_PLUS_NEUTRAL for a vanishing buoyancy flux.
real(kind=wp), public :: m_mass

Canonical melt mass flux (kg/m^2/s of ice), > 0 melting.
integer, public :: n_iter

Outer bisection iterations taken (0 for an explicit law or a destabilising short-circuit).
real(kind=wp), public :: q_ice

Conductive + ice-warming flux interface → ice (W/m^2).
real(kind=wp), public :: q_latent

Latent heat consumed by the phase change (W/m^2).
real(kind=wp), public :: q_ocean

Turbulent heat flux ocean → interface (W/m^2).
real(kind=wp), public :: u_star

Friction velocity the solve was given (m/s).

Functions

public pure function cavity_buoyancy_flux(const, T_w, S_w, T_b, S_b, gamma_t, gamma_s) result(b_flux)

Interfacial buoyancy flux B_b (m^2/s^3), NEGATIVE = stabilising — Yung et al. (2025) eq. (6) p. 5831,

Read more…

Arguments

Type IntentOptional Attributes Name
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle — g, alpha_T, beta_S.
real(kind=wp), intent(in) :: T_w

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

Far-field salinity (g/kg).
real(kind=wp), intent(in) :: T_b

Interface temperature (degC).
real(kind=wp), intent(in) :: S_b

Interface salinity (g/kg).
real(kind=wp), intent(in) :: gamma_t

Heat exchange velocity (m/s).
real(kind=wp), intent(in) :: gamma_s

Salt exchange velocity (m/s).

Return Value real(kind=wp)

public pure function cavity_l_plus_from_state(const, u_star, b_flux) result(l_plus)

Viscous Obukhov scale L+ = L/delta_nu with delta_nu = nu/u*, i.e. L+ = -u*^4/(nu*kappa*B_b) — Yung et al. (2025) eq. (5) p. 5831; the same definition in Vreugdenhil & Taylor (2019) eq. (27) and Rosevear et al. (2022) eqs. (6)+(8) p. 2592. POSITIVE for melting. Returns CAVITY_L_PLUS_NEUTRAL for a vanishing buoyancy flux.

Arguments

Type IntentOptional Attributes Name
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle — nu, kappa_vk.
real(kind=wp), intent(in) :: u_star

Friction velocity (m/s), > 0.
real(kind=wp), intent(in) :: b_flux

Interfacial buoyancy flux (m^2/s^3).

Return Value real(kind=wp)

public pure function cavity_l_plus_is_neutral(l_plus) result(is_neutral)

Is this trial L+ the neutral / unsuppressed limit?

Read more…

Arguments

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

Viscous Obukhov scale, or CAVITY_L_PLUS_NEUTRAL.

Return Value logical

public pure function cavity_m_ice_from_mass(const, m_mass) result(m_ice)

Solid-ice thickness rate (m/s) from the canonical mass flux, m_ice = m_mass/rho_i. REPORTING ONLY — this is Jenkins, Nicholls & Corr (2010)’s a_b convention.

Arguments

Type IntentOptional Attributes Name
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle — rho_i.
real(kind=wp), intent(in) :: m_mass

Melt mass flux (kg/m^2/s).

Return Value real(kind=wp)

public pure function cavity_m_weq_from_mass(const, m_mass) result(m_weq)

Freshwater-equivalent thickness rate (m/s) from the canonical mass flux, m_weq = m_mass/rho_fw. REPORTING ONLY — this is ISOMIP+’s m_w (Asay-Davis et al. (2016) eq. (24) p. 2485), the number their figures are in.

Arguments

Type IntentOptional Attributes Name
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle — rho_fw.
real(kind=wp), intent(in) :: m_mass

Melt mass flux (kg/m^2/s).

Return Value real(kind=wp)

public pure function cavity_obukhov_length(const, u_star, b_flux) result(l_obukhov)

Dimensional Obukhov length L = -u*^3/(kappa*B_b) (m), POSITIVE for a stabilising (melting) buoyancy flux. McPhee, Maykut & Morison (1987) p. 7029; the same scale appears as Yung et al. (2025) eq. (5) p. 5831. (Holland & Jenkins (1999) uses L_O in their eq. (18) but never defines it.)

Read more…

Arguments

Type IntentOptional Attributes Name
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle — kappa_vk.
real(kind=wp), intent(in) :: u_star

Friction velocity (m/s), > 0.
real(kind=wp), intent(in) :: b_flux

Interfacial buoyancy flux (m^2/s^3).

Return Value real(kind=wp)

public pure function parse_cavity_exchange_law(name) result(code)

Translate an exchange-law string into a CAVITY_LAW_* code. RESERVED laws parse successfully — the refusal belongs to cavity_exchange_velocities, which returns CAVITY_MELT_NOT_IMPLEMENTED — so that a typo (CAVITY_LAW_INVALID) and an honest request for unwritten physics stay distinguishable.

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: name

Return Value integer

public pure function parse_cavity_freshwater(name) result(code)

Translate a &ocean_cavity_melt_nml freshwater string into a CAVITY_FW_* code. The accepted set MIRRORS the nml_enum registration in register_ocean_cavity_melt — the two lists move together.

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: name

Return Value integer

public pure function parse_cavity_ice_mode(name) result(code)

Translate an ice-conduction string into a CAVITY_ICE_* code.

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: name

Return Value integer

public pure function parse_cavity_volume_comp(name) result(code)

Translate a &ocean_cavity_melt_nml volume_compensation string into a CAVITY_VC_* code. Mirrors the nml_enum list.

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: name

Return Value integer

private pure function cavity_law_is_implicit(law) result(is_implicit)

Does this law’s (gamma_t, gamma_s) depend on the interfacial buoyancy flux — and therefore on the melt rate it produces? Yung et al. (2025) p. 5833 states the consequence: “Since the transfer coefficients depend on L+, which in turn depends on melt rate via surface buoyancy forcing, iteration is required for convergence of the three-equation parameterisation solution.”

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: law

CAVITY_LAW_* code.

Return Value logical

private pure function cavity_outer_residual(lp_new, x) result(g)

Outer-iteration residual G(x) = ln(L+_new(x)) - x, with the destabilising branch folded in: a non-positive, non-finite or sentinel L+_new means the buoyancy flux at this iterate is destabilising (or zero), which every law treats as L+ = +infinity, so G = +infinity. Mapping it that way keeps the bisection bracket valid instead of taking log() of a negative number.

Arguments

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

L+ re-diagnosed from the state this iterate produced.
real(kind=wp), intent(in) :: x

The trial ln(L+) it was produced at.

Return Value real(kind=wp)

private pure function cavity_t_ice(ice) result(T_ice)

Ice temperature actually used: exactly zero when the mode ignores it, so an unset or stale T_ice cannot leak into an insulating run through L_eff.

Arguments

Type IntentOptional Attributes Name
type(ocean_cavity_ice_t), intent(in) :: ice

Ice-conduction bundle.

Return Value real(kind=wp)


Subroutines

public pure subroutine cavity_exchange_velocities(par, const, u_star, l_plus, T_w, S_w, T_b, S_b, gamma_t, gamma_s, ierr)

Exchange-velocity dispatch with the Coriolis parameter taken from the bundle (par%f_cor). Thin wrapper over cavity_exchange_velocities_f, which holds the dispatch.

Arguments

Type IntentOptional Attributes Name
type(ocean_cavity_exchange_t), intent(in) :: par

Law selector + parameters.
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle.
real(kind=wp), intent(in) :: u_star

Friction velocity (m/s), strictly positive.
real(kind=wp), intent(in) :: l_plus

Trial viscous Obukhov scale, or CAVITY_L_PLUS_NEUTRAL.
real(kind=wp), intent(in) :: T_w

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

Far-field salinity (g/kg) — reserved-law argument.
real(kind=wp), intent(in) :: T_b

Trial interface temperature (degC) — reserved-law argument.
real(kind=wp), intent(in) :: S_b

Trial interface salinity (g/kg) — reserved-law argument.
real(kind=wp), intent(out) :: gamma_t

Heat exchange velocity (m/s). Zero on any non-OK status.
real(kind=wp), intent(out) :: gamma_s

Salt exchange velocity (m/s). Zero on any non-OK status.
integer, intent(out) :: ierr

CAVITY_MELT_* status.

public pure subroutine cavity_heat_fluxes(T_w, T_b, m_mass, gamma_t, ice, const, q_ocean, q_ice, q_latent)

The three heat fluxes of (E2), W/m^2:

Read more…

Arguments

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

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

Interface temperature (degC).
real(kind=wp), intent(in) :: m_mass

Melt mass flux (kg/m^2/s).
real(kind=wp), intent(in) :: gamma_t

Heat exchange velocity (m/s).
type(ocean_cavity_ice_t), intent(in) :: ice

Ice-conduction bundle.
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle.
real(kind=wp), intent(out) :: q_ocean

Turbulent heat flux ocean → interface (W/m^2).
real(kind=wp), intent(out) :: q_ice

Conductive + ice-warming flux interface → ice (W/m^2).
real(kind=wp), intent(out) :: q_latent

Latent heat consumed by the phase change (W/m^2).

public pure subroutine cavity_melt_columns(n, T_w, S_w, p_b, u_star, S_i, par, ice, eos, const, T_b, S_b, m_mass, q_ocean, ierr_col)

Data-parallel driver: solve n independent columns.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n

Number of columns.
real(kind=wp), intent(in) :: T_w(n)

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

Far-field salinity per column (g/kg).
real(kind=wp), intent(in) :: p_b(n)

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

Friction velocity per column (m/s).
real(kind=wp), intent(in) :: S_i(n)

Ice salinity per column (g/kg).
type(ocean_cavity_exchange_t), intent(in) :: par

Exchange-law bundle, shared by every column.
type(ocean_cavity_ice_t), intent(in) :: ice

Ice-conduction bundle, shared by every column.
type(eos_t), intent(in) :: eos

Shared EOS handle — the liquidus. Flat POD, by value.
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle, shared by every column.
real(kind=wp), intent(out) :: T_b(n)

Interface temperature per column (degC).
real(kind=wp), intent(out) :: S_b(n)

Interface salinity per column (g/kg).
real(kind=wp), intent(out) :: m_mass(n)

Melt mass flux per column (kg/m^2/s), > 0 melting.
real(kind=wp), intent(out) :: q_ocean(n)

Ocean -> interface heat flux per column (W/m^2).
integer, intent(out) :: ierr_col(n)

CAVITY_MELT_* status PER COLUMN. A column kernel must not take the run down for one bad column, so the failures are counted by the caller, not raised here.

public pure subroutine cavity_melt_columns_2d(nx, ny, cover, T_w, S_w, p_b, u_far, v_far, S_i, f_cor, cd, u_tide, ustar_min, par, ice, eos, const, u_star, T_b, S_b, m_mass, q_ocean, gamma_t, gamma_s, ierr_col)

Masked 2-D driver: form the friction velocity and solve the interface on every ICE-COVERED column of an (nx, ny) plane, leaving the rest untouched at exactly zero.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx

First dimension (ghosts included — the caller decides).
integer, intent(in) :: ny

Second dimension.
real(kind=wp), intent(in) :: cover(nx,ny)

Ice-cover mask; a column is solved iff cover > 0.5.
real(kind=wp), intent(in) :: T_w(nx,ny)

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

Far-field salinity (g/kg).
real(kind=wp), intent(in) :: p_b(nx,ny)

Interface pressure (Pa) — multilayer_state_t%p_top.
real(kind=wp), intent(in) :: u_far(nx,ny)

Far-field x velocity at the CELL CENTRE (m/s).
real(kind=wp), intent(in) :: v_far(nx,ny)

Far-field y velocity at the CELL CENTRE (m/s).
real(kind=wp), intent(in) :: S_i

Ice salinity (g/kg), >= 0; one scalar for the whole plane.
real(kind=wp), intent(in) :: f_cor(nx,ny)

Coriolis parameter (1/s) per column; read by hj99 only.
real(kind=wp), intent(in) :: cd

Top drag coefficient for the MELT friction velocity.
real(kind=wp), intent(in) :: u_tide

RMS tidal velocity (m/s); melt u* only, never the drag.
real(kind=wp), intent(in) :: ustar_min

Friction-velocity floor (m/s).
type(ocean_cavity_exchange_t), intent(in) :: par

Exchange-law bundle; its f_cor member is OVERRIDDEN per column by the f_cor array above.
type(ocean_cavity_ice_t), intent(in) :: ice

Ice-conduction bundle, shared by every column.
type(eos_t), intent(in) :: eos

Shared EOS handle — the liquidus. Flat POD, by value.
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle, shared by every column.
real(kind=wp), intent(out) :: u_star(nx,ny)

Friction velocity (m/s); 0 where uncovered or refused.
real(kind=wp), intent(out) :: T_b(nx,ny)

Interface temperature (degC); 0 where uncovered.
real(kind=wp), intent(out) :: S_b(nx,ny)

Interface salinity (g/kg); 0 where uncovered.
real(kind=wp), intent(out) :: m_mass(nx,ny)

Melt mass flux (kg/m^2/s), > 0 melting; 0 where uncovered.
real(kind=wp), intent(out) :: q_ocean(nx,ny)

Ocean -> interface heat flux (W/m^2); 0 where uncovered.
real(kind=wp), intent(out) :: gamma_t(nx,ny)

Thermal exchange velocity (m/s) of each column’s converged solve; EXACTLY zero where the column is not solved.
real(kind=wp), intent(out) :: gamma_s(nx,ny)

Haline exchange velocity (m/s), same convention.
integer, intent(out) :: ierr_col(nx,ny)

CAVITY_MELT_* status per column; CAVITY_MELT_OK where uncovered.

public pure subroutine cavity_melt_point(T_w, S_w, p_b, u_star, S_i, par, ice, eos, const, T_b, S_b, m_mass, q_ocean, ierr)

Scalar entry point returning only what a coupling seam consumes: the interface state, the canonical melt mass flux and the ocean -> interface heat flux. A thin wrapper over cavity_solve_melt that keeps the solution BUNDLE inside the callee, so a do concurrent over columns needs no derived-type local(...) clause at all — each iteration writes its own array elements.

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).
real(kind=wp), intent(in) :: p_b

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

Friction velocity (m/s).
real(kind=wp), intent(in) :: S_i

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

Exchange-law bundle.
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(kind=wp), intent(out) :: T_b

Interface temperature (degC).
real(kind=wp), intent(out) :: S_b

Interface salinity (g/kg).
real(kind=wp), intent(out) :: m_mass

Melt mass flux (kg/m^2/s), > 0 melting; EXACTLY zero on any non-OK status.
real(kind=wp), intent(out) :: q_ocean

Turbulent heat flux ocean -> interface (W/m^2).
integer, intent(out) :: ierr

CAVITY_MELT_* status.

public pure subroutine cavity_melt_point_gamma(T_w, S_w, p_b, u_star, S_i, par, ice, eos, const, T_b, S_b, m_mass, q_ocean, gamma_t, gamma_s, ierr)

cavity_melt_point plus the two EXCHANGE VELOCITIES the solve converged on.

Read more…

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).
real(kind=wp), intent(in) :: p_b

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

Friction velocity (m/s).
real(kind=wp), intent(in) :: S_i

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

Exchange-law bundle.
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(kind=wp), intent(out) :: T_b

Interface temperature (degC).
real(kind=wp), intent(out) :: S_b

Interface salinity (g/kg).
real(kind=wp), intent(out) :: m_mass

Melt mass flux (kg/m^2/s), > 0 melting.
real(kind=wp), intent(out) :: q_ocean

Turbulent heat flux ocean -> interface (W/m^2).
real(kind=wp), intent(out) :: gamma_t

Thermal exchange velocity (m/s) of the converged solve.
real(kind=wp), intent(out) :: gamma_s

Haline exchange velocity (m/s) of the converged solve.
integer, intent(out) :: ierr

CAVITY_MELT_* status.

public pure subroutine cavity_salt_fluxes(const, S_w, S_b, m_mass, gamma_s, S_i, f_turb, f_phase)

The two sides of (E3), in (g/kg)*kg/m^2/s:

Read more…

Arguments

Type IntentOptional Attributes Name
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle — rho_w.
real(kind=wp), intent(in) :: S_w

Far-field salinity (g/kg).
real(kind=wp), intent(in) :: S_b

Interface salinity (g/kg).
real(kind=wp), intent(in) :: m_mass

Melt mass flux (kg/m^2/s).
real(kind=wp), intent(in) :: gamma_s

Salt exchange velocity (m/s).
real(kind=wp), intent(in) :: S_i

Ice salinity (g/kg).
real(kind=wp), intent(out) :: f_turb

Turbulent salt flux toward the interface.
real(kind=wp), intent(out) :: f_phase

Phase-change salt flux.

public pure subroutine cavity_solve_melt(T_w, S_w, p_b, u_star, S_i, par, ice, eos, const, sol, ierr)

cavity_solve_melt_f with the Coriolis parameter taken from the bundle (par%f_cor) — the scalar entry point every host caller and the kernel suite use.

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.
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.

public 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.

Read more…

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) — g*rho_i*h_ice, the undiluted ice load.
real(kind=wp), intent(in) :: gamma_t

Heat exchange velocity (m/s), >= 0.
real(kind=wp), intent(in) :: gamma_s

Salt exchange velocity (m/s), > 0.
real(kind=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(kind=wp), intent(out) :: T_b

Interface temperature (degC).
real(kind=wp), intent(out) :: S_b

Interface salinity (g/kg).
real(kind=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).

public 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.

Read more…

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).
real(kind=wp), intent(in) :: p_b

Interface pressure (Pa).
real(kind=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(kind=wp), intent(out) :: T_b

Interface temperature (degC) = T_f(S_w, p_b).
real(kind=wp), intent(out) :: S_b

Interface salinity (g/kg) = S_w exactly.
real(kind=wp), intent(out) :: m_mass

Melt mass flux (kg/m^2/s), > 0 melting.
integer, intent(out) :: ierr

CAVITY_MELT_* status.

public pure subroutine cavity_ustar(u, v, cd, u_tide, ustar_min, u_star, ierr)

Interface friction velocity (m/s),

Read more…

Arguments

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

Far-field velocity component (m/s).
real(kind=wp), intent(in) :: v

Far-field velocity component (m/s).
real(kind=wp), intent(in) :: cd

Top drag coefficient (dimensionless), >= 0.
real(kind=wp), intent(in) :: u_tide

RMS tidal velocity (m/s), >= 0 by squaring.
real(kind=wp), intent(in) :: ustar_min

Friction-velocity floor (m/s), >= 0.
real(kind=wp), intent(out) :: u_star

Friction velocity (m/s). Exactly zero on any non-OK status.
integer, intent(out) :: ierr

CAVITY_MELT_* status.

private pure subroutine cavity_exchange_velocities_f(par, f_cor, const, u_star, l_plus, T_w, S_w, T_b, S_b, gamma_t, gamma_s, ierr)

Exchange-velocity dispatch — ONE argument list for every law, so a later do concurrent kernel dispatches with a single select case and no reshaping. l_plus is the viscous Obukhov scale of the CURRENT iterate: the single scalar that carries the stratification feedback for every implicit law. Pass CAVITY_L_PLUS_NEUTRAL for the neutral / unsuppressed evaluation.

Read more…

Arguments

Type IntentOptional Attributes Name
type(ocean_cavity_exchange_t), intent(in) :: par

Law selector + parameters (its f_cor member is NOT read).
real(kind=wp), intent(in) :: f_cor

Coriolis parameter (1/s) for this column; hj99 only.
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle.
real(kind=wp), intent(in) :: u_star

Friction velocity (m/s), strictly positive.
real(kind=wp), intent(in) :: l_plus

Trial viscous Obukhov scale, or CAVITY_L_PLUS_NEUTRAL.
real(kind=wp), intent(in) :: T_w

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

Far-field salinity (g/kg) — reserved-law argument.
real(kind=wp), intent(in) :: T_b

Trial interface temperature (degC) — reserved-law argument.
real(kind=wp), intent(in) :: S_b

Trial interface salinity (g/kg) — reserved-law argument.
real(kind=wp), intent(out) :: gamma_t

Heat exchange velocity (m/s). Zero on any non-OK status.
real(kind=wp), intent(out) :: gamma_s

Salt exchange velocity (m/s). Zero on any non-OK status.
integer, intent(out) :: ierr

CAVITY_MELT_* status.

private pure subroutine cavity_gamma_hj99(u_star, l_plus, f_cor, const, gamma_t, gamma_s, ierr)

Holland & Jenkins (1999) eqs. (14)-(18) p. 1792:

Read more…

Arguments

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

Friction velocity (m/s), > 0.
real(kind=wp), intent(in) :: l_plus

Trial viscous Obukhov scale, or CAVITY_L_PLUS_NEUTRAL.
real(kind=wp), intent(in) :: f_cor

Coriolis parameter (1/s), used as |f|. A SCALAR, not the exchange bundle: see cavity_exchange_velocities_f.
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle.
real(kind=wp), intent(out) :: gamma_t

Heat exchange velocity (m/s).
real(kind=wp), intent(out) :: gamma_s

Salt exchange velocity (m/s).
integer, intent(out) :: ierr

CAVITY_MELT_* status.

private pure subroutine cavity_gamma_yung25(u_star, l_plus, gamma_t, gamma_s, ierr)

Yung et al. (2025) “StratFeedback”, eqs. (7)-(8) p. 5832:

Read more…

Arguments

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

Friction velocity (m/s), > 0.
real(kind=wp), intent(in) :: l_plus

Trial viscous Obukhov scale, or CAVITY_L_PLUS_NEUTRAL.
real(kind=wp), intent(out) :: gamma_t

Heat exchange velocity (m/s).
real(kind=wp), intent(out) :: gamma_s

Salt exchange velocity (m/s).
integer, intent(out) :: ierr

CAVITY_MELT_* status.

private pure subroutine cavity_ice_terms(ice, melting, const, c_i_eff, kh, ierr)

(c_i_eff, kh) for the requested ice-conduction mode — the two numbers through which every mode enters the quadratic.

Read more…

Arguments

Type IntentOptional Attributes Name
type(ocean_cavity_ice_t), intent(in) :: ice

Ice-conduction bundle.
logical, intent(in) :: melting

Melt/freeze branch, decided from the thermal driving.
type(ocean_cavity_const_t), intent(in) :: const

Constants bundle.
real(kind=wp), intent(out) :: c_i_eff

Effective ice heat capacity (J/kg/K) in L_eff = L_f + c_i_eff*(T_b - T_ice).
real(kind=wp), intent(out) :: kh

Conductance k_ice/h_ice (W/m^2/K). Exactly zero for both shipped modes; carried so the derivation above and the code stay the same expression when the reserved DIFFUSIVE mode lands.
integer, intent(out) :: ierr

CAVITY_MELT_* status.

private pure subroutine cavity_melt_point_gamma_f(T_w, S_w, p_b, u_star, S_i, par, f_cor, ice, eos, const, T_b, S_b, m_mass, q_ocean, gamma_t, gamma_s, ierr)

cavity_melt_point_gamma with the Coriolis parameter passed as a per-column SCALAR — what cavity_melt_columns_2d calls, so no device code has to build a modified ocean_cavity_exchange_t (see cavity_exchange_velocities_f).

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).
real(kind=wp), intent(in) :: p_b

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

Friction velocity (m/s).
real(kind=wp), intent(in) :: S_i

Ice salinity (g/kg).
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.
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(kind=wp), intent(out) :: T_b

Interface temperature (degC).
real(kind=wp), intent(out) :: S_b

Interface salinity (g/kg).
real(kind=wp), intent(out) :: m_mass

Melt mass flux (kg/m^2/s), > 0 melting.
real(kind=wp), intent(out) :: q_ocean

Turbulent heat flux ocean -> interface (W/m^2).
real(kind=wp), intent(out) :: gamma_t

Thermal exchange velocity (m/s) of the converged solve.
real(kind=wp), intent(out) :: gamma_s

Haline exchange velocity (m/s) of the converged solve.
integer, intent(out) :: ierr

CAVITY_MELT_* status.

private pure subroutine cavity_safe_state(eos, S_w, p_b, T_b, S_b, m_mass)

The documented SAFE STATE returned on every non-OK status: zero melt, interface salinity equal to the far field, interface temperature on the liquidus there. m_mass is EXACTLY zero, so a failed column contributes nothing to a heat/salt budget rather than contributing a plausible wrong number.

Read more…

Arguments

Type IntentOptional Attributes Name
type(eos_t), intent(in) :: eos

Shared EOS handle — carries the liquidus coefficient set.
real(kind=wp), intent(in) :: S_w

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

Interface pressure (Pa).
real(kind=wp), intent(out) :: T_b

Interface temperature (degC) = T_f(S_w, p_b).
real(kind=wp), intent(out) :: S_b

Interface salinity (g/kg) = S_w.
real(kind=wp), intent(out) :: m_mass

Melt mass flux (kg/m^2/s) = exactly 0.

private pure subroutine cavity_solution_reset(sol)

Define every component of the solution bundle. Called FIRST by cavity_solve_melt, before any early return, because the type carries no default initialisers (see its docstring: they would bar it from a do concurrent local(...) clause on gfortran).

Arguments

Type IntentOptional Attributes Name
type(ocean_cavity_solution_t), intent(out) :: sol

Solution bundle, zeroed.

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.

Read more…

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.

private pure subroutine cavity_state_at_x(x, 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)

Evaluate the whole interface at a trial x = ln(L+): exchange velocities at that stratification, the closed-form three-equation solve, the buoyancy flux the answer implies and the L+ it re-diagnoses. The fixed point of x -> ln(L+_new) is the solution of the implicit system.

Arguments

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

Trial ln(L+). At or above CAVITY_LP_X_HI the trial L+ is the neutral sentinel.
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.
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(kind=wp), intent(in) :: u_star

Friction velocity (m/s).
real(kind=wp), intent(in) :: T_w

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

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

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

Ice salinity (g/kg).
real(kind=wp), intent(out) :: T_b

Interface temperature (degC).
real(kind=wp), intent(out) :: S_b

Interface salinity (g/kg).
real(kind=wp), intent(out) :: m_mass

Melt mass flux (kg/m^2/s).
real(kind=wp), intent(out) :: gamma_t

Heat exchange velocity (m/s) at this iterate.
real(kind=wp), intent(out) :: gamma_s

Salt exchange velocity (m/s) at this iterate.
real(kind=wp), intent(out) :: b_flux

Interfacial buoyancy flux (m^2/s^3) the answer implies.
real(kind=wp), intent(out) :: lp_new

L+ re-diagnosed from that buoyancy flux.
integer, intent(out) :: ierr

CAVITY_MELT_* status.