[mpaiao]

What and why

This pull request is based on ED2 PR 445. The existing saturation vapour pressure function in ED-MORE used Bolton (1980) over liquid and a Magnus form over ice, blended by the liquid fraction. Both Bolton and Magnus are empirical, so d(e_sat)/dT was not the Clausius-Clapeyron slope of the latent heat of vaporisation the model uses (Kirchhoff, #49), so it is not fully consistent with the second law of thermodynamics.

This PR replaces the function with the formulation by Ambaum (2020). This uses a closed form of the Clausius-Clapeyron equation, built on the model's own latent_heat_vap, for each phase. The slope is e L / (Rv T^2) by construction, and both phases equal es_3ple (611.655 Pa, see Ambaum (2020), eq. 14) at the triple point. This ensures that the function is continuous in temperature and liquid fraction (the old curves differed by 0.06 Pa there). I did not propose this as an option because this formulation makes es consistent with the rest of ED-MORE, and the model is sufficiently new so maybe back-compatibility is still not an issue. The liquid-fraction blend, the "no fliq means liquid" rule and every caller are unchanged. With the model's own heat capacities and latent heat it is within 0.24 % of Murphy and Koop (2005) over liquid from -30 to 40°C (0.5 % at -40°C) and within 0.01 % at -10°C over ice.

The flux-tower tools mirror it: tower_conversions.py holds the Python copy, now with the ice branch, the liquid-fraction blend and slope, so it is equivalent to sat_phase in ed_therm_lib.f90 (liquid is still the default when no liquid fraction is given), and priors.py imports it instead of keeping a second copy. ambaum is added as the first of SATURATION_CURVES, but the bolton option is kept. Science pages, the output attribute text, unit tests (Clausius-Clapeyron against latent_heat_vap, equality at t_3ple, Murphy and Koop values, and Python values pinned to the Fortran) and the changelog are updated.

Refs #49. This pull request builds on #49 (latent heat of vaporisation as function of temperature), which is already merged.

Effect on results (tick exactly one; CI checks)

• Bit-for-bit: output is identical to the base commit's in every run listed under Verification. A new option that is off by default counts when the model is unchanged with it off.
• Round-off: differences at round-off level only, e.g. from reordered sums.
• Answer-changing: results move; give the largest changes below.

Run (site or example)	Variable	Period	Before	After
region test, all cases	mean_le_pa (mean; largest daily change 0.09 W/m2)	Jan-Feb 2021	16.98 W/m2	17.02 W/m2 (+0.22 %)
region test, all cases	mean_h_pa (mean; largest 0.056 W/m2)	Jan-Feb 2021	53.79 W/m2	53.76 W/m2 (-0.04 %)

User-facing changes (tick all that apply; each needs an "Upgrading" note in the changelog)

• None
• Config keys added, renamed or removed, or a value in the shipped config files changed
• Output variables or output files changed
• Input files (forcing, maps, initial conditions) or their format changed
• Building or running changed: command-line or CMake options, dependencies, the Python API
• State files written by the last release no longer load

No config key, output variable or file format changes. The forcing_qair global attribute of ED_ERA5land diagnostic files now names Ambaum (2020) instead of Bolton (1980). In the flux-tower site TOML, variables.VPD.curve still has to be declared; ambaum is now a valid value.

Verification

• Base commit (beta): db5ce8e. It was rebuilt from a clean worktree for the comparison.
• Runs compared with the base commit, and the result: region test daily output, 168 files. 55 of 163 daily variables change; the rest
  are identical. Largest means: latent heat and evapotranspiration +0.22 % (the sublimation and ground-evaporation terms move with them),
  snow water equivalent and depth −0.16 %, surface heat gain −0.20 %; energy and water residuals stay at round-off levels (at most 7e-11W/m2 and 2e-17 kg/m2/s). Everything else moves by 0.06 % or less.
• gfortran only. ifx is not available on my machine; CI covers it.
• The long tests (label long; expected for a science change): not run
• Other builds (production preset, nvfortran): not run
• Conservation warnings in those runs, or "none": not checked
• Full ctest suite (88 tests) and the lint checks: passed, gfortran dev build.
• Python: 138 tests pass, including the ones that load libedmore and run edmore.

Credit and AI use

• Collaborators:
• AI assistance (the tool and what it did), or "none": Claude Code (Sonnet 5.5) read the code, wrote the Fortran, tests, docs and Python changes, and ran the comparison. I specified the method, checked the formulation and es_3ple against Ambaum (2020), and reviewed the results.
• I answer for this pull request: I can explain every change in it, and the verification above was actually run

Not in this PR

• During testing, we spotted some python warnings on dayofyear being deprecated. This was fixed in a separate pull request (refactor(scripts): replace dayofyear with day_of_year (pandas) #52, merged).
• This only changes the formulation of water vapour pressure at saturation, but Claude pointed out that there are instances in which the code assumes liquid water only. These may require revisions and changes, but will be added later, on a separate pull request.

This message is based on draft generated with Claude Code.

[xiangtaoxu]

Marcos, when you got chance to work on this PR, can you rebase this to the newest beta?

[mpaiao]

Yes, I rebased it, but I haven't pushed it yet because I found some inconsistency between some of the python and the fortran code. The python code (tower_conversions.py and priors.py) implicitly assumed liquid-only phase, but I don't think we want this.

[mpaiao]

@xiangtaoxu this one is ready for review, please let me know if anything needs changes. Thanks!

[xiangtaoxu]

Hi, Marcos

Finally got time to go over the PR. Everything looks great.

I ran it for BCI and ET slightly decreased <0.1% because the Aubaum curve is below Bolton in warm temperature.

There could be some minor changes in doc, for example, the doc does not use letter equation number yet so for consistency, we should use (2) not (1b). But we will do an overhaul of science doc before stable release anyway, so I will not worry about it now.

[mpaiao]

I think the equation numbering was defined automatically by Claude, sorry that I overlooked it.