How do I handle sigma coordinate output in Ferret? How do I handle models based on curvilinear coordinates?

Ferret provides both graphical methods and analytical methods for deailing with sigma (and curvilinear) coordinates. The analytical methods involve regridding the sigma coordinates to a grid based on a true depth axis. The results of this method can then be combined in calculations with other depth-gridded fields (The same procedures apply to pressure-gridded fields.) The graphical methods are "faithful" in the sense that they do not involve regridding the data to another grid.

In this FAQ we will begin by:

- creating an artificial sigma coordinate data set. Clearly in your application you would not do this -- you would be working with your own model outputs.
- defining the "depth" function as the vertical integral of layer thickness and produce some reference plots from it. (In your application you may already have a depth variable.)

We will then demonstrate three methods for working with the data:

- purely graphical method using 3-argument SHADE and FILL commands (Ferret v4.9).
- Z axis replacement technique using ZAXREPLACE (simple to use -- Ferret v4.9).
- older technique based on the @WEQ transformation ( pre-v4.9).

(Only needed for demonstration purposes)

Create variables to define a bottom bathymetry on a grid and a vertical sigma coordinate system of 10 layers. We will invent a (non-physical) flow field in a channel.

The variables we create will be

`bathymetry`- the bottom bathymetry for the model

- a channel with a rise along the axis of the channel
`h`- layer thickness

- varies in X, Y, Z, and T in this example
`flow`- the flow field
- faster near surface
- faster at mid-channel than near the edges
- speeds up over the rise
- speed increases with sinusoidal variation in time

See sigma_coordinate_demo.jnl for details of the variable definitions.

The depth of each (x,y,k) grid point is computed by integrating `h`.
We subtract `h/2` because we want the depth of the midpoint of the
layer.

let depth = h[k=@rsum]-h/2 set variable/title="DEPTH function"/unit=meters depth

**Solutions:**

As of Ferret v4.9, you can use the 3 argument SHADE and FILL commands to directly plot data fields on "curvilinear coordinates". To display the field we need only create multidemensional fields specifying the horizontal and vertical positions associated with each data point to be plotted.

The following commands generate the graphic at the top of this page.

! regrid 'Y' to the data grid let ygg = y[g=gg] set variable/title="Y"/unit=kilometers ygg ! use the three argument 'shade' command shade flow[x=0,l=1], ygg, depth[x=0,i=1]

(See the v4.9 release notes for more information on the three-argument SHADE command.)

See sigma_coordinate_demo.jnl for details of the variable definitions.

Note: this demo is graphical but the procedure is valid for analyses, as well

Define flow field on a surface of constant 50 m depth

DEFINE AXIS/Z=0:180:2/UNIT=meters/DEPTH zdepth LET flow_on_depth = ZAXREPLACE(flow,depth, z[gz=zdepth]) LET flow_at_50 = flow_on_depth[Z=50]

Right hand frame: `SHADE/L=1 flow_at_50`

Note: the missing points between the flow and the black mask occur because the deepest layer of flow is considered to be 1/2 grid box above the bottom

`PLOT/X=0/Y=0/L=1:20 flow[k=8], flow_at_50`

**--Cross-channel
sections of flow**

Right hand frame: `SHADE/Z=0:180/Y=0/L=1 flow_on_depth`

**--Along-channel
sections of flow**

Right hand frame: `SHADE/Z=0:180/X=0/L=1 flow_on_depth`

(Compare this last plot with the graphic at the top of this page.)

This presentation will solve precisely the same problems as the previous section, which is based on the ZAXREPLACE() function. See sigma_coordinate_demo_weq.jnl for details of the variable definitions.

Note: this demo is graphical but the procedure is valid for analyses, as well

Define flow field on a surface of constant 50 m depth

let kernel = depth[z=@weq:50] * flow let flow_at_50 = kernel[z=@sum]

Note: the missing points between the flow and the black mask occur because the deepest layer of flow is considered to be 1/2 grid box above the bottom

Note: the procedure is valid for analyses as well as graphics

This procedure (borrowed from depth_to_density_demo.jnl) requires us to "borrow" an axis orientation to serve as the new depth axis. We can only access a single point location along whatever axis we borrow -- though that point can be any point of the axis. In this example we will borrow the TIME axis -- thus we can produce plots only at L=1, L=2, etc.. If we desired to produce a time series plot we would have to borrow the X or the Y axis, instead.

It takes just 5 commands to transform layers to depth in this manner:

! in a new grid replace the (borrowed) time axis with the desired depth axis define axis/t=0:180:2/unit=meters tdepth define grid/like=gg/t=tdepth ggdepth ! define a new variable, r0, with a value of zero wherever depth equals its ! coordinate on the tdepth axis (r0 is a 4-dimensional variable with ! depth in the T axis slot) ! Note that "L=1" here and in the definition of kflow, below, determine ! the fixed value on the TIME axis at which this calculation takes place let r0 = depth[l=1] - t[g=ggdepth] ! define a new variable, kflow, which, when summed along the Z axis, will give ! the (single) value of flow at the location where depth equals its own ! coordinate on the tdepth axis let kflow = r0[z=@weq:0] * flow[l=1] ! sum the variable (integrate) along the Z axis. Since the Z axis reduces ! to a point in this operation the result is 3D - X,Y and depth. let flow_on_depth = kflow[z=@sum]

**--Cross-channel
sections of flow**

**--Along-channel
sections of flow**

(Compare this last plot with the graphic at the top of this page.)

Last modified: Nov 20, 1997