Calculate z from sigma* coordinates

1

Hi,

I am using output from a run with time varying sigma* coordinates (MetOffice AMM7 domain) and want to calculate the z position of the top of each layer.
I tried integrating e3w(:,0:mbathy, : , : ) over depth but the result is offset from the local bathymetry, which I took from hbatf(t,y,x) in a auxiliary file called mesh_mask_co6.nc.

I was expecting the z position of a layer of a specific grid (u,v,w,T) to be given by:

z(lev) = bathymetry - \sum_{level=mbathy}^{lev} layer_thickness(level)

with mbathy being the first level that is not masked, but I I haven’t found the correct variables.

Below is the header for data file YYYYMMDD_grid_shelftmb_grid_W followed by headers for auxiliary files coordinates.nc and mesh_mask_co6.nc

<xarray.Dataset>
Dimensions: (time_counter: 24, depthw: 51, y: 3, x: 3,
axis_nbounds: 2)
Coordinates:

  • depthw (depthw) float32 0.0 10.0 20.0 … 1.452e+04 1.502e+04
    nav_lat (y, x) float32 dask.array<chunksize=(3, 3), meta=np.ndarray>
    nav_lon (y, x) float32 dask.array<chunksize=(3, 3), meta=np.ndarray>
    time_centered (time_counter) datetime64[ns] dask.array<chunksize=(24,), meta=np.ndarray>
  • time_counter (time_counter) datetime64[ns] 2001-01-01T00:30:00 …
    Dimensions without coordinates: y, x, axis_nbounds
    Data variables:
    avm (time_counter, depthw, y, x) float32 dask.array<chunksize=(24, 51, 3, 3), meta=np.ndarray>
    avt (time_counter, depthw, y, x) float32 dask.array<chunksize=(24, 51, 3, 3), meta=np.ndarray>
    depthw_bounds (depthw, axis_nbounds) float32 dask.array<chunksize=(51, 2), meta=np.ndarray>
    e3w (time_counter, depthw, y, x) float32 dask.array<chunksize=(24, 51, 3, 3), meta=np.ndarray>
    time_centered_bounds (time_counter, axis_nbounds) datetime64[ns] dask.array<chunksize=(24, 2), meta=np.ndarray>
    time_counter_bounds (time_counter, axis_nbounds) datetime64[ns] dask.array<chunksize=(24, 2), meta=np.ndarray>
    wo (time_counter, depthw, y, x) float32 dask.array<chunksize=(24, 51, 3, 3), meta=np.ndarray>
    Attributes:
    name: shelftmb_grid_W
    description: TMB ocean W grid variables
    title: TMB ocean W grid variables
    Conventions: CF-1.5
    production: An IPSL model
    timeStamp: 2021-Jan-21 06:56:09 UTC
    history: Mon Mar 27 12:39:05 2023: ncks -d x,57,59 -d y,102,104 2001…
    NCO: netCDF Operators version 4.9.3 (Homepage =

mesh_mask = data_path + ‘coordinates_small.nc’
ds = xr.open_dataset(mesh_mask)
with xr.set_options(display_max_rows=999):
print(ds)

<xarray.Dataset>
Dimensions: (time: 1, y: 3, x: 3)
Dimensions without coordinates: time, y, x
Data variables:
glamf (time, y, x) float64 …
glamt (time, y, x) float64 …
glamu (time, y, x) float64 …
glamv (time, y, x) float64 …
gphif (time, y, x) float64 …
gphit (time, y, x) float64 …
gphiu (time, y, x) float64 …
gphiv (time, y, x) float64 …
Attributes:
history: Mon Mar 27 12:39:38 2023: ncks -d x,57,59 -d y,102,104 coordina…
NCO: netCDF Operators version 4.9.3 (Homepage =

mesh_mask2 = data_path + ‘mesh_mask_co6_small.nc’
ds = xr.open_dataset(mesh_mask2)
with xr.set_options(display_max_rows=999):
print(ds)

<xarray.Dataset>
Dimensions: (t: 1, y: 3, x: 3, z: 51)
Dimensions without coordinates: t, y, x, z
Data variables:
e1f (t, y, x) float64 …
e1t (t, y, x) float64 …
e1u (t, y, x) float64 …
e1v (t, y, x) float64 …
e2f (t, y, x) float64 …
e2t (t, y, x) float64 …
e2u (t, y, x) float64 …
e2v (t, y, x) float64 …
e3t_0 (t, z, y, x) float64 …
e3u_0 (t, z, y, x) float64 …
e3v_0 (t, z, y, x) float64 …
e3w_0 (t, z, y, x) float64 …
esigt (t, z) float64 …
esigw (t, z) float64 …
ff (t, y, x) float64 …
fmask (t, z, y, x) int8 …
fmaskutil (t, y, x) int8 …
gdept_0 (t, z, y, x) float32 …
gdept_1d (t, z) float64 …
gdepw_0 (t, z, y, x) float32 …
gdepw_1d (t, z) float64 …
glamf (t, y, x) float32 …
glamt (t, y, x) float32 …
glamu (t, y, x) float32 …
glamv (t, y, x) float32 …
gphif (t, y, x) float32 …
gphit (t, y, x) float32 …
gphiu (t, y, x) float32 …
gphiv (t, y, x) float32 …
gsi3w (t, z) float64 …
gsigt (t, z) float64 …
gsigw (t, z) float64 …
hbatf (t, y, x) float64 …
hbatt (t, y, x) float64 …
hbatu (t, y, x) float64 …
hbatv (t, y, x) float64 …
isfdraft (t, y, x) float32 …
mbathy (t, y, x) int16 …
misf (t, y, x) int16 …
nav_lat (y, x) float32 …
nav_lev (z) float32 …
nav_lon (y, x) float32 …
rx1 (t, y, x) float64 …
time_counter (t) float64 …
tmask (t, z, y, x) int8 …
tmaskutil (t, y, x) int8 …
umask (t, z, y, x) int8 …
umaskutil (t, y, x) int8 …
vmask (t, z, y, x) int8 …
vmaskutil (t, y, x) int8 …
Attributes:
file_name: mesh_mask.nc
TimeStamp: 08/06/2017 09:17:02 +0000
history: Fri Mar 24 17:08:34 2023: ncks -d x,121,123 -d y,147,149 mesh…
NCO: netCDF Operators version 4.9.3 (Homepage =

2

The time invariant position of cell interfaces is given by the gdepw_0 array.
If you want the time varying position (e.g. dependent on the sea level), i’m afraid you have to compute it offline from gdepw_0 and ssh such as gdepw(i,j,k,t) = gdepw_0(i,j,k) * (1 + ssh(i,j,t)/hbatt(i,j)). It will be referenced to the surface, therefore, if you want the absolute position, substract ssh.

3

Many thanks, unfortunately I am still confused with the meaning of these variables.

  • In sigma coordinates, what is the meaning of gdepw_0 (t, z, y, x)? The NEMO documentation only refers to it in the context of z coordinates.

  • My output files (listed above) don’t include ssh, but the cell thickness in m, e3w, is present. So would this be a correct way to determine the position of the top surface for each layer in relation to the mean sea level?

depth_{MSL}(t,k,i,j) = hbatt(i,j) - \sum_{l=NLevels}^{k} (e3w[t,l,i,j])

4

gdepw_0 refers to the depths of cell interfaces in s-coordinates as in z-coordinates. The only difference being that in z-coordinates, gdepw_0(z,y,x) = gdepw_1d(z) everywhere except at the ocean floor.
If you really need the depth of interfaces (the top of layer), you need e3t arrays, the distance between two consecutive interfaces. I guess you don’t have it.
If your concern is to get the depth of the T-point in the layer, then you can use e3w array since this is the distance between two consecutive points in consecutive layers.
Hence your formula is wrong, since e3t should be used not e3w

5

Thank you jchanut. This is confusing, so e3t is the distance between W points (top surface of layer) and e3w is the distance between T points in the middle of the layers?!

6

Yes it is. Formally speaking, T-points may not be located right “in the middle” of a layer since their locations (depths) are defined according to analytical functions.