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:

[Output Graphic]

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:

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
h
layer thickness
flow
the flow field

See sigma_coordinate_demo.jnl for details of the variable definitions.

Define the "depth" function

Solutions:

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

  2. 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]
          

    a 2-frame plot comparing flow field    in a single layer with flow at fixed depth--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

    a time series plot    making a similar comparison--Time series of flow

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

    comparing cross-section of flow field in a single layer with flow at fixed depth--Cross-channel sections of flow

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

    comparing longitudinal-section of flow field in a single layer with flow at fixed 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.)

  3. 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]
        

    a 2-frame plot comparing flow field    in a single layer with flow at fixed depth--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

    a time series plot    making a similar comparison--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]
        

    comparing cross-section of flow field in a single layer with flow at fixed depth--Cross-channel sections of flow

    comparing longitudinal-section of flow field in a single layer with flow at fixed depth--Along-channel sections of flow

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


Last modified: Nov 20, 1997