Frequently Asked Questions (FAQ)
Questions and Answers from correspondence between POM users.
Please check the list before submitting inquiries to the entire
pom-group. File updated 03-31-03 T.E.
TOPICS:
*** ADVECTION
Q:
> I am attempting to simulate a shallow estuary (average depth
> just over 1m), with salinities varying from nearly fresh to
> hypersaline over the course of each year. Fairly strong salinity
> gradients exist at certain times of the year.
>
> The difficulty is that after about six days, salinity calculations
> start to overshoot: salinity in isolated cells becomes slightly
> negative, then shoots up to 50ppt+.
A:
I don't think you are necessarily doing anything wrong.
Whenever you have strong scalar gradients and use the
CENTRAL advection scheme, you can get short wavelength
overshoots.
For my river plume simulations, I also get
some large salinities in the vicinity of the river mouth.
An interesting situation arises when the river meets an
fairly steady alongshore current, and the CENTRAL
differencing scheme (via the overshoot) generates dense
water that then sinks and flows down the shelf as a
completely fictitious density current.
The only way I know to get rid of this is to use an
advection scheme that doesn't overshoot. Upwind is one,
but it strongly damps. Higher order schemes such as the
Smolarkiewicz schemes or TVD flux limiters are a better option,
but of course more complex (and CPU time-consuming).
--
Rich Signell | rsignell@crusty.er.usgs.gov
U.S. Geological Survey | (508) 457-2229 | FAX (508) 457-2310
384 Woods Hole Road | http://crusty.er.usgs.gov/rsignell.html
Woods Hole, MA 02543-1598 |
--- See the paper by Pietrzak (Mon. Weather Rev., 126, 812-830, 1998)
about testing different advection schemes with POM.
--- Another advection scheme by Lin et al. (1994) has been used in POM
by several users (S. Hakkinen, T. Ezer) it is now available from the
POM ftp site in directory contrib-code, it replaces ADVT subroutine
in POM.
--- Smolarkiewicz (1984) iterative scheme is also available in contrib-code
( ADVT.Sml.sub), read on on questions related to it.
>
> QUESTIONS
> (1) It is necessary to calculate the first order upwind
> advection before calculating anti- diffusion. But, in ADVT
> .Sml.sub, if NItera=1 then the scheme is just the first
> order upwind, if NItera=2 or 3, it calculates only
> anti-diffusion. This is not Smolarkiewicz scheme.
I looked at the Smol. (1984) paper, and the scheme does follow it closely.
For 1 iteration the scheme should be a simple upwind. For NItera>1 XMassFlux
etc. are recalculated each iteration & used in do 140.
>
> (2)In the original scheme, there is no 'switch option' in anti -
> diffusion program. ( determine sw_u sw_v and sw_w )
> Does it depend on user's situation? Must we be cautious
> on determining ' ES' ?
>
The switch option (described in the 84 paper) is quite suspicious to me too
and does not seem to work. The paper did say that it is quite empirical & the
constants are not universal. I tried instead to skip the switch & simply set SW
to constant (should be SW=1, but I get some overshooting so I found that setting
smaller values such as SW=0.5 give smoother solution).
>
> (3) In SMOL_ADIF, in 'Do 110' loop
> There is the following statement.
> ' IF(FF(I,J,K).LT.VALUE_MIN.OR.FF(I-1,J,K).LT.VALUE_MIN---'
>
> It seems this is for no advection on land area. But, FF
> sometimes can approach 0. There is advection even if FF is
> nearly equal to 0.
>
Yes, the values of FF (say temp.) should not be too close to zero (no
TBIAS) , this is a problem of many positive definition advection schemes.
*** BOUNDARY CONDITIONS
For more information concerning open boundary conditions in POM see:
1. POM users guide.
2. Table of possible BC (by P. Holloway, see POM web page)
3. Shulman, I., and J. K. Lewis, Optimization approach to the treatment of open
boundary conditions, J. Phys. Oceanogr., 25, 1006-1011, 1995.
Q:
>
> In the diagram
> illustrating the C-grid that POM uses in subroutine BCOND, it shows that
> the boundary points are at I=1+1/2, and 2 at the west boundary,
> I=IM,IM+1/2 at the east boundary, J=1+1/2, and 2, at the south boundary,
> and J=JM,JM+1/2 at the north boundary. My question is, for a closed
> lateral boundary, where exactly the 'solid wall' presents? for instance,
> if the western boundary is a solid wall, the wall is located at I=1+1/2 or
> 2?
A:
It depends on the variable. For example, for a west wall H(1,J)=1
(FSM(1,J)=0), but U(1,J)=U(2,J)=0 (DUM(2,J)=0). The first boundary points in
the west and south are not really used (U(1,J) & V(I,1)), but for cosmetics
put same values as in point 2.
Q:
This is a general question. When I use POM to do some simple
simulations, I found out that whenever I have a net volume of fresh
water inflow into the model domain, the extra volume of water can not get
out of model domain, no matter what kind of open boundary condition
is.
I tried to nudge the sea level back to zero near the open boundary,
it works in terms of getting rid of the extra water, but the temperature
and salinity values near the boundary are abnormally high. The nudging of
sea level casued the errors in calculation of horizontal mass fluxes, thus
the net advection fluxes of temperature and salinity were wrong.
I couldn't find a decent way to get rid of water, and in the mean
time, do not mess up the physics.
A:
To maintain sea level:
1. If you use UA*D on east/west boundaries and VA*D on north/south boundaries,
the horizontally integrated flow must be zero.
2. If any part of the boundary uses elevation (not its derivative) as the
boundary condition - whether directly or as part of a radiation condition - the
elevation will remain bounded but not necessarilly at its initial condition.
The above of course applies to the external mode.
Q:
> The bottom drag coefficients 'Cz' is usually set
> constant. But, in POM,
> Cz = max ( k^2/(ln((1+sigma(kb-1))H/z0))^2 , 0.0025)
> What is the ground of this eq.? ( Where does this come
> from? ) Could you let me know the reference paper or
> book.
A:
In PROFU and PROFV the (quite complicated) B.C. at the bottom is
KdU/dz= tau0x, KdV/dz= tau0y
The form given in the Users Guide, Eqn (14c,d,e), can be derived
using the law of the wall formulas (see any boundary layer text)
U(z')=(tau0x/k*utau)*ln(z'/z0)
V(z')=(tau0y/k*utau)*ln(z'/z0)
where z' is distance from the nearest stationary surface and
utau**4=(tau0x**2+tau0y**2) and z0 is the roughness parameter
related to the actual roughness of the bottom, presumed known.
Thus, the value of the "drag coefficient" as in (14e) depends on the
distance from the bottom and is generally not a constant. In (14e),
however, a constant is imposed when the velocity grid point
nearest the bottom is large such that the (bottom) boundary layer
is not resolved.
The same methodology, generalized for a moving surface, can be applied
at the surface to find the surface velocity where the stress is
generally prescribed. The question is: what is z0 at the surface
in a wave environment?
Stay tuned.
Prof. George L. Mellor
*** COMPILATION PROBLEMS
Q:
I have been getting error messages while trying to run POM97. I get the
following results when I compile and run the program.
f77 pom97.f -o pom
Note: IEEE floating-point exception flags raised:
Inexact; Underflow;
See the Numerical Computation Guide, ieee_flags(3M)
I am using a SUN Sparc10 workstation with system SunOS 5.4 and the most
current FORTRAN77 compiler for SUN.
A:
This is simply a warning and can safely be ignored. It has
nothing to do with your output. I get rid of it by ending my
program with "call exit(0)" rather than "stop" on the Sun. "stop"
is better on the IBM RS/6000 because flushes the output buffers
while "call exit(0)" will not.
Hmmmmm... While I agree that there is no problem, I am not
so confident in your 'call exit(0)' solution. What happens if
there is a real problem? Does this eliminate the message?
By a real problem, I mean a divide by zero or other infinite
number result, a parity error, or some other serious problem.
I tried the 'call exit(0)' with a piece of code that divides by 0
and discovered that the error messages disappear! Better not to use
this fix for the warning messages unless you check the specifics of
the problem before a 'call exit(0)' statement.
*** CPU RUN TIME ON DIFFERENT COMPUTERS
----- Please see information on this topic on the POM web page ------
*** DENSITY
Q:
> I am trying to get a hold on the reference that explains the
> calculations in subroutine Density in the POM code.
The reference is included in the users guide, but here it is again:
Mellor, G. L., An equation of state for numerical models of ocean and
estuaries, J. Atmos. Oceanic Technol., 8, 609-611, 1991.
Q:
> I'm using POM to model currents field in a deep lake of fresh =
>water, but I have some doubts on right POM setting and particularly
>on right DENS subroutine for fresh water, in fact density and salinity =
>in this case are very different from these of the test-case.
>Is there someone who has already prepared a DENS subroutine, based on =
>limnological equation of state developed by Chen and Millero, to compute =
>density and could help me ?
>With thanks in anticipation!
A: (William O'Connor)
The SUBROUTINE DENS in the POM should be able to calculate freshwater
denisty if you set the salinity to the low value found in fresh water. For
Great Lakes applications the salinity is typically 0.2 psu (0.2 ppt), so
SB(I,J,K)=0.2
You might want to take a look at the application of the POM to Lake Erie
and Lake Michigan, in the articles:
Beletsky, D., O'Connor, W.P., Schwab, D.J., and Dietrich, D.E., 1997.
Numerical Simulation of Internal Kelvin Waves and Coastal Upwelling
Fronts, Journal of Physical Oceanography, 27(7), 1197-1215.
Schwab, D.J., Beletsky, D., O'Connor, W.P., and Dietrich, D.E., 1996.
Numerical Simulation of Internal Kelvin Waves with Z-level and Sigma
Level Models, in: M.L. Spaulding and R.T. Cheng (eds.),Estuarine and
Coastal Modeling, Proceedings of the 4th International Conference,
Oct. 26-28, 1995, San Diego, CA, American Society of Civil Engineers,
New York, NY, 298-312.
Schwab, D.J., O'Connor, W.P., and Mellor, G.L., 1995. On the Net
Cyclonic Circulation in Large Stratified Lakes, Journal of Physical
Oceanography, 25(6), Part II, 1516-1520.
O'Connor, W.P., and Schwab, D.J., 1994. Sensitivity of Great Lakes
Forecasting System Nowcasts to Meteorological Fields and Model Parameters,
in: M.L. Spaulding, K. Bedford, A. Blumberg, R. Cheng, and C. Swanson
(eds.), Estuarine and Coastal Modeling III, Proceedings of the 3rd
International Conference, Sept. 8-10, 1993, Oak Brook, IL, American
Society of Civil Engineers, New York, NY, 149-157.
The GREAT LAKES COASTAL FORECASTING SYSTEM has been set up for Lake Erie
by Ohio State University and the NOAA Great Lakes Environmental Lab, using
the POM.
>
> I am not using potencial T (temperatures). Can this take to very bad
> solutions?
>
Not necesarily. Just keep in mind that when pressure effect in DENS is
neglected (or when using a linear eq. of state) you have to comment out the
pressure-dep. part in PROFQ (the line before 102 CONTINUE).
*** DIAGNOSTIC CALCULATIONS
Q:
> I read in the paper, 'Diagnostic and Prognostic Numerical Circulation
> Studies of the South Atlantic Bight', about running the model in a
> diagnostic mode to obtain the climitological velocity fields. I am not
> sure practically how to run the model in this mode.
Simply set MODE=4 (see POM users guide), it can be done with constant
surface wind forcing. The following paper describes the adjustment process
in such calculations.
Ezer, T. and G. L. Mellor, Diagnostic and prognostic calculations of the
North Atlantic circulation and sea level using a sigma coordinate ocean
model, J. Geophys. Res., 99(C7), 14,159-14,171, 1994.
*** HIGH ORDER SCHEME FOR POM
A high order interpolation scheme used by Chu & Fan (1998) has been
kindly provided to pom users and is now available in the pom ftp
directory: contrib-code/pom_highsub.f (it is basically a new BAROPG
subroutine).
The reference is (more references are included in the code):
Chu, P.C., and C.W. Fan,"A three-point combined compact difference
scheme," Journal of Computational Physics, 140, 370-399, 1998.
Tests I made show significant reduction in pressure gradient errors
(when horizontal resolution is insufficient & RMEAN is not used).
However, the present code does not vectorized well and on a vector
machines like the T90, cpu may be more than 10 times the standard
pom97.f code (on SGI the additional cost is only around 60%). If any
user has a suggestion on how to improve the efficiency of the code,
please let me know.
*** HORIZONTAL DIFFUSION
NOTE: see Ezer & Mellor, Sensitivity studies with the North Atlantic sigma
coordinate Princeton Ocean Model, Dyn. Atm. Oceans, 32, 185-208, 2000.
Q:
I am bewildered by the choice of the Smagorinsky diffusivity for horizontal
diffusion. The followings are the AM of different ahthors:
Value(m2/s) Author
400 Stanev(1990)
50(GCM) Figueroa and Olson(1994)
1000 Cessi(1995)
100-400 Meacham and Berloff(1997)
about 1E+4 G.L.Mellor(1985)
about 1E+3 Masuda(1995)
about 1E+4 Nishigaki(1997)
...
Does the diffusive coefficient depend on the horizontal scale of problems
studied? I do not know why there are so much difference among them.
Then, what is the relationship between AM and grid length? Is the AM
larger(smaller) in the model with high horizontal resolution than one with
low resolution?
A:
In the Smagorinsky formulation AM (and AH=TPRNI*AM) depends on
gridsize and velocity gradients, based on the physical assumption that
AM should be larger near sharp fronts and in coarser resolution grids
where more subgrid scale eddies are unresolved; in POM it also depends
on the coefficient HORCON (with a typical value of 0.1-0.2). In a
submitted paper I presented at the last POM meeting, sensitivity studies
with various values of HORCON and TPRNI did not show a too large effect
on a high resolution North Atlantic model, but they may affect processes
such as eddy shedding in the Gulf of Mexico (see L. Oey's paper on
this).
Q:
What is the order of magnitudes of AM for a basin with
horizontal scale with about 1000km?
A: (T.Ezer)
The answer depends on each application and how small value the numerical
scheme can tolerate. For example, a new general coordinate POM that compares
sigma to z-level grids using a coarse resolution (1x1deg) grid (Mellor et al.
1998) tests AM values of 1E+5, 1E+4 and 1E+3 m2/s. With the highest
viscosity, both models gave similar results, but with the two low
viscosity cases, the z-level model became unstable, while the sigma grid
seems to be able to tolerate the small values (in some cases even zero
viscosity).
Bill O'Connor adds:
The horizontal diffusion coefficient is set by a CELL REYNOLDS NUMBER
argument. The nondimensional Reynolds number is the ration of the advection
to the diffusion terms:
advection U (delta U) / (delta X)
R = --------- = ----------------------------
diffusion nu (delta U)/ (delta X)**2
where R= Reynolds number and nu is eddy diffusivity. Then
nu = U (delta X) / R
We might assume that R=10, at least, and this can be used to see what values
are expected for the viscosity nu (AAM in POM) in m**2/s.
And George Mellor adds:
To add to confusion, I note that a stability analysis (Appendix, Bryan
et al., J. Phys. Oceanogr. 5, 30-46, 1975), on the 1-D advection-diffusion
equation, proclaims that R < 2 lest a 2 DX computational mode be enabled.
POM repeatedly violates this rule; Apparently, the mode is not excited.
We have seen successful cases where R = O(100) and, in fact, in a highly
resolved estuary (Oey et al., J. Phys. Oceanogr., 15, 1676-1692, 1985)
R = infinity.
*** INTERPOLATION AND INITIALIZATION PROBLEMS
Q:
I would like to ramp up my BCs around the Makassar Strait.
The mid-latitude of the grid is ~5N which gives an inertial period of about
11-12 days. Since ramp = time/period, can I change the value of period here
to speed up the ramping ?
A:
The ramp intends for mid latitudes to depress inertial oscillations.
It is obviously not too useful in low latitudes, but you can specify whatever
period you want, say 1-2 days which should depress other sudden variations
associated with the initialization or forcing.
Q:
I want to use POM in orthogonal curvilinear grid. But I don't know how to
alter pom.f to an orthogonal curvilinear coordinate system because the basic
equations and the test case listed in POM user guide is in horizontal
cartesian coordinates.
A:
You do not need to modify the basic pom.f code, it is capable of using
curvilinear grids.
However, you will have to generate your own grid and initial condition
ALON(I,J), ALAT(I,J), H(I,J) etc. before running the model; READ(40) is an
example of grid and initial conditions that are needed and they will
overwrite the cartesian parameters used in the test case of pom.f. You
may use curvigrid.f in the ftp directory (or other more user friendly
utilities indicated on the POM web page) to generate your grid. See also
grid.f in the ftp subdirectory contrib-code for an example of a code to
generate grid and interpolate data into the grid.
Q:
I am going to start numerical similation in a orthogonal curvilinear grid
system, but I am not sure that how to compute dx(i,j) and dy(i,j) from the
coordinates value of the grid system. Could you please give me some help
A:
DX & DY are simply the grid size in meters; look in the pom ftp site in
directory contrib-code for grid.f for an example how to calculate them from
longitude & latitude.
Q:
How far away from orthogonal can one's numerical grid be before
noticeable errors or changes in results occur? I realize I can
formulate the governing equations in a generalized coordinate system
and evaluate the terms missing from a cartesian representation, but
am hoping here that someone might have a "rule-of-thumb" response.
A:
The way POM is set up, scalars should be
conserved - I think. The momentun equations will be in error; how much I don't
know. Conservation can be checked easily for a closed basin.
Alan Blumberg violated orthogonality in a Chesepeake Bay application.
Check the POM pub list.
Of course it would be good to have the horizontal grid generalized; I started
working on this but other problems intervened. (GLM)
More on grids (from Le Ly):
Our studies on orthogonal and nearly-orthogonal grids (NOG) using POM show that
errors are dependent on skew angles. POM can run with very large skew angles
(30- 45 deg) or more which show that the momentum equations do not make large
errors to blow the model up. Our Monterey Bay (MOB) multi-block grid version
runs with curvilinear NOG with a high grid density packed along steep slopes
and Montey Submarine Canyon reduces the sigma coordinate errors by 40 %
compared to orthogonal grid with the same number of grid points. With the same
forcing and number of horizontal and vertical grid point and other paramters
the MOB orthogonal model blows up at around 100 model days. In our MOB case,
only CNO model works!!!. Our experience with CNO and orthogonal grids show that
skew angles less than 20-25 deg are considered as a "slight" degree of
nonorthogonality. Regarding this problem, enginering flow and CFD people
really did very good jobs. Mastin ("Error induced by coordinate system" in
Numerical
Grid Generation, ed. J. Thompson, North-Holland, 1982) studied the source of
truncation error in the numerical solution on curvilinear coordinate systems
and concluded that the truncation error is dependent not only on the higher
order derivatives and the local grid spacing, but also on the rate of change of
the grid space and grid nonorthogoanlity. A "slight" degree of
nonorthogonality has a negligible effect on truncation error.
Q:
I wonder if the equations of POM are still valid for spherical coordinates
or need to be transformed. As I know, various Coriolis coefficient considers
the infulence of curvature. Is it enough?
A: POM should work fine for any curvilinear orthogonal grids (as spherical
coordinates are), and in fact POM has been used for many basin-scale problems
(see publications on web). Just make sure that Coriolis=f(latitude). I belive
that curvature errors are insignifficant (at least compared to other numerical
errors).
For spherical coordinates
just specify DX(I,J)=RE*COS(RAD*LAT(I,J))*DELLON
and DY(I,J)=RE*DELLAT
where RE=radius of earth and RAD=PI/180.. DELLON and DELLAT are
the longitude and lattitude increments.
Q:
I run pom model in intertidal region.
So, I use the subroutine that is moving boundary method to represent the
tidal flat.
In Pom, land is represented 'h(i,j)=1.0' but I wanna represent land
'h(i,j)=0.0' to denote water depth less than 1.0 .
Used above method, pom program don't work.
A:
Land is set in POM as H=1 instead of 0 to allow division by H.
Q:
What is the difference between TMEAN, SMEAN and RMEAN and how do I
get them.
A:
It should be noted first that TMEAN & SMEAN are now called TCLIM & SCLIM
(since the pom96.f version). They are usually climatological data interpolated
into the model grid (in most cases are the same as the initial condition) and
used only in the diffusion terms (in ADVT); they are subtracted from T & S
before diff. terms are calculated in order to reduce the along-sigma diapycnal
mixing along sloping bottoms. This procedure was found to be useful for
better BBL dynamics (see Mellor and Blumberg 1985).
RMEAN is the area averaged climatological density field (density(z))
interpolated into the sigma grid, it is subtracted from RHO to reduce
pressure gradient errors (see Mellor et al. 1994, for more detail).
Note that numerically POM can run with TMEAN=SMEAN=RMEAN=0.
See grid.f in contrib-codes for an example how to calculate these three fields
in your initialization program.
Q:
>
> I learned from page 4 in User's Guide for POM that you can "supply
> simple subroutines that convert constant z-level coordinate system to a
> sigma coordinate system and vice versa". But I can not find them in
> POM-ftpsite. Would you please inform me how I can get those simple
> subroutines? Thank you in advance for your help.
>
A:
In the pom ftp (ftp.gfdl.gov, pub/glm) go to subdirectory
contrib-code. The subroutines you are looking for are ZTOSIG and SIGTOZ.
In the same directory, there are examples of programs that use those
routines such as plot.tmp.f and grid.f.
Q:
I dont know how to use ZTOSIG subroutine.
A:
ZTOSIG is used to interpolate data from z-levels to the model
sigma coordinate, for example, to generate initial conditions from
Levitus climatology (see grid.f in ftp site subdir=contrib-code).
C ------------------------------------------------------------------
SUBROUTINE ZTOSIG(ZS,TB,ZZ,H,T,IM,JM,KS,KB)
C
C --- input ---
C ZS(KS) z-level depths
C ZZ(KB) sig-levels
C H(IM,JM) ocean depth
C TB(IM,JM,KS) data on KS levels (e.g., Levitus climatology)
C --- output ---
C T(IM,JM,KB) data interpolated onto the sigma grid
>
> We seem to find that the model bomb because of the rough topography,
> although this is not sure yet. Is there a simple criterion for how
> smooth the topography needs to be for POM not to bomb? Something like
> a critical ratio of changes in total depth per delta x?
>
Despite possible false near bottom velocities, POM usually does not blow up
because of pressure grad. errors so make sure its not for other reasons.
Attached is a subroutine (SLPMIN) we often use to smooth the topography based on
a criterion for max bound on dH/H of neighboring points. It is also useful to
add a one iteration with a Laplacian smoother (SMOOTH, also attached below).
Q: I am trying to get a hold on the reference that explains the
calculations in subroutine Density in the POM code.
A:
Mellor, G. L., An equation of state for numerical models of ocean and
estuaries, J. Atmos. Oceanic Technol., 8, 609-611, 1991.
Q:
>I need a model that can resolve the 3D circulation pattern in an
>archipelago with a number of fairly short and narrow straits. In the
>strait I'm working with right now, a 10-20m horizontal resolution is
>needed to resolve the topographical features of interest. Since the
>strait is also about 20m deep, I believe this is beyond the limit
>within which the hydrostatic approximation is valid.
A: It is not so much the size of the domain as the physical process you
wish to measure. If you wish to model a non-hydrostatic flow such
as flow over a dam or obstacles, or wakes in the lee of islands where vertical
velocities are important, or the free surface is perturbed in the lee, or
a strong internal bore, then yes, you would need a non-hydrostatic model.
As a comment, for those depths and grid spacing, you will need a very
short time step for the CFL criterion.
*** NETCDF VISUALIZATION
Q:
Is there anyone using, or working with, the NetCDF-ready POM97 code supplied
by Rich Signell/USGS at http://crusty.er.usgs.gov/omviz/; specifically the
binary NetCDF output file generated by the latter?
I downloaded this seamount test case code and ran it "as is". However, when
I try to read, import or NCdump the resulting "pom97.cdf" file, the various
software I've tried (e.g. OpenDX (http://www.opendx.org/) and GrADS
(http://grads.iges.org/grads/)) have difficulty with it. Specifically,
the GradsNC executable returns the following error:
A:
I just downloaded the NetCDF-enabled test case, compiled
it, ran it, and successfully looked at the output in Matlab using our
Matlab/Netcdf toolkit, and everything worked fine. I also looked
at the file using the NetCDF "ncdump" program with no problems.
In general, software that "reads NetCDF files" either reads NetCDF
files with specific conventions, or is a flexible interface that needs to
be programmed to deal with specific conventions. I believe GrADS
is the former, and DX is the latter.
Rich Signell rsignell@usgs.gov
U.S. Geological Survey (508) 457-2229 | FAX (508) 457-2310
384 Woods Hole Road http://crusty.er.usgs.gov/rsignell.html
Woods Hole, MA 02543-1598
*** PARALLEL POM
Q:
> Has POM been paralellized? How can I get the code?
A:
Due to recent requests by many users, I am making parallel (MPI & HPF)
codes available through our anonymous ftp Use at own risk, we can not provide
support for these codes at this time. The codes (tested mostly on sgi)
were converted by a group at U. of Minnesota and provided to us by S.
Piacsek (results presented at the Maine users meeting). A version of the
above MPI code, configured to look like the more familiar standard
seamount pom97.f case is being tested on T3E and is also provided (see
also benchmark comparisons on pom web). See README file for more info.
References to parallelizing pom by different groups (see POM PUBL):
Boukas et al. (1999); Luong et al. (2000); Oberpriller et al. (1999).
*** PARTICLE TRACKING ROUTINES
lmost two months ago I asked the group about the state-of-the
art in particle tracking routines for POM. I receieved a lot
of responses, which I have summarized below. As expected, there
are a lot of pieces of code using many different techniques, but
none are "ready-to-go" codes for POM that include the ability
to handle arbitrary orthogonal curvilinear geometry and also
include random walk in both the horizontal and vertical. It seems
we are pretty close, however, with all the code that exists.
--
Dr. Richard P. Signell | signell@saclantc.nato.int
NATO/SACLANT Undersea Research Centre | Tel: (+39) 0187 527 381
Viale San Bartolomeo 400 | Fax: (+39) 0187 527 331
19138 La Spezia, ITALY --> From USA/CANADA, use: APO AE 09613-5000
Here is the summary:
Eugene Wei (eugene.wei@noaa.gov) reports:
"In response to your inquiry about Lagragian trajectory model, I
would like to contribute the work that Frank Aikman and I did at
NOAA/NOS for Dumpsite 106 porject back in 1994. There was a NOAA
technical report and a journal article from this project. They
are:
Lagrangian Drifter Trajectory Modeling, by Eugene Wei, 1994,
106-Mile Dumpsite Oceanography Project, NOAA Technical Report NOS
OES 004.
Aikman, F. III and E.J. Wei, 1996. A Comparison of
Model-Simulated Trajectories and Observed Drifters in the Middle
Atlantic Bight. Journal of Marine Environmental Engineering,
Vol.2, pp. 123-139.
The computer code is designed particularly for taking POM
curvilinear grid velocity output and compute the trajectory. "
Briefly Wei compares 2 different numerical Lagrangian advection
techniques (Awaji vs. fourth-order Runge-Kutta) with analytical
solutions, and provides the Runge-Kutta computer code in his
report. No random walk component.
**************************************************************
Mark Hadfield (m.hadfield@niwa.cri.nz) reports that a particle
tracking routine for ROMS has just been developed which handles
arbitrary orthogonal curvilinear coordinates and a stretched
vertical coordinate. This would also have to be adapted to
POM, as the vertical coordinate is somewhat different. No random
walk component.
*****************************************************************
Brian Williams (brian.williams@newcastle.edu.au) reports:
I've been playing around with this stuff for some years (since
about 1975 actually). I have a cartesian coordinate version which
reflects particles off coastal boundaries and does a random
diffusive step which can be based on diffusivities determined by
POM. This version does not integrate across time steps ( although
I've done that in the past) and it does not deal with the higher
order terms in the random diffusive step due to steep bottom
slopes and rapid spatial change of the diffusion coefficient.
Order of magnitude analysis in my problems (shallow coastal
lakes) have suggested this is not a problem although the issue
was raised in a J.Comp Phys article a few years ago:
Hunter, J.R., P.D. Craig, and H.E. Phillips. 1993. On the use of
random walk models with spatially variable diffusivity. Journal
of Computational Physics 106(2): p 366-376. The comment which
followed is also useful: Heemink, A.W. and P.A. Blokland. 1995.
On random walk models with space varying diffusivity (Note).
Journal of Computational Physics 119: p 388-389.
*****************************************************************
David Schwab (schwab@glerl.noaa.gov) made some improvements
to the tracer.f code obtained from the CONTRIB area, making these
modifications:
C 1. Assume free surface displacement is not important for trajectory
C calculations
C 1. Fixed error in vertical displacement calculation
C 2. Moved all calling parameters to common block
C 3. Calculate z-coordinate based on (horizontally) interpolated depth
C 4. Fixed error in PARTZ calculation (found by Mike Settles)
C changed: PARTZ = (SIGMA - DEPTHP + DZ(KCELL))/DZ(KCELL)
C to: PARTZ = (DEPTHP + DZ(KCELL) - SIGMA)/DZ(KCELL)
C 5. Implemented Bennett and Clites 2nd order scheme for horizontal
components
Dave said that his version worked in his application, but should be
tested to
make sure it works in generic 3D situations. No random walk
component.
*******************************************************************
Roger Proctor (rp@pol.ac.uk) pointed out that there is a particle
tracking option in COHERENS (COupled Hydrodynamical Ecological
model for REgioNal Shelf seas). You can get a free CD-ROM, with
code, test cases and an incredibly detailed 914 page users manual
(pdf) at: http://www.mumm.ac.be/~patrick/mast/coherens.html
COHERENS uses sigma-levels in the vertical, but is limited to
Cartesian or spherical horizontal grids, and therefore would
have to be adapted to the fully orthogonal curvilinear
coordinates of POM. It does include a random walk component in
the vertical and horizontal.
**********************************************************
Konstantin Korotenko reports:
I used coupled 3D Flow/particle models in different aspects of studying
matter transport in coastal ocean. The models based on POM and other
flow
models and my own particle tracking models designed for oil pollution,
sediments, phytoplankton..etc transport and dispersal.
2001 A couple high-resolution circulation and particle-tracking
model for study oil spill transport in the Black Sea (with
M. J. Bowman, D. E. Dietrich)
2001 Prediction of the Dispersal of Oil Transport in the
Caspian Sea Resulting from a Continuous Release (with R.M.
Mamedov and C. N. K. Mooers) Oil Pollution Bulletin (in Press).
2001 Modeling oil spill transport processes in sea coastal
waters of the Caspian Sea. Oceanology, 41(1), 37-48.
2000 3D Flow/Particle Model for Transport and Dispersal Oil
Spills In Coastal Waters: the Caspian Sea. Physics and Chemistry
of the Earth, EGS (in press)
2000 Matter Transport in Meso-Scale Oceanic Fronts of River
Discharge Type. J. Marine Systems. Vol 24 No1-2. p.p.85-95.
1999 Eddy structure of horizontal water movement in shallow
seas with a special reference to the North Sea (with G. van Dam
and R.Ozmidov). J. Marine Systems. Vol 21, Nos1-4, pp. 207-228
************************************************************
>From Bruce Nairm (Bruce.Nairn@METROKC.GOV):
I have done some slight modifications to the tracer.f routine:
changing the interpolation scheme to one that does linear
interpolation on each velocity component using the eight
surrounding velocity vectors. I also added the transform to
convert the velocities from sigma to z for the calculation. This
has all been done for a Cartesian grid.
****************************************************************
I am also aware of some thesis work by Xue Yong Zhang at MIT in
1995 on particle tracking for ECOM-si, an semi-implicit two time
level version of POM. The nice thing about this work is that he
transformed the Fokker-Planck equation into orthogonal
curvilinear coordinates so that the random walk can be done
easily on the (i,j,k) coordinates. He had some very nice
examples demonstrating Gaussian spreading of a cloud of particles
under uniform diffusivity on highly stretched grids. Eric Adams
(eeadams@mit.edu) said he will see one of his students can put
this into the current version of POM.
*** PC COMPILATION AND RUNNING
Q:
Can POM be compiled & run on PC nachines and what Fortran compiler
should I use?.
A:
Here are results from experiments I did with pom97.f (seamount case)
Running Time of pom97.f on Different Computers
----------------------------------------------
All runs with DTE=6 sec, ISPLIT=30, DAYS=1 (Seamount topography)
Grid Points
Computer Compiler (20x10x12) (65x49x21)
-------- -------- ---------- -----------
MEM=0.3MW MEM=2.8MW
T90 f90 4s 56s
WinNT Digital Visual Fortran 39s (x10) 1336s (x23)
(Dell PentiumII 333MHz)
GNU/LINUX pgf90 52s (x13) 1969s (x35)
(Dell PentiumIII 600MHz)
(-fast option cut 10%
-Mconcur (2 process.) 60 & 50%)
SGI f77 96s (x24) 3080s (x55)
(IRIX5.3 IP22)
(IRIX6.2 IP28 - x3 faster than IRIX5.3)
Please see below comments of different users on their experience of running
pom on PCs and see also on the main web page a separate contribution by Oey
of benchmark tests including ORIGIN200 & DEC Alpha.
-----
From: Steve Brenner
On PC I have run POM with SVS Fortran, Lahey Fortran and Digital's Visual
Fortran all with no problem. I have also tried Linux and the performance
of the standard fortran is rather disappointing.
Steve Brenner
Israel Oceanographic and Limnological Research
PO Box 8030
Haifa 31080 Israel
phone: +972 4 8515202
fax: +972 4 8511911
e-mail: sbrenner@mail.biu.ac.il
steve@ocean.org.il
-----
Use Windows95/98 or WindowsNT as your operation system and Microsoft Fortran
PowerStation as Fortran compiler.
Yue Fang
Dept of Physical Oceanography
Institute of Oceanology, Chinese Academy of Sciences
7 Nanhai Rd., Qingdao 266071, China
Work : 86-532-2879062 ext 5810
Fax : 86-532-2870882
Home : 86-532-5895867
Email: yfang@ms.qdio.ac.cn
-----
From: Bryan Pearce
I've run POM exclusively on a PC. I've had to fix a few things. There
were some zero divides etc. I don't mind sending the code but would
feel sorry for the user as some things are hard wired. I've also changed
to code in places to speed it up a little. So it may be confusing. I've
also focused on the 2D part, thus there may be unexplored problems with
the 3D part. It runs OK but may or may not be correct.
I've used MS and now Digital fortran and liked it a lot.
I suppose this isn't much help but it does seem to work fine on the PC.
I've using a PenPro machine with dual processors. Only one thread in
the code though.
-----
From: Robin Robertson
I am running POM on a PC. I run it either on a Pentium 200 with
Windows 95 or on a Pentium 333 with Windows NT. The 333 is much
faster probably more due to using Windows NT. I use Fortran Power
Station for compiling and debugging. I have a slightly modified
version of POM, because I am running it for the Weddell Sea in a
region which includes ice shelves. Therefore there is a surface
frictional layer in part of the domain. I think it still works
for wind forcing, but it should be tested before I say for sure.
I use tidal forcing. Anyway, the person could contact me, but
I am not sure if my modified version of POM would be the right
thing for them. I also modified outputs and inputs for my uses.
Let me know if I can be of any assistance. The model runs faste
on my 200 PC than on the Unix machine I have access to by about
3X. I think there must be faster Unix machines though.
-----
From: Jerry Miller
I have recently run POM97 on a 266Mhz Pentium laptop with 192 Mbytes of RAM.
No changes were made to the version of the code downloaded from the web
site. It was compiled with Lahey Fortran90 under Windows95.
One day of simulation required 24 minutes of cpu time.
I know of others running POM on fast Pentium PC's with the Linux
operating system and associated compiler. They are able to do useful
runs but I have no information on required cpu time.
-----
I suggest installing Linux on the PC. That's by large the best of the
solutions. I have some 300Mhz Pentiun II running with g77 (Gnu Fortran).
They are much faster than some of the Suns and DEC's I have.
**********************************************************************
* Dr. Edmo J. D. Campos *
* Associate Professor of Physical Oceanography *
* Depto. de Oceanografia Fisica *
* Instituto Oceanografico - Universidade de Sao Paulo *
* Pca. do Oceanografico, 191 * 05508-900 Sao Paulo, SP, Brazil *
* edmo@usp.br - ph: (55) (11) 818-6597 - fax: 818-6597/210-3092 *
**********************************************************************
-----
From: Barbara Robson
I had a (modified, two-dimensional) version of POM working on a PC
about a year ago, using the Lahey Fortran 77 compiler. I don't have
the code to hand, but if you think it would be useful, I can probably
dig it up from my archives.
*** PROFQ (TURB. SCHEME)
Q:
>
> I have been trying to use POM to simulate internal
> waves in the Weddell Sea (near Antarctica). To
> determine background variations, I ran a simulation
> with a constant potential temperature and constant
> salinity. I found that POM determined that there
> was a N ~ .5 cph for a constant salinity and pot.
> temperature. This N results from the pressure
> term in rho, specifically due to the variation
> from eta. If I change the pressure calculation
> in Dens to use just H and not H + eta, the N goes
> away. My question is what will changing this
> to H do to the rest of the code. Will it cause
> problems in other calculations?
>
> By the way, why is the pressure correction for
> c different in Dens than it is in Profq when
> it is corrected for the N squared calculation.
> Or is one a rho correction for pressure and sound
> speed and the other one just a sound speed correction.
>
> Oh, the false N due to eta in the pressure term in
> Dens is significant enough to induce a internal tide
> near the M2 critical latitude.
>
A:
If you calculate N from in situ density, you will get the wrong
result. See the middle of p. 8 and Appendix A in the Users Guide.
N**2 = -BOYGR as calculated in S.R. PROFQ.
S.R. DENS deals with density. PROFQ deals with density gradients.
If, as an approximation, you wish to deal with potential density relative
to the surface, you had better get rid of the pressure corrections
in both DENS and PROFQ.
George Mellor
Q:
I'm curious about the variable 'sef' in profq, used when setting the
surface and bottom boundary conditions on q2.
There is a data statement which sets sef to 1.0. I assume that this isn't
the only valid value for sef, otherwise why include it? Is it just a
tuning parameter, or do sef values ne 1 have meaning? Assuming sef is
adjustable, what is a reasonable range and how is it determined?
A:
When the wind forcing is not resolved; e.g. monthly winds are used,
deepening of the mixed layer or mixing in general will be under-
estimated by the model. SEF is therefore a fudge factor to increase
mixing. Unfortunately, there needs to be research on ways of
estimating SEF as a function of the time over which the winds are
averaged and other relevant parameters.
George Mellor
*** RELAXATION OF T/S
Q:
I'd like to know if folks have used a weak damping of T,S back to
climatology for the deep ocean and if so, where/how is this inserted into
the model.
A:
From: Dong-Shan Ko
We have been using a weak damping of T,S back to climatology for
all our pom models. Here is the code we use (following Leo Oey),
do k = 1, kbm1
do j = 1, jm
do i = 1, im
ZRlx = DRlx*(1.-exp(ZZ(k)*H(i,j)/500.))
ZRlxr= 1.-ZRlx
UF(i,j,k)=ZRlxr*UF(i,j,k)+ZRlx*Tclim(i,j,k)
VF(i,j,k)=ZRlxr*VF(i,j,k)+ZRlx*Sclim(i,j,k)
end do
end do
end do
where DRlx is the damping coefficient (relaxation) defined as
DRlx = dt/T_scale
where dt = 2*DTI as pom uses a leap frog time stepping and T_scale,
damping time scale, is 180 days (NPAC2) or 250 days (NPAC4) in our
experiments.
A vertical weighting function, 1-exp(-z/z_scale),
is applied such that damping is in the deeper ocean only.
The code can be inserted in the subroutine BCOND before do loop 460.
We have some experiences applying damping, or relaxation as we call
it. Recently, we are working on a North Pacific model. First we applied
the relaxation without the vertical weighting (to save some computer time)
and found that although the El Nino is well simulated the longer term
variation is not for obvious reason (damping). With vertical weighting and
a smaller damping coeff. it improved as you can see from attached plot on
the tide gauge - model comparison for the two cases.
George Mellor has been suggested that the way pom does the diffusion
has an effect of damping back to climatology. In the hindsight, explicit
damping may not be needed in this simulation.
*** SEA-ICE COUPLING
Q:
I am interested in climate modelling of Antarctica region. We are
using POM as our ocean model. POM does not include sea-ice component.
So sea-ice model will have to be coupled with POM externally. But, I
do not know how to do this. Can you help me? It will be a great help
if you can providing names of scientists who have coupled POM with
Sea-Ice modeld.
A:
A sea-ice model, which is a very complicated issue, is not publically available
(at least not from Princeton). It has been coupled to POM by several users
(search the POM publications), the most prominent pom user scientist on this issue
is Sirpa Hakkinen of NASA, who published many papers on sea-ice modeling.
See also an early paper by Kantha & Mellor (1989).
*** STREAM FUNCTION
Q:
In the POM-code there
is a subroutine called findpsi. This subroutine can be easily adapted
to the special topography and region of interest. The stream function
have to be set to zero on points at the coast and then the calculation
have to be carried out considering the sweep direction of the loops as
it is done in the POM-code. Then the results of the loop in northern or
southern direction should give the same results compared to the loop
in western or eastern direction. As mentioned in the comment of the
subroutine the elevation field have to be near steady state. What I
miss is a consideration of the islands. I would expect, that the
stream function around an island have to be calculated constraining
that the mass transport around an island is zero. Anyone who thought
about it? Moreover the values of the stream function are very high.
Again any help would be appreciated.
A:
Subroutine FINDPSI can handle islands, just mask them out in your
output. Multiply PSI by E-6 to get units of Sverdrups.
*** SURFACE FORCING
Q:
About the surface flux (both heat flux: WTSURF and Salinity Flux: WSSURF),
if Both the climatological Sea Surface Temperature(SST) and Sea surface Sal.
(SSS) are known and I want to try 'resotring' forcing based on those data,
But how to selct the coefficents?
A1:
Here is a piece of code that I thought users may find useful for relaxation
surface BC, a similar approach can be used for lateral buffer zone BCs.
(Reference:
Ezer & Mellor, Sensitivity studies with the North Atlantic sigma coordinate
Princeton Ocean Model, Dyn. Atm. Ocean, 32, 185-208, 2000).
-------------------------------------------------
C RELAXATION TIME SCALE: 30d at surface to 360d at level 5
S1=DTI2/(30.*24.*3600.)
S2=DTI2/(360.*24.*3600.)
DO 1 K=1,5
COEFF=S1-(S1-S2)*(K-1)/4.
DO 1 I=1,IM
DO 1 J=1,JM
S(I,J,K)=( S(I,J,K)+COEFF*(Sobs(I,J)-S(I,J,K)) )*FSM(I,J)
1 CONTINUE
-------------------------------------------------
A2:
Tal's way to relax the t/s is similar to sarmiento and bryan [1982; also
used by oey and chen, jgr-oceans 1992, p.20,087, see eqn.27 and
discussion therein], which in some way is done already by
the horizontal mixing terms in subroutine advt [and use of t/sclimatology
there], except that it provides more control over the time-scales used
[rather than the variable smagorinsky]. one can also use the relaxation
as surface flux ONLY [see also oey and chen, eqn. 16 & 17] - which i
think is what Zhang asked; i have now modified pom97?? [which i did not
use] version of proft to use oey and chen's
surface flux relaxation - it also explains why/what etc. leo.
(Leo's modified PROFT.sub can be found in the ftp's contrib-codes)
Q:
What is the sign of the windstress and of the currents of water movement?
A:
EASTWARD and NORTHWARD WINDS STRESS have NEGATIV signs,
EASTWARD and NORTHWARD CURRENTS have POSITIV signs.
WUSURF and WVSURF have the same sign as the turbulent fluxes which,
conventionally, are the negatives of the surface stresses.
Q:
This is a question on how to set windstress direction in POM, especially,
for continuous model running (month1---> month12) in Indian ocean where
wind changes direction from NE (NOV-MAR)---> SW (June-SEP).
A:
The wind stress is set in the arrays WUSURF(I,J) and WVSURF(I,J).
WUSURF(I,J) is the X-direction (east-west) momentum flux at the ocean surface.
It is the negative of the wind stress divided by the water desity, and so it is
negative for westerly winds. The winds stress is in newtons/m**2 and the water
desnity is 1024 kg/m**3, so that WUSURF has dimensions of m**2/s**2. A wind
stress of 1 dyne/cm**2 from a westerly wind would be WUSURF(I,j)=- 1.e-4
A wind stress of 1 dyne/cm**2 from a southerly wind is WVSURF(I,J)=-1.e-4
(Bill O'Connor)
When using monthly (or even daily) wind stress fields it is advised to
interpolate in time at each time step to prevent artificial sudden
wind changes. (T.E.)
Q:
> I read about interpolating wind stress in time at each time step to
> prevent artificial sudden wind changes. I am using monthly averaged wind
> stress, and there is only one value for each month. Does it mean that I
> only have to do the interpolation at the last day of each month so that
> the wind stress can change smoothly to the new value for following month?
A:
No, the best way is to do the interpolation at each time step. You assume
that the monthly mean represents a value at the middle of each month and
interpolate between. For example, if at month 1 (say Jan. 15) you have a
value of V1 & at month 2 the value is V2 and you have 100 time steps
t=1,..100 per month (t=50 is about Jan 31), then the value at any point
between Jan 15 and Feb 15 can be estimated by
V(t)=V2+(100-t)*(V1-V2)/100. i.e., V(0)=V1, V(100)=V2, V(50)=(V1+V2)/2
Q:
I have heat flux data in W/m2, how do I convert it to WTSURF?
A:
The model surface heat flux forcing WTSURF (as well as the short wave
radiation SWRAD) are in units of Km/s; so you have to divide the fields by
the factor 4.1876E6. Note also that WTSURF is the flux from the ocean upward
into the atmosphere, so that WTSURF>0 -> cooling; some atmospheric data sets
(e.g., Oberhuber's COADS) give downward fluxes from the atmosphere to the
ocean, so you will need to change the sign too.
Q:
Can anyone tell me where the conversion factor comes from in dealing with
heat flux unit (watt/m^2 ---> K m/s) ?
A:
The surface temperature flux array WTSURF(IM,JM) has dimensions
(deg K)*m/s.
The heat flux is positive when the water is cooling.
The heat flux is negative when the water is warming.
Physically, the term in the equations is in
energy_flux watts/(m**2)
-------------------- = -------------------
density*specific_heat density*specific_heat
and
1 watt = 1 joule/s
density = 1024 kg/m**3
specific_heat=1 cal/(gm*degK) or find MKS units in physics textbook
1 cal = 4.19 joules
The result is
1 watt/(m**2) 1 degK*m
--------------------- = ------- ------
density*specific_heat 4.19E6 s
Q:
How can I use short wave radiation penetration in pom? and how surface
salinity fluxes are used?
A:
Set WTSURF=total heat flux minus short wave radiation & SWRAD=short
wave radiation (always negative) and in CALL PROFT set the parameter NBC=2.
The same PROFT is used also for salinity fluxes putting WSSURF (in units
of psu m/s) instead of WTSURF and set NBC=1.
Please look in the POM users guide (p.23-26, the description of PROFT)
for more detail.
*** SURFACE VELOCITY
Q:
> We are interested in using POM in a remote sensing context. We want
> to compare surface velocities with measurment. Is there a way to calculate
> actual surface velocities or is the uppermost velocity always used as the
> surface quantity?
> Also, how do people calculate vertical velocity shears at the
> surface, given that the relationship between vertical shear and windstress
> depends on the kinematic viscocity?
>
>
> Robert Fusina
> Naval Research Lab, code 7250
>
>
A:
The velocity asymptotes to the "law of the wall' which can be used
to extrapolate to the surface IF ONE KNOWS THE EFFECTIVE SURFACE ROUGHNESS
PARAMETER. The later is a research subject which I and others are working
on. It involves surface waves.
Kinematic viscosity should not be a factor. It is so for a
turbulent flow over a smooth surface. It is not a factor over a rough
surface or under a wavy surface I would think.
George Mellor
*** TIDAL FORCING
Q:
> We are interested in specifying the tides, together with subtidal frequency
> forcing, on the open boundary of our Prince William Sound Nowcast/Forecast
> System. If I recall correctly, you did some interesting work of
> this
> type a few years ago for the East Coast Forecasting System. Is there a ms
> or a technical report available?
A:
This paper was:
Chen, P. and G.L. Mellor, Determination of tidal boundary
forcing using tide station data. In: Coastal Ocean Prediction, Coastal and
Estuarine Studies, Vol. 56, C. N. K. Mooers (Ed.), American Geophysical
Union, Washington, DC, 329-351, 1999.
You can also search the publication list on the POM web page for other
tide-related papers. An example of a model with tides GOMtide.f is provided
in contrib-code ftp.
*** VERTICAL MIXING (M-Y TURBULENCE)
Q:
In the subroutine profq of pom98.f and earlier versions
KQ is calculated as .41 q l SM in pom2k it is .41 q l SH. Neither is
documented that I can find.
A1:
pom.change has an entry for 7 February 2001, which says:
In PROFQ, in the 280 do loop (near the bottom), change SM to SH in
the expression for KQ. This error will hardly be noticed, but it is
a bug.
A2:
Here is the history of SQ. SQ = constant = 0.2 was determined in the
'82 paper as that which gave the best results for neutral flow, both boundary
layer and channel flow. Later, I felt that the constant should be stability
dependent and set SQ = 0.41 SH; at the neutral limit, SH =0.49 and SQ=0.2.
Actually, there does not seem to be much difference in 1-D calculations whether
one uses SQ=0.2 or SQ=.41*SH. The flow properties are dominated by the behavior
of SM and SH stability functions, whereas SQ, like the constants in the length
scale equation, is a tuned constant. How SM instead of SH got used in one or
more versions of pom is a mystery to me.
George
*** VERTICAL VELOCITY
Q:
WE use the POM to calculate the flow field in lake Kinneret. The lake
is 22km long, 13km wide and 41m deep. We come up against difficulties
when calculate the vertical velocity W(i,j,k). The maximum value of
W(i,j,k) is only 2.5E-4 and the minimun value is only -6.5E-5 during the
period of strong west wind 20m/s of 20 days from my calculation. I think
the values are too small. Could you tell me why?
A:
W should be quite small, especially when wind is
homogenious and bottom topography is shallow and smooth.
By simple scaling of the continuity eq. W=U*H/L and current vel. are about
0.01x wind vel = 0.2m/s so W should be about 0.2*40/20000=4E-4 as you have.
Q:
Is there anyone who have written a routine to calculate the "real"
vertical velosity and who would share it? The output w of the POM
model is a transformed velocity - in the user guide it is called a
velocity normal to the sigma surfaces.
A:
The full expression of the transformation from omega (vertical velocity
on sigma coordinate) to vertical velocity W can be found in Ezer & Mellor
1997: Simulations of the Atlantic Ocean with a free surface sigma
coordinate ocean model, JGR, 102(C7), 15,647-15,657. For a steady state
deep flow where surface elevation<< water depth, a good approximation
is W=OMEGA+SIGMA(uDH/Dx+vDH/Dy) which is solved by the following:
C
C --------- CALC. REAL (UPWARD) VERTICAL VELOCITY ----------------
C (assuming EL << H and D(EL)/Dt=0)
C
DO 810 K=1,KB-1
DO 810 J=2,JM-1
DO 810 I=2,IM-1
810 W(I,J,K+1)=W(I,J,K+1) + Z(K+1) * (
1 (U(I,J,K)+U(I,J,K+1))*(H(I,J)-H(I-1,J))/(DX(I,J)+DX(I-1,J))
2 +(V(I,J,K)+V(I,J,K+1))*(H(I,J)-H(I,J-1))/(DY(I,J)+DY(I,J-1)) )
C
Q:
W-field looks very checkered. I have read that this is
a problem of ARAKAWA B grids but about C grids I haven't notice it.
Is there any idea? Perhaps about filtering or so? Any help would be
appreciated.
A:
Yes, W fields are quite noisy, but after filtering with Laplacian smoother
they can be plotted as shown in the paper mentioned in the previous question.
*** WAVE TANK USAGE ***
Q:
I am trying to develop a Hydrodynamic
characterization model of a rectangular wave tank that contains a sloped
beach surface and inlets and outlets for variable tidal flows as well as
computer controlled waves. I was wondering if the POM model would work in
this situation.
A:
POM should work for the wave tank but it is a hydrostatic model
(an enhancenment ot non-hydrostatic applications is a future
possibility) and thus may not apply unless the wave number x
depth is less than, say, 0.5.