Skip to content

Commit efeadac

Browse files
sbrus89sbrus89
authored
Merge pull request #230 from sbrus89/omega/vert-coord-design
Add vertical coordinate design document
2 parents 10d7afb + 2080878 commit efeadac

File tree

2 files changed

+268
-0
lines changed

2 files changed

+268
-0
lines changed

‎components/omega/doc/design/VertCoord.md‎

Lines changed: 267 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,267 @@
1+
(omega-design-vert-coord)=
2+
# Vertical Coordinate
3+
4+
## 1 Overview
5+
The vertical coordinate module will be responsible for computing and storing information related to the vertical mesh in Omega.
6+
Since Omega will be a non-Boussinesq model, the vertical independent variable will be pressure-based (expressed in terms of the pseudo-height) as opposed to $z$-based.
7+
Similar to MPAS-Ocean's support for tilted height ($z^\star$) coordinates, Omega V1 will support a general vertical coordinate, where pseudo thickness $\tilde{h}$ varies in proportion to the total pseudo thickness of the water column.
8+
The prognostic variable in the layered continuity equation is pseudo thickness, which can be used to calculate the total pressure at a given vertical layer interface.
9+
The geometric height of a given layer interface is still necessary in the model to compute sea level, the geopotential, and derivatives with respect to $z$, all of which affect the dynamics.
10+
The geometric thickness of a layer ($\Delta z$) can be found using the pseudo thickness, density, reference density, and known bottom elevation (positive up) relative to the reference geoid.
11+
The vertical coordinate module will also serve as a container for information related to the variable number of active vertical layers for a given ocean column (due to variations in bottom elevation and ice shelf cavities).
12+
13+
14+
## 2 Requirements
15+
16+
### 2.1 Requirement: Compute and store the total pressure at a given layer interface
17+
18+
The total pressure at each layer interface will be computed and stored within the vertical coordinate module.
19+
This involves a surface-down summation of the pseudo thickness variable times the gravitational acceleration, $g$, and the reference density, $\rho_0$, starting with the surface pressure.
20+
The pressure variable will be needed in the computation of the TEOS-10 equation of state and the pressure gradient.
21+
The variable will be registered with the `IOStreams` framework to enable it to be output during runtime.
22+
23+
### 2.2 Requirement: Compute and store the $z$ location of a given layer interface
24+
25+
The $z$ location at each layer interface will be computed and stored within the vertical coordinate module.
26+
In the model, $z$ is defined as being positive up from the geoid.
27+
The $z$ location is calculated in a bottom-up summation of the pseudo thickness times the specific volume and reference density, starting with the known bottom elevation.
28+
The $z$ location of the interfaces is needed to compute the geopotential gradient and the sea level.
29+
The variable will be registered with the `IOStreams` framework to enable it to be output during runtime.
30+
31+
### 2.3 Requirement: Compute and store the vertical layer bounds for each horizontal cell
32+
33+
All tendency terms require a variable vertical loop range to account for variation in bottom elevation and ice shelf cavities.
34+
These bounds will be computed on construction in the vertical coordinate module and used for setting the loop bounds for active layers in calculations over the vertical.
35+
36+
### 2.4 Requirement: Compute and store the geopotential
37+
38+
The geopotential will be computed and stored within the vertical coordinate module.
39+
Initially this will be calculated by summing together the $z$ location times $g$.
40+
In the future, tidal potential and self attraction and loading will also have optional (default off) contributions to the geopotential.
41+
42+
43+
### 2.5 Requirement: Compute and store the desired thickness used in the $p^\star$ vertical coordinate
44+
45+
The $p^\star$ vertical coordinate specifies a desired layer thickness that stretches the initial thicknesses by the bottom pressure anomaly.
46+
This involves adding a vertically-weighted bottom pressure perturbation to the reference (initial) thickness.
47+
This desired thickness will be computed and stored within the vertical coordinate module.
48+
The vertical advection algorithm will use this desired thickness to compute the vertical transport.
49+
50+
### 2.6 Desired: General Arbitrary Lagrangian Eulerian (ALE) coordinate support
51+
52+
In future versions of Omega, the $p^\star$ coordinate will be extended to a more general ALE coordinate
53+
54+
### 2.7 Desired: Vertical Lagrangian remapping (VLR) support
55+
56+
The vertical coordinate module should retain flexibility to support VLR in future versions of Omega.
57+
58+
## 3 Algorithmic Formulation
59+
In the layered non-Boussinesq equations solved in Omega (see [V0 governing equation document](OmegaV1GoverningEqns) for details), the prognostic variable for cell $i$ and layer $k$ is the pseudo thickness, $\tilde{h}_{i,k}$, so that the geometric thickness (in meters) is a diagnostic variable defined as:
60+
61+
$$ \Delta z_{i,k} = \rho_0 \alpha_{i,k} \tilde{h}_{i,k}. $$
62+
63+
The pressure at vertical cell interfaces is the found by summing the pseudo thicknesses:
64+
65+
$$ p_{i,k+1/2} = p_{i}^{surf} + g\rho_0 \sum_{k^\prime=1}^k \tilde{h}_{i,k^\prime}. $$
66+
67+
and at the midpoint by
68+
69+
$$ p_{i,k} = p_{i}^{surf} + g\rho_0 \sum_{k'=1}^{k-1} \tilde{h}_{i,k'} + \frac{1}{2} g\rho_0 \tilde{h}_{i,k} $$
70+
71+
The $z$ location of cell interfaces is found by summing the pseudo thicknesses from the bottom multiplied by the specific volume:
72+
73+
$$ z_{i,k+1/2} = z_i^{floor} + \rho_0 \sum_{k^\prime=k}^{K_{max}} \alpha_{i,k^\prime} \tilde{h}_{i,k^\prime}, $$
74+
75+
where $z_i^{floor}$ is the (positive-up) bottom elevation.
76+
The $z$ location of a layer midpoint is given by:
77+
78+
$$ z_{i,k} = z_i^{floor} + \frac{1}{2} \rho_0\alpha_{i,k} \tilde{h}_{i,k} + \rho_0\sum_{k^\prime= k+1}^{K_{max}} \alpha_{i,k^\prime} \tilde{h}_{i,k^\prime}. $$
79+
80+
Initially, the geopotential is the sum of the $z$ height times $g$.
81+
In the future, it will include contributions from the tidal potential ($\Phi_{tp}$) and self attraction and loading ($\Phi_{SAL}$):
82+
83+
$$ \Phi_{i,k} = \left( gz_{i,k} + \Phi_{tp} + \Phi_{SAL} \right). $$
84+
85+
The desired pseudo thickness used for the $p^\star$ coordinate is:
86+
87+
$$\tilde{h}_k^{p^\star} = \tilde{h}_k^{ref} + \left(\frac{p_B-p_{surf}}{g\rho_0} - \sum_{k^\prime=1}^K \tilde{h}_{k^\prime}^{ref} \right)\frac{W_k\tilde{h}_k^{ref}}{\sum_{k^\prime=1}^K W_{k^\prime}\tilde{h}_{k^\prime}^{ref}}, $$
88+
89+
where $\tilde{h}_k^{ref}$ is the initial reference pseudo thickness and the weights $W_k$ determine how pressure perturbations are distributed amongst the layers.
90+
Setting all $W_k$ values to a constant, corresponds to uniform stretching of the layers.
91+
Different weight values can be chosen such that perturbations are distributed unequally.
92+
93+
94+
## 4 Design
95+
96+
### 4.1 Data types and parameters
97+
98+
The `VertCoord` class will be used to compute pressure, $z$ height, and geopotential at each layer.
99+
It will also contain the vertical loop bounds for cells, edges, and vertices.
100+
```c++
101+
class VertCoord {
102+
public:
103+
I4 NVertLevels;
104+
I4 NVertLevelsP1;
105+
106+
// Variables computed
107+
Array2DReal PressureInterface;
108+
Array2DReal PressureMid;
109+
Array2DReal ZInterface;
110+
Array2DReal ZMid;
111+
Array2DReal GeopotentialMid;
112+
Array2DReal LayerThicknessTarget;
113+
114+
// Vertical loop bounds (computed on construction)
115+
Array1DI4 MinLevelCell;
116+
Array1DI4 MaxLevelCell;
117+
Array1DI4 MinLevelEdgeTop;
118+
Array1DI4 MaxLevelEdgeTop;
119+
Array1DI4 MinLevelEdgeBot;
120+
Array1DI4 MaxLevelEdgeBot;
121+
Array1DI4 MinLevelVertexTop;
122+
Array1DI4 MaxLevelVertexTop;
123+
Array1DI4 MinLevelVertexBot;
124+
Array1DI4 MaxLevelVertexBot;
125+
126+
// Variables read in from vert coord stream
127+
/// p star coordinate variables
128+
Array2DReal VertCoordMovementWeights;
129+
Array2DReal RefLayerThickness;
130+
131+
/// Variables read in from mesh file
132+
Array2DReal BottomDepth;
133+
private:
134+
135+
// Variables from HorzMesh
136+
I4 NCells;
137+
I4 NEdges;
138+
I4 NVertices;
139+
Array2DReal CellsOnEdge;
140+
Array2DReal CellsOnVertex;
141+
142+
static VertCoord *DefaultVertCoord;
143+
static std::map<std::string, std::unique_ptr<VertCoord>> AllVertCoords;
144+
}
145+
```
146+
147+
### 4.2 Methods
148+
149+
```c++
150+
class VertCoord{
151+
public:
152+
static VertCoord *create();
153+
void computePressure();
154+
void computeZHeight();
155+
void computeGeopotential();
156+
void computeTargetThickness();
157+
private:
158+
void minMaxLevelEdge();
159+
void minMaxLevelVertex();
160+
161+
}
162+
```
163+
#### 4.2.1 Creation
164+
165+
The constructor will be responsible for storing any static mesh information as private variables and handling the selection any user-specified options.
166+
It will also compute the vertical loop bounds on edges and vertices.
167+
```c++
168+
VertCoord::VertCoord(const HorzMesh *Mesh,
169+
int NVertLevels,
170+
Config *Options);
171+
```
172+
173+
The create method will take the same arguments as the constructor plus a name.
174+
It calls the constructor to create a new vertical coordinate instance, and put it in the static map of all vertical coordinate objects.
175+
It will return a pointer to the newly created object.
176+
```c++
177+
VertCoord *VertCoord::create(const std::string &Name,
178+
const HorzMesh *Mesh,
179+
int NVertLevels,
180+
Config *Options);
181+
```
182+
A vertical coordinate `IOStream` will be defined to read in the `BottomDepth`, `VertCoordMovementWeights` and `RefLayerThickness` variables.
183+
184+
#### 4.2.2 Initialization
185+
186+
The init method will create the default vertical coordinate and return an error code:
187+
```c++
188+
int VertCoord::init();
189+
```
190+
191+
#### 4.2.3 Retrieval
192+
193+
There will be methods for getting the default and non-default vertical coordinate instances:
194+
```c++
195+
VertCoord *VertCoord::getDefault();
196+
VertCoord *VertCoord::get(const std::string &Name);
197+
```
198+
199+
#### 4.2.4 Computation
200+
201+
The public `computePressure` will sum the pseudo thicknesses times $g$ from the top layer down, starting with the surface pressure.
202+
This will be done with a `parallel_scan` inside a `parallel_for` over the mesh cells using hierarchical parallelism.
203+
```c++
204+
void VertCoord::computePressure(const Array2DReal &PressureInterface,
205+
const Array2DReal &PressureMid,
206+
const Array2DReal &LayerThickness,
207+
const Array2DReal &SurfacePressure) {
208+
209+
}
210+
```
211+
212+
The public `computeZHeight` will sum the pseudo thicknesses times $\alpha$ from the bottom later up, starting with the bottom elevation.
213+
This will be done with a `parallel_scan` inside a `parallel_for` over the mesh cells using hierarchical parallelism.
214+
```c++
215+
void VertCoord::computeZHeight(const Array2DReal &ZInterface,
216+
const Array2DReal &ZMid,
217+
const Array2DReal &LayerThickness,
218+
const Array2DReal &SpecVol,
219+
const Array2DReal &BottomDepth) {
220+
221+
}
222+
```
223+
224+
The public `computeGeopotential` will sum together the $z$ height times $g$, the tidal potential, and self attraction and loading:
225+
```c++
226+
void VertCoord::computeGeopotential(const Array2DReal &GeopotentialMid,
227+
const Array2DReal &ZMid,
228+
const Array2DReal &TidalPotential
229+
const Array2DReal &SelfAttractionLoading) {
230+
231+
}
232+
```
233+
Tidal potential forcing and self attraction and loading will be default-off features.
234+
The will be added to (or excluded from) the geopotential based on config flags.
235+
236+
The public `computeTargetThickness` will determine the desired pseudo thickness used for the $p^\star$ vertical coordinate:
237+
```c++
238+
void VertCoord::computeTargetThickness(const Array2DReal &LayerThicknessPStar,
239+
const Array2DReal &VertCoordMovementWeights,
240+
const Array2DReal &RefLayerThickness) {
241+
242+
}
243+
```
244+
245+
The private `minMaxLevel` will determine the various vertical loop bounds on edges and vertices.
246+
It will compute the class member variables: `MinLevelEdgeTop`, `MaxLevelEdgeTop`, `MinLevelEdgeBot`, `MaxLevelEdgeBot`, `MinLevelVertexTop`, `MaxLevelVertexTop`, `MinLevelVertexBot`, and `MaxLevelVertexBot` from `MinLevelCell`, `MaxLevelCell`, `CellsOnEdge`, and `CellsOnVertex`.
247+
```c++
248+
void VertCoord::minMaxLevel( ) {
249+
250+
}
251+
```
252+
253+
#### 4.2.5 Destruction and removal
254+
255+
No operations are needed in the destructor.
256+
The erase method will remove a named vertical coordinate instance, whereas the clear method will remove all of
257+
them.
258+
Both will call the destructor in the process.
259+
```c++
260+
void VertCoord::erase(const std::string &Name);
261+
void VertCoord::clear();
262+
```
263+
264+
## Verification and Testing
265+
266+
### Test:
267+
Unit tests will be used to test each of the computations (computePressure, computeZHeight, computeGeopotential, computePStarThickness) for a given pseudo thickness array. Comparison will be made to a ''truth'' vertical mesh that includes spatially varying `minLevelCell` and `maxLevelCell`.

‎components/omega/doc/index.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -119,6 +119,7 @@ design/Timers
119119
design/TimeStepping
120120
design/Tracers
121121
design/TridiagonalSolver
122+
design/VertCoord
122123
design/VerticalMixingCoeff
123124
124125
design/Template

0 commit comments

Comments
 (0)