Model help:ROMS: Difference between revisions

From CSDMS
No edit summary
No edit summary
Line 1,183: Line 1,183:
The Forward-Backward equation is a second-order accurate and is stable for longer time steps than many other schemes. It is however unstable for advection terms.
The Forward-Backward equation is a second-order accurate and is stable for longer time steps than many other schemes. It is however unstable for advection terms.


The Forward-Backward Feedback (RK2-FB)equation, if we set β=ε
The Forward-Backward Feedback (RK2-FB)equation is one option for solving equations (56) and (57), if we set β=ε=0 in the equations, it becomes a standard second order Runge-Kutta scheme, which is unstable for a non-dissipative system. Adding the β and ε terms adds Forward-Backward feedback to this algorithm, and allows us to improve both its accuracy and stability. The choice of β and ε terms adds Forward-Backward feedback to this algorithm, and allows us to improve both its accuracy and stability. The choice of β = 1/3 and ε = 2/3 leads to a stable third-order scheme with α<sub>max</sub> = 2.14093.


LF-TR and LF-AM3 with FB Feedback is another possibility for solving equations (56) and (57). In this equation, keeping γ = 1/12 maintains third-order accuracy. The most accurate choices for β and ε are β = 17/20 and ε = 11/20, which yields a fourth-order scheme with low numerical dispersion and diffusion and α<sub>max</sub> = 1.851640.
One drawback of Forward-Backward Feedback and LF-TR and LF-AM3 with FB Feedback is the need to evaluate the right-hand-side terms twice per time step. The original forward-backward scheme manages to achieve α<sub>max</sub> = 2 while only evaluating each right-hand-side term once. The Generalized FB with an AB3-AM4 Step is second-order accurate for any set of values for β, γ, and ε. It is third-order accurate if β = 5/12. However, it has a wider stability range for β = 0.281105. Shchepetkin and McWilliams (2009) choose to use this scheme for the barotropic mode in their model with β = 0.281105, γ = 0.0880, and ε = 0.013, to obtain α<sub>max</sub> = 1.7802, close to the value of 2 for pure FB, but with better stability properties for the combination of waves, advection, and Coriolis terms.


==Examples==
==Examples==

Revision as of 16:44, 26 July 2011

The CSDMS Help System

ROMS

ROMS is a community model shared by a large user group around the world, with applications ranging from the study of entire ocean basins to coastal sub-regions.

Model introduction

ROMS incorporates advanced features and high-order numerics, allowing efficient and robust resolution of mesoscale dynamics in the oceanic and coastal domains. The model solves the free surface, hydrostatic, primitive equations of the fluid dynamics over variable topography using stretched; terrain-following coordinates in the vertical and orthogonal, curvilinear coordinates in the horizontal. This allows enhancement of the spatial resolution in the regions of interest.

Model parameters

Parameter Description Unit
Input directory Path to input directory -
Input File name Name of input file -
Site prefix Site prefix for input/output files -
Case prefix Case prefix for input/output files -
ROMS application title - -
C-preprocessing flag - -
Variable information filename - -
Runinterval timestep multiplied by timestep size (In cases where input is read from a file rather than GUI second
Parameter Description Unit
Number of nested grids - -
Number of nodes in I-direction number of I-direction interior rho-points -
Number of nodes in J-direction number of J-direction interior rho-points -
Number of vertical levels - -
Number of sediment bed layers - -
Number of active tracers usually 2; e.g., temperature and salinity -
Number of inactive passive tracers - -
Number of cohesive (mud) sediment tracers - -
Number of non-cohesive (sand) sediment tracers - -
Number of titles in I-direction used for parallel mode -
Number of titles in J-direction used for parallel mode -
Parameter Description Unit
Number of timesteps - -
Timestep size - seconds
Number of barotropic steps - -
Starting perturbation or iteration - -
Ending perturbation or iteration - -
Max number of 4DVar outer loop iterations - -
Max number of 4DVar inner loop iterations - -
Number of stochastic optimals interval divisions - -
Number of eigenvalues for GST analysis - -
Number of eigenvectors for GST analysis - -
Parameter Description Unit
Model restart flag - -
Switch to recycle restart records True (T) or False (F) -
Number of time-steps between restart records - -
Number of time-steps between stations records - -
Number of time-steps between floats records - -
Number of time-steps between info diagnostics - -
Switch to create (T) or append (F) files True (T) or False (F) -
Number of time-steps between history records - -
Number of time-steps between creation of new history file - -
Starting averages timestep - -
Number of time-steps between averages records - -
Number of time-steps between creation of new averages file - -
Starting diagnostics timestep - -
Number of time-steps between diagnostics records - -
Number of time-steps between creation of new diagnostics file - -
Number to recycle TLM time records True (T) or False (F) -
Number of time-steps between creation of new TLM file - -
Switch to recycle ADM time records True (T) or False (F) -
Number of time-steps between ADM records - -
Number of time-steps between creation of new ADM file - -
Number of time-steps between 4DVAR adjustment of SFF Number of time-steps between 4DVAR adjustment of surface forcing flux -
Number of time-steps between 4DVAR adjustment of OBF Number of time-steps between 4DVAR adjustment of open boundary fields -
Switch for GST restart Generalized Stability Theory restart switch (T or F) -
Maximum number of iterations for GST Maximum number of iterations for Generalized Stability Theory -
Checking point interval for GST Checking point interval for Generalized Stability Theory -
Parameter Description Unit
Relative accuracy of Ritz values in GST Relative accuracy of Ritz values in Generalized Stability Theory -
Harmonic horiz. diffusion coeff. for tracers (I-dir, NLM) Harmonic horizontal diffusion coefficient for tracers in NONLINEAR model in I-direction m2/s
Harmonic horiz. diffusion coeff. for tracers (J-dir, NLM) Harmonic horizontal diffusion coefficient for tracers in NONLINEAR model in J-direction m2/s
Biharmonic horiz. diffusion for tracers (I-dir, NLM) Biharmonic horizontal diffusion coefficient for tracers in NONLINEAR model in I-direction m4/s
Biharmonic horiz. diffusion for tracers (J-dir, NLM) Biharmonic horizontal diffusion coefficient for tracers in NONLINEAR model in J-direction m4/s
Harmonic horiz. diffusion coeff. for tracers (I-dir, ADJ) Harmonic horizontal diffusion coefficient for tracers in ADJOINT model in the I-direction m2/s
Harmonic horiz. diffusion coeff. for tracers (J-dir, ADJ) Harmonic horizontal diffusion coefficient for tracers in ADJOINT model in the I-direction m2/s
Biharmonic horiz. diffusion for tracers (I-dir, ADJ) Biharmonic horizontal diffusion coefficient for tracers in ADJOINT model in I-direction m4/s
Biharmonic horiz. diffusion for tracers (J-dir, NLM) Biharmonic horizontal diffusion coefficient for tracers in ADJOINT model in J-direction m4/s
Harmonic horiz. viscosity coeff. (NLM) Harmonic horizontal viscosity coefficient in NONLINEAR model m2/s
Biharmonic horiz. viscosity coeff. (NLM) Biharmonic horizontal viscosity coefficient in NONLINEAR model m4/s
Harmonic horiz. viscosity coeff. (ADJ) Harmonic horizontal viscosity coefficient in ADJOINT model m2/s
Biharmonic horiz. viscosity coeff. (ADJ) Biharmonic horizontal viscosity coefficient in ADJOINT model m4/s
Background vert. mixing coeff. for active tracers (NLM) Background vertical mixing coefficient for active tracers in NOLINEAR model m2/s
Background vert. mixing coeff. for passive tracers (NLM) Background vertical mixing coefficient for passive tracers in NOLINEAR model m2/s
Background vert. mixing coeff. for active tracers (ADI) Background vertical mixing coefficient for active tracers in ADJOINT model -
Background vert. mixing coeff. for passive tracers (ADI) Background vertical mixing coefficient for passive tracers in ADJOINT model -
Background vert. mixing coeff. for momentum (NLM) Background vertical mixing coefficient for momentum in NOLINEAR model m2/s
Background vert. mixing coeff. for momentum (ADJ) Background vertical mixing coefficient for momentum in ADJOINT model -
Parameter Description Unit
Turbulent closure parameter 1 - m2/s
Turbulent closure parameter 2 - m2/s
Turbulent closure parameter 3 - m2/s
Turbulent closure parameter 3 - m4/s
K-epsilon parameter for GLS K-epsilon parameter (Generic Length Scale closure) -
Turbulent kinetic energy exponent for GLS Turbulent kinetic energy exponent (Generic Length Scale closure) -
Turbulent length scale exponent for GLS Turbulent length scale exponent (Generic Length Scale closure) -
Min value of specific turbulent energy for GLS Min value of specific turbulent energy (Generic Length Scale closure) -
Min value of dissipation for GLS Min value of dissipation (Generic Length Scale closure) -
Stability coefficient (GLS, closure indep.) Stability coefficient (Generic Length Scale closure) -
Shear production coefficient (GLS, closure indep.) Shear production coefficient (Generic Length Scale closure) -
Dissipation coefficient (GLS, closure indep.) Dissipation coefficient (Generic Length Scale closure) -
Buoyancy production coefficient, minus (GLS, closure indep) Buoyancy production coefficient, minus (Generic Length Scale closure) -
Buoyancy production coefficient, plus (GLS, closure indep.) Buoyancy production coefficient, plus (Generic Length Scale closure) -
Constant Schmidt number for turb.KE diffusivity (GLS, closure indep.) Constant schmidt number for turbulent kinetic energy diffusivity (Generic Length Scale closure) -
Constant Schmidt number for turb. "psi" field (GLS, closure indep.) Constant schmidt number for turbulent generic statistical field, "psi" (Generic Length Scale closure) -
Charnok surface roughness - -
Roughness from wave amplitude - -
Roughness from wave dissipation - -
Graig and Banner wave breaking coefficient - -
Parameter Description Unit
Momentum stress constant 1 constant used in momentum stress computation m/s
Momentum stress constant 2 constant used in momentum stress computation -
Momentum stress constant 3 constant used in momentum stress computation m
Momentum stress constant 4 constant used in momentum stress computation m
Height of measurement for air humidity (bulk flux) Height of atmospheric measurement for air humidity (bulk flux) m
Height of measurement for air temperature (bulk flux) Height of atmospheric measurement for air temperature (bulk flux) m
Height of measurement for winds (bulk flux) Height of atmospheric measurement for winds (bulk flux) m
Min depth for wetting and drying - m
Jerlov water type for shortwave radiation depth scale Jerlov water type used to set vertical depth scale for shortwave radiation depth scale -
Deepest level to apply surf. momentum stress as a body force - -
Shallowest level to apply surf. momentum stress as a body force - -
Mean water density - kg/m3
Brunt-Vaisala frequency - 1/s2
Time-stamp for model initialization - days
Reference time origin for tidal forcing - days
Model reference time for output NetCDF units attribute - yyyymmdd.dd
Nudging/relaxation time scale 1 - days
Nudging/relaxation time scale 2 - days
Nudging/relaxation time scale 3 - days
Nudging/relaxation time scale 4 - days
Nudging/relaxation time scale 5 - days
Factor between passive and active open BCs Factor between passive (outflow) and active (inflow) open boundary conditions -
Linear equation of state density parameter - kg/m3
Linear equation of state temperature parameter - Celsius
Linear equation of state salinity parameter - PSU
Linear equation of state temperature coeff. - 1/Celsius
Linear equation of state salinity coeff. - 1/PSU
Slipperiness parameter: 1.0 (free slip) or -1.0 (no slip) - -
Switch to use point sources/sinks for temperature True (T) or False (F) -
Switch to use point sources/sinks for salinity True (T) or False (F) -
Switch to use point sources/sinks for intert tracer 1 True (T) or False (F) -
Parameter Description Unit
Choice of vertical transformation equation Choice of vertical transformation equation for terrain-following coordinates (1 or 2) -
Choice of vertical stretching function Choice of vertical stretching function for terrain-following coordinates (1,2 or 3) -
Surface control parameter for terrain-following coords Surface control/stretching parameter for terrain-following coords -
Bottom control parameter for terrain-following coords Bottom control/stretching parameter for terrain-following coords -
Width of surface or bottom layer which requires higher-res stretch - m
Parameter Description Unit
Model name Name of the model -
Author name name of the model author -

Uses ports

This will be something that the CSDMS facility will add

Provides ports

This will be something that the CSDMS facility will add

Main equations

1) Equations of motion

a) Momentum balance in the x-directions

[math]\displaystyle{ {\frac{\partial u}{\partial t}} + \vec{v} \nabla u - fv = -{\frac{\partial \phi}{\partial x}} - {\frac{\partial}{\partial z}} \left ( \overline{u'w'} - \nu {\frac{\partial u}{\partial z}} \right ) + F_{u} + D_{u} }[/math] (1)

b) Momentum balance in the y-directions

[math]\displaystyle{ {\frac{\partial v}{\partial t}} + \vec{v} \nabla v + fu = -{\frac{\partial \phi}{\partial y}} - {\frac{\partial}{\partial z}} \left ( \overline{v'w'} - \nu {\frac{\partial v}{\partial z}} \right )+ F_{v} + D_{v} }[/math] (2)

c) Advective-diffusive equations for a scalar concentration field C(x,y,z,t) (e.g. salinity, temperature, or nutrients)

[math]\displaystyle{ {\frac{\partial C}{\partial t}} + \vec{v} \nabla C = {\frac{\partial}{\partial z}} \left ( \overline{C'w'} - \nu_{\theta} {\frac{\partial C}{\partial z}} \right )+ F_{C} + D_{C} }[/math] (3)

e) Equation of state

[math]\displaystyle{ \rho = \rho \left (T,S,P \right ) }[/math] (4)

f) Vertical momentum equation (assume that the vertical pressure gradient balances the buoyancy force)

[math]\displaystyle{ {\frac{\partial \phi}{\partial z}} = {\frac{- \rho g}{\rho_{o}}} }[/math] (5)

g) continuity equation for an incompressible fluid

[math]\displaystyle{ {\frac{\partial u}{\partial x}} + {\frac{\partial v}{\partial y}} + {\frac{\partial w}{\partial z}} = 0 }[/math] (6)

h) Reynolds stresses (an overbar represents a time average and a prime represents a fluctuation above the mean)

[math]\displaystyle{ \overline{u'w'} = - K_{M} {\frac{\partial u}{\partial z}} }[/math] (7)
[math]\displaystyle{ \overline{v'w'} = - K_{M} {\frac{\partial v}{\partial z}} }[/math] (8)
[math]\displaystyle{ \overline{C'w'} = - K_{C} {\frac{\partial C}{\partial z}} }[/math] (9)

2) Vertical boundary conditions a) top boundary condition ( z = ζ (x,y,t ))

[math]\displaystyle{ K_{m} {\frac{\partial u}{\partial z}} = \tau_{S}^x \left (x,y,t \right ) }[/math] (10)
[math]\displaystyle{ K_{m} {\frac{\partial v}{\partial z}} = \tau_{S}^y \left (x,y,t \right ) }[/math] (11)
[math]\displaystyle{ K_{C} {\frac{\partial C}{\partial z}} = {\frac{Q_{C}}{\rho_{o} c_{P}}} }[/math] (12)
[math]\displaystyle{ w = {\frac{\partial \zeta}{\partial t}} }[/math] (13)

b) bottom boundary condition (z = -h(x,y))

[math]\displaystyle{ K_{m}{\frac{\partial u}{\partial z}} = \tau _{b}^x \left (x,y,t \right ) }[/math] (14)
[math]\displaystyle{ K_{m}{\frac{\partial v}{\partial z}} = \tau _{b}^y \left (x,y,t \right ) }[/math] (15)
[math]\displaystyle{ K_{C}{\frac{\partial C}{\partial z}} = 0 }[/math] (16)
[math]\displaystyle{ -w + \vec{v} \nabla h = 0 }[/math] (17)

c) Horizontal boundary conditions (Eastern and western boundary)(considering biharmonic friction)

[math]\displaystyle{ {\frac{\partial}{\partial x}} \left ( {\frac{h \nu}{mn}}{\frac{\partial^2 u}{\partial x^2}} \right ) = 0 }[/math] (18)

Horizontal boundary conditions (Northern and southern boundary)

[math]\displaystyle{ {\frac{\partial}{\partial y}} \left ( {\frac{h \nu}{mn}}{\frac{\partial^2 u}{\partial y^2}} \right ) = 0 }[/math] (19)

3) Vertical transformation equations

a) vertical coordinate transformation (been available in ROMS since 1999)

[math]\displaystyle{ z \left (x,y,\sigma, t \right ) = S \left (x,y,\sigma \right )+ \zeta \left (x,y,t \right ) [1 + {\frac{S \left (x,y,\sigma \right )}{h \left (x,y\right )}}] }[/math] (20)
[math]\displaystyle{ S \left (x,y,\sigma \right ) = h_{c} \sigma + [h \left (x,y\right ) - h_{c}] C \left (\sigma \right) }[/math] (21)
[math]\displaystyle{ S \left (x,y,\sigma \right ) =\left\{\begin{matrix} 0 & \sigma = 0, C \left (\sigma \right ) = 0 \\ - h \left (x,y \right ) & \sigma = -1, C \left (\sigma \right ) = -1 \end{matrix}\right. }[/math] (22)

b) vertical coordinate transformation (been available in ROMS since 2005)

[math]\displaystyle{ z \left (x,y,\sigma,t \right ) = \zeta \left (x,y,t\right ) + [\zeta \left (x,y,t\right ) + h \left (x,y \right )] S \left (x,y,\sigma \right ) }[/math] (23)
[math]\displaystyle{ S \left (x,y,\delta \right ) = {\frac{h_{c} \delta + h \left (x,y \right ) C \left ( \delta \right )}{h_{c} + h \left (x,y\right )}} }[/math] (24)
[math]\displaystyle{ S \left (x,y,\sigma \right ) =\left\{\begin{matrix} 0 & \sigma = 0, C \left (\sigma \right ) = 0 \\ -1 & \sigma = -1, C \left (\sigma \right ) = -1 \end{matrix}\right. }[/math] (25)

unperturbed depth:

[math]\displaystyle{ \hat{z} \left (x,y,\sigma \right ) \equiv h \left (x,y \right ) S \left (x,y,\sigma \right ) = h \left (x,y\right )[{\frac{h_{c}\sigma + h \left (x,y\right ) C \left (\sigma \right )}{h_{c} +h \left (x,y \right )}}] }[/math] (26)
[math]\displaystyle{ d \hat{z} = d \sigma h \left (x,y\right ) [{\frac{h_{c}}{h_{c} + h \left (x,y\right )}}] }[/math] (27)

4) Vertical stretching functions

a) Vstretching = 1 (Song and Haidvogel (1994))

[math]\displaystyle{ C \left (\sigma \right ) = \left (1-\theta_{B}\right ) {\frac{sinh \left (\theta_{S}\sigma \right )}{sinh \theta_{S}}} + \theta_{B}[{\frac{tanh[\theta_{S}\left (\sigma +{\frac{1}{2}}\right )]}{2 tanh \left ({\frac{1}{2}} \theta_{S}\right )}} - {\frac{1}{2}}] }[/math] (28)

b) Vstretching = 2 (A.Shchepetkin (2005))

[math]\displaystyle{ C \left (\sigma \right ) = \mu C_{sur}\left (\sigma\right ) + \left (1-\mu \right ) C_{bot} \left (\sigma \right ) }[/math] (29)
[math]\displaystyle{ C_{sur} \left (\sigma \right ) = {\frac{1 - cosh \left (\theta_{S} \sigma \right )}{cosh \left (\theta_{S} \right ) - 1}} }[/math] (30)
[math]\displaystyle{ C_{bot} \left (\sigma \right ) = {\frac{sinh[\theta_{B}\left (\sigma + 1 \right )]}{sinh \left (\theta_{B}\right )}} - 1 }[/math] (31)
[math]\displaystyle{ \mu = \left (\sigma +1 \right )^\alpha [1 + {\frac{\alpha}{\beta}} \left ( 1 - \left ( \sigma + 1 \right )^\beta \right )] }[/math] (32)

c) Vstretching = 3 (R.Geyer function for high bottom layer resolution in relatively shallow applications)

[math]\displaystyle{ C \left (\sigma \right ) = \mu C_{bot} \left (\sigma \right ) + \left (1 - \mu \right ) C_{sur} \left (\sigma \right ) }[/math] (33)
[math]\displaystyle{ C_{sur} \left (\sigma \right ) = - {\frac{log[cosh \left (\gamma \right ) abs \left (\sigma \right ) ^\left (\theta_{S}\right )]}{log [cosh \left (\gamma \right )]}} }[/math] (34)
[math]\displaystyle{ C_{bot} \left (\sigma \right ) = {\frac{log [cosh \left (\gamma \right ) \left (\sigma + 1 \right )^ \left (\theta_{B}\right )]}{log [cosh \left (\gamma \right )]}} - 1 }[/math] (35)
[math]\displaystyle{ \mu = {\frac{1}{2}} [1 - tanh \left (\gamma \left (\sigma + {\frac{1}{2}}\right )\right )] }[/math] (36)

d) Vstretching = 4 (A. Shchepetkin (2010) UCLA-ROMS current function)

Surface refinement function

[math]\displaystyle{ \left\{\begin{matrix} C \left ( \sigma \right ) = {\frac{1 - cosh \left (\theta_{S}\sigma \right )}{cosh \left (\theta_{S}\right ) - 1}} & \theta_{S} \gt 0 \\ - \sigma ^2 & \theta_{S} \lt = 0 \end{matrix}\right. }[/math] (37)

Bottom refinement function (θB > 0)

[math]\displaystyle{ C \left (\sigma \right ) = {\frac{exp \left (\theta_{B} C\left (\sigma \right )\right ) - 1}{1 - exp \left (- \theta_{B}\right )}} }[/math] (38)

5) Terrain-following coordinate transformation (the dynamic equations and boundary conditions are transformed based on the following equations)

[math]\displaystyle{ \hat{x} = x }[/math] (39)
[math]\displaystyle{ \hat{y} = y }[/math] (40)
[math]\displaystyle{ \sigma = \sigma \left (x,y,z \right ) }[/math] (41)
[math]\displaystyle{ z = z \left (x,y,\sigma \right ) }[/math] (42)
[math]\displaystyle{ \hat{t} = t }[/math] (43)
[math]\displaystyle{ \left ({\frac{\partial}{\partial x}}\right )_{z} = \left ({\frac{\partial}{\partial x}}\right ) - \left ({\frac{1}{H_{z}}}\right ) \left ({\frac{\partial z}{\partial x}}\right )_{\sigma} {\frac{\partial}{\partial \sigma}} }[/math] (44)
[math]\displaystyle{ \left ({\frac{\partial}{\partial y}} \right )_{z} = \left ({\frac{\partial}{\partial y}}\right )_{\sigma} - \left ({\frac{1}{H_{z}}}\right ) \left ({\frac{\partial z}{\partial y}}\right )_{\sigma} {\frac{\partial}{\partial \sigma}} }[/math] (45)
[math]\displaystyle{ {\frac{\partial}{\partial z}} = \left ({\frac{\partial \sigma}{\partial z}}\right ) {\frac{\partial}{\partial \sigma}} = {\frac{1}{H_{z}}} {\frac{\partial}{\partial \sigma}} }[/math] (46)

6) Curvilinear coordinate transformation

relationship of horizontal arc length to the differential distance

[math]\displaystyle{ \left (ds \right )_{\epsilon} = \left ({\frac{1}{m}} \right ) d \epsilon }[/math] (47)
[math]\displaystyle{ \left (ds \right )_{\eta} = \left ({\frac{1}{n}}\right ) d \eta }[/math] (48)

velocity components in this coordinate system (equations of motion can be re-written with the following equaitons)

[math]\displaystyle{ \overrightarrow{v} \hat{\epsilon} = u }[/math] (49)
[math]\displaystyle{ \overrightarrow{v} \hat{\eta} = v }[/math] (50)

7) Time-stepping schemes

Assuming that multiple equations for any number of variables are time stepped synchronously

[math]\displaystyle{ {\frac{\partial \phi \left ( t \right )}{\partial t}} = F \left ( t \right ) }[/math] (51)

a) Euler time step

[math]\displaystyle{ {\frac{\phi \left( t + \Delta t \right ) - \phi \left (t\right )}{\Delta t}} = F \left (t\right ) }[/math] (52)

b) Leapfrog

[math]\displaystyle{ {\frac{\phi \left ( t + \Delta t \right ) - \phi \left ( t - \Delta t \right )}{2 \Delta t}} = F \left (t \right ) }[/math] (53)
[math]\displaystyle{ {\frac{\phi \left ( t + \Delta t \right ) - \phi \left ( t \right )}{\Delta t}} = {\frac{1}{2}} [F \left ( t \right ) + F ^* \left ( t + \Delta t \right )] }[/math] (54)

c) Third-order Adams-Bashforth (AB3)

[math]\displaystyle{ {\frac{\phi \left ( t + \Delta t \right ) - \phi \left ( t \right )}{\Delta t}} = \alpha F \left ( t \right ) + \beta F \left ( t - \Delta t \right ) + \gamma F \left ( t - 2 \Delta t \right ) }[/math] (55)

Assume time stepping asynchronously

d) Forward-Backward

[math]\displaystyle{ {\frac{\partial \zeta}{\partial t}} = F \left (u \right ) }[/math] (56)
[math]\displaystyle{ {\frac{\partial u}{\partial t}} = G \left (\zeta \right ) }[/math] (57)
[math]\displaystyle{ \zeta^ \left (n + 1 \right ) = \zeta ^n + F \left ( u^n \right ) \Delta t }[/math] (58)
[math]\displaystyle{ u^ \left (n + 1 \right ) = u^n + G \left (\zeta^ \left (n + 1 \right ) \right ) \Delta t }[/math] (59)

e) Forward-Backward Feedback (RK2-FB) (one option for solving equation (56), (57))

[math]\displaystyle{ \zeta^ \left (n+1,* \right ) = \zeta^n + F \left (u^n \right ) \Delta t }[/math] (60)
[math]\displaystyle{ u^ \left (n + 1, * \right ) = u^n + [\beta G \left ( \zeta^ \left (n + 1, * \right )\right ) + \left ( 1 - \beta \right ) G \left (\zeta^n \right )] \Delta t }[/math] (61)
[math]\displaystyle{ \zeta^ \left (n+1 \right ) = \zeta^n + {\frac{1}{2}}[F \left (u^ \left (n+1, * \right ) \right ) + F \left ( u^n \right )] \Delta t }[/math] (62)
[math]\displaystyle{ u^ \left (n+1 \right ) = u^n + {\frac{1}{2}}[\epsilon G \left ( \zeta ^ \left (n+1 \right )\right ) + \left ( 1 - \epsilon \right ) G \left ( \zeta^ \left (n+1, *\right ) \right ) + G \left ( \zeta^n \right )] \Delta t }[/math] (63)

f) LF-TR and LF-AM3 with FB Feedback (another option to solve equation (56), (57))

[math]\displaystyle{ \zeta^ \left (n+1,* \right ) = \zeta^ \left (n-1\right ) + 2 F \left (u^n \right ) \Delta t }[/math] (64)
[math]\displaystyle{ u^ \left (n+1,*\right ) = u^\left (n-1\right ) + 2 \{ \left (1 - 2 \beta \right ) G \left (\zeta^n \right ) + \beta [G \left (\zeta^ \left (n+1,*\right )\right ) + G \left (\zeta^\left (n-1\right )\right )]\} \Delta t }[/math] (65)
[math]\displaystyle{ \zeta^\left (n+1\right ) = \zeta^n + [\left ({\frac{1}{2}} - \gamma \right ) F \left (u^ \left (n+1,*\right )\right ) + \left ({\frac{1}{2}} + 2 \gamma \right ) F \left (u^n \right ) - \gamma F \left (u^ \left (n-1\right )\right )] \Delta t }[/math] (66)
[math]\displaystyle{ u^ \left (n+1\right ) = u^n + \{ \left ({\frac{1}{2}} - \gamma \right ) [\epsilon G \left ( \zeta^ \left (n+1\right )\right ) + \left (1 - \epsilon \right ) G \left (\zeta^ \left (n+1,* \right ) \right )] + \left ({\frac{1}{2}} + 2 \gamma \right ) G \left (\zeta^n \right ) - \gamma G \left (\zeta^ \left (n-1\right )\right ) \} \Delta t }[/math] (67)

g) Generalized FB with an AB3-AM4 Step

[math]\displaystyle{ \zeta^ \left (n+1\right ) = \zeta^n + [\left ({\frac{3}{2}} + \beta \right ) F \left (u^n\right ) - \left ({\frac{1}{2}} + 2\beta \right ) F \left (u^ \left (n-1\right ) \right ) + \beta F \left (u^ \left (n-2\right )\right )] \Delta t }[/math] (68)
[math]\displaystyle{ u^ \left (n+1\right ) = u^n + [\left ({\frac{1}{2}} + \gamma + 2\epsilon \right ) G \left ( \zeta^ \left (n+1\right )\right ) + \left ( {\frac{1}{2}} - 2 \gamma - 3 \epsilon \right ) G \left (\zeta^n \right ) + \gamma G \left (\zeta^ \left (n-1\right )\right ) + \epsilon G \left (\zeta^ \left (n-2\right ) \right )] \Delta t }[/math] (69)


Notes

1 In the Boussinesq approximation, density variations are neglected in the momentum equations except in their contribution to the buoyancy force in the vertical momentum equation.

2 In the vertical transformation equation, if equation (a) is used, the value of hc and other input vertical transformation parameters must be strictly defined. Since the beginning of ROMS, hc = min (hmin,Tcline) where hmin = minval(h) and Tcline is the stretching parameter. It is not possible to build the vertical stretching coordinates when hc > hmin and using transformation (a) because it yields Δz / Δδ < 0. This has been a legacy issue in ROMS for years, and you either consistent with the previous set-up of your application or change the value of Tcline in all input NetCDF files. Therefore, in new applications you need to be sure that Tcline <= hmin when using the transformation (a). However, if transformation (b) is used, hc can be any positive value. It works for both hc <= hmin, or hc > hmin.

Compared to equation (a), Shchepetkin (personal communication) points out that equation (b) has several advantages:

  • Regardless of the design of C(σ), it behaves like equally-spaced sigma-coordinates in shallow regions, where h(x,y) << hc. This is advantageous because it avoids excessive resolution and associated CFL limitation is such areas.
  • Near-surface refinement behaves more or less like geopotential coordinates in deep regions (level thicknesses, Hz, do not depend or weakly depend on bathymetry), while near-bottom like sigma coordinates (Hz is roughly proportional to depth). This reduces the extreme r-factors near the bottom and reduces pressure gradient errors.
  • The true sigma-coordinate system is recovered as hc is close to infinity. This is useful when configuring applications with flat bathymetry and uniform level thickness. Practically, you can achieve this by setting Tcline to 1.0d+16 in ocean.in. This will set hc=1.0*1016. Although not necessary, we also recommend that you set θS=0 and θB=0.

In an undisturbed ocean state, corresponding to zero free-surface, z=S(x,y,σ). Shchepetkin and McWilliams (2005) denotes transformation (a) as an unperturbed coordinate system since all the depths are not affected by the displacements of the free-surface. This ensures that the vertical mass fluxes generated by a purely barotropic motion will vanish at every interface. On the other hand, for transformation equation (b), in an undistrubed ocean state, it yields unperturbed depths as equation (26) and (27), as a consequence, the uppermost grid box retains very little dependency from bathymetry in deep areas, where h_c << h(x,y). For example, if h_c = 250m and h(x,y) changes from 2000 to 6000m, the uppermost grid box changes only by a factor of 1.08 (less than 10%).

Now the transformation equation (b) is the default value for ROMS.

The vertical stretching functions have the following properties:

  • C(σ) is a dimensionless, nonlinear, monotonic function;
  • C(σ) is a continuous differentiable function, or a differentiable piecewise function with smooth transition;
  • C(σ) must be discritized in terms of the fractional stretched vertical coordinate σ;
  • C(σ) must be constrained by -1 <= C(σ) <= 0, that is C(σ)= 0 for σ = 0 at the free-surface, C(σ) = -1 for σ = -1 at the ocean bottom.

The 4 available stretching functions has been shown in equations (28)~(38), equation (28) shows the condition for Vstretching = 1, which has the following features:

  • 0 <= θS <= 20, 0 <= θB <= 1, respectively.
  • It is infinitely differentiable in σ.
  • The larger the value of θS, the more resolution is kept above hc.
  • For θB = 0, the resolution all goes to the surface as θS is increased.
  • For θB = 1, the resolution goes to both the surface and the bottom equally as θS is increased.
  • For θS is not equal to 0, there is a subtle mismatch in the discretization of the model equations, for instance in the horizontal viscosity term. We recommend that you stick with reasonable values of θS, say θS <= 8.
  • Some applications turn out to be sensitive to the value of θS used.

Equation (29)~ (32) shows the condition for Vstretching = 2, which is similar in meaning to equation (28), but the hyperbolic function in Csur is cosh instead of sinh and that when σ is close to 0, ΔC/Δσ is also close to 0.

Equations (33)~(36) shows the condition for Vstretching = 3, which is designed for high bottom boundary layer resolution in relatively shallow applications, and the typical values for the control parameters are:

  • γ = 3;
  • θS = 0.65 for minimal increase of surface resolution;
  • θS = 1.0 for significant surface amplification;
  • θB = 0.58 for no bottom amplification;
  • θB = 1.0 for significant bottom amplification;
  • θB = 1.5 for good bottom resolution for gravity flows;
  • θB = 3.0 for super-high bottom resolution.

Equations (37), (38) shows the condition for Vstretching = 4, where C(σ) is defined as a continuous, double stretching function. equation (38) is the second stretching of an already stretched transform (37). The resulting stretching function is continuous with respect to θS and θB as their values approach zero. The range of meaningful values for θS and θB are: 0 <= θS <= 10; 0 <= θB <= 4. However, users need to pay attention to extreme r-factor values near the bottom.

The terrain-following coordinate transformation is used as from the point of view of the computational model, it is highly convenient to introduce a stretched vertical coordinate system which essentially "flattens out" the variable bottom at z = -h(x,y).

The Curvilinear coordinate transformation is used in many applications of interest (e.g., flow adjacent to a coastal boundary), the fluid may be confined horizontally within an irregular region. In such problems, a horizontal coordinate system which conforms to the irregular lateral boundaries is advantageous. It is often also true in many geophysical problems that the simulated flow fields have regions of enhanced structure (e.g., boundary currents or fronts) which occupy a relatively small fraction of the physical/computational domain. In these problems, added efficiency can be gained by placing more computational resolution in such region.

The Euler time step is the simplest approximation, this method is accurate to first order in Δt, however, it is unconditionally unstable with respect to advection.

The Leapfrog time step is more accurate, but it is unconditionally unstable with respect to diffusion. Also, the even and odd time steps tend to diverge in a computational mode. This computational mode can be damped by taking correction steps. SCRUM's time step on the depth-integrated equations was a leapfrog step with a trapezoidal correction (LF-TR) on every step, which uses a leapfrog step to obtain an initial guess of Φ (t + Δt). This leapfrog-trapezoidal time step is stable with respect to diffusion and it strongly damps the computational mode. However, the right-hand-side terms are computed twice per time step.

The third-order Adams-Bashforth step is used for SCRUM's full 3-D fields. The coefficients for the equations are: α = 23/12, β = -4/3, γ = 5/12. This method requires one time level for the physical fields and three time levels of the right-hand-side information and requires special treatment on startup.

The Forward-Backward equation is a second-order accurate and is stable for longer time steps than many other schemes. It is however unstable for advection terms.

The Forward-Backward Feedback (RK2-FB)equation is one option for solving equations (56) and (57), if we set β=ε=0 in the equations, it becomes a standard second order Runge-Kutta scheme, which is unstable for a non-dissipative system. Adding the β and ε terms adds Forward-Backward feedback to this algorithm, and allows us to improve both its accuracy and stability. The choice of β and ε terms adds Forward-Backward feedback to this algorithm, and allows us to improve both its accuracy and stability. The choice of β = 1/3 and ε = 2/3 leads to a stable third-order scheme with αmax = 2.14093.

LF-TR and LF-AM3 with FB Feedback is another possibility for solving equations (56) and (57). In this equation, keeping γ = 1/12 maintains third-order accuracy. The most accurate choices for β and ε are β = 17/20 and ε = 11/20, which yields a fourth-order scheme with low numerical dispersion and diffusion and αmax = 1.851640.

One drawback of Forward-Backward Feedback and LF-TR and LF-AM3 with FB Feedback is the need to evaluate the right-hand-side terms twice per time step. The original forward-backward scheme manages to achieve αmax = 2 while only evaluating each right-hand-side term once. The Generalized FB with an AB3-AM4 Step is second-order accurate for any set of values for β, γ, and ε. It is third-order accurate if β = 5/12. However, it has a wider stability range for β = 0.281105. Shchepetkin and McWilliams (2009) choose to use this scheme for the barotropic mode in their model with β = 0.281105, γ = 0.0880, and ε = 0.013, to obtain αmax = 1.7802, close to the value of 2 for pure FB, but with better stability properties for the combination of waves, advection, and Coriolis terms.

Examples

An example run with input parameters, BLD files, as well as a figure / movie of the output

Follow the next steps to include images / movies of simulations:

See also: Help:Images or Help:Movies

Developer(s)

Name of the module developer(s)

References

Key papers

Links

Any link, eg. to the model questionnaire, etc.