National Oceanic and
Atmospheric Administration
United States Department of Commerce

Using sigma and curvilinear model coordinates

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 dealingwith sigma (and curvilinear) coordinates. The analytical methods involveregridding the sigma coordinates to a grid based on a true depth axis.The results of this method can then be combined in calculations with otherdepth-gridded fields (The same procedures apply to pressure-gridded fields.)The graphical methods are "faithful" in the sense that they donot involve regridding the data to another grid.

In this FAQ we will show two examples:

I) Open a dataset containing model data that uses a sigma-coordinate "stretching variable". From the datset define variables to use for two techniques, a purely graphical one, and a second which regrids the data to a fixed Z axis.

Start by defining the "depth" functionas the product of a bathymetry variable and the sigma-z stretching variable.

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

  1. purely graphical method using 3-argumentSHADE and FILL commands
  2. Z axis replacement technique using ZAXREPLACE

II) begin by:creating an artificial sigma coordinate dataset.

define the "depth" functionas the vertical integral of layer thickness and produce some referenceplots 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-argumentSHADE and FILL commands
  2. Z axis replacement technique using ZAXREPLACE
  3. older technique based on the @WEQ transformation (pre-v4.9).


Open a sigma coordinate data set:

Your model data might look like this:

yes? use "http://address/thredds/dodsC/directory/roms_forec/ROMS_Ocean_Model_best.ncd"
yes? show dat
currently SET data sets:
1> http://address/thredds/dodsC/directory/roms_forec/ROMS_Ocean_Model_best.ncd (default)
name title I J K L M N
HC S-coordinate parameter, critica ... ... ... ... ... ...
CS_R S-coordinate stretching curves ... ... 1:30 ... ... ...
H bathymetry at RHO-points 1:294 1:194 ... ... ... ...
SALT salinity 1:294 1:194 1:30 1:6626 ... ...
TIME_RUN run times for coordinate = time ... ... ... 1:6626 ... ...

Have a look at the attributes of the salinity variable

yes? show attr salt
attributes for dataset: http://address/thredds/dodsC/directory/roms_forec/ROMS_Ocean_Model_best.ncd salt.units = 1e-3 salt.long_name = salinity salt.time = ocean_time salt.coordinates = time_run time s_rho lat_rho lon_rho salt.field = salinity, scalar, series salt._FillValue = NaN salt.standard_name = sea_water_salinity

The coordinates attribute tells us that the Z axis is to be given by the s_rho variable.


Define the "depth" function, and plot with a purely graphical method

The variable H is the bathymetry on the X-Y grid of our variable. We can define depth using the bathymetry and the sigma-coordinate stretching variable. Depth is the depth in meters, at the X-Y grid of SALT, and at each unitless z level of the grid of SALT.

Now use the 3 argument SHADE and FILL commandsthat directly plot data fields on "curvilinear coordinates". Todisplay the field we need only create multidimensional fields specifyingthe horizontal and vertical positions associated with each data point tobe plotted.


yes? let depth = cs_r* h

yes? shade/y=100/L=1 salt[l=1], lon_rho, depth

X-Z sigma plot, 3-arg plot command


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

Define a Z axis and use ZAXREPLACE to regrid the data to this axis. Note: this demo is graphical but the procedure is valid for analysis as well.

This is the z axis that is used in the Levitus Climatology dataset. You can define any Z axis, with depth in meters matching the units of variable h.

yes? define axis/z/units=meters/depth zaxlevitr = {0,10,20,30,50,75,100,150,200,300,400,600,800,1000,1200,1500,2000,3000,4000,5000}

Variable cs_r has negative values between -1 and 0, and our z axis refers to depth as positive down, so we'll change our definition of the depth variable, making its values positive depths as well.

yes? list cs_r VARIABLE : S-coordinate stretching curves at RHO-points DATA SET : Regional Ocean Modeling System (ROMS): FILENAME : ROMS_Ocean_Model_best.ncd FILEPATH : http://address/thredds/dodsC/directory/roms_forec/ SUBSET : 30 points (Z) -0.983 / 1: -0.9802 -0.95 / 2: -0.9260 -0.917 / 3: -0.8593 -0.883 / 4: -0.7862 -0.85 / 5: -0.7109 -0.817 / 6: -0.6366...

yes? let depth = -1* cs_r* h
yes? let/title="`salt,return=title`"/units="`salt,return=units`" salt_on_depth = zaxreplace (salt, depth, z[gz=zaxlevitr])
yes? fill/Y=100/L=1 salt_on_depth ! The image below

yes? shade/z=10/L=1 salt_on_depth, lon_rho, lat_rho ! a single depth in meters, curvilinear in longitude/latitude

X-Z sigma plot, zaxreplace


OR, Create an artificial sigma coordinate data set:

(Only needed for demonstration purposes)

Create variables to define a bottom bathymetry on a grid and a verticalsigma coordinate system of 10 layers. We will invent a (non-physical) flowfield 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 the Ferret script 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 thelayer.

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


Solutions:

  1. purely graphical method

  2. As of Ferret v4.9, you can use the 3 argument SHADE and FILL commandsto directly plot data fields on "curvilinear coordinates". Todisplay the field we need only create multidimensional fields specifyingthe horizontal and vertical positions associated with each data point tobe 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 releasenotes for more information on the three-argument SHADE command.)

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

  4. See sigma_coordinate_demo.jnlfor details of the variable definitions.

    Note: this demo is graphical but the procedure is valid for analysis, 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 becausethe 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

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

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 previoussection, which is based on the ZAXREPLACE() function. See the sigma_coordinate_demo_weq.jnlfor 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 analysis,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 becausethe deepest layer of flow is considered to be 1/2 grid box above the bottom.

    sigma-xy plots

    Time series of flow: Visualize flow where the "view" of the data requires a rangeof depths

    sig_timeseries

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

    This procedure (borrowed from depth_to_density_demo.jnl) requires usto "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 willborrow 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 theX 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]

    Again, Cross-channel section of the flow; and Along-channel sections of flow:

    sig_xz

    sig_yz