Adds the thermo coupling by alicebarthel · Pull Request #461 · E3SM-Project/Omega

alicebarthel · 2026-07-06T22:49:28Z

This PR adds the thermo coupling - from forcing terms to thickness and tracer tendencies.
The energy of mass fluxes and phase changes are hard-coded into the tendency conversions, like it was in MPAS-O. Hopefully the documentation (inline or otherwise) clarifies the meaning of each term.

This PR is 3/3 to make #418 more digestible.

Checklist

vanroekel · 2026-07-07T03:38:05Z

+- `SensibleHeatFlux`: sensible heat transfer
+- `LongWaveHeatFluxUp`: upward longwave radiation
+- `LongWaveHeatFluxDown`: downward longwave radiation
+- `SeaIceHeatFlux`: heat from sea-ice interaction


I don't understand what you mean by sea ice interaction, can you clarify?

I'll clarify. This includes the direct heat fluxes (conduction, etc.) and the enthalpy of liquid water from sea ice melt.

@vanroekel Is this comment ready to be resolved?

vanroekel · 2026-07-10T19:42:06Z

+         const Real PTopDb = PressureMid(ICell, KTop) * Pa2Db;
+         const Real SaTop  = SaltIndex >= 0
+                                 ? TracerCell(SaltIndex, ICell, KTop)
+                                 : 0.0_Real; // not sure we want zero here?


I think 0 is a fine limit. Does this form of code preclude logging when we go below this threshold?

@vanroekel note that this conditional is checking the validity of the tracer index (SaltIndex >= 0), not the value of salinity. For a check on salinity, it would be preferable to use Kokkos::clamp().
This was to make sure we are not using an invalid index, but we could do a clean fail rather than use 0 salinity...

I'm not sure that it's the job of this function to make sure that the tracer indices are initialized correctly. I'd be inclined to check for that elsewhere (and maybe not every timestep or every time forcing is updated) and remove the conditional. Thoughts?

Yea, I would agree. Are you thinking something like OceanInit, @cbegeman?

Yes, maybe Tracers::Init which is called during OceanInit?

vanroekel · 2026-07-10T19:59:55Z

overall this looks quite good @alicebarthel I only had pretty minor comments

alicebarthel · 2026-07-13T22:53:30Z

One key test to do is to verify that omega_pr is BFB with develop baseline since we do not currently have a polaris test using the surface tracer forcing. @katsmith133

Here is my current test report:

Polaris `omega_pr` suite

Baseline workdir: /global/homes/a/abarthel/pscratch/polaris-scratch/baseline-current-develop-cpu
Baseline build: /global/homes/a/abarthel/pscratch/omega_scratch/develop-cpu
PR build: /global/homes/a/abarthel/pscratch/omega_scratch/add-thermo-forcing-cpu-gnu
PR workdir: /global/homes/a/abarthel/pscratch/polaris-scratch/add-thermo-cpu_vsbaseline
Machine: pm-cpu
Compiler: gnu
Build type: Debug
Log: not found
Result:
- Diffs (1 of 12):
  - ocean/planar/barotropic_gyre/munk/free-slip

alicebarthel · 2026-07-14T00:50:46Z

the log is at
/pscratch/sd/a/abarthel/polaris-scratch/add-thermo-cpu_vsbaseline/case_outputs/ocean_planar_barotropic_gyre_munk_free-slip.log
if that can be helpful

cbegeman · 2026-07-15T17:09:05Z

@alicebarthel Is this ready to test with standalone Omega (I'm thinking column cases) with the flux fields provided in the forcing stream?

alicebarthel · 2026-07-15T17:23:56Z

Yes @cbegeman additional testing with single-column polaris cases (and analytical equilibrated state) or more realistic cases could be performed at this stage. It may be good to resolve omega_pr fail above before then, though? (I will not be able to debug this myself, sorry).
If it can be helpful for the additional testing, I generated forcing files with realistic values (taken from an existing simulation) and the necessary variables for various meshes, stored on pm:
in /global/homes/a/abarthel/pscratch/omega_scratch/
forcing_nCells_7153.nc [forcing.nc]
forcing_nCells_2304.nc [forcingPlanar.nc]
forcing_nCells_2562.nc [forcingSphere.nc]
You may prefer to start with idealized forcing but thought I'd let you know.

katsmith133 · 2026-07-16T18:30:51Z

Pulled this PR down onto PM this morning and verified it passes all CTests on PM-CPU and GPU. Also, the polaris omega_pr suite passed BFB on PM-CPU for me so perhaps @alicebarthel you had something wonky with your set up.

Polaris `omega_pr` suite

Baseline workdir: /pscratch/sd/k/katsmith/polaris_PR461//baseline_omega_pr
Baseline build: /pscratch/sd/k/katsmith/polaris-main/omega_build
PR build: /pscratch/sd/k/katsmith/polaris-PR461/omega_build
PR workdir: /pscratch/sd/k/katsmith/polaris_PR461/pr461_omega_pr
Machine: pm-cpu
Compiler: gnu
Build type: Release
Log: not found
Result: All tests passed

@cbegeman, I think it should be safe to test with single-column now if you want

katsmith133 · 2026-07-16T18:47:53Z

Pretty minor comment for the PR. Seems like overall its pretty good, but I'd also like to see the results of the single-column test. @cbegeman can I help in any way with this?

katsmith133 · 2026-07-16T19:03:37Z

Was just checking the compile logs and saw this warning for PM-GPU:

TracerForcingVars.cpp(198): warning #20011-D: calling a __host__ function("OMEGA::Eos::calcPtFromCt(const double &, const double &) const") from \
a __host__ __device__ function("OMEGA::TracerForcingVars::computeSurfInsituTemp(const  ::Kokkos::View<double ***,  ::Kokkos::LayoutRight,  ::Kokkos::CudaSpace >  &, const  ::OMEGA::VertCoord *, const  ::OMEGA::Eos *) c\
onst::[lambda(int) (instance 1)]::operator () const") is not allowed

working on how to resolve this issue...

katsmith133 · 2026-07-16T22:50:38Z

Was just checking the compile logs and saw this warning for PM-GPU:

TracerForcingVars.cpp(198): warning #20011-D: calling a __host__ function("OMEGA::Eos::calcPtFromCt(const double &, const double &) const") from \
a __host__ __device__ function("OMEGA::TracerForcingVars::computeSurfInsituTemp(const  ::Kokkos::View<double ***,  ::Kokkos::LayoutRight,  ::Kokkos::CudaSpace >  &, const  ::OMEGA::VertCoord *, const  ::OMEGA::Eos *) c\
onst::[lambda(int) (instance 1)]::operator () const") is not allowed

working on how to resolve this issue...

Ok, I figured out a solution to this that removes the compile warning, but I am now adding in Linear and Constant EOS support to this PR, which primarily is just adding non-TEOS10 options to the calcCtFreezing function. I believe I have that finished, but am working on testing. Note: Alice suggests (on L355-356 of Eos.cpp) that we use the UNESCO 1983 polynomial for any of the non-TEOS10 options, but I think this would be inconsistent? So for now I've implemented a linear option and a constant option. I would like feedback on that once I have it tested and pushed to this PR.

…iscussion" This reverts commit 70b0ca2.

…arning

This reverts commit d32ee4d.

cbegeman · 2026-08-11T20:24:42Z

+             SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) +
+             (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * CtTop +
+             (SnowFlux(ICell) + IceRunoffFlux(ICell)) *
+                 (Cp0Sw * CtFrz - LatIce);


What's the justification for including an enthalpy flux at CtFrz instead of MPAS-O's zero-temperature enthalpy? I remember considering whether we should change this for MPAS-O and the thought was that it was inconsistent with cross-component assumptions, e.g., E3SM-Project#7217

I assumed that CtFrz made more sense on the ocean side of things, since I was conceptually melting it in the ocean.

Since we are only using LatIce as the energy of phase change, we could simply shift the assumption to enthalpy at 0 C. That is equivalent to assuming that the frozen runoff arrives at exactly 0 C (all terms are freshwater so it should be ok), gets melted using energy from the ocean (in this implementation only the top ocean layer provides energy), then gets mixed into the ocean layer. For me, that works too.

@cbegeman In your prior discussions, was there concern re. the enthalpy reference point (i.e. 0C or triple point) being inconsistent across components?

Yes, exactly. The intention was that we preserve the assumption that frozen runoff at least arrives at 0 degC and thus the enthalpy flux is 0 as well. I think that our decision for snow flux and ice runoff could be considered separately since they originate from different components. However, for now, I think we should preserve MPAS-O behavior and add a note in the documentation of the potential for an implementation that provides these fluxes at the local freezing point.

Let's just say that the temperature is 0 deg C. Whether that means that the enthalpy is 0 depends on the enthalpy reference. I prefer if we keep the enthalpy formulation in there (which happens to be close to zero) so that we explicitly acknowledge we account for the enthalpy associated with the mass flux.

@alicebarthel Can you propose code or pseudo-code? I'm not sure what you mean here.

Are you saying that all enthalpy terms should be a function of the T, S, p assumed by the other components, i.e., gsw_enthalpy(0 g/kg salt, 0 degC, p_sfc)`?

We just chatted about this. I am just going to add a comment in above this calculation that notes that that term is technically there, but we've made the assumption that that term is 0.

Sorry for being unclear! There are several things that got brought up by this comment, which I am still pondering.

One thing to consider is that assuming CT = 0.00 deg C in omega is not exactly the same as assuming that PotT = 0.00 deg C in mapso, because CT != PotT for freshwater. It is ok to set the total energy change to Lf, because it is by far the biggest contributor but technically the 0.01 K reference means that we are missing ~25 J/kg.

The "enthalpy is zero" comment made me uncomfortable because it depends if we are talking about absolute enthalpy, potential enthalpy etc.

cbegeman · 2026-08-12T22:23:05Z

Testing

Verification of each term and benchmarked against MPAS-O reported here E3SM-Project/polaris#663 (comment)

katsmith133 · 2026-08-12T22:39:09Z

Retesting on Frontier after enthalpy calculation changes

CTest unit tests:

Machine: Frontier
Compiler: craygnu, craygnu-mphipcc
Build type: Release
Result: All tests passed

Polaris `omega_pr` suite

Baseline workdir: /lustre/orion/cli115/scratch/kmsmith/polaris_testing_pr461//main_omega_pr
Baseline build: /ccs/home/kmsmith/polaris-main/omega_build
PR build: /ccs/home/kmsmith/polaris-pr461/omega_build
PR workdir: /lustre/orion/cli115/scratch/kmsmith/polaris_testing_pr461/pr461_omega_pr
Machine: frontier
Partition: batch
Compiler: craygnu-mphipcc
Build type: Release
Log: not found
Result:
- Failures (1 of 22):
  - ocean/column/ekman

alicebarthel · 2026-08-13T15:04:56Z

+         // Heat tendencies are due to direct heat fluxes + enthalpy fluxes
+         // The enthalpy of liquid water is assumed to be:
+         // - local SST for liquid mass fluxes (rain, rivers)
+         // - local freezing point for solid --> liq mass fluxes (snow, frozen
+         // runoff)
+         // - solid mass fluxes are locally melted by the ocean (constant Lat
+         // heat of fusion)
+         const Real HeatFlux =
+             LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) +
+             LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) +
+             SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) +
+             (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * CtTop +
+             (SnowFlux(ICell) + IceRunoffFlux(ICell)) *
+                 (Cp0Sw * CtFrz - LatIce);


Suggested change

// Heat tendencies are due to direct heat fluxes + enthalpy fluxes

// The enthalpy of liquid water is assumed to be:

// - local SST for liquid mass fluxes (rain, rivers)

// - local freezing point for solid --> liq mass fluxes (snow, frozen

// runoff)

// - solid mass fluxes are locally melted by the ocean (constant Lat

// heat of fusion)

const Real HeatFlux =

LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) +

LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) +

SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) +

(RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw * CtTop +

(SnowFlux(ICell) + IceRunoffFlux(ICell)) *

(Cp0Sw * CtFrz - LatIce);

// Potential enthalpy tendencies are due to direct heat fluxes + enthalpy fluxes

// Each mass flux has an associated potential enthalpy flux.

// Levels of simplification can be done here. For now:

// - We approximate PotEnthalpyIce(Tinsitu, P=0) ~ -LatIce (constant);

// A better expression is pot_enthalpy_ice(T, P). At 0C, hI0 = −333360 J/kg

// a 0.020% / 66 J/kg difference with the LatIce value from pcd.

PotEnthalpyIce = - LatIce;

// - We assume dry snow and use PotEnthalpyIce for snow.

// - We assume liquid water comes in with the same specific enthalpy as the top ocean layer;

// Another possibility would be to add it at the same in situ temp as ocean, i.e. CT(Sa=0, T)

PotEnthalpyFw = Cp0Sw * CtTop;

const Real HeatFlux =

LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) +

LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) +

SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) +

(RainFlux(ICell) + RiverRunoffFlux(ICell)) * PotEnthalpyFw +

(SnowFlux(ICell) + IceRunoffFlux(ICell)) * PotEnthalpyIce ;

I think you commented on an outdated version of the code. I updated this section yesterday.

If you want to look at what is there now and comment on that, I can change it.

alicebarthel · 2026-08-13T15:12:02Z

+         // Heat tendencies are due to direct heat fluxes + enthalpy fluxes
+         // The enthalpy of liquid water is assumed to be:
+         // - local SST for liquid mass fluxes (rain, rivers)
+         // - zero degrees for solid --> liq mass fluxes (snow, frozen
+         // runoff)
+         // - solid mass fluxes are locally melted by the ocean (constant Lat
+         // heat of fusion)
+         // Technically there is an enthalpy flux associated with solid water
+         // associated with the energy needed to bring the solid water to the
+         // freezing point, captured by a Cp0Sw * CtFrz term. We make the
+         // assumption that the incoming solid fresh water from snow and ice
+         // is at 0C and that this term is 0. Thus, the enthalpy flux associated
+         // with solid water is only captured by the latent heat of fusion term,
+         // LatIce.
+         const Real HeatFlux =
+             LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) +
+             LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) +
+             SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) +
+             (RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw *
+                 Kokkos::max(0.0_Real, CtTop) -
+             (SnowFlux(ICell) + IceRunoffFlux(ICell)) * LatIce;


Suggested change

// Heat tendencies are due to direct heat fluxes + enthalpy fluxes

// The enthalpy of liquid water is assumed to be:

// - local SST for liquid mass fluxes (rain, rivers)

// - zero degrees for solid --> liq mass fluxes (snow, frozen

// runoff)

// - solid mass fluxes are locally melted by the ocean (constant Lat

// heat of fusion)

// Technically there is an enthalpy flux associated with solid water

// associated with the energy needed to bring the solid water to the

// freezing point, captured by a Cp0Sw * CtFrz term. We make the

// assumption that the incoming solid fresh water from snow and ice

// is at 0C and that this term is 0. Thus, the enthalpy flux associated

// with solid water is only captured by the latent heat of fusion term,

// LatIce.

const Real HeatFlux =

LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) +

LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) +

SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) +

(RainFlux(ICell) + RiverRunoffFlux(ICell)) * Cp0Sw *

Kokkos::max(0.0_Real, CtTop) -

(SnowFlux(ICell) + IceRunoffFlux(ICell)) * LatIce;

// Potential enthalpy tendencies are due to direct heat fluxes + enthalpy fluxes

// Each mass flux has an associated potential enthalpy flux.

// Levels of simplification can be done here. For now:

// - We approximate PotEnthalpyIce(Tinsitu, P=0) ~ -LatIce (constant);

// A better expression is pot_enthalpy_ice(T, P). At 0C, hI0 = −333360 J/kg

// a 0.020% / 66 J/kg difference with the LatIce value from pcd.

const Real PotEnthalpyIce = - LatIce;

// - We assume dry snow and use PotEnthalpyIce for snow.

// - We assume liquid water comes in at the specific enthalpy as the top ocean layer;

// if top layer CT is negative (lower than TEOS ref state), the enthalpy is capped by a lower bound.

// Another possibility would be to add it at the same in situ temp as ocean, i.e. CT(Sa=0, max(0, T))

//- We assume evaporation removes the same specific enthalpy as the top ocean layer;

// this is not capped by 0 to avoid an artificial cooling tendency to CT

const Real PotEnthalpyFwIn = Cp0Sw * Kokkos::max(0.0_Real, CtTop);

const Real PotEnthalpyFwout = Cp0Sw * CtTop;

const Real HeatFlux =

LatentHeatFlux(ICell) + SensibleHeatFlux(ICell) +

LongWaveHeatFluxUp(ICell) + LongWaveHeatFluxDown(ICell) +

SeaIceHeatFlux(ICell) + ShortWaveHeatFlux(ICell) +

(RainFlux(ICell) + RiverRunoffFlux(ICell)) * PotEnthalpyFwIn +

EvaporationFlux(ICell) * PotEnthalpyFwOut

(SnowFlux(ICell) + IceRunoffFlux(ICell)) * PotEnthalpyIce ;

@cbegeman note that this implementation is slightly different from Mpas-O since the max(0, CtTop) is applied to both runoff and rain (while it was only applied to runoff before).
I also have some qualms about applying the 0 limiter to CtTop instead of in situ Temperature, I think it doesn't make much physical sense but I documented it as is. What you meant to do physically would be to use CT(Sa=0, max(0, SST)).

alicebarthel requested a review from vanroekel July 6, 2026 22:49

vanroekel reviewed Jul 7, 2026

View reviewed changes

Comment thread components/omega/doc/devGuide/Forcing.md Outdated

vanroekel reviewed Jul 7, 2026

View reviewed changes

Comment thread components/omega/doc/userGuide/Forcing.md Outdated

vanroekel reviewed Jul 7, 2026

View reviewed changes

Comment thread components/omega/doc/userGuide/Forcing.md

vanroekel reviewed Jul 7, 2026

View reviewed changes

Comment thread components/omega/src/ocn/forcingVars/TracerForcingVars.cpp Outdated

vanroekel reviewed Jul 7, 2026

View reviewed changes

Comment thread components/omega/src/ocn/Forcing.cpp Outdated