The schematic work-flow diagrams in Fig. 2.1 shown the step involved in the SURF-NEMO numerical platform.
The steps in the most recent version release are grouped as follows:
Initialization: the user specifies the values of the input simulation parameters for the ocean model in the configuration file (horizontal and vertical grids, subgrid scale parameterizations, etc.) for the specific experiment selected.
Access and download of the input datasets: this is an automated step in which the input datasets for the selected simulation period are downloaded from remote or local data repositories, as specified in the configuration file. The input data are the bathymetry, the coastline, the atmospheric forcing and the coarse resolution parent ocean model for the initial and lateral boundary condition datasets.
Spatial numerical grid generation: this is an automated step that generates the horizontal and vertical grid for the nested model.
Input data regridding: this is an automated step that generates the bottom topography, surface forcing, initial and open lateral boundary conditions datasets on the child grid.
Forecast: another automated step in which the NEMO code is exectuted to produce the final outputs.
Post-processing: in this step the visualization and analysis procedures of the final outputs are considered. Different options of postoprocessing are available, i.e. comparing parent/child fields, comparing the simulation results with in-situ or satellite datasets and converting the datasets into what?
Work-flow of the Relocatable SURF-NEMO platform.
The graphical calling function flow shown in Figure 2.2 represent all paths traversed through a program during its execution and shows how the program is completed from start to finish, step-by-step. The six macro-tasks identified are: (1) child meshmask generation; (2) atmospheric data regridding; (3) ocean IC data regridding; (4) ocean BC data regridding and OBC data extraction; (5) ocean model simulation; and (6) visualization and data analysis.
Graphical calling function flow of the Relocatable SURF-NEMO platform.
The computational tasks are not independent of each other. The dependency flow graph of the macrotasks is given in Figure 2.3. Each node (from A to F) represents a macro-task and the solid edges represent data dependencies among these macro-tasks. From node A, tree edge can take us to node B, node C and D. Node E can begin once all the lower nodes are completed.
Dependency flow graph of macro-tasks.
The first task executed by each macro-task consists of setting up the values of the input model parameters, which are obtained from the configuration files setParFree.json and setParFixed.json. The configuration file structure and the input parameters are described in Chapter 3 and Appendix A.
In this phase, the read_inJsonFree and read_inJsonFixed procedures are executed to respectively define the user-free and fixed input parameters required to execute the NEMO model and all of the pre- and post-processing tasks. The procedures set_pathData and set_fileData are also called to define all of the paths and file names used by the program for the specific numerical simulation.
The child grid is generated after the model configuration phase. The NEMO model uses the Arakawa C grid for spatial discretization, with the state variables defined on the staggered grid illustrated in Figure 2.4. In the C grid the scalar quantities (temperature T, salinity S, pressure p, density \(\rho\)) are defined at the center of each grid volume. The velocity field components (zonal u, meridional v and vertical w) are shifted by half a grid width in their respective direction so that they are defined at the edges of the grid volumes. The procedures executed in this phase are:
Generation of the child 2D-mesh.
Interpolation of the source bathymetric dataset on the generated child grid.
Generation of the child 3D-meshmask.
The staggered Arakawa C-grid used by NEMO ocean model.
The horizontal grid generation is managed by the NEMO-MESH code. SURF uses ONLY a rectangular (or latitude-longitude) grid in a spherical coordinate system \(\lambda,\varphi\).
The horizontal grid (expressed in degrees) is generated by specifying the number of points
\(n_{\lambda}\) and \(n_{\varphi}\),respectively, in zonal and meridional directions, and the respective grid sizes \(\Delta\lambda\) and \(\Delta\varphi\) (in degrees) and the longitude and latitude \((\lambda,\varphi)_{1,1}\) of the first row and first column of the T grid. On the \(\lambda\varphi\) plane, the location of the T points of the grid are:
$$
\begin{equation}
\begin{array}{ll}
\lambda_{i,j} = \lambda_{11} + (i-1) \Delta \lambda \hspace{0.5cm} \mbox{with} \hspace{0.2cm} i=1.....n_\lambda
\\
\varphi_{i,j} = \varphi_{11} + (j-1) \Delta \varphi \hspace{0.5cm} \mbox{with} \hspace{0.2cm} j=1.....n_\varphi
\end{array}
\end{equation}$$
The u, v points of the grid are shifted by half a grid width in the zonal e/o meridional direction, as indicated in Fig. 2.4.
In this phase the bathymetric dataset is interpolated on the child grid, which is required to generate the 3D meshmask. The procedures in this phase are:
Accessing and downloading the source bathymetry and the coastline datasets.
Manipulating the bathymetry dataset.
The spatial interpolation of the source bathymetric on the child grid.
A procedure for checking if the necessary input datasets are present in the experiment directory $PATH_IDEXP/data/data00/indata/ is executed. If some of the requested data are not present then the procedures downlCoastlineInfile and downlBathyInfile are automatically executed to download the coastline and the bathymetry, respectively, from remote or local data repositories as specified in the configuration file setParFree.json.
Before conducting the spatial interpolation on the child grid, the source bathymetry data product can be
manipulated if specified in the configuration file setParFree.json. Several methods can be used:
Add a constant value to the surface elevation for the whole nested region
(e.g., an inland body of water with water level below the global ocean level such as the Caspian Sea).
Set maximum and minimum values that are different from the source value
(e.g., if a minimum depth of 5 metres or a maximum depth of less than the actual depth is required).
Define the land/sea interface grid points according to the input coastline.
Set maximums and minimums inside sub-regions (e.g., to mask a specific area)
Smooth out bathymetry variations with the Shapiro filter.
This is a high order horizontal filter that efficiently removes small scale grid noise without affecting the physical structures of a field. A Shapiro filter of a 2N order of accuracy is applied to a variable based on the expression:
$$\begin{equation}
\label{eq:shapiro_filter}
\tilde{w_i} = F^{2N}(w_i) = \left[ I + (-1)^{N-1} \frac{\delta^{2N}}{2^{2N}} \right] (w_i) = w_i + (-1)^{N-1} \frac{\delta^{2N} w_i}{2^{2N}}
\end{equation}$$
where \(\tilde{w_i}\) is the filtered value of variable \(w\) at point \(x_i\), \(I\) is the identity operator and \(\delta^{2N}\) is the even composition of the standard difference operator \(\delta\) (Richtmyer, 1957).
This filter is a discrete symmetric operator with a (2N + 1) point stencil. It acts as a low-pass filter that preserves the low-frequency content (i.e., largest wavelengths) and completely dissipates the high-frequency content (i.e., shortest wavelengths) from the original field.
The spatial interpolation of the source bathymetric dataset on the child grid can be conducted after the manipulation of the source bathymetry phase. The procedures are based on the Spherical Coordinate Remapping and Interpolation Package (SCRIP) code.
The interpolation methods available are listed and described in Section 2.3.4.
The vertical grid generation is managed by the NEMO-MESH code. The type of vertical grid used in SURF corresponds to geopotential z-coordinate levels with partial bottom cell representation of the bathymetry. After the bathymetry \(z = H(\lambda,\varphi)\) and the number of levels \(n_{z}\) have been specified, the vertical location of w- and t-levels (expressed in metres) are defined, except in the bottom layer, from the following analytic expression:
$$
\begin{equation}
z(k) = h_{sur} - h_{0} k - h_{1} log [cosh (( k - h_{th}) h_{cr})]
\end{equation}$$
where the coefficients \(h_{sur}\), \(h_0\), \(h_1\), \(h_{th}\) and \(h_{cr}\) are the parameters to be specified.
\(h_{cr}\) represents the stretching factor of the grid and
\(h_{th}\) is the approximate model level at which maximum stretching occurs.
This expression enables stretched z-coordinate vertical levels to be defined,
which are smoothly distributed along the water column, with appropriate thinning designed to better resolve the surface and intermediate layers.
Through partial cell parameterization, the thickness of the bottom layer can vary as a function of the geographical
location \((x,y)_{i,j}\) which allows a better representation of the real bathymetry.
(A) Horizontal grid.(B) Vertical T-grid.
Example of horizontal (left) and vertical (right) numerical grid with, respectively, grid sizes of \(\Delta\lambda\) and \(\Delta\varphi\) in horizontal and \(\Delta z\) in vertical direction.
Regridding, also called remapping, is the process of changing the grid (from a source to a destination grid)
underlying the field data values while preserving the qualities of the original data. In this section we describe the spatial extrapolation and interpolation procedure used in in SURF to remap the input fields on the child grid.
This phase will generate the surface forcing, initial and open lateral boundary condition datasets on the child grid. The procedures in this phase are:
Accessing and downloading the input datasets
Rotation of the vector fields (if required)
Extrapolation of the input datasets
Spatial interpolation of the source dataset on the child grid
Lateral Open Boundary Condition dataset generation
A procedure for checking if the necessary input datasets are in the experiment directory $PATH_IDEXP/data/ is executed. If any of the requested data are not present then the procedures downlAtmSrc, downlOceICSrc and downlOceBCSrc are automatically executed to respectively download the atmospheric forcing and the initial and lateral boundary condition datasets for the selected period of the simulation from remote or local data repositories, as specified in the configuration file setParFree.json.
When the parent coarse resolution model is defined on a rotated or a curvilinear grid
(e.g., the global tripolar grid, Fig. 2.6(a))
one more step is required to interpolate the horizontal velocity fields on the child grid.
In an ocean model with a “distorted” grid, the velocity vectors are obtained according to the direction of the grid lines. In a staggered Arakawa C grid system, the velocity field components are defined at the cell edges (the gray arrows in Fig. 2.6(b)).
A rotation in the latitudinal and longitudinal directions of the velocity components must be applied to change the vectors from the local system \((x,y)\) to a geographical system \((x^{'},y^{'})\), so that U gives the zonal component (W-E direction) and V the meridional component (S-N
direction) of the velocity vector. Therefore, to transform the \((x,y)\) coordinates to the \((x^{'},y^{'})\) coordinates, the vectors must be rotated according to vectors need to be rotated according to
$$
\begin{equation}
\begin{array}{ll}
U^{'}(x_{t}^{'},y_{t}^{'}) = U(x_{t},y_{t}) * cos(\alpha_{t}) - V(x_{t},y_{t})*sin(\alpha_{t})
\\
V^{'}(x_{t}^{'},y_{t}^{'}) = U(x_{t},y_{t}) * sin(\alpha_{t}) + V(x_{t},y_{t})*cos(\alpha_{t})
\end{array}
\end{equation}$$
For parent models with rotated rectangle grids, the angle is almost constant.
For those with curvilinear tripolar grids (as shown in Fig. 2.6(a)) the angle will vary within each grid cell.
(A)(B)
Example of curvilinear grid. Panel A gives an example of a tripolar grid.
Panel B gives the horizontal velocity components defined on the source curvilinear grid (black arrows) and on the destination rectilinear lat/lon grid (red arrows) after the rotation.
The extrapolation procedure used in SURF is referred to as the sea-over-land (SOL) procedure, and provides the ocean field values for areas near the coastline where the parent model solutions are not defined. The SOL procedure iteratively extrapolates the ocean quantities on the land grid-points, so they can be interpolated on the child grid. This also applies to several atmospheric fields and takes into account the atmospheric land-sea mask, to avoid land contaminations near the land-sea boundaries.
The SOL procedure is applied to the coarse resolution ocean fields to extrapolate the salinity, temperature, sea surface height and current fields to the land points. The ocean fields in the nested-grid points near to the coast can then be defined (by interpolation). The procedure applied for the atmospheric forcing fields takes into account the atmospheric land-sea mask to ensure there is no contamination from the land atmospheric fields to the sea points. The source code is provided in the directory $PATH_SURFNEMO/utilities/extrapol/seaoverland.
The extrapolation procedure described in the previous section provides the input data for the interpolator. The procedures are based on the Spherical Coordinate Remapping and Interpolation Package (SCRIP)
code. SCRIP is a software package that computes addresses and weights for remapping and interpolating fields between grids in spherical coordinates. The package should work for any grid on the surface of a sphere. SCRIP currently supports five remapping options:
Conservative remapping: First- and second-order conservative remapping, as described by Jones (1999, Monthly Weather Review, 127, 2204-2210).
Bilinear interpolation: Slightly generalized to use a local bilinear approximation (only logically rectangular grids).
Distance-weighted averaging: The inverse-distance-weighted average of a user-specified number of nearest neighbour values.
Particle remapping: A conservative particle (Monte-Carlo-like) remapping scheme
The source code can be found in the directory $PATH_SURFNEMO/nemo/NEMOGCM/TOOLS/WEIGHTS.
Regridding can be separated into two stages. The first stage is the generation of an interpolation weight matrix that describes how points in the source grid contribute to points in the destination grid. The second stage is the multiplication of values on the source grid by the interpolation weight matrix to produce the appropriate values on the destination grid.
The SCRIP spatial interpolation procedure is applied to the input fields for the generated bathymetry, atmospheric forcing and initial and boundary conditions files that are necessary to run the NEMO code. The generated files are stored in the directory $PATH_EXP/IDEXP/data/.
The lateral open boundary condition for the selected nested-domain is implemented using the BDY module of NEMO. Two numerical algorithms are used to treat open boundary conditions depending on the prognostic simulated variables. The Flather scheme (Oddo and Pinardi, 2008) is used for barotropic velocities, while the flow relaxation scheme (Engerdahl, 1995) is considered for baroclinic velocities, active tracers and sea surface height.
In our formulation, we provide external data along straight open boundary lines, and the relaxation area is equal to one internal grid point. As the parent coarse resolution ocean model provides only the total velocity field, the interpolated total velocity field in the child grid is split into barotropic and baroclinic components. An integral constraint method is imposed to preserve the total transport after the interpolation.
This process involves the following steps: (1) defining the open boundary geometry (for each of the T, U and V grids) and physical fields (active tracers, sea-surface height, barotropic and baroclinic velocities) at the open boundary points using the geometry_bdy and fields_bdy procedures, respectively; (2) writing these data arrays to the files that are necessary to run the NEMO code. The algorithms used for the different fields are the Flather radiation scheme for the barotropic velocities and the sea surface height and the Flow relaxation scheme for the baroclinic velocities and active tracers.
The generated files are stored in the directory $PATH_EXP/IDEXP/data/.
The downscaling is designed to ensure that the volume transport across the open boundary (OB) of the child model matches that across the corresponding section of the parent model. At the eastern/western boundaries U-Points are imposed using the following conditions
$$
\begin{equation}
\begin{array}{ll}
\int_{y_2}^{y_1} \int_{-H_{child}}^{\eta_{child}} U_{child} dz dy = \int_{y_2}^{y_1} \int_{-H_{parent}}^{\eta_{parent}} U_{parent} dz dy
\end{array}
\end{equation}$$
where \(y_1, y_2\) are the extremes of the open boundary section;
\(\eta_{child}, H_{child}\) are the surface elevation and the bathymetry of the child model at the boundary, respectively;
\(\eta_{parent}, H_{parent}\) are the surface elevation and the bathymetry of the parent model at the boundary, respectively; and
\(U_{parent},U_{child}\) are the parent/child total zonal velocities (normal velocity to the W/E boundaries).
The corrected velocity component normal to the boundary \(V_{child}\) is given
(see N. Pinardi et al., 2003) by:
$$
\begin{equation}
\begin{array}{ll}
U_{child} (x,y,z,t) = U_{interp} - U_{correction}
\end{array}
\end{equation}$$
where \(U_{interp}\) is the \(U_{parent}\)
interpolated on the child open boundary points and the velocity correction is given by
$$
\begin{equation}
\begin{array}{ll}
U_{correction} = \frac{M_{interp} - M_{parent}}{S}
\end{array}
\end{equation}$$
where \( M_{interp} = \int_{y_2}^{y_1} \int_{-H_{child}}^{\eta_{child}} U_{interp} dz dy \) is the volume transport across the OB,
the \( M_{parent} = \int_{y_2}^{y_1} \int_{-H_{parent}}^{\eta_{parent}} U_{parent} dz dy \) is the volume transport across the corresponding OB
and \( S = \int_{y_2}^{y_1} \int_{-H_{child}}^{\eta_{child}} dz dy \) is the area of the section.
These conditions are similarly imposed for the meridional velocity at the northern/southern boundaries (V-Points).
The Integral Constraint procedure ensures that the interpolation does not modify the net transport across the child model lateral open boundary.
Finally, the NEMO code is numerically integrated by the SURF platform. The output files are continuously updated, given a fixed output frequency, throughout the execution of the program. The logfile is a text file NEMO produces that contains data on how far the run has advanced by providing the time step. After the model run is completed, the output files are stored in the experiment directory $PATH_IDEXP/data/data00/outdata/ (see appendix B).
In this step the visualization and analysis procedures of the forecast are considered.
These can be activated after the execution of each macro-task. The user can visualize the input, regridded and output datasets, compare parent/child fields and convert the model output datasets.
