Sigma coordinates

["Donald S. Dunbar" <dsd@xxxxxxx>]

I am working with 4D data sets (X,Y,Z,T) where the vertical axis Z is
defined as Z = (Z'+E)/(H+E); -E<=Z<=H; E,H = elevation,depth.  Z' is the
depth below LLW and Z: 0<=Z<=1 is the sigma coordinate.  I want to
convert data from Z to Z'.  To date I have been doing the following (at
10 m, for example) for a 4D variable Q and a vector of SIGMA values:

let ZP = (H+E)*SIGMA - E   ..... define the depth below LLW
let MASK = ZP[Z=@WEQ:10]   ..... find where ZP = 10 m
let TMP = Q*MASK           ..... A matrix of weighted Q values
let Q_10 = TMP[Z=@SUM]     ..... Q at 10 m

Then I repeat the procedure to get Q_20, Q_30, etc.  All of this is in a
GO file generated by a PERL script (so it is automated somewhat). I then
would like to combine all of the resulting horizontal slabs back into a
4D array which will have a fixed vertical spacing (10 m in this
example).

My question is: Is all of this necessary, and isn't there an easier way?
Sigma coordinate systems are common enough that I would think some nice
transformation would be built in.  I don't see any reference in my
manual, however.

If there is no other way, what is the easiest way to stack the slabs
into a 4D array?  Say I have a new grid GQ given by:

DEFINE GRID/LIKE=Q/Z=10:100:10 GQ

How do I define a variable on GQ which consists of the slabs Q_10, etc.
determined above?

Thanks for any pointers.
-- 
------------------------------------------------------
D.S. Dunbar (dsd@oar.com)  Ocean Applied Research Ltd.
1864 Duchess Ave.  West Vancouver BC  V7V 1R1  CANADA
------------------------------------------------------