diagnosing vertical coordinate from history/grid files

[jpringle]

Hi all --

I have a simple question; I fear perhaps the answer is not so simple. I am involved in a couple of efforts to read in ROMS history and grid files, and either plot them or compute other products from them (e.g. Lagrangian particle paths). I have several questions

1) What is the best way to figure out what version of ROMS made the output?

2) For what versions of the ROMS code is Vtransform and Vstretching specified somewhere in the output, and how? I see them as variables in some history files.

3) Same question as (2), but for average files. I have on hand several average files that seem to be very complete, but do not have the Vtransform flag.

4) If the Vtransform flag is missing from the history and average files, can I assume Vtransform=1?

5) in ROMS 2.2, hc and other vertical grid parameters are attributes in the output; in 3 they are variables. Is this true for all versions of ROMS 3?

I can guess some of these from the code that I have on hand; but it is easy to over generalize from a few examples, and I would prefer to know the "official" answers. File sniffing is easy to do poorly. I have to deal with ROMS 2.2* and 3* output, and can't just ask the folks providing me data to re-run the models with up to date code.

Thanks everyone in advance.
Jamie Pringle
University of New Hampshire

[kate]

1) What is the best way to figure out what version of ROMS made the output?

Maybe not so easy to even know when running it if you are thinking of the 3.x kind of versions. That value is certainly not in the ROMS/Version file.

2) For what versions of the ROMS code is Vtransform and Vstretching specified somewhere in the output, and how? I see them as variables in some history files.

3) Same question as (2), but for average files. I have on hand several average files that seem to be very complete, but do not have the Vtransform flag.

4) If the Vtransform flag is missing from the history and average files, can I assume Vtransform=1?

yes

5) in ROMS 2.2, hc and other vertical grid parameters are attributes in the output; in 3 they are variables. Is this true for all versions of ROMS 3?

I can't remember back to ROMS 2.

I can guess some of these from the code that I have on hand; but it is easy to over generalize from a few examples, and I would prefer to know the "official" answers. File sniffing is easy to do poorly. I have to deal with ROMS 2.2* and 3* output, and can't just ask the folks providing me data to re-run the models with up to date code.

Thanks everyone in advance.
Jamie Pringle
University of New Hampshire

Maybe Hernan can pipe up with the "official" answer. We can ask that ROMS 4 has the 4.x version number in ROMS/Version and that that value gets put into the output (for those of us not getting svn numbers into our output).

[shchepet]

Jamie Pringle: ...I would prefer to know the "official" answers. File sniffing is easy to do poorly...

Actually it looks like you already provided the official answer yourself: it is what you call
"file sniffing".

The whole story of "new" coordinate emerged from the dissatisfaction with the "old" coordinate,
specifically inability to have good behavior of surface boundary layer, excessive water-mass mixing
and few other things which were cured or mitigated. Also to encourage people to have open mind and
propose/develop new ideas about vertical coordinates, depending of someone's specific needs, e.g,
coastal/sediment/estuary people may have different preferences than open-ocean KPP-minded people.

The code itself is designed to be mathematically coordinate free in sense that one can plugin any.

When I proposed it to Hernan, the idea was initially screwed up by keep imposing the restriction of
hc < hmin, which kills the point: the point is that it is the hc < hmin restriction itself which was the
primary drawback of the "old" coordinate. [The hc < hmin condition for new coordinate was
erroneously posted on the relevant ROMS Wiki page
https://www.myroms.org/wiki/index.php?t ... ntable=yes
and was kept there for more than a year until Brian Powell pointed out the mistake in web suite
(not in the code!) back in April 2009.

I also proposed a taxonomy standard to avoid confusion, but that proposal was rejected by
Hernan, and things like Vtransform and Vstretching=1 or 2 appeared instead in the official ROMS.
It is my experience that negotiating convention standards is just a plain waste of time: these
negotiations always fail -- it is just someone is not willing to change his Matlab script, no
matter what. It is as simple as that.

Yes, my preferred way is
(1) Between using netCDF global attributes vs. variable: if a functionality can be achieved
by using an attribute, one should use attribute, not a variable. NetCDF variables require two
stages in creating: first nf_def_var(...), then, after call to nf_enndef(...), put a meaningful
value into that variable; handling an attribute require just a single call. In practice I saw
too many netCDF files where one creates a variable and then never assigns a value to it.
This would never happen to an attribute.

(2) Use meaningful TEXT label as an identifier instead of numbers, for example

Code: Select all

 VertCoordType='SH94' 

is preferred over

Code: Select all

 Vtransform=1 

because we already have people who came late enough to that the only transform they knew is the
"new", and we do not have to explain legacies, i.e., what was first and what was second, and why
the numbers are assigned in certain order. Also to prevent a spurious illusion of the existence of
some kind of "default" transform.

(3) Keep the specification in netCDF file at the absolute mathematical minimum of the information,
i.e., for the "new" coordinate is is sufficient to store only 4 attributes: say

Code: Select all

VertCoordType='new'
hc=400
Cs_w=-1.,-0.9654, .......-0.0032, 0.  <-- an array of values
Cs_r= ......                             <-- array of values

Everything else is mathematically redundant.
Parameters theta and theta_b are useful only as a descriptive information for user, but not to be
used by the software processing output (ROMS plotting package; Matlab or whatever) -- all software
must read Cs_w,r and use the provided values rather than attempting to generate them from theta and
theta_b parameters. [It is expected that one can modify the functions generating Cs_w,r, so theta
and theta_b may end up having a different meaning, but this should be transparent to the post-
processing software. For this reason it is actually non necessary to have something like Vstretching=1
or 2 or 3, whatever.]
Storing variables or attributes for s_w and s_r (primary or non-stretched sigma-coordinate) is
unnecessary because these are always just equally spaced numbers, s_w(k)=(k-N)/N and
s_r(k)=(k-N+1/2)/N, so there is no point of storing them.

Note that all of the three principles above (not just some of them) were rejected as the outcome
of negotiations for a convention standard. In reality it is just an inconvenience for Matlab users, but
overall is not a big deal. Also helps to enforce some sort of code loyalty -- if someone start using
a branch of ROMS, he tends to stick with it.

[arango]

If the variables Vtransform and Vstretching are not in the output NetCDF files, it implies that you are using a version of ROMS that it is old. In this case, the vertical transformation equations are that for Vtransform=1 and Vstretching=1. This is the original transformation equations.

We think that CF compliance is extremely important. This will allow us to plot ROMS results correctly with several generic NetCDF plotting packages. Please see the WikiROMS for detailed information. We always try to keep the correct information in WikiROMS and appreciate input from users. As in any kind of web publications, sometimes there are typos in the provided information. We always correct them when reported.

[shchepet]

I do not object adhering with CF conventions in principle -- there is nothing wrong with them, but
I do not see how CF compliance changes this discussion, in sense that that I do not see what it
has to do with the way how the information about the specific type of vertical coordinate is stored
in a netCDF file, nor it has to say anything about preference of variable vs. attribute, nor am I
aware about any generic plotting package (other than the several ones developed within the ROMS
community) which can handle ROMS output, especially if we decide to introduce a new vertical
coordinate within a month from now.

What I did notice just now is that Eq. (6) on Wiki page corresponds to a long-obsolete version in
its refinement-toward-bottom part and this should be corrected/updated. I recall that it was coded
correctly in ROMS v.3.4.xxx back in 2009, but, paradoxically, I see the older version in v.3.4.535.

The relevant piece of code should be

Code: Select all

      function CSF (sc, theta_s,theta_b)     ! limits of CSF,csrf for
      implicit none                          ! theta_s, theta_b --> 0
      real*8 CSF, sc, theta_s,theta_b,csrf   ! match that under "else"
                                             ! logical branches.
      if (theta_s.gt.0.D0) then
        csrf=(1.D0-cosh(theta_s*sc))/(cosh(theta_s)-1.D0)
      else
        csrf=-sc**2
      endif
      if (theta_b.gt.0.D0) then
        CSF=(exp(theta_b*csrf)-1.D0)/(1.D0-exp(-theta_b))
      else
        CSF=csrf
      endif
      return
      end

This corresponds exactly to Eq. (2.4) from the 2009 Correction Note published in JCP.

The difference from the one in ROMS Wiki is only within the bottom part -- the two are
equivalent if theta_b=0.

For basin-scale typical parameter values are theta=8.0...10.0, theta_b=2.0, hc=400.0, which is
a bit more aggressive than in the ~6.5, 0, and 300, respectively found in the 2009 Note, ...but life
is moving. Also note that the larger setting of theta requires some theta_b>0 even if you do not
care about bottom boundary layer, because otherwise the bottom grid boxes are too coarse.
The above choice of 10, 2, 400 allows to match quite closely the stretching function from Fig. 3
of Griffies et al., Ocean Science, 1, 45–79, 2005,
http://www.gfdl.noaa.gov/bibliography/r ... mg0501.pdf
resulting in a nearly-uniform vertical resolution in the upper 200m.

[jpringle]

Kate, Sasha & Hernan--

Thanks for your replies, they have been very helpful.

Cheers,
Jamie

[kate]

What I did notice just now is that Eq. (6) on Wiki page corresponds to a long-obsolete version in
its refinement-toward-bottom part and this should be corrected/updated. I recall that it was coded
correctly in ROMS v.3.4.xxx back in 2009, but, paradoxically, I see the older version in v.3.4.535.
...

This corresponds exactly to Eq. (2.4) from the 2009 Correction Note published in JCP.

The difference from the one in ROMS Wiki is only within the bottom part -- the two are
equivalent if theta_b=0.

So the distributed code had one version and now has a different, older version? It changed without warning? And now we've built all these tools with the old version? Gotta love it.

Hernan - what are your future plans with regard to this issue?

[shchepet]

So the distributed code had one version and now has a different, older version? It changed
without warning? And now we've built all these tools with the old version? Gotta love it.

Kate, I understand your frustration, which is, obviously an outcome of a communication
screw up on our side. But to say that it changed without warning is a bit overstatement:
Did you receive the following e-mail (check you mailbox around Wednesday,
May 20, 2009 10:52 PM; the e-mail an attachment which is version of "set_scoord.F"):

Date: Wednesday, May 20, 2009 10:52 PM
Re: vertical stretching with bottom refinement
From: "Alexander Shchepetkin" <old_galaxy@yahoo.com>
To: "Hernan G. Arango" <arango@marine.rutgers.edu>, "John CWarner"
<jcwarner@usgs.gov>, "Brian Powell" <powellb@hawaii.edu>,
"Kate Hedstrom" <kate@arsc.edu>, "Alexander Shchepetkin" <old_galaxy@yahoo.com>
Message contains attachments
1 File (5KB)

I think I found a good functional form which allows bottom refinement in a well-controllable
manner, please check it out, say code it in matlab and see how it behaves. See the attached
file set_scoord.F. This is to be used in conjunction with so-called new transform, that it

z^(0)= h * [hc*s + h*C(s) ]/[hc+h]

and the near-surface property of C(s) is that its derivative --> as as s-->0, which is nothing
new at this moment. What is new is how it handles refinement near the bottom. It is no
longer "blending" of two functions, but is rather works as the second stretching of already
stretched transform: note that in the code below "csrf" computed in the upper part is used
as argument in the lower.


Also note that the transfor is a continuous function with respect to both theta_s and
theta_b when both --> 0: mathematical limits of the upper branches of "if" statements
in the code below are equal to the expressions in the lower branches.

Therefore, the behavior with respect to vatiation of theta_s , while theta_b=0 is as it was
before, while the bottom part of the curve is smoothly deformed when theta_b departs
smoothly from zero.

Just program it with Matlab, plot it and, especially its derivative with respect to s (this is
what setts grid-box heights) and play with its parameters theta_s and theta_b it a little bit.

The range of meaningful values for both theta_s and theta_b is from 0 to 6.5 or so,
although, as usual, one should pay attention to extreme rx1 values near the bottom.



function CSF (sc, theta_s,theta_b)
implicit none
real*8 CSF, sc, theta_s,theta_b,csrf
if (theta_s.gt.0.D0) then
csrf=(1.D0-cosh(theta_s*sc))/(cosh(theta_s)-1.D0)
else
csrf=-sc**2 !<-- note continuity between the two if-branches
endif
if (theta_b.gt.0.D0) then
CSF=(exp(theta_b*(csrf+1.D0))-1.D0)/(exp(theta_b)-1.D0) -1.D0
else
CSF=csrf
endif
return
end

The upper part was settled few years ago, and this is well documented in Hernan's
own power-point presentations.

The bottom part -- the formula which uses function of function,

Code: Select all

S(s)=upper stretching function, e.g. [1-cosh(theta*s)]/[cosh(theta)-1]
C[S(s)]= bottom refinement function, e.g., [exp(theta_b*S)-1]/[1-exp(-theta_b)] 

[note that S in exp(theta_b*S) in the second formula in the uppercase S.] -- is from
early spring 2009, an by early summer we have acquired sufficient experience to
build confidence that it works better that the "alpha-beta blending" of surface
and bottom stretching which was there before. [Also note that the line of code

Code: Select all

       CSF=(exp(theta_b*(csrf+1.D0))-1.D0)/(exp(theta_b)-1.D0) -1.D0

appearing in the e-mail above and in the code attached therein is
mathematically equivalent to

Code: Select all

        CSF=(exp(theta_b*csrf)-1.D0)/(1.D0-exp(-theta_b))

as it appears in the today's code.]

It was never my intend to make too big noise out of it, not too keep it under
wraps. It was briefly mentioned among features of ROMS in my talk June 15 2009
ONR meeting in Chicago, and ended up in 2009 JCP paper. I checked dates: the
JCP paper was finalized on June 19, 2009, that is within a weak after Chicago.
There were multiple communications between me and Hernan during that
period, but this particular point somewhat did not make it to the "official"
svn branch.

I have checked my e-mail and there were two series of communications related
to this matter, one in throughout December, 2008 initiated by Rocky Geyer and
Rich Signell; the other one in April-May, 2009 initiated by Brian Powell, and
involving, Hernan, myself, John Warner, and some (but not all) of these e-mails
was copied to you. However, much of the emphasis in these e-mails was mainly
about taxonomy, and how to label different options withing netCDF files.
Also I understand that back in December 2008, Rich and Rocky were going to
some sort of CSTMS steering committee meeting focusing on bbl and wanted
some input.

Did you receive this e-mail?

Date: Sunday, May 31, 2009 10:59 AM
Re: vertical stretching with bottom refinement
From: "Alexander Shchepetkin" <old_galaxy@yahoo.com>
To: "Brian Powell" <powellb@hawaii.edu>
Cc: "Hernan G. Arango" <arango@marine.rutgers.edu>, "John CWarner" <jcwarner@usgs.gov>,
"Kate Hedstrom" <kate@arsc.edu>, "Alexander Shchepetkin" <old_galaxy@yahoo.com>

--- On Sat, 5/30/09, Brian Powell <powellb@hawaii.edu> wrote:

>
> Hernan, I have added this as Vstretching=4 into my set_scoord.F, I recommend
> we incorporate it into the public version.
>


Brian and Hernan,

You should avoid things like introducing

Vstretching = some number


to avoid future confusion because I expect new functions will appear at some point;
and we should come up with a smart taxonomy befor its too late.

I do not consider the latest version of my C(s) as "new", since if theta_b=0, it
is identically equivalent to the previous one (which you call Vstretching = 2), and
too date all usage of that form with non-zero theta_b was basically limited and is
already in process of being phased out (there is a handful of solutions which we
are recomuting any way).


To summarize, as of right now we have only two versions of vertical coordinate,

z^(0) = hc*s + [h(x,y)-hc]*C(s)
and
z^(0) = h(x,y) * [hc*s + h(x,y)*C(s)] / [ [h(x,y)+hc]

which I call VertCoordType = Legacy and VertCoordType = NEW
respectively. They can be used in combination with any stretching function C(s).

Both of the types above yield straight POM-style sigma, if hc=0.


As for C(s) we have three different formulas [Legacy SH94; my newest one; and
Rocky Geyer] so in theory we already have a total of six permutations.


[ For some time in the past (early 2005) I was using VertCoordType = NEW in combination
with Legacy SH94 stretching, and almost universally setting theta_b=0, then replaced
sinh(theta_s*s) with 1-cosh(theta_s*s), but the rest of the formula was still deactivated
by setting theta_b=0: I happened to stay away from bottom boundary layer business until
very recently....

...my 1-cosh(theta_s*s) formula does not yield uniformly spaced C(s)=s
vertical grid, if theta_s --> 0, instead it asumptotes to s^2. I do not consider it as an issue
in open-ocean problem, because we always have upper boundary layer, we need to have
a region of quasi-uniform vertical grid at the top. My rationale that this should be entirely
controlled by setting "hc" with as little as possible interference from other parameters [and it
more or less does so, since (1) hc is liberated from hc<hmin restriction, and (2) C(s) has
zero derivative at surface, hence theta_s does not affect the hight of the uppermost grid
box). In contrast, the role of theta_s is reserved to set resolution in the main thermocline
area.

Formally speaking the uniformly-stretched sigma coordinate (same as setting
theta_b=theta_s=hc=0 in SH94 formula) is achieved by setting infinitely large hc. However,
setting large "hc" also makes this formula be insensitive to theta_b: in fact, it cannot achieve
bottom-only refinement. I do not consider it as a drawback for any open or coastal ocean
problem, but I envision river people who will ask one day, that we do not want surface
boundary layer because there is very little stratification and almost no wind, but the river
flows any way. So they need bottom-only refinement, and one will be tempted to replace

if (theta_s.gt.0.D0) then
csrf=(1.D0-cosh(theta_s*sc))/(cosh(theta_s)-1.D0)
else
csrf=-sc**2
endif

with

if (theta_s.gt.0.D0) then
csrf = (exp( - theta_s*sc)-1.D0) / (1.0 - exp(theta_s) )
else
csrf=sc
endif

while the bottom if-block of my CSF function remains unchanged. (Note that it is kind
of bottom-surface symmetric in terms of function types.)

Now the overall CSF has smooth transition toward uniform C(s)=s when theta_s=0, and
one can achieve bottom-only refinement by setting theta_s=0, theta_b>0, and to avoid
over-resolution near bottom in the shallowest places sets a reasonable value for "hc"
LARGER than the minimum depth.

Again, its near-surface behavior is not the best idea for open ocean because it somewhat
lost its property of flattening levels in areas where here, but for a river problem -- why not.

Finally, why not sinh(theta_s*sc)/sinh(theta_s) kind of function?

SURFACE: unlike (exp( - theta_s*sc)-1.D0) / (1.0 - exp(theta_s) ), sinh
yields nearly uniform resolution at sc --> 0. Although this may seem like an
advantage, this is redundant with the functionality of "hc". Recall, that in SH94
formula "hc" takes over control of grid spacing when s-->0 because it is
presumed that |C(s)| << |s| when s --> 0. The problem is that in practice
SH94 must restrict hc by minimum depth, and that makes it less useful.
So the use of "sinh" has some merit to cope with that. Once you get rid of
hc<hmin limitation, there is no point of having "sinh".

BOTTOM: Actually, at first I considered using sinh-type in the bottom part of CSF,
but it turns out that its vanishing second derivative is not a good idea there, because
it must "over bend" surface-refined csrf which has non-vanishing second derivative
as it approaches bottom. As the result, there always "kink" in resolution, and one of few
bottom-most grid boxes end up being larger than above.

Recall that SH94 never had bottom-bound refinement, but instead it actually yields
the placements of relatively finer resolution somewhere in the bottom-half of s-space,
when one sets theta_b>0, but NOT near the bottom. ]





For the purpose of post-processing there is no need to introduce any identifier for
stretching functions C(s): they must be stored in netCDF files and all post-processing
software must read C(s) curves rather than attempt to reconstruct them from parameters.
I strongly encourage (even force) people around here and there to modify their
matlab/python scripts to adhere with this principle (as opposite to reconstruction C(s) from
theta_s, theta_b using specific formula). Then no changes are needed in the software
if C(s) formula is modified.


Thus, the mathematically minimal information to characterize all what we have thus far is

1. VertCoordType = Lgacy/NEW
2. hc =
3. Cs_w =
4. Cs_r =

I no longer store arrays "sc_r" and "sc_w" in netCDF files -- there are trivial to reproduce.

theta_s and theta_b are still stored in netCDF, but no software relies on them; just for human
convenience.

ALL OF THE ABOVE, including hc,Cs_w, Cs_r are stored as netCDF global attributes
(They are NOT VARIABLES).

[kate]

I'm sorry I let my snarkiness and my senility show. I do have that email in my archive.

We still need to generate the Cs arrays from theta_s, etc. in our tools if we are to build the ROMS boundary/initial conditions before ROMS is run. One could suggest that ROMS would be the one to read the Cs arrays rather than compute them - or compute them in the ana_xxx routines when nothing is read.