Helper for outputting a NetCDF "grid file"

[glwagner]

After learning a bit about how people work with NetCDF (😺 ) and discussing with @NoraLoose it seems there may be scope for a helper function that writes out a grid file. Something like:

write_grid_to_netcdf!(filename, grid)

Also internally, this functionality could be added to NetCDFOutputWriter, so that users could write something like

ow = NetCDFOutputWriter(model, grid; filename, with_grid_metrics=true)

and include all the same grid metrics that are saved when using write_grid_to_netcdf!.

Under the hood, I think this is accomplished by creating a special NetCDFOutputWriter and treating the grid metrics like output variables. One possible difference is that unlike other kinds of output the grid metrics do not have a time dimension. So we also need to generalize the scope of outputs for NetCDFOutputWriter to outputs that don't have a time dimension? My NetCDF knowledge is a little hazy so please let me know if this isn't right.

@ali-ramadhan @tomchor

[glwagner]

For extending output utilities to variables that are not time-dependent I think this is just a matter of extending define_output_variable!, which currently assumes that time is a dimension here:

Oceananigans.jl/src/OutputWriters/netcdf_output_writer.jl

Line 487 in 8878541

defVar(dataset, name, FT, (dims..., "time"); deflatelevel, attrib)

So perhaps just a kwarg time_dependent=true will do.

save_output! may also have to be generalized:

Oceananigans.jl/src/OutputWriters/netcdf_output_writer.jl

Line 507 in 8878541

ds[name][colons..., time_index:time_index] = data

[glwagner]

Each grid will need to define a function that returns a list of all the possible metrics too. Because they are different for each grid. Or do we need to save every location even if they are redundant for some grids?

[NoraLoose]

Thanks for opening this issue! I don't think the metrics need to be different for different grids. Here is a list of metrics that would be helpful to have in the grid file, assuming we are on a C-grid:

• x-spacing:
  • centered at T-points
  • centered at u-points
  • centered at v-points
  • centered at vorticity points
• y-spacing:
  • centered at T-points
  • centered at u-points
  • centered at v-points
  • centered at vorticity points
• z-spacing (even though the user can easily infer this from zC and zF in case of a z-coordinate system)
  • dzC
  • dzF
• grid cell area of
  • T-cells
  • u-cells
  • v-cells
  • vorticity cells
• grid cell volumes
• masks

The grid cell area could be skipped if we simply have

grid area of P-cell = x-spacing at P-points * y-spacing at P-points

Is the above equation true for all grid types?

Is it also worth outputting the grid cell volume if there are partial cells?

Am I forgetting any useful metrics @tomchor?

[glwagner]

Is the above equation true for all grid types?

I don't think that is true for a curvilinear grid, for example for LatitudeLongitudeGrid:

Oceananigans.jl/src/Grids/latitude_longitude_grid.jl

Line 437 in 8878541

@inline Azᶜᶜᵃ(i, j, k, grid::LatitudeLongitudeGrid) = @inbounds grid.radius^2 * deg2rad(grid.Δλᶜᵃᵃ[i]) * (hack_sind(grid.φᵃᶠᵃ[j+1]) - hack_sind(grid.φᵃᶠᵃ[j]))

Azᶜᶜᵃ is our notation for the cell area with a z-normal vector, at the staggered location "Center, Center, Any" in "x, y, z". The third direction is "Any" because the z-areas do not depend on the vertical position, because we make the thin shell approximation.

[glwagner]

Is it also worth outputting the grid cell volume if there are partial cells?

Yes, although in principle because we make the thin shell approximation, one can derive the volume from area * dz. Still though I think it is better to output everything rather than miss a few things and create headaches / complicate user scripts that do analysis.

[tomchor]

Is it also worth outputting the grid cell volume if there are partial cells?

I think so. It's simple and it makes users' lives easier :)

Am I forgetting any useful metrics @tomchor?

Not that I can tell. Although I'd standardize and output the same variables for all directions. Oceananigans is very general, so I don't think we should assume a priori that any simulation treats the vertical direction any differently from the others. There'll be some redundancy in the output for lots of simulations, but given how cheap it is to output these metrics, that's okay I think.

[glwagner]

Hmm so the missing metrics are the w locations perhaps then.