**Quote:**

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"):

**Quote:**

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:**

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:**

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:**

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?

**Quote:**

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).