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 05-31-01 TE. TOPICS: *** ADVECTION *** BOUNDARY CONDITIONS *** COMPILATION PROBLEMS *** CPU RUN TIME ON DIFFERENT COMPUTERS *** DENSITY *** DIAGNOSTIC CALCULATIONS *** HIGH ORDER SCHEME FOR POM *** HORIZONTAL DIFFUSION *** INTERPOLATION AND INITIALIZATION PROBLEMS *** NETCDF VISUALIZATION *** PARALLEL POM *** PC COMPILATION AND RUNNING *** PROFQ (TURB. SCHEME) *** RELAXATION OF T/S *** STREAM FUNCTION *** SURFACE FORCING *** SURFACE VELOCITY *** TIDAL FORCING *** VERTICAL VELOCITY *** WAVE TANK USAGE *** 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 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 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.gfdl.gov, dir=perm/POM/MPI (note its not the usual /pub/glm). 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 last 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). However, a word of warning, on the T3E I was able to get identical results to the serial pom97.f only for a constant grid, the source of the problem is still unknown. Any advice from expert parallel users who are willing to test the code will be greatly appreciated. T.E. (10/05/00) References to parallelizing pom by different groups (see POM PUBL): Boukas et al. (1999); Luong et al. (2000); Oberpriller et al. (1999). *** 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. *** 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. *** 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<