rdb_ice_enthalpy.F90 Source File

!! Sea-ice enthalpy library (SIS2 port, PR 2).
module rdb_ice_enthalpy
   !! Specific-enthalpy <-> temperature inversions for the sea-ice
   !! thermodynamic column model.
   !!
   !! * The thermodynamic prognostic of the ice model is SPECIFIC
   !!   ENTHALPY (J/kg), never temperature; these are the exact T<->E
   !!   inversions every later rung (PR 3 Winton column, PR 4
   !!   transport re-layering) consumes.
   !! * Source: ../SIS2/src/SIS2_ice_thm.F90 (Apache-2.0) — T_Freeze
   !!   (:1622), enth_from_TS (:1647), enthalpy_liquid_freeze (:1679),
   !!   enthalpy_liquid (:1691), Temp_from_En_S (:1830); constant
   !!   defaults from ice_thermo_init (:1534-1617).
   !! * THE SIMPLIFICATION (ratified in PLAN_SEA_ICE.md): Cp_brine ==
   !!   Cp_ice. SIS2's default Cp_brine = Cp_water makes Temp_from_En_S
   !!   transcendental (a T_fr*log(T_fr/T) term) and needs the
   !!   20-iteration Newton/false-position solve at SIS2_ice_thm.F90:
   !!   1876-1933; with Cp_brine = Cp_ice the log term vanishes and the
   !!   T<->E map is a closed-form quadratic, invertible EXACTLY. Only
   !!   the Cp_ice == Cp_brine branches are ported; the else branches
   !!   are deliberately NOT.
   !! * WHY THE ROUND-TRIP IS EXACT: for S > 0, T < t_fr, substituting
   !!   ice_enth_from_ts's En into ice_temp_from_en_s gives
   !!   BB = 0.5*(ICE_LAT_FUS*t_fr/T + ICE_CP_ICE*T), and the quadratic
   !!   ICE_CP_ICE*T**2 - 2*BB*T + ICE_LAT_FUS*t_fr = 0 has T as its
   !!   smaller root by construction — so
   !!   ice_temp_from_en_s(ice_enth_from_ts(T,S), S) == T to machine
   !!   precision, not merely to a solver tolerance. The fresh-water
   !!   and melted branches invert trivially.
   !! * Everything pure elemental + `!$acc routine seq`: host scalars,
   !!   whole-array elemental calls, or per-point inside a
   !!   do concurrent device kernel.
   use rdb_constants, only: wp
   implicit none
   private

   public :: ice_t_freeze, ice_enth_from_ts, ice_enthalpy_liquid_freeze, &
             ice_enthalpy_liquid, ice_temp_from_en_s

   ! SIS2 ice-thermodynamics constants (ice_thermo_init defaults,
   ! SIS2_ice_thm.F90:1534-1617), unscaled SI. Co-located with their
   ! kernel per the rdb_eos TFR_S_COEFF precedent.
   real(wp), parameter, public :: ICE_LAT_FUS = 3.34e5_wp
      !! Latent heat of fusion (J/kg) — SIS2 `LATENT_HEAT_FUSION`.
   real(wp), parameter, public :: ICE_CP_ICE = 2100.0_wp
      !! Heat capacity of fresh ice (J/kg/K) — SIS2 `CP_ICE`.
   real(wp), parameter, public :: ICE_CP_WATER = 4200.0_wp
      !! Heat capacity of seawater as carried by the ICE model
      !! (J/kg/K) — SIS2 `CP_SEAWATER`. Deliberately DISTINCT from the
      !! ocean's `SEAWATER_CP = 3992` (`rdb_ocean_surface_flux`): SIS2's
      !! ice model owns its own Cp_water = 4200 as the enthalpy
      !! zero-of-reference, and reusing 3992 here would break bit-parity
      !! with SIS2's enthalpy and mis-close the PR-3 frazil energy
      !! handshake. Do not "unify" them.
   real(wp), parameter, public :: ICE_CP_BRINE = ICE_CP_ICE
      !! Heat capacity of brine pockets (J/kg/K). THE simplification:
      !! set equal to `ICE_CP_ICE` (SIS2's documented computational-
      !! convenience option for `CP_BRINE`, not its default Cp_water) so
      !! the T<->E map closes in quadratic form with no Newton loop.
   real(wp), parameter, public :: ICE_DTF_DS = -0.054_wp
      !! Liquidus slope dT_f/dS (degC/PSU) — SIS2 `DTFREEZE_DS`. Same
      !! value as the ocean-side `TFR_S_COEFF` (rdb_eos) by
      !! construction — SIS2 parity on both sides of the seam.
      !!
      !! DELIBERATELY NOT switched by `&ocean_eos_nml tfreeze_set`.
      !! That knob selects the OCEAN-side liquidus that
      !! `eos_freezing_point` returns (the sea-surface freezing
      !! temperature the frazil / basal-flux seam works against); THIS
      !! constant is internal to the SIS2 enthalpy relation, where it
      !! fixes the brine-pocket melting temperature inside the ice and
      !! is baked into the closed-form quadratic T<->E map (and into
      !! SIS2 bit-parity). Retuning it is an ice-thermodynamics change,
      !! not an ocean-liquidus one. Consequence, documented rather than
      !! papered over: under `tfreeze_set = "isomip"` the ocean surface
      !! freezing point and the ice-internal brine liquidus disagree by
      !! ~0.03 degC. The ISOMIP+ set exists for ICE-SHELF-CAVITY work,
      !! where the sea-ice column model is normally off; running both at
      !! once is legal but means accepting that offset.
   real(wp), parameter, public :: ICE_ENTH_LIQ_0 = 0.0_wp
      !! Enthalpy of liquid fresh water at 0 degC (J/kg) — SIS2
      !! `ENTHALPY_LIQUID_0`. Zero-point of the enthalpy scale; kept in
      !! the algebra (even though it is 0) so a future non-zero
      !! reference is a one-line change.
   real(wp), parameter, public :: ICE_LIQ_LIM = 0.99_wp
      !! `LIQUID_LIMIT` (PR 3a) — max liquid fraction before the
      !! excess-enthalpy clamp (SIS2_ice_thm.F90:500-524) and the
      !! bottom-freeze `min_dEnth_freeze` floor (:1164-1188). Co-located
      !! here (not `rdb_ice_column`) so both `rdb_ice_column` and
      !! `rdb_ice_mass` can depend on this leaf module without a cycle.
   integer, parameter, public :: ICE_NK_MAX = 8
      !! Compile-time cap on per-column local-array size (kappa-shear
      !! `NZ_STACK_MAX` precedent). `nk_ice <= ICE_NK_MAX` is asserted
      !! at slot init (`ocean_sea_ice_init`, `rdb_ice_state`). Co-located
      !! on this leaf (not `rdb_ice_column`) — like `ICE_LIQ_LIM` — so
      !! `rdb_ice_mass` can size its fixed device-stack locals by it
      !! without a `rdb_ice_column` → `rdb_ice_mass` → `rdb_ice_column`
      !! cycle. `rdb_ice_column` re-exports it for existing importers.

contains

   pure elemental function ice_t_freeze(s) result(t_fr)
      !! Freezing temperature (degC) of the ice brine at bulk salinity
      !! `s` — SIS2 `T_Freeze` (SIS2_ice_thm.F90:1622), linear liquidus.
      !! RAW `s` (no max(0,s) clamp) — faithful to SIS2, which clamps
      !! only in enth_from_TS. Do not unify with ice_enth_from_ts's
      !! internal t_fr.
      !$acc routine seq
      real(wp), intent(in) :: s
         !! Ice bulk salinity (PSU).
      real(wp) :: t_fr

      t_fr = ICE_DTF_DS*s
   end function ice_t_freeze

   pure elemental function ice_enthalpy_liquid_freeze(s) result(enth)
      !! Enthalpy (J/kg) of liquid water at the freezing point for
      !! salinity `s` — SIS2 `enthalpy_liquid_freeze` (SIS2_ice_thm.F90:1679).
      !$acc routine seq
      real(wp), intent(in) :: s
         !! Ice bulk salinity (PSU).
      real(wp) :: enth

      enth = ICE_CP_WATER*(ICE_DTF_DS*s) + ICE_ENTH_LIQ_0
   end function ice_enthalpy_liquid_freeze

   pure elemental function ice_enthalpy_liquid(t, s) result(enth)
      !! Enthalpy (J/kg) of liquid water at temperature `t` — SIS2
      !! `enthalpy_liquid` (SIS2_ice_thm.F90:1691). `s` is unused in
      !! this linear form; the argument is kept for SIS2 call-site
      !! parity (`enthalpy_liquid(T, S, ITV)`).
      !$acc routine seq
      real(wp), intent(in) :: t
         !! Water temperature (degC).
      real(wp), intent(in) :: s
         !! Ice bulk salinity (PSU) — unused, kept for call-site parity.
      real(wp) :: enth

      enth = ICE_ENTH_LIQ_0 + ICE_CP_WATER*t
   end function ice_enthalpy_liquid

   pure elemental function ice_enth_from_ts(t, s) result(enth)
      !! Ice specific enthalpy (J/kg) from temperature + bulk salinity —
      !! SIS2 `enth_from_TS` (SIS2_ice_thm.F90:1647). NOTE: here (and
      !! only here) the freezing point uses max(0, s), per SIS2.
      !! The SIS2 `else` branch at :1667-1671 (Cp_brine /= Cp_ice, the
      !! T_fr*log(T_fr/T) form) is deliberately NOT ported — with
      !! ICE_CP_BRINE == ICE_CP_ICE it is unreachable, and dropping it
      !! keeps the map closed-form-invertible (see module docstring).
      !$acc routine seq
      real(wp), intent(in) :: t
         !! Ice temperature (degC).
      real(wp), intent(in) :: s
         !! Ice bulk salinity (PSU).
      real(wp) :: enth

      real(wp) :: t_fr

      t_fr = ICE_DTF_DS*max(0.0_wp, s)

      if (s == 0.0_wp .and. t <= 0.0_wp) then
         ! Fresh water at/below freezing is assumed all ice, due to the
         ! degeneracy in inverting temperature for enthalpy (SIS2 note).
         enth = (ICE_ENTH_LIQ_0 - ICE_LAT_FUS) + ICE_CP_ICE*t
      else if (t >= t_fr) then
         ! Already melted: just the sensible heat relative to 0 degC.
         enth = ICE_ENTH_LIQ_0 + ICE_CP_WATER*t
      else
         ! Cp_ice == Cp_brine closed form (SIS2 :1664-1666).
         enth = (ICE_ENTH_LIQ_0 - ICE_LAT_FUS*(1.0_wp - t_fr/t)) + &
                (ICE_CP_ICE*t + (ICE_CP_WATER - ICE_CP_ICE)*t_fr)
      end if
   end function ice_enth_from_ts

   pure elemental function ice_temp_from_en_s(en, s) result(t)
      !! Ice temperature (degC) from specific enthalpy + bulk salinity —
      !! SIS2 `Temp_from_En_S` (SIS2_ice_thm.F90:1830), Cp_ice ==
      !! Cp_brine path only: the quadratic
      !!   ICE_CP_ICE*T**2 - 2*BB*T + ICE_LAT_FUS*t_fr = 0
      !! solved for its smaller root. The SIS2 Newton/false-position
      !! refinement (:1876-1933) is only needed when Cp_brine /= Cp_ice
      !! and is deliberately NOT ported. Uses RAW `s` for t_fr (no max),
      !! faithful to SIS2.
      !$acc routine seq
      real(wp), intent(in) :: en
         !! Ice specific enthalpy (J/kg).
      real(wp), intent(in) :: s
         !! Ice bulk salinity (PSU).
      real(wp) :: t

      real(wp) :: t_fr, en_j, bb

      t_fr = ICE_DTF_DS*s
      en_j = en - ICE_ENTH_LIQ_0

      if (s <= 0.0_wp) then
         ! Step function for fresh water: liquid / mushy plateau / solid.
         if (en_j >= 0.0_wp) then
            t = en_j/ICE_CP_WATER
         else if (en_j >= -ICE_LAT_FUS) then
            t = 0.0_wp
         else
            t = (en_j + ICE_LAT_FUS)/ICE_CP_ICE
         end if
      else if (en_j >= t_fr*ICE_CP_WATER) then
         ! Completely melted layer.
         t = en_j/ICE_CP_WATER
      else
         ! Closed-form quadratic root (SIS2 :1868-1874).
         bb = 0.5_wp*((en_j - t_fr*(ICE_CP_WATER - ICE_CP_ICE)) + ICE_LAT_FUS)
         t = (bb - sqrt(bb*bb - t_fr*ICE_CP_ICE*ICE_LAT_FUS))/ICE_CP_ICE
      end if
   end function ice_temp_from_en_s

end module rdb_ice_enthalpy