Re: ZAXREPLACE and Sigma to Z depth

[Ansley Manke <ansley@xxxxxxxxxxxxx>]

Hello Steve,
Take another look at the documentation for the ZAXREPLACE
function.  (Ch3 Sec2.3.24. of the on-line
Users' Guide).  I'll quote the
text about the arguments to the function below.  This section
was recently
revised in the Users' Guide to clarify it.   Note that all
the arguments to
ZAXREPLACE  must have Z axes, and the first two arguments share
the same Z axis.  So you don't want the variable "depth" as the
second
argument.   That's just the bathymetry and has no z axis.

You need a variable having the same Z axis as "temp", whose values
give the conversion to the desired Z axis in the third argument.  In your
case this is the depth of each sigma layer.   How do the sigma-depths
zpos = 0, -0.015625, etc. relate to the depths over the range of the output
z axis, Z=0:5000:10?

You're right in general about the kind of interpolation that the function does.
It just needs the translation.

Here's the paragraph about the arguments to ZAXREPLACE:

> The ZAXREPLACE function takes three arguments. The first
> argument, V, is the field of data values, say temperature or salinity.
> This variable is available on what we will refer to as the "source"
> Z-axis -- say in terms of layer number. The second argument, ZVALS,
> contains the values of the desired destination Z axis defined on the
> source Z axis -- for example, it may contain the depth values associated
> with each vertical layer. It should always share the Z axis from the
> first
> argument. The third argument, ZAX, is defined on the destination Z
> axis. 
> Only the Z axis of this variable is relevant -- the values of the variable,
> itself, and its structure in X, Y, and T are ignored. Often "Z[gz=zaxis_name]"
> is used for the third argument.
I hope this helps!

Ansley Manke

Steve Cousins wrote:

I'm trying to do the conversion from Sigma layers to Z-depth for our
Princeton Ocean Model (POM) data that is produced by Rich Signells NetCDF
putcdf() function.

The 15 Sigma layers are fractions of depth from 0.0 to -1.0.  There is a
2D depth variable which has the true depth at each grid point.  My first
goal is to get X,Y plots at specified depths in meters. My second goal is
to get X,Z and Y,Z plots that use depth in meters as opposed to sigma
layers.

Here is the relevant NetCDF header information:

netcdf pom.24day {
dimensions:
        xpos = 151 ;
        ypos = 121 ;
        zpos = 15 ;
        time = UNLIMITED ; // (46 currently)
variables:
        float xpos(xpos) ;
                xpos:long_name = "Columns" ;
                xpos:units = "column" ;
        float ypos(ypos) ;
                ypos:long_name = "Rows" ;
                ypos:units = "row" ;
        float zpos(zpos) ;
                zpos:long_name = "Layer" ;
                zpos:units = "sigma_level" ;
        float time(time) ;
                time:units = "days since 30-MAR-1998 00:00:00" ;
                time:time_origin = "30-MAR-1998 00:00:0" ;
                time:modulo = " " ;
        float depth(ypos, xpos) ;
                depth:valid_range = 0.f, 5000.f ;
                depth:long_name = "Bathymetry" ;
                depth:units = "meters" ;
                depth:missing_value = -99999.f ;
        short temp(time, zpos, ypos, xpos) ;
                temp:valid_range = -5.f, 30.f ;
                temp:long_name = "Temperature" ;
                temp:units = "Celsius" ;
                temp:scale_factor = 0.0005340739f ;
                temp:add_offset = 12.5f ;
                temp:missing_value = -99999.f ;

The sigma layer fractions look like:

 zpos = 0, -0.015625, -0.03125, -0.0625, -0.125, -0.25, -0.375, -0.5, -0.625,
    -0.75, -0.875, -0.9375, -0.96875, -0.984375, -1 ;

I've been looking at the User's Guide and the sigma_coordinate_demo.html
but am not getting it to work.  If I try:

     DEFINE AXIS/Z=0:5000:10/UNIT=meters/DEPTH zdepth
     LET temp_on_depth = ZAXREPLACE(temp,depth, z[gz=zdepth])
     LET temp_at_50 = temp_on_depth[Z=50]
     SHADE/L=1 temp_at_50, 'X','Y'

ferret produces:

**ERROR: illegal limits: Z limits of data and depth fields must match
SHADE/L=1 temp_at_50, 'X','Y'

I think I understand this to mean that since the Z dimension of the
temperature variable is only 15, I can't ask for temp_on_depth[Z=50].

I'm obviously missing something fairly basic.  I would think that an
interpolation would need to be done based on the Sigma Layers and the
Depth variable.  At a given grid point, you could find the temperature 50
meters down by doing something like (linear interpolation):

c    Make sure the total depth is at least 50 meters and calculate
c    the fraction of depth that 50.0 meters is:
     IF (depth(x,y) .gt. 50.0)
        fraction = 0.0 - (50.0 / depth(x,y))
     else
        temp_at_50 = -99999.0
        GOTO 10
     endif

     do k=2,kb  ! kb = number of sigma layers
        if (fraction .ge. zpos(k)) then
c          Sandwich the depth between two layers:
           d_top = ABS(zpos(k-1) * depth(x,y))
           d_bot = ABS(zpos(k) * depth(x,y))
           d_diff = d_bot - d_top
           d_frac = (50.0 - d_top) / d_diff

c          Sandwich the temp between two layers:
           t_top = temp(time,k-1,y,x))
           t_bot = temp(time,k,y,x)
           t_diff = t_bot - t_top
           temp_at_50 = (d_frac * t_diff) + t_top
           GOTO 10
        endif
     enddo
 10  continue

     print*,'Temperature at X,Y,50.0': temp_at_50

I know this isn't bullet-proof (possible divide by zero) but I just
wanted to show what I was thinking.  Is ZAXREPLACE supposed to be able to
do something like this?

My ultimate goal is to get LAS to be able to do this.  Any information,
tips, suggestions are very welcome.

Thanks,

Steve
_____________________________________________________________
 Steve Cousins                 Email: cousins@umit.maine.edu
 Research Associate            Phone: (207) 581-4302
 Ocean Modeling Group
 School of Marine Sciences     208 Libby Hall
 University of Maine           Orono, Maine 04469

--
Ansley Manke  Pacific Marine Environmental Laboratory  Seattle WA  (206)526-6246