eos_freezing_point Function

public pure elemental function eos_freezing_point(eos, S, p) result(T_f)

Seawater freezing point T_f (degC) at salinity S and pressure p — the ocean-side prerequisite for the sea-ice port (PLAN_SEA_ICE.md, PR 1). Same point-function style as eos_density_point (flat-POD eos_t by value, device- callable), plus elemental so callers can evaluate whole salinity arrays in one reference.

The form is LINEAR for every eos%variant:

T_f = λ1·S + λ2 + λ3·p

with the three coefficients carried ON THE HANDLE (eos%tfr_s / eos%tfr_0 / eos%tfr_p) and selected as a named set by &ocean_eos_nml tfreeze_set:

  • "seaice" (DEFAULT) — SIS2/MOM6, λ = (−0.054, 0, −7.53e-8); T_f(35 PSU, 0 Pa) = −1.89 °C.
  • "isomip" — ISOMIP+ / Asay-Davis et al. (2016) Table 4, λ = (−0.0573, 0.0832, −7.53e-8); T_f(34.5, 0) = −1.89365 °C.

Keeping the linear form under the Wright / Roquet density branches is MOM6 parity, not a shortcut: MOM6’s TFREEZE_FORM = "LINEAR" is its default under any density branch.

VARIANT / FORM DISPATCH SEAM. A NONLINEAR liquidus — MOM6’s TFREEZE_FORM = "MILLERO_78" (Millero 1978, UNESCO TP28) or a TEOS-10 t_freezing(SA, p) polynomial — is a different FUNCTIONAL FORM, not another coefficient triple, so it does NOT belong in tfreeze_set. It slots in HERE, as a leading branch

if (eos%tfreeze_form == TFREEZE_FORM_MILLERO78) then … else

mirroring the eos%variant dispatch in eos_density_point and leaving the linear expression below untouched. Not implemented: the Millero (1978) coefficients are UNVERIFIED here (the primary document could not be obtained — see the prototype’s ice_shelf_melt/CITATIONS.md §1 “Millero (1978)”), and this repository does not ship unverified constants.

BIT-IDENTITY, and why the parentheses are load-bearing. The pre-knob expression was TFR_S_COEFF*S + TFR_P_COEFF*p, i.e. (λ1·S) + (λ3·p) by Fortran’s left-to-right evaluation. The expression below keeps EXACTLY that pair together in its own parenthesised subexpression and adds the intercept LAST, which is the order that makes the default set a no-op: parentheses are binding in Fortran, so a reassociating compiler (-fast / -ffast-math without -Kieee) may not fold λ2 into a different sum. Writing it as λ2 + λ1·S + λ3·p would have put the intercept INSIDE the pair and changed the legacy grouping.

At every production call site the result is bitwise unchanged. All four callers — ice_frazil_accumulate, ice_frazil_uptake{,_multicat}_impl, ice_compute_basal_flux_impl — pass p = 0.0_wp, and there λ3·p is a signed zero, so the sum is λ1·S regardless of whether the toolchain contracts the two products into an FMA: fma(λ1, S, ±0) = round(λ1·S) is the same value as round(λ1·S) + (±0) for every nonzero product. Adding λ2 ≡ TFR_0_COEFF ≡ +0.0 is then the IEEE-754 x + 0.0 identity — exact for every finite x EXCEPT x = −0.0. MEASURED (gfortran 15.1, -O3 -march=native): bitwise identical at p = 0 over 400 001 salinities spanning [0, 40], with the single exception below.

SIGNED ZERO, decided and documented: T_f is a zero at all only when λ1·S and λ3·p are both zero, i.e. only at S = 0 AND p = 0, and there the legacy expression returned −0.0 while this one returns +0.0. That difference is ACCEPTED, because −0.0 == +0.0 is .true., no consumer divides by T_f or forms 1/T_f, no consumer branches on sign(T_f), and every call site uses T_f only inside the difference T − T_f, where the sign of a zero cannot survive.

OFF the production envelope (p /= 0, which nothing passes yet — wiring the cavity pressure in is a later PR) the answer may move by at most 1 ulp from the pre-knob expression, and only on a toolchain whose FMA contraction is sensitive to whether the coefficients are compile-time parameters or runtime handle members. gfortran 15.1 at -O3 -march=native is such a toolchain: it contracts eos%tfr_s*S + eos%tfr_p*p into a vfmadd but did NOT contract the old all-constant form (measured: 33 825 of 160 040 (S, p /= 0) samples differ, max gap exactly 1 ulp). That is a rounding-mode difference in a more accurate direction, not a change of formula. test_default_set_bit_identical pins both arms.

Arguments

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

Shared EOS handle, by value — carries the liquidus coefficient set (tfr_s/tfr_0/tfr_p) written at configure by eos_apply_tfreeze_set, plus the variant tag the future nonlinear-form branch will read.
real(kind=wp), intent(in) :: S

Salinity (PSU / g/kg).
real(kind=wp), intent(in) :: p

Pressure (Pa), hydrostatic surface-relative (0 at the surface — the frazil kernel’s use case).

Return Value real(kind=wp)


Called by

proc~~eos_freezing_point~~CalledByGraph proc~eos_freezing_point eos_freezing_point proc~cavity_safe_state cavity_safe_state proc~cavity_safe_state->proc~eos_freezing_point proc~cavity_solve_melt_f cavity_solve_melt_f proc~cavity_solve_melt_f->proc~eos_freezing_point proc~cavity_solve_melt_f->proc~cavity_safe_state proc~cavity_state_at_x cavity_state_at_x proc~cavity_solve_melt_f->proc~cavity_state_at_x proc~cavity_three_equation cavity_three_equation proc~cavity_three_equation->proc~eos_freezing_point proc~cavity_three_equation->proc~cavity_safe_state proc~ice_compute_basal_flux_impl ice_compute_basal_flux_impl proc~ice_compute_basal_flux_impl->proc~eos_freezing_point proc~ice_frazil_accumulate_impl ice_frazil_accumulate_impl proc~ice_frazil_accumulate_impl->proc~eos_freezing_point proc~ice_frazil_uptake_impl ice_frazil_uptake_impl proc~ice_frazil_uptake_impl->proc~eos_freezing_point proc~ice_frazil_uptake_multicat_impl ice_frazil_uptake_multicat_impl proc~ice_frazil_uptake_multicat_impl->proc~eos_freezing_point proc~tfreeze_ib_impl tfreeze_ib_impl proc~tfreeze_ib_impl->proc~eos_freezing_point proc~thermal_driving_impl thermal_driving_impl proc~thermal_driving_impl->proc~eos_freezing_point 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_state_at_x->proc~cavity_safe_state proc~cavity_state_at_x->proc~cavity_three_equation proc~cavity_two_equation cavity_two_equation proc~cavity_two_equation->proc~cavity_safe_state proc~fill_tfreeze_ib fill_tfreeze_ib proc~fill_tfreeze_ib->proc~tfreeze_ib_impl proc~fill_thermal_driving fill_thermal_driving proc~fill_thermal_driving->proc~thermal_driving_impl proc~ice_compute_basal_flux ice_compute_basal_flux proc~ice_compute_basal_flux->proc~ice_compute_basal_flux_impl proc~ice_frazil_accumulate ice_frazil_accumulate proc~ice_frazil_accumulate->proc~ice_frazil_accumulate_impl proc~ice_frazil_uptake ice_frazil_uptake proc~ice_frazil_uptake->proc~ice_frazil_uptake_impl proc~ice_frazil_uptake->proc~ice_frazil_uptake_multicat_impl 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~engine_step_ice engine_step_ice proc~engine_step_ice->proc~ice_compute_basal_flux proc~engine_step_ice->proc~ice_frazil_accumulate proc~engine_step_ice->proc~ice_frazil_uptake proc~cavity_melt_columns cavity_melt_columns proc~cavity_melt_columns->proc~cavity_melt_point proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_step_ice proc~ocean_cavity_flux_step ocean_cavity_flux_step proc~ocean_cavity_flux_step->proc~cavity_melt_columns_2d proc~rdb_ocean_step rdb_ocean_step proc~rdb_ocean_step->proc~engine_step_ice proc~driver_run driver_run proc~driver_run->proc~driver_run_ocean proc~engine_step_finalize engine_step_finalize proc~engine_step_finalize->proc~ocean_cavity_flux_step

Source Code

   pure elemental function eos_freezing_point(eos, S, p) result(T_f)
      !! Seawater freezing point T_f (degC) at salinity `S` and pressure
      !! `p` — the ocean-side prerequisite for the sea-ice port
      !! (PLAN_SEA_ICE.md, PR 1).  Same point-function style as
      !! `eos_density_point` (flat-POD `eos_t` by value, device-
      !! callable), plus `elemental` so callers can evaluate whole
      !! salinity arrays in one reference.
      !!
      !! The form is LINEAR for every `eos%variant`:
      !!
      !!   T_f = λ1·S + λ2 + λ3·p
      !!
      !! with the three coefficients carried ON THE HANDLE
      !! (`eos%tfr_s` / `eos%tfr_0` / `eos%tfr_p`) and selected as a
      !! named set by `&ocean_eos_nml tfreeze_set`:
      !!
      !!   * `"seaice"` (DEFAULT) — SIS2/MOM6, λ = (−0.054, 0, −7.53e-8);
      !!     T_f(35 PSU, 0 Pa) = −1.89 °C.
      !!   * `"isomip"` — ISOMIP+ / Asay-Davis et al. (2016) Table 4,
      !!     λ = (−0.0573, 0.0832, −7.53e-8); T_f(34.5, 0) = −1.89365 °C.
      !!
      !! Keeping the linear form under the Wright / Roquet density
      !! branches is MOM6 parity, not a shortcut: MOM6's
      !! `TFREEZE_FORM = "LINEAR"` is its default under any density
      !! branch.
      !!
      !! VARIANT / FORM DISPATCH SEAM.  A NONLINEAR liquidus — MOM6's
      !! `TFREEZE_FORM = "MILLERO_78"` (Millero 1978, UNESCO TP28) or a
      !! TEOS-10 `t_freezing(SA, p)` polynomial — is a different
      !! FUNCTIONAL FORM, not another coefficient triple, so it does NOT
      !! belong in `tfreeze_set`.  It slots in HERE, as a leading branch
      !!
      !!   if (eos%tfreeze_form == TFREEZE_FORM_MILLERO78) then ... else
      !!
      !! mirroring the `eos%variant` dispatch in `eos_density_point` and
      !! leaving the linear expression below untouched.  Not implemented:
      !! the Millero (1978) coefficients are UNVERIFIED here (the primary
      !! document could not be obtained — see the prototype's
      !! `ice_shelf_melt/CITATIONS.md` §1 "Millero (1978)"), and this
      !! repository does not ship unverified constants.
      !!
      !! BIT-IDENTITY, and why the parentheses are load-bearing.  The
      !! pre-knob expression was `TFR_S_COEFF*S + TFR_P_COEFF*p`, i.e.
      !! `(λ1·S) + (λ3·p)` by Fortran's left-to-right evaluation.  The
      !! expression below keeps EXACTLY that pair together in its own
      !! parenthesised subexpression and adds the intercept LAST, which is
      !! the order that makes the default set a no-op: parentheses are
      !! binding in Fortran, so a reassociating compiler (`-fast` /
      !! `-ffast-math` without `-Kieee`) may not fold λ2 into a different
      !! sum.  Writing it as `λ2 + λ1·S + λ3·p` would have put the
      !! intercept INSIDE the pair and changed the legacy grouping.
      !!
      !! **At every production call site the result is bitwise unchanged.**
      !! All four callers — `ice_frazil_accumulate`,
      !! `ice_frazil_uptake{,_multicat}_impl`, `ice_compute_basal_flux_impl`
      !! — pass `p = 0.0_wp`, and there `λ3·p` is a signed zero, so the
      !! sum is `λ1·S` regardless of whether the toolchain contracts the
      !! two products into an FMA: `fma(λ1, S, ±0) = round(λ1·S)` is the
      !! same value as `round(λ1·S) + (±0)` for every nonzero product.
      !! Adding `λ2 ≡ TFR_0_COEFF ≡ +0.0` is then the IEEE-754 `x + 0.0`
      !! identity — exact for every finite `x` EXCEPT `x = −0.0`.
      !! MEASURED (gfortran 15.1, `-O3 -march=native`): bitwise identical
      !! at `p = 0` over 400 001 salinities spanning [0, 40], with the
      !! single exception below.
      !!
      !! SIGNED ZERO, decided and documented: `T_f` is a zero at all only
      !! when `λ1·S` and `λ3·p` are both zero, i.e. only at `S = 0` AND
      !! `p = 0`, and there the legacy expression returned `−0.0` while
      !! this one returns `+0.0`.  That difference is ACCEPTED, because
      !! `−0.0 == +0.0` is `.true.`, no consumer divides by `T_f` or
      !! forms `1/T_f`, no consumer branches on `sign(T_f)`, and every
      !! call site uses `T_f` only inside the difference `T − T_f`, where
      !! the sign of a zero cannot survive.
      !!
      !! OFF the production envelope (`p /= 0`, which nothing passes yet —
      !! wiring the cavity pressure in is a later PR) the answer may move
      !! by **at most 1 ulp** from the pre-knob expression, and only on a
      !! toolchain whose FMA contraction is sensitive to whether the
      !! coefficients are compile-time `parameter`s or runtime handle
      !! members.  gfortran 15.1 at `-O3 -march=native` is such a
      !! toolchain: it contracts `eos%tfr_s*S + eos%tfr_p*p` into a
      !! `vfmadd` but did NOT contract the old all-constant form (measured:
      !! 33 825 of 160 040 `(S, p /= 0)` samples differ, max gap exactly
      !! 1 ulp).  That is a rounding-mode difference in a more accurate
      !! direction, not a change of formula.
      !! `test_default_set_bit_identical` pins both arms.
      !$acc routine seq
      type(eos_t), intent(in) :: eos
         !! Shared EOS handle, by value — carries the liquidus
         !! coefficient set (`tfr_s`/`tfr_0`/`tfr_p`) written at configure
         !! by `eos_apply_tfreeze_set`, plus the variant tag the future
         !! nonlinear-form branch will read.
      real(wp), intent(in) :: S
         !! Salinity (PSU / g/kg).
      real(wp), intent(in) :: p
         !! Pressure (Pa), hydrostatic surface-relative (0 at the
         !! surface — the frazil kernel's use case).
      real(wp) :: T_f

      T_f = (eos%tfr_s*S + eos%tfr_p*p) + eos%tfr_0
   end function eos_freezing_point