CliMA · tapios · Mar 9, 2026 · Mar 7, 2026 · Mar 7, 2026 · Mar 7, 2026
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
@@ -31,7 +31,7 @@ steps:
 
   - label: "CPU tests"
     command:
-      - "THERMODYNAMICS_INCLUDE_DEPRECATED=false julia --project=test/ --color=yes test/runtests.jl"
+      - "julia --project=test/ --color=yes test/runtests.jl"
     agents:
       slurm_ntasks: 1
     timeout_in_minutes: 60

diff --git a/NEWS.md b/NEWS.md
@@ -1,10 +1,13 @@
 Thermodynamics.jl Release Notes
 ========================
 
-v1.xxx  (TODO: release once refactoring is done)
+v1.0.0
 --------
 
-- ![][[badge-✨feature/enhancement]] Complete rewrite of documentation
+- ![][[badge-💥breaking]] Removal of the deprecated object-oriented and state-based API (e.g., `ThermodynamicState`, `PhasePartition`, `PhaseEquil`). The package now relies exclusively on a stateless, completely functional API.
+- ![][[badge-🚀performance]] Enforced zero-allocations across core routines (`air_temperature`, `air_pressure`, `saturation_adjustment` with fixed iterations) backed by explicit testing.
+- ![][[badge-✨feature/enhancement]] Upgraded testing infrastructure: `Documenter.doctest` continuously checks all mathematical examples in the docstrings.
+- ![][[badge-✨feature/enhancement]] Rewrite of documentation
   - Restructured documentation for better organization and clarity
   - Updated all function documentation to be consistent and comprehensive
   - Improved code examples and usage patterns throughout documentation
@@ -17,10 +20,6 @@ v1.xxx  (TODO: release once refactoring is done)
   - Added comprehensive tests and physical consistency validation for these functions
   PR [259](https://github.com/CliMA/Thermodynamics.jl/pull/259)
   PR [263](https://github.com/CliMA/Thermodynamics.jl/pull/263)
-
-main
-----
-
 - Renamed `specific_enthalpy*` to `enthalpy*`.
 - Renamed `specific_entropy*` to `entropy*`.
 - Renamed `latent_heat_liq_ice` to `humidity_weighted_latent_heat`.

diff --git a/Project.toml b/Project.toml
@@ -1,11 +1,10 @@
 name = "Thermodynamics"
 uuid = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 authors = ["Climate Modeling Alliance"]
-version = "0.15.8"
+version = "1.0.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 
@@ -18,7 +17,6 @@ CreateParametersExt = "ClimaParams"
 [compat]
 ClimaParams = "1"
 ForwardDiff = "1.0.1"
-
 Random = "1"
 RootSolvers = "0.4, 1"
 julia = "1.6"
diff --git a/README.md b/README.md
@@ -55,7 +55,7 @@ params = TD.Parameters.ThermodynamicsParameters(Float64)
 
 # 2. Define your thermodynamic variables
 ρ     = 1.1        # Density [kg/m³]
-e_int = 200000.0   # Internal energy [J/kg]
+e_int = -36000.0   # Internal energy [J/kg, can be negative]
 q_tot = 0.015      # Total specific humidity [kg/kg]
 q_liq = 0.005      # Liquid specific humidity [kg/kg]
 q_ice = 0.001      # Ice specific humidity [kg/kg]

diff --git a/docs/src/API.md b/docs/src/API.md
@@ -156,76 +156,3 @@ TemperatureProfiles.DryAdiabaticProfile
 Thermodynamics.DataCollection
 ```
 
-## Deprecated Functions
-
-These functions are deprecated and will be removed in a future release.
-
-### Backward Compatibility Wrappers
-
-These wrappers exist for backward compatibility with older versions of the package.
-
-```@docs
-specific_enthalpy
-specific_enthalpy_dry
-specific_enthalpy_vapor
-specific_enthalpy_liquid
-specific_enthalpy_ice
-dry_pottemp
-total_specific_enthalpy
-q_vap_saturation_generic
-latent_heat_liq_ice
-```
-
-### Other Deprecated Functions
-
-```@docs
-air_temperature_given_hq
-air_temperature_given_pρq
-air_temperature_given_pθq
-air_temperature_given_ρθq
-air_temperature_given_ρθq_nonlinear
-saturated
-total_specific_humidity
-liquid_specific_humidity
-ice_specific_humidity
-mixing_ratios
-specific_volume
-q_vap_from_RH_liquid
-temperature_and_humidity_given_TᵥρRH
-liquid_ice_pottemp_sat
-PhasePartition_equil
-PhasePartition_equil_given_p
-```
-
-## Thermodynamic State Constructors (Deprecated)
-
-```@docs
-ThermodynamicState
-PhasePartition
-PhaseDry
-PhaseDry_ρe
-PhaseDry_pT
-PhaseDry_pθ
-PhaseDry_pe
-PhaseDry_ph
-PhaseDry_ρθ
-PhaseDry_ρT
-PhaseDry_ρp
-PhaseEquil
-PhaseEquil_ρeq
-PhaseEquil_ρTq
-PhaseEquil_pTq
-PhaseEquil_pθq
-PhaseEquil_peq
-PhaseEquil_phq
-PhaseEquil_ρθq
-PhaseEquil_ρpq
-PhaseNonEquil
-PhaseNonEquil_ρTq
-PhaseNonEquil_pTq
-PhaseNonEquil_ρθq
-PhaseNonEquil_pθq
-PhaseNonEquil_peq
-PhaseNonEquil_phq
-PhaseNonEquil_ρpq
-```
diff --git a/docs/src/Formulation.md b/docs/src/Formulation.md
@@ -47,7 +47,7 @@ Specific thermodynamic formulations often vary in how they approximate the relev
 2. **Consistency**: The formulation is thermodynamically consistent (e.g., it conserves energy and satisfies the Clausius-Clapeyron relation).
 3. **Simplicity and Efficiency**: The formulation leads to closed-form expressions that can be evaluated efficiently.
 
-The implementation follows the thermodynamic formulation of the CliMA Earth System Model ([Yatunin2026](@cite)). It relies on the **Rankine-Kirchhoff approximations**, which provide a consistent framework for moist thermodynamics ([Romps2021](@cite)).
+The implementation follows the thermodynamic formulation of the CliMA Earth System Model [Yatunin2026](@cite). It relies on the **Rankine-Kirchhoff approximations**, which provide a consistent framework for moist thermodynamics [Romps2021](@cite).
 
 Specific choices of thermodynamic constants have been made in [ClimaParams.jl](https://github.com/CliMA/ClimaParams.jl) to maximize accuracy given the Rankine-Kirchhoff approximations ([Ambaum2020](@cite), [Yatunin2026](@cite)).
 
@@ -412,7 +412,7 @@ The sum of the specific enthalpy of moist air and the specific gravitational pot
 \end{equation}
 ```
 
-Moist static energy arises naturally as the static energy component that is transported in moist air. "Static" here refers to the fact that the (small) kinetic energy contribution to the total energy is neglected. The global integral of moist static energy is approximately conserved in adiabatic processes, even in the presence of reversible phase transitions and latent heat release. It is also approximately materially conserved ([Romps2015](@cite)).
+Moist static energy arises naturally as the static energy component that is transported in moist air. "Static" here refers to the fact that the (small) kinetic energy contribution to the total energy is neglected. The global integral of moist static energy is approximately conserved in adiabatic processes, even in the presence of reversible phase transitions and latent heat release. It is also approximately materially conserved [Romps2015](@cite).
 
 ## 9. Saturation Vapor Pressure
 
@@ -427,7 +427,7 @@ The Clausius-Clapeyron relation describes how the saturation vapor pressure $p_v
 
 Here, $L$ is the specific latent heat of the phase transition, which is $L_v$ for the saturation vapor pressure over liquid, or $L_s$ for the saturation vapor pressure over ice.
 
-Substituting the linear relation \eqref{e:LHTemperature} between latent heat and temperature, and taking $p_\mathrm{tr}$ to be the vapor pressure at the triple point (by definition equal to the saturation vapor pressures both over liquid and ice), the Clausius-Clapeyron relation can be integrated to give a closed-form expression for the saturation vapor pressure. This closed-form expression is known as the **Rankine-Kirchhoff approximation** ([Romps2021](@cite)), named after its independent discoverers. It is the unique solution consistent with our thermodynamic assumptions (constant heat capacities and ideal gas law):
+Substituting the linear relation \eqref{e:LHTemperature} between latent heat and temperature, and taking $p_\mathrm{tr}$ to be the vapor pressure at the triple point (by definition equal to the saturation vapor pressures both over liquid and ice), the Clausius-Clapeyron relation can be integrated to give a closed-form expression for the saturation vapor pressure. This closed-form expression is known as the **Rankine-Kirchhoff approximation** [Romps2021](@cite), named after its independent discoverers. It is the unique solution consistent with our thermodynamic assumptions (constant heat capacities and ideal gas law):
 
 ```math
 \begin{equation}
@@ -441,7 +441,7 @@ Substituting the linear relation \eqref{e:LHTemperature} between latent heat and
 !!! tip "Implementation Note"
     The saturation vapor pressure is implemented in the [`saturation_vapor_pressure`](@ref) function. The closed-form expression enables efficient computation without numerical integration.
 
-With $L_0 = L_{v,0}$ or $L_0 = L_{s,0}$ and the corresponding heat capacity difference $\Delta c_p$, this gives saturation vapor pressures over liquid or ice that are accurate within 3% for temperatures between 200K and 330K ([Ambaum2020](@cite)). The accuracy of this approximation depends on the choice of thermodynamic constants; the values used in `Thermodynamics.jl` (specified in [ClimaParams.jl](https://github.com/CliMA/ClimaParams.jl)) are chosen to minimize errors in the Rankine-Kirchhoff approximation ([Yatunin2026](@cite)).
+With $L_0 = L_{v,0}$ or $L_0 = L_{s,0}$ and the corresponding heat capacity difference $\Delta c_p$, this gives saturation vapor pressures over liquid or ice that are accurate within 3% for temperatures between 200K and 330K [Ambaum2020](@cite). The accuracy of this approximation depends on the choice of thermodynamic constants; the values used in `Thermodynamics.jl` (specified in [ClimaParams.jl](https://github.com/CliMA/ClimaParams.jl)) are chosen to minimize errors in the Rankine-Kirchhoff approximation [Yatunin2026](@cite).
 
 !!! example "Typical Values"
     At $T = 300$ K:
@@ -459,7 +459,7 @@ With $L_0 = L_{v,0}$ or $L_0 = L_{s,0}$ and the corresponding heat capacity diff
 
     The ratio of liquid to ice saturation vapor pressure at 300 K is approximately 12.4, reflecting the higher energy required for sublimation compared to vaporization.
 
-To obtain the saturation vapor pressure over a mixture of liquid and ice (e.g., in mixed-phase clouds), using a weighted average of the relevant specific latent heats in the vapor pressure \eqref{e:SatVaporPressure} leads to a thermodynamically consistent formulation ([Pressel2015](@cite)). That is, if a fraction $\lambda_p$ of the condensate is liquid and the complement $1-\lambda_p$ is ice, calculating the saturation vapor pressure with a specific latent heat $\lambda_p L_v + (1-\lambda_p)L_s$ gives a thermodynamically consistent saturation vapor pressure over the mixture.
+To obtain the saturation vapor pressure over a mixture of liquid and ice (e.g., in mixed-phase clouds), using a weighted average of the relevant specific latent heats in the vapor pressure \eqref{e:SatVaporPressure} leads to a thermodynamically consistent formulation [Pressel2015](@cite). That is, if a fraction $\lambda_p$ of the condensate is liquid and the complement $1-\lambda_p$ is ice, calculating the saturation vapor pressure with a specific latent heat $\lambda_p L_v + (1-\lambda_p)L_s$ gives a thermodynamically consistent saturation vapor pressure over the mixture.
 
 !!! note "Physical Interpretation"
     The weighted latent heat approach ensures thermodynamic consistency when computing saturation vapor pressure over mixed-phase conditions. This is important for modeling mixed-phase clouds where both liquid and ice coexist.
@@ -704,7 +704,7 @@ The speed of sound in (moist) unstratified air is
 \end{equation}
 ```
 
-with the appropriate gas constants for moist air. In the presence of stratification, additional terms arise ([Durran1999](@cite)).
+with the appropriate gas constants for moist air. In the presence of stratification, additional terms arise [Durran1999](@cite).
 
 !!! note "Physical Interpretation"
     The speed of sound represents the rate at which pressure perturbations propagate through the fluid. The expression accounts for the effect of moisture on the gas constant and specific heat capacities of the air mixture.

diff --git a/docs/src/TestedProfiles.md b/docs/src/TestedProfiles.md
@@ -41,7 +41,7 @@ import .TestedProfiles
 FT = Float64
 param_set = TP.ThermodynamicsParameters(FT)
 
-profiles = TestedProfiles.PhaseDryProfiles(param_set, Array{FT});
+profiles = TestedProfiles.DryProfiles(param_set, Array{FT});
 (;T, ρ, z) = profiles
 p1 = Plots.scatter(ρ, z./10^3, xlabel="Density [kg/m^3]", ylabel="z [km]", title="Density");
 p2 = Plots.scatter(T, z./10^3, xlabel="Temperature [K]", ylabel="z [km]", title="Temperature");
@@ -73,7 +73,7 @@ import .TestedProfiles
 FT = Float64
 param_set = TP.ThermodynamicsParameters(FT)
 
-profiles = TestedProfiles.PhaseEquilProfiles(param_set, Array{FT});
+profiles = TestedProfiles.EquilMoistProfiles(param_set, Array{FT});
 (;T, ρ, q_tot, z) = profiles
 p1 = Plots.scatter(ρ, z./10^3, xlabel="Density [kg/m^3]", ylabel="z [km]", title="Density");
 p2 = Plots.scatter(T, z./10^3, xlabel="Temperature [K]", ylabel="z [km]", title="Temperature");

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -1,6 +1,6 @@
 # Thermodynamics.jl
 
-A Julia package for atmospheric thermodynamics, providing consistent thermodynamic functions for moist air including all phases of water (vapor, liquid, and ice). `Thermodynamics.jl` implements the thermodynamic formulation of the CliMA Earth System Model ([Yatunin2026](@cite)), providing a framework for moist thermodynamics based on the **Rankine-Kirchhoff approximations** ([Romps2021](@cite)).
+A Julia package for atmospheric thermodynamics, providing consistent thermodynamic functions for moist air including all phases of water (vapor, liquid, and ice). `Thermodynamics.jl` implements the thermodynamic formulation of the CliMA Earth System Model [Yatunin2026](@cite), providing a framework for moist thermodynamics based on the **Rankine-Kirchhoff approximations** [Romps2021](@cite).
 
 ## Table of Contents
 

diff --git a/perf/common.jl b/perf/common.jl
@@ -17,7 +17,7 @@ using .TestedProfiles
     functional_inputs(param_set, ::Type{FT}; n=256)
 
 Construct equilibrium-consistent scalar inputs for functional benchmarks, derived from
-`test/TestedProfiles.jl` (no deprecated PhasePartition/state types).
+`test/TestedProfiles.jl`.
 """
 function functional_inputs(param_set, ::Type{FT}; n = 256) where {FT}
     ps = TestedProfiles.PhaseEquilProfiles(param_set, Array{FT})

diff --git a/perf/microbenchmarks.jl b/perf/microbenchmarks.jl
@@ -1,6 +1,6 @@
 include("common_micro_bm.jl")
 
-function benchmark_thermo_states(::Type{FT}) where {FT}
+function benchmark_functional_api(::Type{FT}) where {FT}
     summary = OrderedCollections.OrderedDict()
     ArrayType = Array{FT}
     param_set = TP.ThermodynamicsParameters(FT)
@@ -25,5 +25,5 @@ function benchmark_thermo_states(::Type{FT}) where {FT}
     tabulate_summary(summary)
 end
 
-benchmark_thermo_states(Float64)
-benchmark_thermo_states(Float32)
+benchmark_functional_api(Float64)
+benchmark_functional_api(Float32)
diff --git a/src/Parameters.jl b/src/Parameters.jl
@@ -84,7 +84,8 @@ end
 "Internal energy of vaporization at reference temperature `T_0`."
 @inline e_int_v0(ps::ATP) = LH_v0(ps) - R_v(ps) * T_0(ps)
 
-"Internal energy of fusion at reference temperature `T_0`."
+"Internal energy of fusion at reference temperature `T_0`. Equal to the latent heat of fusion
+because for incompressible condensed phases the `p Δv` term is negligible (`p Δv ≪ L_f`)."
 @inline e_int_i0(ps::ATP) = LH_f0(ps)
 
 "Ratio of dry air gas constant to isobaric specific heat `R_d / cp_d`."

diff --git a/src/Thermodynamics.jl b/src/Thermodynamics.jl
@@ -60,8 +60,6 @@ const APS = TP.AbstractThermodynamicsParameters
 
 include("ThermoTypes.jl")
 
-include("depr_PhasePartitionTypes.jl")
-
 @inline solution_type() = RS.CompactSolution()
 include("DataCollection.jl")
 import .DataCollection
@@ -80,13 +78,6 @@ include("air_entropies.jl")
 include("air_dry_adiabatic.jl")
 include("TemperatureProfiles.jl")
 
-# Soon to be removed
-include("depr_air_temperatures.jl")
-include("depr_saturation_adjustment.jl")
-include("depr_air_states.jl")
-include("depr_state_methods.jl")
-include("depr_phase_partition_methods.jl")
-
 Base.broadcastable(dap::DryAdiabaticProcess) = tuple(dap)
 Base.broadcastable(phase::Phase) = tuple(phase)
 

diff --git a/src/air_entropies.jl b/src/air_entropies.jl
@@ -53,6 +53,8 @@ The specific entropy of dry air at its partial pressure.
  - `s_d`: specific entropy of dry air [J/(kg·K)]
 
 In the dry limit (`q_tot = q_liq = q_ice = 0`, the default), the dry-air partial pressure equals the total pressure.
+Note: `entropy_dry` diverges logarithmically as `q_tot → 1` (since `p_d → 0`). See also
+the analogous warning in [`entropy_vapor`](@ref).
 """
 @inline function entropy_dry(
     param_set::APS,
@@ -89,6 +91,9 @@ The specific entropy of water vapor at its partial pressure.
  - `s_v`: specific entropy of water vapor [J/(kg·K)]
 
 Note: the entropy of water vapor diverges logarithmically as `q_tot → 0` (since `p_v → 0`).
+The analogous divergence occurs in [`entropy_dry`](@ref) as `q_tot → 1` (since `p_d → 0`).
+A small numerical regularization `ϵ_numerics(FT)` is added to the pressure before taking the
+logarithm to avoid floating-point exceptions, but results should not be trusted in these limits.
 """
 @inline function entropy_vapor(
     param_set::APS,

diff --git a/src/air_humidities.jl b/src/air_humidities.jl
@@ -6,13 +6,12 @@ export partial_pressure_vapor
 export specific_humidity_to_mixing_ratio
 export q_vap_from_p_vap
 export q_vap_from_RH
-export q_vap_from_RH_liquid
 export relative_humidity
 
 """
     vapor_specific_humidity(q_tot=0, q_liq=0, q_ice=0)
 
-The vapor specific humidity (clamped to be non-negative).
+The vapor specific humidity.
 
 # Arguments
  - `q_tot`: total specific humidity [kg/kg]
@@ -24,6 +23,8 @@ The vapor specific humidity (clamped to be non-negative).
 
 If the specific humidities are not given, the result is zero.
 The formula is `q_vap = q_tot - q_liq - q_ice`.
+
+Note: callers must ensure `q_liq + q_ice ≤ q_tot`; no clamping is applied.
 """
 @inline function vapor_specific_humidity(q_tot = 0, q_liq = 0, q_ice = 0)
     return q_tot - q_liq - q_ice
@@ -185,6 +186,8 @@ Compute the vapor specific humidity from the relative humidity.
 
 # Returns
  - `q_vap`: vapor specific humidity [kg/kg]
+
+Note: this function assumes all water is in vapor form (no liquid or ice condensate).
 """
 @inline function q_vap_from_RH(param_set::APS, p, T, RH, phase::Phase)
     p_vap_sat = saturation_vapor_pressure(param_set, T, phase)
@@ -193,26 +196,6 @@ Compute the vapor specific humidity from the relative humidity.
     return p_vap / Rv_over_Rd / (p - (1 - 1 / Rv_over_Rd) * p_vap)
 end
 
-"""
-    q_vap_from_RH_liquid(param_set, p, T, RH)
-
-Compute the vapor specific humidity from the relative humidity over liquid.
-
-# Arguments
- - `param_set`: thermodynamics parameter set, see [`Thermodynamics`](@ref)
- - `p`: pressure [Pa]
- - `T`: temperature [K]
- - `RH`: relative humidity [dimensionless], 0 ≤ RH ≤ 1
-
-# Returns
- - `q_vap`: vapor specific humidity [kg/kg]
-
-This function is deprecated. Use `q_vap_from_RH` with `Liquid()` instead.
-"""
-@inline function q_vap_from_RH_liquid(param_set::APS, p, T, RH)
-    return q_vap_from_RH(param_set, p, T, RH, Liquid())
-end
-
 """
     relative_humidity(param_set, T, p, q_tot=0, q_liq=0, q_ice=0)
 

diff --git a/src/air_properties.jl b/src/air_properties.jl
@@ -86,6 +86,7 @@ The specific heat capacities are assumed to be constant (calorically perfect air
     cv_v = TP.cv_v(param_set)
     cv_l = TP.cv_l(param_set)
     cv_i = TP.cv_i(param_set)
+    # rearranged formula for cv_m to avoid computation of vapor specific humidity
     return cv_d +
            (cv_v - cv_d) * q_tot +
            (cv_l - cv_v) * q_liq +