Properties on a level surface

[cwinant]

Using Roms output I am trying to make maps of properties (in this case
rho) on level surfaces. I have encountered the following problem. The
basin is a rectangle(10000 km x 2500 km) of depth 4000 m surrounded by
a continental shelf of nominal width 150 km and depth 200m. The depth
varies as a tanh between the shelf and deep values, the width of the
transition region is about 100 km.

The case starts from an initial condition where rho varies in the
vertical as (-tanh alpha(z+hbar)) where hbar is the thermocline depth,
and alpha is a moderate slope. I can take the rho values given in
the output at the initial time, and use set_depth.m and stretching.m
to compute the depth corresponding to each vertical bin (there are
40). A blow up of the resulting profiles at 5 points is in the
attached figure. Each line corresponds to a different value of x_rho,
ranging from the edge of the basin (blue) to twice the shelf depth
(green).

Even though the density profile should be exactly the same at all
points, the sigma-coordinates result in uneven sampling depths, and as
a consequence, the density profiles cross each other. The surface
corresponding to a density of say, -5 is highly contorted. The density corresponding to each bin and the bin depth are illustrated in the second plot.

Finding properties along level surfaces is such a common thing to do
that I am certain someone has encountered and solved this problem
before. I have not found reference to the subject in this forum, although that is probably because I havent look hard enough. Perhaps someone would share how best to handle the issue.

Clinton Winant