(design-ocean-analysis-initial)= # Initial Omega Analysis Capabilities date: 2026/08/11 Contributors: Xylar Asay-Davis, Claude ## Summary The E3SM Ocean Team has committed to delivering a first set of analysis capabilities for Omega's initial coupled runs by **September 15, 2026**, as a stepping stone to a more complete zppy/MPAS-Analysis-style workflow later in the calendar year. This document designs that deliverable. The broader context is in {ref}`design-ocean-analysis`. The deliverable is a Polaris suite, `omega_analysis`, that is pointed at a completed Omega simulation through a user-supplied config file and produces: 1. **Map-view climatologies** (monthly, seasonal, and annual) of sea surface temperature, sea surface salinity, zonal and meridional velocity, mixed-layer depth, and vertically integrated ocean heat content. Fields with a vertical dimension are plotted at a configurable set of elevations. 2. **Global time series** of the quantities in Omega's `GlobalStats` output. 3. **A global time series of ocean heat content** integrated over configurable elevation ranges. 4. **A latitude-elevation plot of the global meridional overturning circulation (MOC)** from the MOC diagnostic that Omega computes in situ. Every plot is accompanied by a netCDF file containing the data plotted, and the expensive intermediate products (climatologies, reduced monthly heat content) are written to netCDF as well. The same simulation can be analyzed repeatedly over different date ranges. Results accumulate in a range-keyed staging tree, and re-analyzing a new range inherits the reduced monthly values earlier ranges already computed instead of recomputing them. The staging tree is published with a thumbnail for every plot and a generated gallery over them, as static HTML that a web server serves unchanged. The gallery follows MPAS-Analysis's familiar layout, and is designed around the fact that the portal hosting it throttles: what a page costs to load is a requirement here, not a detail of presentation. MPAS-Analysis is the scientific reference for what each of these diagnostics means, but the implementation is written from scratch with Polaris and Omega in mind rather than ported. The reasoning is in {ref}`design-ocean-analysis`. Four things are deliberately **out of scope** for this deliverable: - **Analyzing MPAS-Ocean output.** The analysis locates a simulation's output by reading the simulation's own Omega configuration, and MPAS-Ocean describes its output with namelists and streams files instead. Supporting it means writing a translator from those into the same form, which is separate work. This costs the deliverable nothing, since the design already develops against Omega output for the reasons given under `omega-monthly-means`. - **Integration with zppy.** For September 15, the analysis is run by hand by members of the Ocean Team, who make the results available to the coupled group within several days of each simulation period completing. MPAS-Seaice analysis continues to be delivered through zppy's MPAS-Analysis. - **Comparison with observations.** When observational comparison arrives, the intent is to remap the observations onto the MPAS mesh and compare there, not to remap the model onto a comparison grid. Everything in this document therefore stays on the native mesh, and nothing here should be read as a step toward interpolating model output to a lat-lon grid. - **A global MOC time series.** This was in the original proposal but has been dropped: MOC strength integrated globally is not a standard metric. The standard metric --- maximum Atlantic MOC near 26.5°N --- is deferred to a later delivery along with the rest of the regional analysis. The deliverable depends on Omega work that is outside this design: monthly-mean output of full model fields, including the geometric vertical coordinate, and a mixed-layer-depth diagnostic. These are described in the `omega-monthly-means` requirement below so that the dependencies are explicit, along with the fallback for mixed-layer depth if it cannot be delivered in time. Success is that a member of the Ocean Team can, with a config file and two Polaris commands, produce the full set of plots above for an Omega coupled run on a machine where that run's output lives. ### Conventions Two conventions are used throughout this document and throughout the code it describes. Both are Polaris-wide conventions rather than choices made here. **Names are MPAS-Ocean names.** MPAS-Ocean variable and dimension names are Polaris's internal standard. Analysis code refers to `temperature`, `salinity`, `zMid`, `zInterface`, `minLevelCell`, `maxLevelCell`, `bottomDepth`, and `areaCell`, and to the dimensions `nCells`, `nEdges`, `nVertLevels`, `nVertLevelsP1`, and `Time`, regardless of which model produced the data. Omega's names --- `Temperature`, `GeomZMid`, `GeomZInterface`, `NCells`, `NVertLayers`, and so on --- are translated to the Polaris standard automatically when a dataset is opened with `OceanIOStep.open_model_dataset`, using the mapping in `polaris/ocean/model/mpaso_to_omega.yaml`. Analysis steps therefore never branch on the model to get a field name, and config options that name fields use the MPAS-Ocean names. Where a field is new to Omega and has no MPAS-Ocean counterpart --- and is not expected to gain one --- it keeps its Omega name, and there is no entry to add: a mapping exists to reconcile two spellings of the same quantity, not to assign names. **Pseudo-thickness is not translated.** There is one deliberate exception to the rule above, and it matters enough to state up front. Omega is non-Boussinesq and prognoses **pseudo-height**, $\tilde{z}$, a normalized pressure with units of meters, defined in [Omega's governing equations](https://github.com/E3SM-Project/Omega/blob/develop/components/omega/doc/design/OmegaV1GoverningEqns.md) as $$ \tilde{z} = -\frac{p}{\rho_0 g}, \qquad\text{so that}\qquad d\tilde{z} = \frac{\rho}{\rho_0} \, dz $$ with $\rho_0$ a constant reference density used purely as a normalization --- its presence does not make the model Boussinesq. `PseudoThickness` is $\tilde{h} = \Delta \tilde{z}$, in meters, and $\tilde{h} \approx h$ only to the extent that $\rho \approx \rho_0$. Two consequences follow, and analysis code needs both: - **Reference density times pseudo-thickness equals full density times geometric thickness**, $$ \rho_0 \tilde{h} = \rho h = \frac{\Delta p}{g}, $$ which is the **mass per unit area** of the layer, exactly, by hydrostatic balance. It needs no equation of state. - The **geometric** thickness is the derived quantity: $h = \rho_0 \, \alpha \, \tilde{h}$, where $\alpha = 1/\rho$ is specific volume, so recovering it requires the equation of state. This is the relation Omega uses to build `GeomZInterface` and `GeomZMid`. MPAS-Ocean's `layerThickness` is the other way around: it is the geometric thickness $h$, and because MPAS-Ocean is Boussinesq, $\rho_0 h$ is *its* mass per unit area. The two variables therefore coincide in what they mean for a mass-weighted integral and differ in what they mean for a geometric one, so renaming one to the other would hide exactly the distinction that matters. `polaris/ocean/model/mpaso_to_omega.yaml` already declines to map them; this design records why. Analysis therefore never asks for "the thickness". It asks for one of two things by name: - **Geometric positions and thicknesses** come from `zInterface` and `zMid`, which both models write (Omega as `GeomZInterface` and `GeomZMid`) and which are translated as usual. A geometric layer thickness, where one is needed, is a difference of interface elevations, not a separate variable --- which is also the only way to get it offline, since $\alpha$ is not written. - **Mass per unit area** comes from the model's own mass-like thickness variable --- `PseudoThickness` for Omega, `layerThickness` for MPAS-Ocean --- read under its native name and multiplied by $\rho_0$. A single helper, `polaris.ocean.model.get_layer_mass`, returns $\rho_0 \tilde{h}$ for Omega and $\rho_0 h$ for MPAS-Ocean, so this is the one place in the analysis code that knows which model wrote the file. This is what makes the heat content integral below a mass integral rather than a volume integral scaled by a reference density, and it is the same convention used for evaluating conservation in Omega. **The vertical coordinate is elevation, positive up.** All vertical positions in config options, algorithms, and output are elevations $z$ in meters with $z = 0$ at the resting sea surface and $z$ increasing upward, so that positions within the ocean are negative. A map "at 100 m below the surface" is requested as `-100.0`, and the ocean heat content range conventionally called "0 to 700 m" is written `top:-700.0` --- `top` being the free surface of each column, which is what "0 m" means in that phrase, rather than the resting sea surface at $z = 0$. This matches `zMid` and `zInterface` as the models write them and avoids sign flips scattered through the code. Where the text uses the word "depth" it is describing a quantity that is positive down, such as `bottomDepth`, and says so. **Attributes are the CF ones.** Omega writes each variable attribute twice: once in its own capitalized spelling --- `Units`, `Name`, `StdName`, `Description`, `ValidMin`, `ValidMax` --- and once in the CF spelling --- `units`, `name`, `standard_name`, `long_name`, `valid_min`, `valid_max`. Analysis reads **only** the CF form and ignores the capitalized one entirely. This is not a stylistic preference. The CF attributes are the ones every tool in the stack already understands: `xarray` ingests them, uses `units` and `calendar` to decode times, and writes them back out. The capitalized duplicates are understood by nothing, so code that reads them works on files Omega wrote and fails on every other file --- including one that has merely been through `xarray`, which moves `units` into `.encoding` on decoding while leaving `Units` sitting in `.attrs`. A dependence on the capitalized form is therefore invisible until it is tested against a file we wrote ourselves, which is exactly when it is most expensive to find. `polaris.ocean.model.time.get_days_since_start` has such a dependence today, in `ds['Time'].Units`. It predates this design and is used across the ocean component, so fixing it is framework work rather than analysis work, but no new analysis code should follow it. ## Requirements ### Requirement: analysis-suite Date last modified: 2026/08/27 Contributors: Xylar Asay-Davis, Claude Polaris shall provide a suite that computes all of the analysis in this document for a single Omega simulation, run standalone or within E3SM. Analyzing MPAS-Ocean output is out of scope for this deliverable. The user shall supply a config file that provides: - the path to the simulation's own Omega configuration file, from which the mesh, the vertical coordinate and the output files shall be found without the user restating any of them; - the start and end dates of the climatology; - the fields for which climatology maps should be produced (if different from the defaults); - the elevations at which climatology maps should be plotted, for fields that have a vertical dimension (if different from the defaults); - the start and end dates of the time series; - the fields from Omega's global statistics for which time series should be plotted (if different from the defaults); - the elevation ranges over which ocean heat content should be integrated (if different from the defaults). Setting up and running the suite shall not require an Omega or MPAS-Ocean build, since no model is run. The suite shall be usable on the machine where the simulation output resides, without copying that output. Each analysis product shall be a separate task within the suite, so that a user can run a single product without running the rest, and the expensive climatology computation shall be shared between the tasks that need it. ### Requirement: data-products Date last modified: 2026/08/25 Contributors: Xylar Asay-Davis, Claude Every plot the suite produces shall be accompanied by a netCDF file containing exactly the data that were plotted, so that the values can be inspected, compared against other tools, or re-plotted without recomputation. Every plot and its data shall also be described in a machine-readable manifest naming the facets that identify it --- at minimum the field, the season, the vertical reduction, and the date range --- so that a reader or an index can find a product without knowing how the work was divided into steps. Intermediate products that are expensive to compute --- climatologies and reduced monthly ocean heat content in particular --- shall be written to netCDF, and a step that finds an intermediate product from a previous run shall be able to reuse it rather than recomputing it. Reuse shall be conditional on the product having been computed for the same simulation, by the same kernel, under the config options that govern it, and shall be reported rather than silent. ### Requirement: repeated-analysis Date last modified: 2026/08/29 Contributors: Xylar Asay-Davis, Claude Polaris shall support analyzing the same simulation repeatedly with different climatology and time-series date ranges. Analyzing a new range shall not overwrite or remove the results of a range analyzed earlier, and ranges analyzed at different times shall appear together in the published output described under `publication` below. Re-running the analysis with a changed range shall recompute the products that depend on that range. It shall not be necessary to delete the work directory or to pass a flag to force recomputation, and it shall not be possible to obtain a plot labeled with one range whose contents were computed for another. Re-running with a changed range shall reuse the intermediate results that do not depend on the range. Extending a time series from twenty years to forty shall cost twenty years of work, not forty. Re-running with an unchanged range shall recompute nothing by default, and it shall be possible to ask for the plots to be redrawn --- after a change to colormaps or other styling --- without recomputing the intermediate results behind them. ### Requirement: publication Date last modified: 2026/08/29 Contributors: Xylar Asay-Davis, Claude Results shall be published to a single location that a reader can browse without knowing anything about Polaris work directories, and that a web server can serve unchanged. That location shall carry a generated index over the results rather than requiring the reader to navigate directories. The index shall be static HTML, generated from the merged manifest, and shall require no server-side code, no build step, and no network access of its own, so that it works equally from a local filesystem and from a web portal. **Every plot shall have a thumbnail.** With a few hundred plots, thumbnails are what let a reader decide which full images to open, and they frequently answer the question without any full image being opened at all. A thumbnail shall be a separate, small, lossy file rather than the full image scaled down by the browser, which would cost the full image's bytes. **A page shall stay within what a throttled link can deliver.** The analysis is hosted on the LCRC public web portal, which throttles in a way that stalls a page asking for too much at once, and this is the constraint that the presentation has to be designed around rather than discovering later. Three things shall bound what a page costs: thumbnails shall be small, no page shall carry the whole result set, and images the reader has not scrolled to shall not be fetched. The parameters controlling this shall be config options, since the binding constraint is a property of the host rather than of the analysis. **The published output shall not need to be regenerated to be extended.** The index shall be derived from the manifest alone, so that facets added later --- region, observational reference, a second simulation to compare against --- and richer presentation --- per-product pages carrying provenance and the code that reproduces a figure, filtering, search --- can be added without changing any step, any manifest fragment already written, or any published path. The broader intent this serves is in {ref}`design-ocean-analysis`. What this deliverable does *not* owe is a considered visual design. A gallery following MPAS-Analysis's familiar layout, and no more, is what Phase 1 ships. ### Requirement: omega-monthly-means Date last modified: 2026/08/25 Contributors: Xylar Asay-Davis, Claude *This requirement describes work in Omega, not in Polaris. It is stated here because everything else in this document depends on it.* Omega shall be able to write monthly means of a configurable list of model fields. The monthly-mean output shall: - cover, at minimum, the fields needed by this analysis: conservative temperature, absolute salinity, pseudo-thickness, sea surface height, reconstructed zonal and meridional velocity at cell centers, and mixed-layer depth. Reconstructed velocities are **required** rather than preferred: the machinery to write them exists in Omega, two cell-centered fields are about two thirds the size of one edge field, and `normalVelocity` is not itself plotted by anything here, so there is no reason for a simulation to write the larger field and no reason for Polaris to carry a reconstruction it would never use; - include the **geometric vertical coordinate**, `GeomZMid` and `GeomZInterface`, so that Polaris does not have to reconstruct it (see the vertical-geometry algorithm design for why it cannot); - carry CF-compliant time metadata: a `time` coordinate with `units`, `calendar` and a `bounds` attribute, and a **`time_bnds` variable holding the start and end of the averaging period**. This is not a formality --- it is the specific thing `ncclimo` needs in order to read Omega output at all, and it has to hold the right values, not merely be present. See the implementation section, where this is now confirmed against real output rather than assumed; - use a file-name convention that encodes the year and month, or that groups whole years into one file, and that is stable across a simulation. Omega's Analysis module already provides temporal reduction with a configurable `ReductionPeriod`, including `1Month`, for the `GlobalStats` group. The work required is to make the same reduction available for full model fields. Checkpointing the monthly reduction so that it resumes correctly across a restart is a longer-term goal and is not required for this deliverable: the simulations to be analyzed are not expected to restart more often than once a month, so a reduction period that divides the restart interval is sufficient. #### Mixed-layer depth Omega should provide a **mixed-layer depth** diagnostic computed in situ from the instantaneous state, using a density-threshold criterion, and make it available for monthly averaging. This is the scientifically correct source: mixed-layer depth is a strongly nonlinear function of the temperature and salinity profiles, so a mixed-layer depth computed from monthly-mean temperature and salinity is not the monthly mean of the mixed-layer depth. Whether this can be delivered in Omega by September 15 is a question for the team rather than something this design can settle. **If it cannot**, Polaris shall fall back to computing mixed-layer depth offline from the monthly-mean conservative temperature and absolute salinity, using the same density threshold criterion, and shall label the resulting plots to make clear that they are computed from monthly means. The fallback is described in the mixed-layer-depth algorithm design, along with what it does and does not capture. ### Requirement: climatology Date last modified: 2026/08/11 Contributors: Xylar Asay-Davis, Claude Polaris shall compute monthly, seasonal, and annual climatologies from a simulation's monthly means, over a start and end year given by the user. The seasons shall include, at minimum, the annual mean (`ANN`) and the four standard three-month seasons (`DJF`, `MAM`, `JJA`, `SON`), and the user shall be able to request additional seasons. The twelve monthly climatologies shall also be available. Climatologies shall be computed only for the fields that are needed, so that the cost of the computation scales with what is being analyzed rather than with the full contents of the monthly-mean files. ### Requirement: climatology-maps Date last modified: 2026/08/25 Contributors: Xylar Asay-Davis, Claude Polaris shall produce global maps of climatological fields on the native MPAS mesh, for each requested field and each requested season. A field with a vertical dimension has to be reduced to a horizontal map before it can be plotted, and there is more than one useful way to do that. Polaris shall provide a general **vertical reduction** for this purpose, of which the user shall be able to request any combination: - **the sea surface** --- the topmost valid layer of each column; - **a fixed geometric elevation** --- a given elevation $z$ (positive up, so negative within the ocean), obtained by linear interpolation in the vertical; - **a fixed layer index** --- a given vertical index, common to all columns; - **the seafloor** --- the bottommost valid layer of each column; - **an integral over an elevation range** --- the mass-weighted integral of the field between two elevations, which is how ocean heat content maps are produced. The "topmost" and "bottommost" valid layers shall respect `minLevelCell` and `maxLevelCell`, so that columns under ice-shelf cavities and columns with partial bottom cells are handled correctly. Where a requested elevation falls outside a column --- below the seafloor, or under land --- the map shall be masked rather than showing an extrapolated value. Sea surface temperature and sea surface salinity are obtained by requesting the sea surface for the temperature and salinity fields; they are not separate fields. Where a model does not output reconstructed zonal and meridional velocity, Polaris shall reconstruct them from the normal velocity on edges. ### Requirement: ocean-heat-content-maps Date last modified: 2026/08/25 Contributors: Xylar Asay-Davis, Claude Polaris shall compute ocean heat content integrated over elevation ranges from a climatology of conservative temperature, and shall produce a global map for each elevation range and each requested season. This is the elevation-range case of the vertical reduction required above, and is delivered as a field of the climatology maps rather than as a product of its own. A heat content map is a climatology map of a field that happens to be derived, and separating the two would mean two code paths, two step trees, and two config conventions for the same operation. The elevation ranges shall be set by config options. The defaults shall be the whole ocean, the surface to $-700$ m, $-700$ m to $-2000$ m, and $-2000$ m to the seafloor. A range boundary that falls in the interior of a model layer shall contribute that layer in proportion to the fraction of the layer within the range, and a range boundary below the seafloor of a given column shall be truncated at the seafloor. ### Requirement: global-stats-time-series Date last modified: 2026/08/11 Contributors: Xylar Asay-Davis, Claude Polaris shall plot time series of the global statistics that Omega's `GlobalStats` analysis group writes, over a start and end year given by the user, for a list of fields given by the user. Polaris shall ship a config file defining a default list of fields and statistics, so that a user who has not thought about which quantities to plot gets a useful set. A field or statistic in that list that the simulation did not write shall be **skipped with a message in the log, not treated as an error**: the defaults describe what we would like to see, and any given simulation will have written some subset of it. A field the user has asked for explicitly is treated the same way, since the user has no more control over what the completed simulation wrote than we do. For each field, the plot shall show the global mean together with the global minimum, maximum, and standard deviation, and shall show the change relative to the beginning of the time series as well as the absolute values, since drift is usually what the reader is looking for. Statistics that are absent are simply omitted from that field's plot. The time axis shall be labeled in simulation years. ### Requirement: ocean-heat-content-time-series Date last modified: 2026/08/11 Contributors: Xylar Asay-Davis, Claude Polaris shall compute a time series of globally integrated ocean heat content, over the same elevation ranges as the ocean heat content maps, from each monthly-mean conservative temperature field over a start and end year given by the user, and shall plot that time series. The plot shall show both the absolute heat content and the anomaly relative to the start of the time series, since the anomaly is the quantity of interest for drift and for the planetary energy budget while the absolute value is dominated by the mean state. The computation shall stream over the monthly files rather than loading the full record into memory, and shall cache the reduced monthly values so that extending the time series with additional simulation years does not require reprocessing the years already covered. ### Requirement: moc-plot Date last modified: 2026/08/11 Contributors: Xylar Asay-Davis, Claude Polaris shall produce a latitude-elevation plot of the global meridional overturning streamfunction from the MOC diagnostic that Omega computes in situ, averaged over the climatology period. The plot shall be in Sverdrups, with a diverging color map centered on zero, contour lines at a configurable interval, and elevation on the vertical axis so that the sea surface is at the top. Omega shall provide, alongside the streamfunction, the mean geometric elevation of each layer interface over the same period, so that the vertical axis is meaningful. Polaris shall not reconstruct it. Polaris shall not compute the MOC itself. The overturning streamfunction requires the full three-dimensional velocity field at every time step to be computed correctly, and Omega computes it in situ for exactly that reason. ### Requirement: regression-test Date last modified: 2026/08/25 Contributors: Xylar Asay-Davis, Claude Polaris shall provide a regression test that runs a short Omega simulation and then analyzes its output, so that the analysis capability is exercised end-to-end by something we run ourselves. Everything else in this document consumes output from a simulation Polaris did not run, on a machine where that simulation happens to live. That is the point of the capability, and it is also why nothing in it would otherwise be covered by a suite: a regression that broke the climatology, the vertical reduction, or the accumulator would be found by a person analyzing a coupled run, which is the most expensive place to find it. The test shall run a coarse-resolution Omega forward run configured to write the monthly-mean output this analysis reads, and shall then run the analysis products against it. It shall be a suite of its own rather than an addition to `omega_pr` or `omega_nightly`, since it costs a forward run and it is blocked on Omega capabilities that the PR suite must not be blocked on. Adding it to `omega_nightly` once it is stable and its cost is known is the expected follow-up. The test verifies that the products are produced and are self-consistent, not that they are scientifically meaningful: a simulation short enough to run in a test is too short for the diagnostics to say anything. ## Algorithm Design ### Algorithm Design: climatology Date last modified: 2026/08/11 Contributors: Xylar Asay-Davis, Claude Climatologies are computed with `ncclimo` from the NCO package, which Polaris already depends on. `ncclimo` is used rather than an `xarray` implementation because it is the tool the rest of the E3SM post-processing workflow uses, it is substantially faster than a naive `xarray` implementation on large files, it handles the season-weighting conventions correctly, and using it keeps our climatologies comparable with those produced by zppy. The monthly climatology for month $m$ is the unweighted mean over the requested years of the monthly means for that month: $$ \overline{\phi}_m = \frac{1}{N_{yr}} \sum_{y=y_0}^{y_1} \phi_{y,m} $$ The seasonal climatology weights each month by the number of days in that month: $$ \overline{\phi}_s = \frac{\sum_{m \in s} d_m \overline{\phi}_m} {\sum_{m \in s} d_m} $$ and the annual mean is the same expression with $s$ running over all twelve months. `ncclimo` implements both, using the calendar of the input files. Two conventions matter and are worth stating explicitly: - **December handling.** `DJF` needs a December, and the December that is contemporaneous with a given January and February belongs to the previous year. `ncclimo`'s `-a sdd` ("seasonally discontinuous December") option takes December from the same year as January and February, which means every year in the requested range contributes exactly one December and no data outside the range are needed. This is the convention MPAS-Analysis uses, and we adopt it for consistency. - **Averaging is over the fields written, not over derived quantities.** A climatology of layer thickness and a climatology of temperature are not the same as a climatology of their product. This matters for heat content and is discussed under the heat content algorithm below. ### Algorithm Design: repeated-analysis Date last modified: 2026/08/25 Contributors: Xylar Asay-Davis, Claude The whole approach rests on separating what depends on the requested range from what does not: | Product | Depends on the range? | Cost | | --- | --- | --- | | Monthly means (model output) | no | read-only input | | Reduced monthly ocean heat content | no, keyed by month | expensive | | Offline monthly mixed-layer depth | no, keyed by month | expensive | | `ncclimo` climatologies | yes | expensive | | Climatology and heat content maps | yes | cheap, from the climatology | | Global stats time series | yes | cheap | | MOC time average | yes | cheap | The rows in the middle are the ones that matter. A given month's vertically integrated heat content is the same quantity no matter which range asked for it, so it should be computed once and reused forever, while everything else is either cheap to redo from those monthly values or is the climatology itself. Every step is keyed by what the user asked for --- a date range --- and the expensive range-independent work is made incremental *inside* a step by the **seeded accumulator** of principle 6 in {ref}`design-ocean-analysis`. The accumulator finds the cache files left by earlier runs of the same product, inherits the months they cover, computes only the rest, and writes a complete cache for its own range. An earlier draft of this design instead gave every simulation year its own shared step at a year-keyed subdirectory, so that reuse fell out of Polaris's completion markers with no machinery of our own. That was appealing but wrong on two counts. It paid a step's overhead --- a directory, a pickle, a config copy, a log file --- for a chunk no user ever asked for, and it turned the directory tree into a bucket named after a mechanism. Worse, the completion marker was then the *only* validity check: change the heat content kernel or a constant and every one of those directories would still report itself complete. The accumulator is cheaper and, with the provenance stamp described below, safer. The steps that remain are **shared steps** in the sense of [Shared steps](shared_steps.md), created with `Component.get_or_create_shared_step()` at a subdirectory computed from the range, so that the climatology for a range is built once no matter how many products read it. The implementation section works through the details. #### How many steps, and how much do they do? Principle 9 in {ref}`design-ocean-analysis` asks for step counts in the low hundreds and for steps that do at least tens of seconds of work at production resolution. For a typical analysis --- a 60-year record with a 20-year climatology --- this design produces: | Steps | Count | | --- | --- | | `ncclimo` climatology | 1 | | Climatology maps, one per field group | 6 | | Heat content series accumulator | 1 | | Offline mixed-layer depth, accumulator plus its climatology (fallback only) | 2 | | Global stats, MOC | 2 | | **Total** | **12** | The count is **independent of the length of the record**, which is the property worth having. It grows with the number of field groups and products, and --- in a later phase, when accumulators are split for a scheduler --- with how much concurrency the machine can use. Neither is something the length of a simulation can run away with. An earlier draft, with one shared step per simulation year, gave about 130 steps for the same analysis and grew linearly, so a century-long record would have given about 230. That was never going to break Polaris --- a few hundred steps is ordinary, and the existing `ocean` component builds several hundred across all its tasks --- but the count grew with something the user has no control over, which is the wrong shape even when the numbers are small. Step *size* is the other half. On a first run the heat content accumulator reads the whole record, which at 6to18 km over several decades is hours; on a re-run over an overlapping range it reads only the new months. A climatology map step plots the seasons and reductions of one field group, which is seconds to minutes per plot at production resolution. Neither is close to the regime where a step's overhead --- a work directory, a pickle, a config copy, a log file --- would be a meaningful fraction of its runtime. Principle 9's other target, that no step be a large fraction of the suite, is knowingly missed: the heat content accumulator is the bulk of a first run. That is acceptable only because Phase 1 is serial, so there is nothing the imbalance could have been traded against. It is the first thing to fix when there is a scheduler, and it is fixed by splitting the accumulator, which costs nothing in reuse. The count grows with the number of field groups and, in later phases, with whatever new products are added --- not with the record, the seasons, the elevations, or the regions. That is what principle 4 is protecting: the dimensions that multiply are loops inside steps, so adding regional analysis multiplies the *plots* without multiplying the *steps*. The climatology is the one expensive computation that genuinely depends on the range, and it is recomputed for each new range. In principle a climatology over a new range could be assembled incrementally from per-year seasonal partial sums, but `ncclimo` has no such mode, and writing our own incremental climatology to save a rerun is not a trade we should make for this deliverable. Two ranges' climatologies coexist without special handling, because `ncclimo` already encodes the range in its output file names (`____climo.nc`). ### Algorithm Design: climatology-maps Date last modified: 2026/08/30 Contributors: Xylar Asay-Davis, Claude #### Vertical geometry Everything in this document that involves a vertical position needs the geometric elevation of layer midpoints, $z^{mid}_{k}$, and of layer interfaces, $z^{int}_{k}$, for each column. These are read directly from `zMid` and `zInterface` --- Omega's `GeomZMid` and `GeomZInterface` --- which is why the `omega-monthly-means` requirement asks for them as monthly-mean output. Polaris cannot reconstruct them from the monthly means of the other fields. Omega builds the geometric coordinate by accumulating upward from $-\mathrm{BottomGeomDepth}$ using layer thicknesses $h = \rho_0 \, \mathrm{SpecVol} \times \mathrm{PseudoThickness}$, and it does so the same way regardless of which vertical coordinate the simulation uses --- the choice of z-star, p-star, or sigma determines how `PseudoThickness` is initialized and how it evolves, not how geometric elevation is computed from it. Reconstructing $z$ offline therefore requires specific volume, which means evaluating the TEOS-10 equation of state on the monthly-mean state. That is not the monthly mean of $\mathrm{SpecVol} \times \mathrm{PseudoThickness}$, so the reconstruction would introduce an error that has nothing to do with the diagnostic being computed. Having Omega write the geometric coordinate removes the problem entirely, since the monthly mean of $z$ is exactly the mean layer geometry we want. Note the direction of the dependence, which is the reason for the conventions stated at the top of this document: geometric thickness is *derived* from pseudo-thickness and specific volume, and it is the derived quantity that cannot be recovered offline. Mass per unit area, $\rho_0 \tilde{h} = \rho h$, needs no equation of state at all. That is why the heat content integral is written in terms of it and why only quantities that genuinely need the geometry --- elevation slices, and the partial layers at a heat content range boundary --- read `zMid` and `zInterface`. Because the input to the map steps is a climatology, the $z^{mid}_{k}$ they use is the climatological-mean layer geometry. This does mean that a $-100$ m map is a map on the time-mean position of the $-100$ m surface rather than the time mean of maps on the instantaneous $-100$ m surface. The difference is small away from regions with a large seasonal cycle in layer thickness. #### Elevation selection Let $f_k$ be the field in a given column, $k_{min}$ = `minLevelCell` and $k_{max}$ = `maxLevelCell` be the zero-based indices of the topmost and bottommost valid layers, and let the requested elevation specification be one of `top`, `bottom`, `k`, or an elevation $z$ in meters (negative within the ocean). - `top`: $f = f_{k_{min}}$. - `bottom`: $f = f_{k_{max}}$. - `k`: $f = f_n$, masked in columns where $n < k_{min}$ or $n > k_{max}$. - elevation $z$: find $k_1$, the largest valid index with $z^{mid}_{k_1} \ge z$, and interpolate linearly between $k_1$ and $k_1 + 1$: $$ w = \frac{z - z^{mid}_{k_1+1}}{z^{mid}_{k_1} - z^{mid}_{k_1+1}}, \qquad f = w f_{k_1} + (1 - w) f_{k_1+1} $$ If $z$ is above the midpoint of the topmost layer, $f = f_{k_{min}}$; if it is below the midpoint of the bottommost layer but above the seafloor, $f = f_{k_{max}}$. If $z < z^{int}_{k_{max}+1}$ --- below the seafloor --- or the column is land, the result is masked. Clamping rather than masking between the top layer midpoint and the sea surface matters in practice: a request for $0$ m or $-5$ m would otherwise be masked everywhere, which is not what a user asking for a near-surface map intends. The search for $k_1$ is vectorized over columns as the count of valid layer midpoints at or above $z$, which avoids a Python loop over cells: ```python k1 = k_min + (z_mid >= z).sum(dim='nVertLevels') - 1 ``` with `z_mid` set to `NaN` outside the valid range so invalid layers do not contribute, followed by a clip into `[k_min, k_max - 1]` and a mask for the out-of-column cases. The $k_{min}$ term is easy to drop and worth keeping in view. Because `z_mid` is `NaN` above the valid range as well as below it, the count is of *valid* midpoints at or above $z$, which is $k_1 - k_{min} + 1$ rather than $k_1 + 1$. Without the term the index lands $k_{min}$ layers too high, which is invisible wherever $k_{min} = 0$ and wrong in every ice-shelf cavity: for $k_{min} = 3$, $k_{max} = 5$ and a $z$ below every valid midpoint, the count is 3, so the index clips to $k_{min}$ and clamps to $f_3$ where the rule above says $f_{k_{max}} = f_5$. Two details of the implementation follow from the rules above rather than adding to them. The two clamping cases come out of clipping the weight $w$ into $[0, 1]$ after the index is clipped, which reproduces them exactly without a branch. A column with a single valid layer, $k_{min} = k_{max}$, does need its own branch: the clip range `[k_min, k_max - 1]` is empty, and the answer is $f_{k_{min}}$ rather than an interpolation against a layer below the seafloor. #### Velocity The map steps plot zonal and meridional velocity at cell centers, read directly from the monthly means. Polaris does not reconstruct them, and Phase 1 has no offline reconstruction path at all. An earlier draft had one, reconstructing from `normalVelocity` on edges with the least-squares weights designed in [Vector Reconstruction](vector_reconstruction.md), so that this product would not block on Omega work. It is not needed: the Omega side is in progress in [Omega #525](https://github.com/E3SM-Project/Omega/pull/525), which adds velocity-component reconstruction for I/O, so reconstructed velocities become a required output rather than a preferred one, and a fallback for a case that will not arise is code we would write, test and maintain for nothing. That work is still a draft, so this is the one place the design depends on something not yet landed upstream. The consequence is contained: until it does, the mock-up files carry `NormalVelocity` and no components, so the velocity maps are the one product that reports missing fields and skips, in the way described under `omega-monthly-means`. Nothing else waits on it. Writing the offline path against that gap would cost more than the wait, and would leave us maintaining two ways to obtain the same field. Polaris's reconstruction itself is not going away and is being fixed independently --- [Polaris #721](https://github.com/E3SM-Project/polaris/pull/721) corrects vector reconstruction on planar meshes, found while doing the Omega work. It stays available for tasks that need it; this analysis simply does not read edge velocities. The accuracy question that would otherwise decide this does not arise either. Reconstruction is linear, so reconstructing from a climatology of normal velocity gives exactly the climatology of the reconstructed velocity --- doing it in the model costs nothing in accuracy, and it saves writing and reading an edge field nothing else here wants. ### Algorithm Design: mixed-layer depth (fallback only) Date last modified: 2026/08/11 Contributors: Xylar Asay-Davis, Claude *This algorithm applies only if Omega cannot deliver an in-situ mixed-layer depth diagnostic in time. If Omega provides one, mixed-layer depth is an ordinary monthly-mean field and flows through the climatology and map steps like any other, and none of this section applies.* The fallback computes mixed-layer depth from each monthly-mean profile of conservative temperature and absolute salinity using the same density-threshold criterion Omega plans to use: the mixed layer extends to the elevation at which the potential density referenced to $10$ m exceeds the density at $10$ m by $\Delta \rho = 0.03$ kg m⁻³, with linear interpolation between the bounding layers. Density is evaluated with `gsw`, which Polaris already depends on and which implements the same TEOS-10 formulation Omega uses. The mixed-layer depth is computed for each month and then averaged over the climatology period, rather than being computed from the climatology. Computing it from the seasonal or annual climatology would be a second, much larger approximation on top of the one described next. What this fallback does *not* capture, and what should be said on the plots and in the netCDF metadata: - Monthly-mean profiles are smoother than instantaneous profiles. Averaging over a month blends the stratification before, during, and after a mixing event, so the mixed-layer depth derived from the mean profile is not the mean of the mixed-layer depths. - The effect is largest exactly where mixed-layer depth matters most: winter deep-convection regions, where a few days of deep mixing set the monthly mean of the true mixed-layer depth but leave only a muted signature in the monthly-mean profile. - Monthly maximum mixed-layer depth, which MPAS-Analysis reports and which is the more useful deep-convection diagnostic, cannot be recovered at all from monthly means. The fallback is therefore adequate for a first look at the seasonal cycle of the mixed layer in a coupled run and is not adequate as a deep-convection diagnostic. Replacing it with Omega's in-situ diagnostic should be the first follow-up after September 15 if the fallback is what ships. ### Algorithm Design: ocean-heat-content Date last modified: 2026/08/25 Contributors: Xylar Asay-Davis, Claude *This algorithm is shared by the heat content maps and the heat content time series; the two differ in what they integrate over and in what they read, not in the kernel.* Ocean heat content per unit area is a *mass*-weighted integral of conservative temperature. Over an elevation range $[z_{bot}, z_{top}]$ with $z_{bot} < z_{top}$, where either boundary may be a fixed elevation, the free surface, or the seafloor, $$ Q(z_{bot}, z_{top}) = c_p^0 \int_{z_{bot}}^{z_{top}} \rho \, \Theta \, dz = \rho_0 c_p^0 \int_{\tilde{z}_{bot}}^{\tilde{z}_{top}} \Theta \, d\tilde{z} \approx \rho_0 c_p^0 \sum_k \Theta_k \, \tilde{w}_k $$ where $\Theta$ is conservative temperature, $\rho$ is in-situ density, and $\tilde{w}_k$ is the pseudo-thickness of the overlap between layer $k$ and the requested range. The middle step is the change of variable $\rho \, dz = \rho_0 \, d\tilde{z}$ from the conventions, and it is an identity, not an approximation: a mass-weighted integral in $z$ is a plain integral in $\tilde{z}$ scaled by $\rho_0$. The final step is the layer quadrature, so the only error is the usual one of treating $\Theta$ as uniform within a layer. In particular, the reference density in the discrete sum is **not** a Boussinesq approximation. Since $\rho_0 \tilde{h}_k = \rho_k h_k$ is the mass per unit area of layer $k$ exactly --- for Omega by the definition of pseudo-height, for MPAS-Ocean because it is Boussinesq --- a range covering whole layers gives the mass-weighted integral with no reference-density error at all. This is the substantive difference from the MPAS-Analysis formulation, which weights $\Theta$ by a *geometric* thickness and multiplies by a reference density, and so carries an in-situ-versus-reference density error of a few tenths of a percent. The geometric coordinate enters only through the partial layers at the range boundaries, because the range is specified in $z$ while the integral is in $\tilde{z}$. That the range is specified in geometric elevation is a deliberate choice rather than an oversight. "Ocean heat content, 0 to 700 m" means 700 *geometric* meters in MPAS-Analysis, in the observational products this diagnostic is compared against, and in the literature, and the same is true of a map "at 100 m". Specifying ranges in pseudo-depth instead would be more natural for a mass-conserving model --- a fixed $\tilde{z}$ range is a fixed pressure range and therefore exactly a fixed mass per unit area, which is a cleaner control volume for a heat budget --- and it would remove the geometric coordinate from this algorithm entirely. We do not do it, because it would silently redefine a number everyone else reports geometrically, for a difference that is a fraction of a percent in the upper ocean and one to two percent at abyssal depths. Config options are documented as geometric elevations, and the conversion happens here, in one place. Let $$ w_k = \max\left(0, \; \min\left(z^{int}_{k}, z_{top}\right) - \max\left(z^{int}_{k+1}, z_{bot}\right)\right) $$ be the geometric thickness of the overlap, for layers within `[minLevelCell, maxLevelCell]` and $w_k = 0$ elsewhere, and let $h_k = z^{int}_{k} - z^{int}_{k+1}$ be the layer's geometric thickness. Then $$ \tilde{w}_k = \frac{w_k}{h_k} \, \tilde{h}_k $$ is the pseudo-thickness of the overlap, and $\tilde{w}_k = \tilde{h}_k$ whenever the layer lies entirely within the range. Splitting a layer by a geometric fraction rather than by a pseudo-height fraction is exact with respect to the model's own discretization, not a further approximation: within a layer Omega uses a single specific volume $\alpha_{i,k}$, so $h = \rho_0 \alpha \tilde{h}$ makes $z$ linear in $\tilde{z}$ across the layer and the two fractions are equal. For MPAS-Ocean the question does not arise, since `layerThickness` is the geometric thickness and $\tilde{w}_k$ reduces to $w_k$. Where these weights are formed from a *climatology* rather than from a single month, the ratio $w_k/h_k$ is a ratio of monthly means rather than a mean of ratios. This is second order --- it affects only the two boundary layers of a range, and only through the seasonal cycle of layer thickness --- and it is one face of the covariance term discussed below. The $w_k$ expression handles all of the cases the requirement calls for without special casing: a range boundary in the interior of a layer contributes a partial thickness; a range extending below the seafloor is truncated because $z^{int}_{k_{max}+1} = -H$; a `bottom` boundary is expressed as $z_{bot} = -\infty$; a `top` boundary as $z_{top} = +\infty$, which resolves per column to the free surface $z^{int}_{k_{min}}$; and a column whose seafloor lies above $z_{top}$ contributes zero. Expressing the upper boundary as `top` rather than as $0.0$ matters more than it looks. A range written `0.0:-700.0` would exclude the water between the resting sea surface and the free surface, and would exclude a different amount of it in each column and in each season. `top` is what "0 to 700 m" means everywhere it is reported, and it makes the whole-column range `top:bottom` cover every valid layer, so that every layer is whole and the geometric coordinate drops out of that answer entirely. The globally integrated heat content used for the time series is the area-weighted sum $$ Q_{tot} = \sum_i A_i \, Q_i $$ over cells $i$ with area `areaCell`. #### Choice of constants Omega uses TEOS-10, in which conservative temperature is *defined* so that potential enthalpy is $c_p^0 \Theta$ with the exact constant $c_p^0 = 3991.86795711963$ J kg⁻¹ K⁻¹. Using that constant with conservative temperature is therefore not an approximation but the definition of heat content. Polaris's Physical Constants Dictionary provides `seawater_specific_heat_capacity_reference` = 3996.0 J kg⁻¹ K⁻¹ and `seawater_density_reference` = 1026.0 kg m⁻³, which are the values MPAS-Analysis uses. The specific heat capacity differs from the TEOS-10 constant by 0.1%, which is small compared to other uncertainties but is a systematic offset in any comparison. **Decision:** Phase 1 uses the PCD value for $c_p^0$. The TEOS-10 constant is not in the PCD today, and adding a constant to it is not something we can complete by September 15, so using the PCD is what keeps Polaris consistent with the constants the rest of E3SM is using in the meantime. It is exposed as a config option so that a user can experiment with the TEOS-10 value without a code change. In the long run we do want $c_p^0$: because Omega carries conservative temperature, $c_p^0 \Theta$ is potential enthalpy per unit mass by definition, so the mass-weighted integral above is the heat content rather than an approximation to it. Adding $c_p^0$ to the PCD and switching to it is recorded as a deferred item in {ref}`design-ocean-analysis`. **$\rho_0$ is not a free parameter.** Unlike $c_p^0$, the reference density in the discrete sum is not a modeling choice we are making --- it is the constant that *defines* pseudo-height in Omega, and the Boussinesq reference density in MPAS-Ocean. Using any other value would make $\rho_0 \tilde{h}$ something other than the layer's mass per unit area. It must therefore be the same $\rho_0$ the model used, and it is not exposed as a config option. Polaris already takes that value from the PCD's `seawater_density_reference` when it builds Omega's p-star vertical coordinate in `polaris/ocean/vertical/pstar.py`, so the analysis reads it from the same place; if Omega's `RhoSw` ever diverges from the PCD value, that is a bug in Polaris's vertical coordinate before it is a bug in this diagnostic. The practical consequence of the mass-weighted form, noted above, is that the in-situ-versus-reference density error in the MPAS-Analysis formulation --- a few tenths of a percent, nearly uniform, so it biased the absolute heat content much more than the anomaly --- is gone at no cost, since the model already writes the mass-like thickness we need. #### Heat content from a climatology versus from monthly means The heat content maps are computed from the climatology of $\Theta$ and of the layer thicknesses, per the requirement. Because heat content is a product of those two, this omits the covariance term: $$ \overline{\Theta \tilde{h}} = \overline{\Theta}\,\overline{\tilde{h}} + \overline{\Theta' \tilde{h}'} $$ This term is already dropped at every timescale shorter than a month, the moment the analysis works from monthly means rather than from model time steps. Neglecting it again from month to season is the same approximation applied over a longer averaging period, not a new one, so there is little sense in being scrupulous about the second while accepting the first in silence. Its size settles the question. Over a fixed elevation range the term enters only through the partial layers at the range boundaries and through the free surface, and for the $0$ to $-700$ m range a seasonal sea surface height of order $0.1$ m against a near-surface seasonal temperature anomaly of a few kelvin gives order $10^6$ J m⁻², against a total near $2.9 \times 10^{10}$ J m⁻². That is order $10^{-4}$ of the signal --- an order-of-magnitude sketch rather than a bound, but several orders below anything that would change a conclusion drawn from these maps. The sub-monthly version is of similar magnitude and partly cancels, being eddy-driven rather than systematic. We therefore compute the maps from the climatology and do not plan to revisit it. The alternative --- computing per-month vertically integrated maps and averaging those --- would cost a full pass over the three-dimensional monthly output for every season plotted, which is a large price for $10^{-4}$. What this does require is that the climatology include both `PseudoThickness` and the geometric interfaces: the former sets the mass weight, the latter the partial-layer fraction at the range boundaries. The heat content time series does not have the issue at all, because it integrates each monthly mean separately and averages afterward. ### Algorithm Design: moc-plot Date last modified: 2026/08/11 Contributors: Xylar Asay-Davis, Claude Omega's MOC analysis group computes the overturning streamfunction on latitude bins and layer interfaces, in Sverdrups, and writes it with the same temporal-reduction machinery as `GlobalStats`. Polaris averages the reductions over the requested climatology years and plots the result. Only the **global** streamfunction is plotted. Omega's MOC group can compute the streamfunction for named regions, but regional analysis is out of scope for this deliverable, so no region masks, region names, or region config options appear here. Regional overturning --- the Atlantic MOC in particular --- comes with the rest of the regional analysis in a later delivery, and the plotting primitive introduced below is written so that it will not need changing then. Two details need care: - **The vertical coordinate of the plot.** The streamfunction lives on layer interfaces (`nVertLevelsP1`), which are not at a fixed elevation. Omega shall provide the mean interface elevation alongside the streamfunction; the natural way to produce it is a climatology of `GeomZInterface` over the same period as the MOC, reduced to a mean per layer interface. Polaris plots what Omega provides and does not reconstruct it: a Polaris-side reconstruction would either repeat that averaging on a different period or fall back on resting thicknesses, and in both cases the plot's vertical axis would quietly disagree with the diagnostic it is labeling. - **Averaging over the requested period.** Omega writes per-period reductions (e.g. monthly or annual means). The streamfunction is linear in the velocity, so a mean of the reductions weighted by the length of each reduction period is the streamfunction of the mean flow over the period. The exact variable, dimension, and coordinate names in Omega's MOC output must be confirmed against the implementation before this step is written; see the open questions. ## Implementation ### Implementation: analysis-suite Date last modified: 2026/08/29 Contributors: Xylar Asay-Davis, Claude #### Tasks, steps, and the suite New tasks live in `polaris/tasks/ocean/analysis/`, added to the ocean component by `add_analysis_tasks(component)` in `polaris/tasks/ocean/add_tasks.py`. The work-directory layout is: ```none ocean/analysis/ ├── climatology/0021-0040/ (shared: ncclimo) ├── climatology_maps/0021-0040/ │ ├── temperature/ salinity/ velocity/ (one step per field group) │ ├── ssh/ mixed_layer_depth/ │ └── heat_content/ ├── mixed_layer_depth/0021-0040/ (only if computed offline) │ ├── monthly/ (accumulator) │ └── climatology/ (second ncclimo call) ├── heat_content_series/0001-0060/ ├── global_stats/0001-0060/ └── moc/0021-0040/ ``` The layout follows one rule --- **`//[]`** --- and the third level exists only for the one product that is chunked. This is the work tree, so it is built for predictability and for finding a step's log, not for browsing; browsing is the staging tree's job. See principles 1 and 5 in {ref}`design-ocean-analysis`. Two levels from an earlier draft are gone, and it is worth saying why, since the reasoning generalizes: - **`years/`**, which held one shared step per simulation year, is gone because those steps are gone --- the expensive per-year work is now a seeded accumulator inside a single step. Renaming the level would not have helped: a directory whose best name describes how the work was chunked is a level that should not exist. - **`maps/`, `time_series/`, and `plot/`**, the single-step levels between a product and its period, are gone because a level with one child and a name that repeats its parent earns nothing and costs a level of depth on every path. Every step is a shared step in the sense of [Shared steps](shared_steps.md), created at `ocean/analysis`, the highest level at or below which all of the tasks that use them live. This matters most for the climatology, which every field group of `climatology_maps` reads and which therefore runs once no matter how many of them are run. Tasks are a thin grouping over these steps and no longer introduce directory levels of their own: `climatology_maps` is one task with a step per field group, and `heat_content_series`, `global_stats`, and `moc` are one task with one step each. Running a subset is `polaris serial --steps`. The suite is `polaris/suites/ocean/omega_analysis.txt`, named to match the existing `omega_pr` and `omega_nightly` suites: ```none ocean/analysis/climatology_maps ocean/analysis/heat_content_series ocean/analysis/global_stats ocean/analysis/moc ``` A user analyzes a simulation with: ```bash polaris suite -c ocean -t omega_analysis -w -f analysis.cfg \ --model omega polaris serial ``` No `-p`/`--component_path` is given, because no model executable is needed. The `--model` flag (or an `[ocean] model` option in the config file) tells Polaris which model's naming conventions to expect on read. Analyzing two simulations means two work directories. The simulation name is not part of the work-directory path, because task subdirectories are fixed when the component is constructed, before the user's config file has been read. #### Config options Defaults ship in `polaris/tasks/ocean/analysis/analysis.cfg` and the user overrides the ones that describe their simulation. The sections below are the proposed starting point: ```ini [ocean_analysis] # The absolute path to the simulation's Omega configuration file. This is the # analysis' only description of where the simulation's output lives: Polaris # reads the mesh, the output streams and their file-name templates from it, so # that none of them have to be restated here. It is required, since an Omega # run always has one. MPAS-Ocean output is not supported; reading it would # need a translator from its namelists and streams into the same form. omega_config_filename = # The absolute path to the directory containing the simulation's output. # Defaults to the directory containing the Omega configuration file, which is # what its relative file names are resolved against. simulation_path = # A short name for the simulation, used in plot titles and file names simulation_name = omega # Where to publish plots, their netCDF files, thumbnails, and the generated # gallery. Defaults to /analysis_output. Point this somewhere # web-servable if you want to share the results. The thumbnail options that # go with it are given under `publication`. output_path = # The horizontal mesh file, absolute or relative to simulation_path. Defaults # to the mesh the Omega config names. mesh_filename = # The vertical-coordinate file (Omega only), absolute or relative to # simulation_path vert_coord_filename = [ocean_analysis_climatology] # The first and last year of the climatology, inclusive start_year = 1 end_year = 10 # The seasons to compute, in addition to the 12 monthly climatologies seasons = ANN, DJF, MAM, JJA, SON # The seasons to plot; may include the monthly climatologies JAN through DEC plot_seasons = ANN, DJF, JJA # The fields for which climatology maps are produced, using MPAS-Ocean # (Polaris standard) names fields = temperature, salinity, velocityZonal, velocityMeridional, ssh, mixedLayerDepth # The elevations at which fields with a vertical dimension are plotted. # Elevations are in m, positive up, so values within the ocean are negative: # top the topmost valid layer of each column (the sea surface) # bottom the bottommost valid layer of each column (the seafloor) # an elevation in m, linearly interpolated # k a fixed, zero-based vertical index elevations = top, -100.0, -500.0, -2000.0, bottom # Whether to compute mixed-layer depth offline from monthly-mean temperature # and salinity, for simulations whose output does not include it compute_mixed_layer_depth = False # The density threshold used when computing mixed-layer depth offline, in # kg m-3, relative to a reference elevation of -10 m mixed_layer_depth_threshold = 0.03 [ocean_analysis_ohc] # The elevation ranges over which heat content is integrated, given as # : in m, positive up. "bottom" means the seafloor. These are # geometric elevations, matching the convention used by MPAS-Analysis and by # observational heat content products; the integral itself is mass-weighted. elevation_ranges = top:-700.0, -700.0:-2000.0, -2000.0:bottom, top:bottom # The specific heat capacity used to convert conservative temperature to heat # content. By default, this comes from the Physical Constants Dictionary. # The reference density is not a config option; see the algorithm design. #seawater_specific_heat_capacity = 3996.0 [ocean_analysis_time_series] # The first and last year of the time series, inclusive start_year = 1 end_year = 10 # The fields from the model's global statistics output to plot. Fields the # simulation did not write are skipped with a message, not an error. fields = temperature, salinity, normalVelocity, kineticEnergyCell, ssh # The statistics to plot for each field. As with fields, missing statistics # are skipped. stats = mean, min, max, std [ocean_analysis_moc] # The contour interval in Sv contour_interval = 2.0 # The maximum absolute value of the color map in Sv; the color map is # symmetric about zero #max_streamfunction = 30.0 ``` Per-field plotting options follow the existing Polaris convention of a section per field, as in `realistic_global.cfg`: ```ini [ocean_analysis_map_temperature] colormap_name = cmo.thermal norm_type = linear norm_args = {'vmin': -2., 'vmax': 32.} ``` Sections for the fields we expect to plot ship with defaults; fields without a section fall back to the defaults in `polaris.viz.get_viz_defaults`. The section name is not the field name pasted onto a prefix. Field names are the models' own and are therefore camel case, while Polaris config sections are lower case with underscores everywhere else in the codebase, and the sections this design adds are no exception. A field's section is `ocean_analysis_map_`, so `velocityZonal` gets `[ocean_analysis_map_velocity_zonal]` and `mixedLayerDepth` gets `[ocean_analysis_map_mixed_layer_depth]`. That conversion is mechanical, so no step spells a section out. A small helper module, `polaris/tasks/ocean/analysis/config_sections.py`, provides `camel_to_snake(name)` and `map_section(field)`, and every caller goes through `map_section`, so exactly one place knows the prefix and the spelling rule. It is a leaf module with no Polaris imports, for the same reason `sim_files.py` is one: it can be unit tested on its own, and any step can use it without pulling in a step. #### Field and dimension naming Per the conventions stated in the summary, config options name fields using MPAS-Ocean names, and analysis code uses MPAS-Ocean variable and dimension names throughout. `OceanIOStep.open_model_dataset` performs the translation from Omega names on read, driven by `polaris/ocean/model/mpaso_to_omega.yaml`, which already maps the vertical geometry this design depends on: ```yaml zMid: GeomZMid zInterface: GeomZInterface ``` **The mapping reconciles two spellings; it is not a naming authority.** An entry exists when both models write the same quantity under different names, and only then. Three cases follow from that, and they cover everything this design needs: - **Both models have the field.** It is mapped, and analysis uses the MPAS-Ocean name --- `temperature`, `zMid`, `areaCell`. - **Only Omega has the field.** A mixed-layer depth diagnostic, for instance. There is nothing to reconcile, so there is no entry and analysis uses the Omega name as written. Inventing an MPAS-Ocean-styled synonym would put a name in the codebase that no model ever writes, and the entry recording it would be a rename in appearance only. - **Both models have a similar name for different quantities.** `layerThickness` and `PseudoThickness` are the case in point. There is no entry, deliberately, for the reasons given under the conventions. The practical consequence of the middle case is that Omega's spelling appears in analysis code, and in the values of config options that name fields, for fields Omega alone provides, which is a small inconsistency of style and an accurate one: those fields come from Omega and from nowhere else. Config section names are the exception, because they are Polaris' own rather than either model's: `map_section` converts the field name, so no model's spelling reaches a section header. If such a field later gains an MPAS-Ocean counterpart, or Omega renames it, that is the point at which an entry earns its keep --- and adding one then is a one-line change, so there is nothing to be gained by adding it in advance. No analysis step branches on `config.get('ocean', 'model')` to choose a field or dimension name, except `get_layer_mass`, which exists precisely to be that one branch. Anywhere else, a branch on the model means either that the mapping is missing an entry, or that a field one model does not have is being read unconditionally. #### Locating input files A small helper module, `polaris/tasks/ocean/analysis/sim_files.py`, expands the file-name templates over a year range into lists of files and checks that they exist, reporting the missing years clearly. It is shared by every step that reads simulation output. This is deliberately a separate module rather than a method on a step, so that it can be unit tested and reused. Where the templates come from is the point. The user's input is `omega_config_filename`, the path to the simulation's own Omega configuration, and the same module reads the output streams from it: each stream gives a file-name template, a reduction period and a directory, which is enough to identify the monthly means, the `GlobalStats` output and the MOC output without the user restating any of it. The mesh and the vertical coordinate are read from the same place. **The templates are not config options, and the Omega configuration is required.** An earlier draft had `monthly_mean_template` and its siblings as config options overriding what the Omega configuration says, which cannot work: Polaris config files use `ExtendedInterpolation`, so a bare `$Y` in a value raises `invalid interpolation syntax` and takes down the whole config combine, not merely that option. Escaping it as `$$Y` would work and would be a trap, since a template pasted out of `omega.yml` would then fail confusingly. Removing the override layer rather than escaping it is the better trade in any case. An Omega run always writes an `omega.yml`, so there is no case where the configuration is genuinely unavailable and the fallback would earn its keep; making it required removes the precedence rules, the reporting of which source won, and the one path by which a `$` could reach a config value. What remains in the config file are plain paths with no templating in them --- `simulation_path`, `mesh_filename` and `vert_coord_filename` --- for output or a mesh that has moved since the run. Each may be absolute or relative to `simulation_path`, and the step reports which of the two sources each path came from, so that a surprising file list can be diagnosed without guessing. **MPAS-Ocean output is therefore not supported, and says so.** Reading it needs a translator from its namelists and streams into the form this module reads, which is separate work and is not in scope here. This costs the deliverable nothing: the design already develops against Omega output rather than MPAS-Ocean output, for the reasons given under `omega-monthly-means`. A step set up against a simulation whose `[ocean] model` is `mpas-ocean` reports that rather than failing obscurely later. Reading Omega's configuration is done defensively --- a missing stream, a stream that names no file, and an analysis group that is turned off are told apart and reported as such --- since its schema is Omega's to change and this is the one place Polaris depends on its shape rather than on its output. The names of an analysis group's output streams have to be reconstructed the way Omega's analysis manager builds them, as `_