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])