Add pseudothickness to write_model_dataset for Omega

[cbegeman]

Add the pseudo-thickness to model datasets for Omega initial conditions.

Surface pressure (constant value from config file), temperature, salinity and layer thickness are used to compute the pressure at layer interfaces. One iteration is currently used, pending testing of how many iterations is needed for convergence to desired accuracy.

The pressure at interfaces and a reference density (from config file) are then used to compute the pseudo-thickness.

The pseudo-thickness is automatically added to any write_model_dataset call when Omega is the ocean model AND the optional config argument is provided. (Note: checks for presence of config argument are still needed)

Fixes #461

Depends on E3SM-Project/Omega#327 to rename LayerThickness to PseudoThickness
OR we can map variable name PseudoThickness (added to MPAS-O dataset before variable mapping) to LayerThickness

Checklist

• [N/A] User's Guide has been updated
• Developer's Guide has been updated
• API documentation in the Developer's Guide (api.md) has any new or modified class, method and/or functions listed
• Documentation has been built locally and changes look as expected
• Testing comment in the PR documents testing used to verify the changes

[cbegeman]

Testing

Polaris omega_pr suite

• Baseline workdir: /pscratch/sd/c/cbegeman/polaris-scratch/main-omega-develop-20260425
• Baseline build: /pscratch/sd/c/cbegeman/polaris-scratch/ztilde-omega-develop-20260425/build
• PR build: /pscratch/sd/c/cbegeman/polaris-scratch/ztilde-omega-develop-20260425/build
• PR workdir: /pscratch/sd/c/cbegeman/polaris-scratch/ztilde-omega-develop-20260425
• Machine: pm-cpu
• Compiler: gnu
• Build type: Release
• Log: not found
• Result:
  • Diffs (3 of 10):
    • ocean/planar/manufactured_solution/convergence_both/default
    • ocean/planar/manufactured_solution/convergence_both/del2
    • ocean/planar/manufactured_solution/convergence_both/del4

All 3 tests that fail have machine-precision level diffs in the forward/200km step

I think this is acceptable because we are now using the constant EOS and may have a machine-precision-level diff in density.

[cbegeman]

Noting that this will need a rebase onto #489

[cbegeman]

@xylar This is basically done so it would be good to get your review before you leave. I've tested the pr suite so all that's left on my list is rebase and clean up the docs.

[xylar]

Great, I'll review this in an hour or so.

[xylar]

Other than this rho0 issue, I think this looks great. Once rho0 is handled as in #489, I'll do some quick tests and hopefully approve.

[cbegeman]

@xylar I was wrong that this was ready because I think we still need some PR to actually call Eos to compute SpecVol. I think the soonest one to be merged would be E3SM-Project/Omega#353.

I also realized we might want to assume SpecVol=1/RhoSw so that PseudoThickness = LayerThickness for initial conditions where T,S don't need to be defined.

[cbegeman]

@hyungyukang I just rebased it and included a dependence on constant EOS for some tests. I'll retest now but this should be ready to build on for developing the baroclinic channel case. I would recommend first following this line

polaris/polaris/tasks/ocean/barotropic_gyre/init.py

Line 145 in 0ea9278

self.write_model_dataset(ds, 'init.nc', config)

for the baroclinc_channel case, adding config to the write_model_dataset call which will use the eos options in the config file to infer pseudothickness from the state. I would first try just running the default case with this polaris-generated IC and the eos options in the config set consistently with the Omega namelist and see if you can reproduce your results with the temporary code. And of course you should feel free to review the code changes in this PR and let me know if you spot anything that you think should be changed.

[hyungyukang]

@cbegeman , Thanks for the guide.

I've tested this PR with the time stepping PR (E3SM-Project/Omega#392) using the overflow test (ocean/planar/overflow/default).

Here're a list of changes I made to support Omega:

1. Added config to this line in tasks/ocean/overflow/init.py:

133         self.write_model_dataset(ds, 'init.nc', config)

2. Changed the tasks/ocean/overflow/forward.yaml . I added Omega part after MPAS-O part:

Omega part

Omega:
  TimeIntegration:
    TimeStepper: RungeKutta4
    TimeStep: 0000_00:00:15
    StopTime: none
    RunDuration: 0000_06:00:00
  Advection:
    FluxTracerType: Upwind
  Tendencies:
    ThicknessFluxTendencyEnable: true
    PVTendencyEnable: true
    KETendencyEnable: true
    SSHTendencyEnable: false
    VelDiffTendencyEnable: true
    ViscDel2: 1.0e3
    VelHyperDiffTendencyEnable: false
    TracerHorzAdvTendencyEnable: true
    TracerDiffTendencyEnable: false
    TracerHyperDiffTendencyEnable: false
    ThicknessVertAdvTendencyEnable: true
    VelocityVertAdvTendencyEnable: true
    TracerVertAdvTendencyEnable: true
    PressureGradTendencyEnable: true
  Tracers:
    Base: [Temperature, Salinity]
  IOStreams:
    InitialVertCoord:
      Filename: initial_state.nc
    InitialState:
      Filename: initial_state.nc
      Contents:
        - State
        - Tracers
    History:
      Filename: output.nc
      Freq: 1
      FreqUnits: hours
      IfExists: append
      # effectively never
      FileFreq: 9999
      FileFreqUnits: years
      Contents:
        - State
        - Tracers
        - KineticEnergyCell

3. Link OmegaMesh.nc  from initial_state.nc after create the case

And then I ran three simulations:

• Polaris_PseudoThick : Use the initial condition directly from Polaris (pseudo thickness created by Polaris with this PR)
• Omega_PseudoThick : Use geometric thickness from Polaris with my workaround code (see Enable Omega time stepping with diagnostic computations Omega#392 (comment)) in Omega
• Omega_GeometricThick : Use geometric thickness from Polaris without the workaround code (not the correct setup for Omega, but I ran it intentionally to demonstrate the behavior)

So, main comparison should be Polaris_PseudoThick and Omega_PseudoThick.

This figure shows time series of the domain integrated kinetic energy.

The domain integrated KE between Polaris_PseudoThick (red) and Omega_PseudoThick (blue) is not BFB as expected, but they differ by around 1.0E-6 which is reasonably small:

Total KE at Day 10

Polaris_PseudoThick
97996305150432.22

Omega_PseudoThick
97996301948114.11

Based on those simulations, I was able to confirm that this Polaris PR works as intended. It could serve as a long-term solution for generating Omega initial conditions from Polaris.

As shown in the Omega_GeometricThick experiment, Omega suffers from numerical noise and instability if PseudoThickness is not used in the initial condition.

[cbegeman]

@hyungyukang Thank you so much for this testing! Just to clarify, this case used a linear EOS, correct?

[hyungyukang]

@hyungyukang Thank you so much for this testing! Just to clarify, this case used a linear EOS, correct?

Yes, I used the linear EOS. I was also wondering if this PR supports teos10 too.

[cbegeman]

@hyungyukang The pseudothickness conversion feature should support TEOS10. I don't see any reason why the overflow test case wouldn't support it but I haven't tested it. If you have time to try it out, that would be great!

[cbegeman]

@xylar This is ready for your re-review. According to a comment above, you were interested in doing some testing.

[xylar]

I don't think I'm in a good position to test this on a truly layered case. I appreciate your testing, @hyungyukang!

I'm just going to double check that omega_pr runs fine for me on Frontier and approve.

[xylar]

My test is stuck in the queue on Frontier. It is late here and I don't think this test will add value beyond testing already done. I'm happy to approve. If the test shows any surprises tomorrow when I take a look, I'll be happy to investigate and follow up.

[xylar]

Thanks for seeing this through, @cbegeman. It's been in the works for a long time!

[xylar]

For what it's worth, my tests did run. As expected in manufactured_solution/default/forward/200km_300s, I'm seeing:

$ ncdump -h init.nc
netcdf init {
dimensions:
	NCells = 2500 ;
	MaxEdges = 6 ;
	NEdges = 7500 ;
	MaxCellsOnEdge = 2 ;
	NVertices = 5000 ;
	VertexDegree = 3 ;
	MaxEdges2 = 12 ;
	time = 1 ;
	NVertLayers = 1 ;
	NVertLayersP1 = 2 ;
variables:
...
	double LayerThickness(time, NCells, NVertLayers) ;
		LayerThickness:units = "m" ;
		LayerThickness:long_name = "pseudo-thickness" ;
		LayerThickness:note = "h_tilde = -dp / (RhoSw * g)" ;
...
	double PseudoThickness(time, NCells, NVertLayers) ;
		PseudoThickness:units = "m" ;
		PseudoThickness:long_name = "pseudo-thickness" ;
		PseudoThickness:note = "h_tilde = -dp / (RhoSw * g)" ;
...

Looks great!

[xylar]

Shoot, I spoke too soon. I'm seeing some analysis errors:

Exception raised while running the steps of the task
Traceback (most recent call last):
  File "/autofs/nccs-svm1_home1/xylar/e3sm_work/polaris/omega/add-pseudothickness-to-write-model-dataset/polaris/run/serial.py", line 395, in _log_and_run_task
    baselines_passed = _run_task(task, available_resources)
  File "/autofs/nccs-svm1_home1/xylar/e3sm_work/polaris/omega/add-pseudothickness-to-write-model-dataset/polaris/run/serial.py", line 559, in _run_task
    _run_step(
    ~~~~~~~~~^
        task,
        ^^^^^
    ...<3 lines>...
        step_log_filename,
        ^^^^^^^^^^^^^^^^^^
    )
    ^
  File "/autofs/nccs-svm1_home1/xylar/e3sm_work/polaris/omega/add-pseudothickness-to-write-model-dataset/polaris/run/serial.py", line 689, in _run_step
    step.run()
    ~~~~~~~~^^
  File "/autofs/nccs-svm1_home1/xylar/e3sm_work/polaris/omega/add-pseudothickness-to-write-model-dataset/polaris/ocean/convergence/analysis.py", line 175, in run
    convergence_failed = self.plot_convergence(
        variable_name=var['name'], title=var['title'], zidx=var['zidx']
    )
  File "/autofs/nccs-svm1_home1/xylar/e3sm_work/polaris/omega/add-pseudothickness-to-write-model-dataset/polaris/ocean/convergence/analysis.py", line 228, in plot_convergence
    error_res = self.compute_error(
        refinement_factor=refinement_factor,
    ...<2 lines>...
        error_type=error_type,
    )
  File "/autofs/nccs-svm1_home1/xylar/e3sm_work/polaris/omega/add-pseudothickness-to-write-model-dataset/polaris/ocean/convergence/analysis.py", line 374, in compute_error
    field_exact = self.exact_solution(
        refinement_factor,
    ...<2 lines>...
        zidx=zidx,
    )
  File "/autofs/nccs-svm1_home1/xylar/e3sm_work/polaris/omega/add-pseudothickness-to-write-model-dataset/polaris/tasks/ocean/manufactured_solution/analysis.py", line 70, in exact_solution
    init = self.open_model_dataset(
        f'init_r{refinement_factor:02g}.nc', self.config
    )
  File "/autofs/nccs-svm1_home1/xylar/e3sm_work/polaris/omega/add-pseudothickness-to-write-model-dataset/polaris/ocean/model/ocean_io_step.py", line 86, in open_model_dataset
    return self.component.open_model_dataset(
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
        filename, config=config, **kwargs
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "/autofs/nccs-svm1_home1/xylar/e3sm_work/polaris/omega/add-pseudothickness-to-write-model-dataset/polaris/tasks/ocean/__init__.py", line 297, in open_model_dataset
    ds['LayerThickness'] = geom_thickness_from_ds(ds, config=config)
                           ~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^
  File "/autofs/nccs-svm1_home1/xylar/e3sm_work/polaris/omega/add-pseudothickness-to-write-model-dataset/polaris/ocean/vertical/diagnostics.py", line 42, in geom_thickness_from_ds
    raise ValueError(
    ...<3 lines>...
    )
ValueError: Geometric layerThickness is not present in the initial condition and SpecVol is not present to compute it

Similar things are showing up in manufactured solution, rotation_2d and cosine_bell.

[xylar]

Looking into these. I'll post and issue if I can't come up with a quick fix.

[cbegeman]

Thanks. Did you confirm that the initial condition contains LayerThickness? The error suggests that it doesn't. I checked my test on pm and it did not hit that error.

[xylar]

The initial condition has LayerThickness, not layerThickness (as it seems it should). It also has PseudoThickness but not SpecVol. And I think this isn't quite right:

    elif 'SpecVol' in ds.keys() and 'PseudoThickness' in ds.keys():
        return RhoSw * ds['SpecVol'] * ds['layerThickness']

ds['layerThickness'] should be ds['PseudoThickness']. Obviously, I missed that in review.

[xylar]

I'm really sorry about this. I shouldn't have rushed the review. I'm too tired to figure out what's supposed to be happening here but I'll try to pick it up in the morning.

[cbegeman]

@xylar I think I know how to fix it. I should be able to get the bugfix PR open this evening.