Implement P-star Vertical Coordinate by brian-oneill · Pull Request #401 · E3SM-Project/Omega

brian-oneill · 2026-05-05T21:34:55Z

Changes include:

Read NVertLayers from InitVertCoord stream
Rename RefLayerThickness to RefPseudoThickness to avoid conflation with MPAS field refLayerThickness
Add Field metadata for RefPseudoThickness to InitVertCoord stream
Needed updates to IOStream class
Add computeTargetThickness method to computeMomVertAux
Add ALE contribution to computeVerticalVelocity

Checklist

Fixes #391

brian-oneill · 2026-05-05T21:48:31Z

Testing

Machines: pm-cpu, pm-gpu, frontier
Compilers: gnu, gnugpu, craycray, craycray-mphipcc
Build type: Release
Result: All tests passed

Polaris `omega_pr` suite

Baseline workdir: /lustre/orion/proj-shared/cli115/boneill/polaris/testing/baseline_omega_pr
Baseline build: /lustre/orion/proj-shared/cli115/boneill/polaris/polaris/e3sm_submodules/Omega/_build
PR build: /lustre/orion/proj-shared/cli115/boneill/myFork/prCPU
PR workdir: /lustre/orion/proj-shared/cli115/boneill/polaris/testing/my_branch_omega_pr
Machine: frontier
Partition: batch
Compiler: craygnu
Build type: Release
Log: /lustre/orion/proj-shared/cli115/boneill/polaris/testing/my_branch_omega_pr/polaris_omega_pr.o4531781
Result: All tests passed

xylar · 2026-05-06T05:57:11Z

@brian-oneill, thanks for making these changes! I'll give this a look.

But, in response to your question, yes, we want to read in VertCoordMovementWeights as part of the vertical initial condition stream. See the vertical coordinate design doc:
https://docs.e3sm.org/Omega/Omega/design/VertCoord.html#creation
The design is missing MinLayerCell and MaxLayerCell but we have those already in the stream.

xylar

@brian-oneill, this looks excellent to me. I just have two comments based on reviewing the code itself, both quite minor. I will go on to revising my Polaris branch with your (excellent) variable name change and make sure I see passing tests, the same as you.

xylar · 2026-05-06T07:28:34Z

+   if (FieldName == OmegaStr) {
+      OldFieldName = "restingThickness";
+   }


restingThickness is a geometric thickness, whereas RefPseudoThickness is obviously a pseudo-thickness. So I don't think we should allow any fallback if RefPseudoThickness is not present.

In other words, we can just take out this fallback. Since refPseudoThickness will obviously not be found, we should be good.

The immediate problem with removing the fallback is that the stream read will throw an error if one or the other isn't found. The mesh we use as ocean_test_mesh.ncin our ctests (https://web.lcrc.anl.gov/public/e3sm/inputdata/ocn/mpas-o/oQU240/ocean.QU.240km.151209.nc) doesn't have RefPseudoThickness. I could make the initialization just ignore the error if not found, but perhaps it's a good opportunity to just set up a dedicated Omega input mesh to point to in the quick start, as opposed to continuing reusing an old MPAS one.

I'm working on making new files with today's date as a date stamp. I think we need these files to be fixed, rather than relying on the fallbacks.

I just added new CTest meshes:

'ocean.QU.240km.omega_vars.260506.nc', 'PlanarPeriodic48x48.omega_vars.260506.nc', 'cosine_bell_icos480.omega_vars.260506.nc',

But those are only appropriate after the rename.

Let's leave this in for now and I'll take it out in #327

xylar · 2026-05-06T08:01:48Z

    int TracerTimeLevel,            ///< [in] Time level
-    TimeInstant Time                ///< [in] Time
+    TimeInstant Time,               ///< [in] Time
+    TimeInterval ProjCoeff          ///< [in] Time interval


The name ProjCoeff isn't particularly clear. Maybe documentation can take care of this but it would be helpful to me if you add a comment here and possibly also elsewhere where ProjCoeff has been added as a method argument indicating that it is a fractional time step dependent on the time stepper and used to compute movement of layer interfaces.

@sbrus89 and I have found ProjCoeff to be a confusing name. Can we call it something related to time step or Dt?

Yeah, wasn't quite sure what to name this one. Maybe AleDt, but if the interval ends up being useful for some other projected state calculation, I'd lean against using an ALE-specific name. Is ProjDt clear enough, provided I add a description like this to the argument list?

Suggested change

TimeInterval ProjCoeff ///< [in] Time interval

TimeInterval ProjDt ///< [in] Time interval for projection over the current time stepper stage

Yep, ProjDt together with your comment would do the tick for me.

xylar · 2026-05-06T10:22:43Z

Excellent! I'm also seeing that the omega_pr suite passes with these changes and E3SM-Project/polaris#573:

Polaris `omega_pr` suite

Baseline workdir: /lustre/orion/cli115/scratch/xylar/polaris_1.0/frontier/test_20260504/omega-pr-baseline
Baseline build: /lustre/orion/cli115/scratch/xylar/polaris_1.0/frontier/test_20260504/omega-pr-baseline/build
PR build: /lustre/orion/cli115/scratch/xylar/polaris_1.0/frontier/test_20260506/omega-pr-p-star/build
PR workdir: /lustre/orion/cli115/scratch/xylar/polaris_1.0/frontier/test_20260506/omega-pr-p-star
Machine: frontier
Partition: batch
Compiler: craygnu
Build type: Release
Log: /lustre/orion/cli115/scratch/xylar/polaris_1.0/frontier/test_20260506/omega-pr-p-star/polaris_omega_pr.o4533646
Result: All tests passed

Once my two small concerns above get addressed, I'll be happy to approve and merge.

xylar · 2026-05-06T10:23:22Z

Oh, and after the docs get updated. Presumably the changes there won't be large.

xylar · 2026-05-06T11:47:02Z

@brian-oneill, when you update the docs, could you take out or replace the references to LayerThicknessPStar in the vertical coordinate docs? That variable doesn't seem to exist.

xylar · 2026-05-06T14:34:23Z

@hyungyukang, can you check if this works with the overflow test case as part of a review here?

xylar · 2026-05-06T14:34:58Z

@sbrus89, if it's feasible, we should test this with the baroclinic channel.

xylar · 2026-05-06T14:38:01Z

+   RKProj[0] = 1. / 2;
+   RKProj[1] = 1. / 2;
+   RKProj[2] = 1.;
+   RKProj[3] = 1.;


Could you comment (preferably in the code, but maybe here in the PR) on how these were derived?

@xylar , I don’t see it in the current develop branch, and RKProj does not appear to be used anywhere in the RK4 time stepper. I’m not sure where it came from either. My guess is that it may have been introduced during a rebase from an older version?

It is used in computing the vertical motion of the ALE coordinate. It is not in the time stepper, it is in the vertical advection (as Dt there).

@xylar , Sorry, I confused myself there. I understand now. Thanks!
@brian-oneill , As Xylar requested, adding a comment there would help make the code easier to understand.

These are the RKA coefficients shifted by 1 RK stage, so essentially the time interval over which the layer evolves from the old thickness to the projected thickness, which contributes to the transport velocity:
(projected thickness - old thickness) / ProjDt

MPAS-O uses the same coefficients in the RK4 time stepper for computing the ALE contribution to the vertical transport velocity.

Sounds good. Can you add that comment in the code?

cbegeman · 2026-05-06T18:35:38Z

                    const I4 K = KStart + KVec;
                    LocLayerThickTarget(ICell, K) =
-                        LocRefLayerThick(ICell, K) *
+                        LocRefPseudoThick(ICell, K) *


Does this imply that LocLayerThickTarget is actually in pseudo-thickness, not geometric thickness? If so, should we rename it LocPseudoThickTarget and should AleTerm above be using (LocThickTarget - PseudoThickness) or convert LocLayerThickTarget to geometric thickness?

This is not the purpose of this PR. I will rename this in #327

If we rename things other than *RefLayer* to *RefPseudo*, there will be essentially no end in what needs fixing.

I took the initiative to rename RefLayer to RefPseudo specifically because the reference thickness is added to the initialization stream in this PR and I just felt uneasy keeping the name as is, since refLayerThickness exists as a distinct 1D field in the MPAS meshes we've been using for tests.

As for the ALE term, all these quantities are pseudo-thicknesses, and the vertical velocities we use in Omega are pseudo-velocities $\tilde{w} \equiv d\tilde{z} / dt$, so no conversion is necessary.

Thanks for clarifying. Good to know that it's all pseudothickness here.

brian-oneill · 2026-05-07T19:23:44Z

But, in response to your question, yes, we want to read in VertCoordMovementWeights as part of the vertical initial condition stream.

@xylar (and anyone else): it looks like MPAS reads the movement weights from the mesh, but also the config_vert_coord_movement option can override the values read from the mesh. Do we want to preserve this behavior in Omega, allowing the weights to be switched via the MovementWeightType config option, or should I remove that option and we rely solely on the weights read from the stream?

xylar · 2026-05-07T19:31:47Z

@mark-petersen, you might be best equipped to answer @brian-oneill's question. My feeling would be that we typically want to change VertCoordMovementWeights in idealized tests and that it would be easier and more consistent with how we do testing if this is something Polaris handles, not something Omega overrides. Do you foresee any circumstances where an Omega override would be important?

cbegeman · 2026-05-07T21:11:09Z

@brian-oneill For what it's worth, I'm comfortable with you removing this option and only reading it in from the mesh.

xylar · 2026-05-08T11:44:02Z

@sbrus89 and @hyungyukang, any update on how your testing is going on this? We need to keep this moving.

cbegeman · 2026-05-08T14:29:08Z

@sbrus89 and @hyungyukang, any update on how your testing is going on this? We need to keep this moving.

@hyungyukang Let me know if you'd like me to test the overflow case.

hyungyukang · 2026-05-08T14:32:18Z

@xylar, @cbegeman , I'm working on it now and will post the results soon!

hyungyukang · 2026-05-08T15:18:37Z

Here's the overflow test result:

It looks the same as before. Is that expected, or should I update something in omega.yml?

sbrus89 · 2026-05-08T15:45:39Z

I'm just finalizing the baroclinic_channel branch in polaris. I'll report back soon.

xylar · 2026-05-08T16:59:05Z

It looks the same as before. Is that expected, or should I update something in omega.yml?

@hyungyukang, my understanding is that something should have changed. We are now calling:
https://github.com/brian-oneill/E3SM/blob/8e98997e4317f0d93787dc2528bd5c88b7666dbc/components/omega/src/ocn/AuxiliaryState.cpp#L102-L103
which we weren't before and we are now computing vertical advection accounting for coordinate movement:
https://github.com/brian-oneill/E3SM/blob/8e98997e4317f0d93787dc2528bd5c88b7666dbc/components/omega/src/ocn/VertAdv.cpp#L442-L458
which, again, we weren't before. Together, these should mean that the vertical coordinate and vertical advection behave differently from before.

I don't see anything in your overflow configuration (https://github.com/E3SM-Project/polaris/pull/572/changes#diff-3b61ea3e604ca93614415174a1dc7c8844506c2d3052e119f26a56ac301e37aa) to explain why you wouldn't be seeing any differences.

@brian-oneill, any guess why this isn't changing answers?

xylar · 2026-05-08T17:00:32Z

@hyungyukang, the differences might be small. Have you verified that you're not seeing changes, or just that things look similar by eye?

hyungyukang · 2026-05-08T17:34:37Z

@hyungyukang, the differences might be small. Have you verified that you're not seeing changes, or just that things look similar by eye?

@xylar , thanks a lot for the explanations. It was just an eyeball comaprison. I will check the field differences of the overflow test between p* and z* and post here later.

xylar · 2026-05-08T17:37:43Z

between p* and z* and post here later

How are you getting z*? I don't think that's what Omega produces before this change, but I could be mistaken.

Or do you mean z* from MPAS-Ocean?

brian-oneill · 2026-05-08T20:26:35Z

Some meshes like the QU240 mesh (including both the old one being used for ctests and the new ocean.QU.240km.omega_vars.260506.nc mesh with the newer omega variable names) have all movement weights equal 0. This leads to division by 0 in computeTargetThickness.

This doesn't seem intentional here, but I'm wondering what the desired behavior should be when this is encountered. Should a mesh with all zero movement weights be treated as a flag for an impermeable interfaces configuration and disable vertical transport? Or should initialization default to p* when invalid weights like this are read? Or just abort with an error?

brian-oneill added 13 commits May 4, 2026 13:03

Add getFilename function to IOStream

6a875e8

Read NVertLayers from InitVertCoord stream

dfa11af

Create calendar if needed in IOStream::init(void)

188d104

Add IOStream::init to tests where needed

a996d36

Rename RefLayerThickness to RefPseudoThickness

d09370c

Add RefPseudoThickness field

479f153

Exchange halo for RefPseudoThickness

87f9049

Add ALE contribution to vertical transport velocity

307c363

Add computeTargetThickness to computeMomVertAux

7c7b289

Propagate ProjDt through call tree to computeVerticalVelocity

4e4f39f

Add ALE projection time intervals to time steppers

0281bbd

Update ctests

ee31325

Correct TotalVerticalVelocity calculation

8e98997

brian-oneill mentioned this pull request May 6, 2026

VertCoord getting NVertLayers from wrong stream and not reading RefLayerThickness #391

Open

xylar reviewed May 6, 2026

View reviewed changes

xylar self-assigned this May 6, 2026

xylar requested a review from hyungyukang May 6, 2026 14:33

xylar requested a review from sbrus89 May 6, 2026 14:34

xylar reviewed May 6, 2026

View reviewed changes

cbegeman reviewed May 6, 2026

View reviewed changes

	TimeInterval ProjCoeff ///< [in] Time interval
	TimeInterval ProjDt ///< [in] Time interval for projection over the current time stepper stage

Conversation

brian-oneill commented May 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

brian-oneill commented May 5, 2026

Testing

Polaris omega_pr suite

Uh oh!

xylar commented May 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

xylar left a comment

Choose a reason for hiding this comment

Uh oh!

xylar May 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

brian-oneill May 7, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

xylar commented May 6, 2026

Polaris omega_pr suite

Uh oh!

xylar commented May 6, 2026

Uh oh!

xylar commented May 6, 2026

Uh oh!

xylar commented May 6, 2026

Uh oh!

xylar commented May 6, 2026

Uh oh!

xylar May 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

cbegeman May 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

brian-oneill commented May 7, 2026

Uh oh!

xylar commented May 7, 2026

Uh oh!

cbegeman commented May 7, 2026

Uh oh!

xylar commented May 8, 2026

brian-oneill commented May 5, 2026 •

edited

Loading

Polaris `omega_pr` suite

xylar commented May 6, 2026 •

edited

Loading

xylar May 6, 2026 •

edited

Loading

brian-oneill May 7, 2026 •

edited

Loading

Polaris `omega_pr` suite

xylar May 6, 2026 •

edited

Loading

cbegeman May 6, 2026 •

edited

Loading

xylar commented May 8, 2026 •

edited

Loading

xylar commented May 8, 2026 •

edited

Loading