Return to Ferret FAQ

Using sigma and curvilinear model coordinates

Question:

How do I handle sigma coordinate output in Ferret?

How do I handle models based on curvilinear coordinates?

Example:

Explanation:

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:

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

Create an artificial sigma coordinate data set:

(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.

Define the "depth" function

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
    

--Bathymetry and DEPTH(layer)

Solutions:

purely graphical method

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

Z axis replacement technique (using the ZAXREPLACE() function)

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]
    

--Plan views of flow

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

--Time series of flow

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

pre-V4.9 technique (based on @WEQ transformation)

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.

Visualize flow where z is a single fixed depth

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]
    

--Plan views of flow

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

--Time series of flow

Visualize flow where the "view" of the data requires a range of depths

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