Ocean Modeling with a quasiLagrangian flowfollowing vertical coordinate

1
Ocean Modeling with a quasi-Lagrangian
(flow-following) vertical coordinate

• Rainer Bleck
• Shan Sun
• NASA Goddard Institute for Space Studies
• New York
• Qingdao training workshop, June 2008

2
Vertical grid considerations

• Ocean is shallow but still rich in vertical
  structure. Inadvertent vertical mixing must be
  avoided.
• Strong flows often occur near boundaries (top,
  bottom, side). Grid should provide good
  resolution there and make it easy to apply
  boundary conditions. (Sigma coordinate )
• Grid points that follow vertical motion
  (Lagrangian grid) can prevent numerical
  dispersion during wave-induced vertical
  transport. (Isopycnic coordinate )
• Sloping coordinate surfaces can make it difficult
  to compute the horizontal pressure gradient.
  (Cartesian coordinate )
• Fluids tend to form discontinuities (fronts).
  High resolution near fronts would be desirable.
  (Isopycnic coordinate )

3
(No Transcript)
4
Vertical grid considerations

• Ocean is shallow but still rich in vertical
  structure. Inadvertent vertical mixing must be
  avoided.
• Strong flows often occur near boundaries (top,
  bottom, side). Grid should provide good
  resolution there and make it easy to apply
  boundary conditions. (Sigma coordinate )
• Grid points that follow vertical motion
  (Lagrangian grid) can prevent numerical
  dispersion during wave-induced vertical
  transport. (Isopycnic coordinate )
• Sloping coordinate surfaces can make it difficult
  to compute the horizontal pressure gradient.
  (Cartesian coordinate )
• Fluids tend to form discontinuities (fronts).
  High resolution near fronts would be desirable.
  (Isopycnic coordinate )

5
General recipe
In particular
(pressure gradient in u,v equations)
Special nonsolenoidal cases
(popular in meteorology)
(a) s p
(popular in oceanography)
(b) s r
6
Vertical grid considerations

• Ocean is shallow but still rich in vertical
  structure. Inadvertent vertical mixing must be
  avoided.
• Strong flows often occur near boundaries (top,
  bottom, side). Grid should provide good
  resolution there and make it easy to apply
  boundary conditions. (Sigma coordinate )
• Grid points that follow vertical motion
  (Lagrangian grid) can prevent numerical
  dispersion during wave-induced vertical
  transport. (Isopycnic coordinate )
• Sloping coordinate surfaces can make it difficult
  to compute the horizontal pressure gradient.
  (Cartesian coordinate )
• Fluids tend to form discontinuities (fronts).
  High resolution near fronts would be desirable.
  (Isopycnic coordinate )

7

• Principal design element of isopycnal models
  Depth and (potential) density trade places as
  dependent / independent variables
• - same number of unknowns, same number of
  (prognostic) equations, but very different
  numerical properties
• The driving force for isopycnal model
  develop-ment is genetic diversity. Isopycnal
  models are not inherently better they are just
  different.

8
Why more than one ocean model?

• Physics laws are stated in terms of differential
  equations, but computers solve algebraic
  equations (gt truncation errors).
• Subgrid-scale processes must be formulated in
  terms of resolved-scale processes (gt closure
  errors).
• Both truncation and closure errors lead to model
  errors.
• Model development benefits if the effects of
  truncation and closure errors can be separated.
• Model intercomparison is one of the few available
  tools.

9

Main benefits of isopycnic coordinate

• explicit potential vorticity and potential
  enstrophy conservation (C grid only)
• reduction of numerically induced diapycnal mixing
  during advection diffusion.
• Dispersive effects of finite difference operators
  are hidden behind smoke screen of naturally
  occurring isopycnic subgridscale mixing

10

Main pitfalls of isopycnic coordinate

• degeneracy in unstratified water columns
• 2-term horizontal PGF is error-prone in steeply
  inclined layers (reduction to 1 term possible at
  the price of approximating state eqn)
• layer outcropping (gt "massless" layers)
• strongly varying layer thickness requires
  sophisticated advection schemes

11
General recipe
In particular
(pressure gradient in u,v equations)
Special nonsolenoidal cases
(popular in meteorology)
(a) s p
(popular in oceanography)
(b) s r
12
Pressure force problems in HYCOM

• In idealized isopycnic models that disregard
  separate effects of T and S on compressibility,
  the pressure gradient is a single-term expression
  (involving M Fpa).
• Thermobaricity adds a second term to the pressure
  gradient expression. The added term can become
  large in steeply inclined coordinate layers.
• The magnitude of the added term depends on an
  arbitrarily chosen reference T/S profile.
• The choice of reference profile affects the
  modeled circulation.

13
Grid degeneracy is main reason for introducing
hybrid vertical coordinate "Hybrid" means
different things to different people - linear
combination of 2 or more conventional coordinate
s (examples zsigma, zrho, zrhosigma) -
ALE (Arbitrary Lagrangian-Eulerian)
coordinate ALE maximizes size of isopycnic
subdomain.
14
ALE Arbitrary Lagrangian-Eulerian coordinate

• Original concept (Hirt et al., 1974) maintain
  Lagrangian character of coordinate but re-grid
  intermittently to keep grid points from fusing.
• In HYCOM, we apply ALE in the vertical only and
  re-grid for 2 reasons
• (1) to maintain minimum layer thickness
• (2) to nudge an entropy-related thermo-
  dynamic variable toward a prescribed
  layer-specific target value by importing fluid
  from above or below.
• Process (2) renders the grid quasi-isopycnic

15
Idealized vertical-meridional section through the
world ocean
warm / light
cold / dense
cold / dense
equator
Pole
Pole
16
Schematic coordinate layout in MICOM
coordinate layer 1
coordinate layer 2
coordinate layer 3
equator
Pole
Pole
17
MICOMs hybrid coordinate cousin HYCOM
c o o r d i n a t e l a y e r 1
c o o r -
d i n a t e
l a y e r 2
c o o r d i n a t e l a y e r 3
equator
Pole
Pole
18
Montevideo
south
Vertical section through HYCOM solution. Heavy
black lines coordinate surfaces. Shaded
contours potential density
19
The HYCOM grid generator

• Adjustable parameters
• minimum thickness of each layer
• target density of each layer

20
Continuity equation in generalized (s)
coordinates
(zero in fixed grids)
(zero in material coord.)
(known)
21
Layer 1
Layer 2
Layer 3
Stairstep profile of potential density r versus
depth
Layer 4
r1
r2
r3
r4
22
Layer 1
Layer 2
Layer 3
Blue arrows indicate some diabatic process
Layer 4
r1
r2
r3
r4
23
Altered profile
r1
r2
r3
r4
24
equal areas
The regridding step find new interface pressure
equal areas
r1
r2
r3
r4
25
Layer 1
Layer 2
Layer 3
Final outcome diabatic process translated into
interface movement
Layer 4
r1
r2
r3
r4
26
The prototype HYCOM re-gridder or grid
generator

• Design Principles
• T/S conservative
• Monotonicity-preserving (no new T/S extrema
  during re-gridding)
• Layer too dense gt entrain lighter water from
  above
• Layer too light gt entrain denser water from
  below
• Maintain finite layer thickness in upper ocean
  but allow massless layers on sea floor
• Minimize seasonal vertical migration of
  coor-dinate layers by keeping non-isopycnic
  layers near top of water column.

27

The prototype HYCOM grid generator (contd)

• Determine how much water from the neighboring
  layer (source layer) would be needed to
  restore target density.
• The amount needed, Dpneed, may exceed the amount
  available, Dpavail, in source layer.
• The amount ultimately transferred is min(Dpneed
  ,Dpavail - Dpmin).
• The minimum thickness Dpmin is prescribed.

28

The prototype HYCOM grid generator (contd)

• The condition Dpneed gt Dpavail typically occurs
  under the following conditions
• receiving layer is too dense
• restoration to target density requires more water
  from source layer than is available.
• The likelihood for this to happen is greatest at
  high latitudes immediately below the surface
• gt high-latitude near-surface layers are more
  likely to end up with constant thickness than
  layers elsewhere.

29

The prototype HYCOM grid generator (contd)

• Major challenge achieve smooth lateral
  transition between fixed-depth and isopycnic
  segments of a coordinate layer.
• Goal avoid sideways-looking algorithms, i.e.,
  accomplish transition through clever vertical
  re-gridding alone.
• Solution employ a cushion function. Details of
  the algorithm are as follows .

30

The prototype HYCOM grid generator (contd)

• The cushion function, which sets the final
  thickness of the source layer,
• leaves large positive Dp values unchanged
  cush(Dp)Dp
• returns a (small) constant value if Dp is large
  negative cush(Dp) const.
• links the two cases above by a smoothly varying
  function for intermediate values of Dp.

31
Loop through interfaces k1,2,3,(top-down)
Is layer k-1/2 lighter than target ?
Is layer k1/2 denser than target ?
no
no
yes
yes
D1 amount of layer k1/2 water needed to
restore layer k-1/2 to target
D1 amount of layer k-1/2 water needed to
restore layer k1/2 to target
D1 0
D1 0
D2 downward interface k displacement needed to
inflate layer k-1/2 to minimum thickness
D2 actual thickness minus minimum thickness of
layer k-1/2
D3 thickness of layer k1/2
Move interface k up by max(0,minD1,D2)
Move interface k down by minD3,max(D1,D2)
32
Practice Session

• Goal explore tradeoffs in vertical coordinate
  layout, e.g. .
• Eulerian versus Lagrangian grid
• choice of target densities
• choice of layer minimum thickness
• Tool toy numerical circulation model BOXOCEAN,
  written in Fortran 90
• Fast enough to allow simultaneous experimentation
  by several participants

33
BOXOCEAN

• Rectangular basin, beta plane
• 5 hybrid-isopycnic layers
• Flat bottom
• Wind- and buoyancy-forced
• Constant salinity
• Extremely primitive mixed layer
• ALE (Arbitrary Lagrangian-Eulerian) vertical
  coordinate
• Parallelized with OpenMP

34
Prognostic variables

• Horizontal velocity u,v
• Layer thickness dp
• Potential density th

Primary diagnostic variables

• Montgomery potential gz p/r
• Barotropic stream function
• Generalized vertical velocity

35
Parameters read in via namelist

• Initial layer thickness (x,y - independent)
• Minimum layer thickness
• Target potential density in each layer
• Initial potential density in each layer
• Double gyre wind forcing (yes/no)
• Thermal forcing (yes/no)
• North-south atmospheric temperature contrast
• Total run time (yrs)

36
Changeable parameters in module control

• Basin dimensions idm,jdm,kdm
• Mesh size (m) meshsz
• Time step (sec) delt
• Parameters controlling lateral mixing of
• momentum
• potential density
• layer interface pressure
• Sidewall slip boundary conditions
  (free-slip/no-slip)

37
Output options

• To generate graphics output using NCAR graphics,
  compile code with FC ncargf90 and set
  precompiler option DNCAR
• To omit graphics, i.e., generate line printer
  output only (in stdout), compile with FC ifort
  or other Fortran 90 compiler.
• stdout will show poor-mans contour line plots
  where characters are used to create zebra stripe
  effects (as in 1950s).
• Interface depth in cross sections is plotted
  using digits that indicate coordinate index.