:math:` `
The TRACMASS documentation starts with the basic theory of the analytical trajectory solutions and the Lagrangian stream functions. Followed by a description of how the code is constructed with it corresponding subroutines. There is at the end a short manual of how to set up the TRACMASS code. This documentation is far from complete but is continuously upgraded in order to help the users all over the world.
The code was originally written in Fortran 77 for the FRAM ocean model at the Institute of Oceanographic Sciences, Deacon Laboratory (IOSDL) in Wormley, UK in the early 90’s. The name TRACMASS comes from the EU project with the same name, where it served together with the similar trajectory code Ariane by The present code is written in Fortran 95 and should be possible to apply to most GCMs based on finite differences. The TRACMASS code is continuously upgraded and adapted. The code can be downloaded from http://doos.misu.su.se/tracmass/.
There is, however, no user support for the code but all constructive comments or advice from the users are welcome. Furthermore the user must be familiar, in order to be able to use the TRACMASS code, with 1) the equation of motion for the ocean-atmosphere circulation, 2) the finite differences of these equations, 3) The TRACMASS theory, 4) Unix and 5) Fortran. There is, however, no user support for the code but all constructive comments or advice from the users are welcome. This manual has been developed with the support from the BONUS project BalticWay. Stockholm University, * *Kristofer Döös, Bror Jönsson, Joakim Kjellsson and Hanna Corell
The TRACMASS model calculates Lagrangian trajectories from Eulerian velocity fields, which have been simulated by ocean or atmosphere general circulation models (GCM). The trajectories are calculated off-line, i.e. after the GCM has been integrated and the velocity fields have been stored. This makes it possible to calculate many more trajectories than would be possible on-line (i.e. simultaneously with the GCM run). TRACMASS has been applied to many different general circulation models, both for the ocean and the atmosphere. The original feature of the method is that it solves the trajectory path through each grid cell with an analytical solution of a differential equation which depends on the velocities on the grid-box walls. The scheme was originally developed by and for stationary velocity fields and hereafter further developed by for time-dependent fields by solving a linear interpolation of the velocity field both in time and in space over each grid box, this in contrast to the Runge-Kutta methods where the trajectories are iterated forward in time with short time steps. The TRACMASS code has been further developed over the years and used in many studies in the atmosphere and ocean for global large scales studies and regional ones for the and Mediterranean and Baltic Seas .
[h]
[h]
image
[h]
image
[h]
This section is here only for pedagogical reasons, since it is only valid for rectangular grids. The TRACMASS code is written on a more general way in order to enable non-rectangular grids and is presented in the next section.
The velocity on a C-grid box (Fig. [fig:BCgrid]) is given by u_{i,j,k}, where i,j,k denote the discretised longitude, latitude and model level, respectively; u is the zonal velocity. Meridional and vertical velocities are defined analogously. It is possible to define the velocity inside a grid box by interpolating linearly between the discretised velocity values of the opposite walls. For the zonal direction, for example one obtains
\label{eq:ux}
u(x) = u_{i-1,j,k} + \frac{(x-x_{i-1})}{ \Delta x} (u_{i,j,k}-u_{i-1,j,k})
Local zonal velocity and position are related by u=dx/dt. The approximation in equation ([eq:ux]) can now be written in terms of the following differential equation:
\label{eq:dxdt}
\frac{dx}{dt} + \beta \, x + \delta = 0
with \beta \equiv \left(u_{i-1,j,k} - u_{i,j,k} \right)/ \Delta x delta equiv - u_{i-1,j,k} - beta , x_{i-1}`. Using the initial condition x(t_0)=x_0 , the zonal displacement of the trajectory inside the considered grid box can be solved analytically and is given by
\label{eq:xt}
x(t) = \left(x_0 + \frac{\delta}{\beta} \right) e^{- \beta (t-t_0)} - \frac{\delta}{\beta}
The time t_1 when the trajectory reaches a zonal wall can be determined explicitly:
\label{eq:t1}
t_1 = t_0 - \frac{1}{\beta} log \left[ \frac{x_1+\delta / \beta}{x_0+\delta / \beta} \right]
where x_1=x(t_1) is given by either x_{i-1} or x_i. For a trajectory reaching the wall x=x_i, for instance, the velocity u_i must necessarily be positive, so in order for equation ([eq:t1]) to have a solution, the velocity u_{i-1} must then be positive also. If this is not the case, then the trajectory either reaches the other wall at x_{i-1} or the signs of the transports are such that there is a zero zonal transport somewhere inside the grid box that is reached exponentially slow. For the meridional and vertical directions, similar calculations of t_1 are performed determining the meridional and vertical displacements of the trajectory, respectively, inside the considered grid box. The smallest transit time t_1-t_0 and the corresponding x_1 denote at which wall of the grid box the trajectory will exit and move into the adjacent one. The exact displacements in the other two directions are then computed using the smallest t_1 in the corresponding equation ([eq:xt]). The resulting trajectory through the grid box is illustrated by Fig. [fig:grid:sub:xytraj] and Fig. [fig:grid:sub:xztraj].
The disadvantage with the scheme presented in previous section is that it requires rectangular grid cells and GCMs generally use some sort of spherical grids as in the case of Fig. [fig:conbelt], where two spherical grids have been used for the world ocean.
The longitudinal size (\Delta x_{i,j}) on the northern and southern walls will typically be different on a spherical grid and on a curvilinear grid both \Delta x_{i,j} and \Delta y_{i,j} (latitudinal size) can be a function of its horizontal position.
The Equations ([eq:ux]-[eq:t1]) are not valid for a non-rectangular grid. This can however be solved by replacing the the velocities by volume transports. The transport through the eastern wall of the {i,j,k} grid box is given by
\label{eq:U}
U_{i,j,k} = u_{i,j,k} \Delta y_{i,j} \Delta z_k
The distance is non-dimensionalised by using r=x/\Delta x, and the linear interpolation of the velocity (Eq. ([eq:ux])) is replaced by
\label{eq:Ur}
U(r) = U_{i-1,j,k} + (r-r_{i-1})(U_{i,j,k}-U_{i-1,j,k})
The local transport and position are now related by U=dr/ds, where the scaled time variable s \equiv t/(\Delta x_{i,j} \Delta y_{i,j} \Delta z_k), where the denominator is the volume of the particular grid box. The differential equation ([eq:dxdt]) is replaced by
\label{eq:drds}
\frac{dr}{ds} + \beta \, r + \delta = 0
with \beta \equiv U_{i-1,j,k} - U_{i,j,k} and \delta \equiv - U_{i-1,j,k} - \beta \, r_{i-1}. Using the initial condition r(s_0)=r_0 , the zonal displacement of the trajectory is now given by
\label{eq:rs}
r(s) = \left(r_0 + \frac{\delta}{\beta} \right) e^{- \beta (s-s_0)} - \frac{\delta}{\beta}
The scaled time s_1 becomes
\label{eq:s1}
s_1 = s_0 - \frac{1}{\beta} log \left[ \frac{r_1+\delta / \beta}{r_0+\delta / \beta} \right]
where r_1=r(s_1) is given by either r_{i-1} or r_i. With the use of equation ([eq:U]), the logarithmic factor can be expressed as log[ U(r_1)/U(r_0)]. For a trajectory reaching the wall r=r_i, for instance, the transport U(r_1) must necessarily be positive, so in order for equation ([eq:s1]) to have a solution, the transport U(r_0) must then be positive also. If this is not the case, then the trajectory either reaches the other wall at r_{i-1} or the signs of the transports are such that there is a zero zonal transport somewhere inside the grid box that is reached exponentially slow.
The meridional solution is calculated similarly as the zonal one but using the meridional transport defined as
\label{eq:W}
V_{i,j,k} = v_{i,j,k}\Delta x \Delta z_k\\
The vertical solution is calculated similarly as the zonal one but using the meridional transport defined as
\label{eq:V}
W_{i,j,k} = w_{i,j,k}\Delta x \Delta y
The scheme is mass conserving since the vertical transport is directly calculated from the continuity equation in the same way as in the ocean GCM, which is due to the incompressibility in the ocean
\label{eq:contz}
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0
that is discretised with finite differences on a C-grid into
\label{eq:contzdisc}
\frac{u_{i,j,k} - u_{i-1,j,k}}{\Delta x_{i,j}} + \frac{v_{i,j,k} - v_{i,j-1,k}}{\Delta y_{i,j}} + \frac{w_{i,j,k} - w_{i,j,k-1}}{\Delta z_k} = 0
Equation ([eq:contzdisc]) can be explained by that the sum of all the volume fluxes in or out of the grid box is zero. We rewrite in flux form so that the continuity equation becomes
\label{eq:cont}
W_{i,j,k}=W_{i,j,k-1} -(U_{i,j,k}-U_{i-1,j,k}+V_{i,j,k}-V_{i,j-1,k})
and is integrated from the bottom and upwards with the boundary condition W=0. Since the trajectory solutions are exact and that the continuity equation is respected the TRACMASS trajectories will therefore never hit the coast or the sea floor.
The calculations of s_1 are performed determining the zonal, meridional and vertical displacements of the trajectory, respectively, inside the considered grid box. The smallest transit time s_1-s_0 and the corresponding r_1 denote at which wall of the grid box the trajectory will exit and move into the adjacent one. The exact displacements in the other two directions are then computed using the smallest s_1 in the corresponding equation ([eq:rs]).
The trajectories can be both calculated forward and backward in time and hence it is possible to trace origins of water or air masses. The differential equation ([eq:dxdt]) is strictly only valid for stationary velocity fields. developed an extension of the TRACMASS code for time dependent velocities. It is however possible to use the present stationary code with negligible loss of accuracy by simply changing the velocity and pressure fields at regular time intervals, which are the same as the output data from the used GCM is stored.
[!h]
image
The atmospheric version of TRACMASS uses however the conservation of mass instead of volume. Most Atmospheric GCMs today use terrain-following vertical coordinates. Following the atmosphere is divided into NLEV layers, which are defined by the pressures at the interfaces between them and these pressures are given by p_{k+1/2} = A_{k+1/2} + B_{k+1/2} \, p_{S} for k=0,1,..,NLEV, with k=0 at the top of the atmosphere and k=NLEV at the Earth’s surface. The A_{k+1/2} and B_{k+1/2} are constants, whose values effectively define the vertical coordinate and p_{S} is the surface pressure. The dependent variables, which are the zonal wind (u), the meridional wind (v), the temperature (T) and the specific humidity (q) are defined in the middle of the layers, where the pressure is defined by p_{k} = \frac{1}{2} (p_{k-1/2} + p_{k+1/2}), for k=1,2,..,NLEV. The vertical coordinate is \eta=\eta(p,p_S) and has the boundary value \eta(0,p_S)=0 at the top of the atmosphere and \eta(p_S,p_S)=1 at the Earth’s surface.
For the ocean, in the previous sections, we used volume transport because of the incompressibility approximation. In the atmosphere we need instead to use mass transport so Equation ([eq:U]) is now replaced by the zonal and meridional mass transports in the model layers:
\label{eq:Um}
U_{i,j,k} = u_{i,j,k} \frac{ \Delta y \, \Delta p_k}{g} \;;\; V_{i,j,k} = v_{i,j,k} \frac{ \Delta x_{j} \, \Delta p_k}{g}
where :math:` Delta p_k = Delta A_k + Delta B_k p_{S , i,j,k} , :math: Delta A_k = A_{k+1/2} - A_{k-1/2} ` and :math:` Delta B_k = B_{k+1/2} - B_{k-1/2} ` .
The vertical mass transport through the model layer interfaces W is also required by the trajectory caclulations but cannot be as simply solved as for the ocean with its continuity equation. The mass of the grid box is, due to hydrostaticity,
\label{eq:M}
M_{i,j,k} = \frac{\Delta p_{i,j,k}}{g} \Delta x_{i,j} \Delta y_{i,j} = \frac{ \Delta A_k + \Delta B_k p_{S \, i,j}}{g} \Delta x_{i,j} \Delta y_{i,j}
The mass conservation of a grid box as illustrated in Fig. [fig:atm-grid] yields d M_{i,j,k}/d t = 0. The rate of mass change is hence equal to the mass transports through the grid box’s walls:
\label{eq:Masscons}
\frac{\partial M_{i,j,k}}{\partial t} = - \left( U_{i,j,k} - U_{i-1,j,k} + V_{i,j,k} - V_{i,j-1,k} +W_{i,j,k} - W_{i,j,k-1} \right)
Eleminating M_{i,j,k} between Eqs. ([eq:M])-([eq:Masscons]):
\label{eq:M2}
\frac{\partial p_{S \, i,j}}{\partial t} \frac{\Delta B_k \Delta x_{i,j} \Delta y_{i,j} }{g} = - \left( U_{i,j,k} - U_{i-1,j,k} + V_{i,j,k} - V_{i,j-1,k} +W_{i,j,k} - W_{i,j,k-1} \right)
The rate of change of surface pressure ({\partial p_{S}}/{\partial t}) can be obtained by integrating Eq. ([eq:M2]) over the entire air column from the Earth’s surface to the top of the atmosphere so that W is eliminated by using the boundary conditions W_{k=0}=W_{k=NLEV}=0 so that
\label{eq:intvert}
\frac{\partial p_{S \, i,j}}{\partial t} = - \frac{g }{ \Delta x_{i,j} \Delta y_{i,j} } \sum_{k=1}^{NLEV} \left( U_{i,j,k} - U_{i-1,j,k} + V_{i,j,k} - V_{i,j-1,k} \right)
Eq. [eq:M2] can now be succesively integrated upwards form the surface with the boundary condition W_{k=0}=0 to each k-level so that
\label{eq:W}
W_{i,j,k} = W_{i,j,k-1} - U_{i,j,k} + U_{i-1,j,k} - V_{i,j,k} + V_{i,j-1,k} - \frac{\partial p_{S \, i,j}}{\partial t} \frac{\Delta B_k \Delta x_{i,j} \Delta y_{i,j} }{g}
The trajectory differential Eq. ([eq:drds]) and its solutions Eqs. ([eq:rs] - [eq:s1]) remain the same but now feeded with mass transports on atmospheric terrain-following vertical coordinates and the scaled time is now s \equiv t g /(\Delta x_{i,j} \Delta y_{i,j} \Delta p_{k} ).
[sec:timesteping]
The schemes in the previous sections are strictly only valid for steady state velocity fields. The scheme can however be used with time dependent fields by assuming that the velocity fields are in steady state during a chosen time step dtmin.
In the present section, we will present a truly time dependent scheme which is solved analytically, which was introduced by . This scheme is however not as robust as the stationary one.
Given a set of velocities \mathbf{V}_n for each model point, where n represents a discretised time, a bi-linear interpolation of transports in space as well as in time leads to
\begin{gathered}
\label{eq:Frs}
F(r,s) = F_{i-1,n-1} + (r-r_{i-1})(F_{i,n-1}-F_{i-1,n-1}) + \\
\frac{s-s_{n-1}}{\Delta s} \left[ F_{i-1,n} - F_{i-1,n-1} + (r-r_{i-1})(F_{i,n}-F_{i-1,n} - F_{i,n-1} + F_{i-1,n-1}) \right]\end{gathered}
which is the genral expression for the three directions where i signifies either a longitudinal, meridional, or vertical direction. The transport F=(U,V,W) and as before r=(x,y,z), r=(x/\Delta x, y/\Delta y, z/\Delta z) , where the denominator is the volume of the particular grid box. \Delta s is the scaled time step between two data sets
\label{eq:ds}
\Delta s = s_n-s_{n-1} = (t_n-t_{n-1})/(\Delta x \Delta y \Delta z) = \Delta t /(\Delta x \Delta y \Delta z)
where dt the time step between two data sets in true time dimension (seconds).
Connecting the local transport to the time derivative of the position with F=dr/ds, we get the differential equation
\label{eq:drds_time}
\frac{dr}{ds} + \alpha \, r \, s + \beta \, r + \gamma \, s + \delta= 0
where the coefficeints are defined by
\label{eq:alpha}
\alpha \equiv - \frac{1}{\Delta s}\, (F_{i,n}-F_{i-1,n} - F_{i,n-1} + F_{i-1,n-1})
\label{eq:beta}
\beta \equiv F_{i-1,n-1} -F_{i,n-1} -\alpha \, s_{n-1}
\label{eq:gamma}
\gamma \equiv - \frac{1}{\Delta s} \, (F_{i-1,n}- F_{i-1,n-1}) - \, \alpha \, r_{i-1}
\label{eq:delta}
\delta \equiv F_{i-1,n-1} + r_{i-1}(F_{i,n-1}-F_{i-1,n-1}) -\gamma \, s_{n-1}
Analytical solutions can be obtained for the following three cases: \alpha > 0, \alpha < 0 and \alpha = 0, Note that inside the grid box, the acceleration, d^2r/ds^2 = -\alpha r - \gamma dependent term proportional to \alpha. For \alpha > 0, the latter term implies a varying deceleration across the grid box. If \alpha > 0, we define the timelike variable \xi = (\beta + \alpha \, s)/\sqrt{2 \alpha} and get
\label{eq:rs_time}
r(s) = \left(r_0 + \frac{\gamma}{\alpha} \right) e^{\xi^2_0- \xi^2} - \frac{\gamma}{\alpha} + \frac{\beta \gamma -\alpha \delta}{\alpha} \sqrt{\frac{2}{\alpha}} \left[ D(\xi) - e^{\xi^2_0- \xi^2} D(\xi_0) \right]
using Dawson’s integral
\label{eq:dawson}
D(\xi) \equiv e^{- \xi^2} \int_0^\xi e^{x^2} dx
and the initial condition r(s_0)=r_0. If \alpha < 0, \xi becomes imaginary. By defining \zeta \equiv i \xi = (\beta + \alpha s)/\sqrt{-2 \alpha}, Eq. [eq:rs:sub:time] can be re-expressed as
\label{eq:rs_time2}
r(s) = \left(r_0 + \frac{\gamma}{\alpha} \right) e^{\zeta^2- \zeta^2_0} - \frac{\gamma}{\alpha} -
\frac{\beta \gamma -\alpha \delta}{\alpha} \sqrt{\frac{\pi}{-2\alpha}} \, e^{\zeta^2} \left[ erf(\zeta) - erf(\zeta_0) \right],
where erf(\zeta)=(2/\sqrt(\pi)\int_0^\zeta e^{-x^2} dx. The case \alpha = 0 will occur occasionally in practice. Its corresponding solution of The differential Eq. ([eq:drds:sub:time]) is
\label{eq:rs_time3}
r(s) = \left(r_0 + \frac{\delta}{\beta} \right) e^{-\beta(s-s_0)} - \frac{\delta}{\beta} +
\frac{ \gamma }{\beta^2} \left[ 1 - \beta s +(\beta s_0-1) e^{-\beta(s-s_0)} \right]
Compare this with the solution ([eq:rs]) for time-independent velocity fields. A major difference compared to the time-independent case is that now the transit times s_1 - s_0 cannot in general be obtained explicitly. Using the solutions ([eq:rs:sub:time])-([eq:rstime3]), the relevant root s_1 of
\label{eq:rs1}
r(s_1) -r_1 = 0
has to be computed numerically for each direction. In the following subsection, we describe how the roots s_1 and the corresponding exiting wall r_1 can be determined. The displacement of the trajectory inside the considered grid box then proceeds as discussed previously for stationary velocity fields.
Here we consider how to determine the roots s1 of Equation ([eq:rs1]) and the corresponding r_1 needed to compute trajectories inside a grid box. In the following, s_{n-1} \leqslant s_0 < s_n and the relevant roots s_1 are to obey s_0 < s_1 \leqslant s_n . We also focus on the cases \alpha > 0 and \alpha < 0, since the considerations below can easily be adapted for \alpha = 0. For numerical purposes, we use
\label{eq:betagamma}
\frac{\beta \gamma -\alpha \delta}{\alpha} = \frac{ F_{i,n} F_{i-1,n-1} - F_{i,n-1} F_{i-1,n} }{F_{i,n}-F_{i-1,n} - F_{i,n-1} + F_{i-1,n-1}}
\label{eq:gammalpha}
\frac{\gamma}{\alpha} = \frac{ F_{i-1,n} - F_{i-1,n-1} }{F_{i,n}-F_{i-1,n} - F_{i,n-1} + F_{i-1,n-1}} - r_{i-1}
\label{eq:xi}
\xi = \frac{ F_{i-1,n-1} - F_{i,n-1} + \alpha (s-s_{n-1})}{ \sqrt{2\alpha}}
\label{eq:zet}
\zeta = \frac{ F_{i-1,n-1} - F_{i,n-1} + \alpha (s-s_{n-1})}{ \sqrt{-2\alpha}}
The coefficients [eq:betagamma] appearing in the solutions [eq:rs:sub:time] and [eq:rs:sub:time2] is exactly zero when either the r_{i-1} or r_i wall represents land, the transports F_i or F_{i-1} being zero for all n, respectively. In these instances, the opposite wall fixes r_1 , and the root s_1 > s_0 can then be computed analytically. If there is no solution, we take s_1 = s_n. When all three transit times equal s_n, the trajectory will not move into an adjacent grid box but will remain inside the original one. Its new position is subsequently computed, and the next time inter val is considered.
If :math:` (beta gamma -alpha delta ) / alpha neq 0`, the computation of the roots of Eq. [eq:rs1] can only be done numerically. This is also true for locating the extrema of the solutions [eq:rs:sub:time] and [eq:rs:sub:time2]. Alternatively, one can consider F(r, s) = 0 using Eq. [eq:Frs] to analyze where possible extrema are located. It follows that in the sr plane, extrema lie on a hyperbola of the form r = (as + b)/(c + ds). Of course, only the parts defined by :math:` s_{n-1} leq s leq s_n` and :math:` r_{i-1} leq r leq r_i` are relevant. Depending on which parts of the hyperbola, if any, lie in this "box" and on the initial condition r(s_0) = r_0 , the trajectory r(s) exhibits none, one, or at most two extrema. In the latter case, the trajectory will not cross either the wall at r_{i-1} or the one at r_i (see Fig. [fig:time:sub:analyt] for an example). Hence, those trajectories r(s) determining the transit time s_1 - s0 will have at most one extremum, that is, there is at most one change of sign in the local transport F.
An efficient way to proceed then is as follows. First, consider the wall at r_i. For a trajectory to reach this wall, the local transport must be nonnegative, which depends on the signs of the transports F_{i-1,n} and F_{i,n}. Four distinct configurations may arise:
- F(r_i,s) > 0 for :math:` s_{n-1} < s < s_n`
- sign of F(r_i,s) changes from positive to negative at :math:` s = tilde{s} < s_n`
- sign of F(r_i,s) changes from negative to positive at :math:` s = hat{s} < s_n`
- F(r_i,s) < 0 for :math:` s_{n-1} < s < s_n`
For case 1, evaluate r(s_n) using the appropriate analytical solution. If r(s_n) \geq r_i, the trajectory has crossed the grid-box wall for s_1 \leq s_n. If the initial transport F(r_0,s_0) < 0, the trajector y may have crossed the opposite wall at an earlier time. The latter is only possible if case 3 applies for the wall at r i_1 and \hat{s} > 0, in which case one checks whether r(\hat{s}) \leq r_{i-1}. If this is not so, then there is a solution to r(s_1) - r_1 = 0 for r_1 = r_i and s_0 < s_1 \leq s_n. Subsequently, this root can be calculated numerically using a root-solving algorithm. On the other hand, if r(s_n) < r_i or, if applicable, r(\hat{s} ) \leq r_{i-1} continue with considering the other wall. The arguments for the wall at r_{i-1} are similar to those relating to r.
If case 2 applies and s_0 < \tilde{s}, follow the considerations given for case 1 using \tilde{s} instead of s_n. If there is a root for r_1 = r_i , then :math:` s_0 < s_1 leq hat{s} ` .
For case 3, follow the considerations given for case 1. If there is a root for r_1=r_i, then :math:` hat{s} < s_1 leq s_n` .
For case 4, no solution of Eq. [eq:rs1] is possible for r_1 = r_i. Turn attention to the other wall.
All these considerations is applied to each direction.
The trajectory solutions in the previous sections only include the implicit large scale diffusion due to along-trajectory changes of temperature and salinity/humidity, and by the models parameterisation of turbulent mixinging in the momentum equations. These trajectories do not, however, explicitly represent sub-gridscale turbulence. There are two ways to incorporate a sub-grid scale turbulence in TRACMASS.
There is, however, no general consesus wether it is physical to include these sub-grid parameterisations since we are adding scales that are not present in the GCM.
This scheme, which was introduced by , adds a random horizontal turbulent velocity u', v' to the horizontal velocity U, V from the GCM and is illustrated by Fig. [fig:turbu]. This is added to each trajectory and each horizontal grid wall every time step. This enabled us to use the TRACMASS code as it is but with a velocity field that is somewhat shaken, resulting in stirred trajectory particles. The new velocity from which the transport is calculated in equation ([eq:U]) is now u=U+u'. The amplitude of the random turbulent velocity is set to the same size as the circulation model velocity U, so that u'=RU, where R is a random number uniformly distributed between -a and a. The best results were obtained for a=1 in the study by .
The effect of this superimposed sub-grid turbulence is clearly visible in Fig. [fig:trajturb], where a particle cluster is traced with and without this sub-grid parameterisation. The turbulence smoothes the trajectory positions and spreads them more evenly. The stirred particles in Fig. [fig:trajturb-b] fill visibly regions where no particles were present without sub-grid turbulence in Fig. [fig:trajturb-a].
This sub-grid parameterisation was tested in and will however need more tests and evaluation in the future.
[h]
image
[h]
image
[h]
This subroutine adds a random displacement to the trajectory position in order to incorporate a sub-grid parameterisation of the non resolved scales. The scheme was introduced in TRACMASS by .
The horizontal Laplacian diffusion equation, which is
\label{eq:laplace}
\frac{\partial P}{\partial t} = A_H \left( \frac{\partial^{2}P}{\partial x^{2}} + \frac{\partial^{2}P}{\partial y^{2}} \right)
is replaced by a random walk using a Gaussian distribution
\label{eq:2dgaus}
P(x_d,y_d, \Delta t) = \frac{1}{4 \pi A_H \Delta t} \exp \left( \frac{-x_d^2 -y_d^2}{4 A_H \Delta t} \right)
Figure [fig:diffusion] illustrates the added displacements to the original position of the particle after each time-step of length \Delta t, which is:
\label{eq:diff_ disp_x}
x_d = \sqrt{ - \, 4 \, A_H \, \Delta t \, \log(1-q_1) } \, \cos{(2 \, \pi \, q_2)}
\label{eq:diff_ disp_y}
y_d = \sqrt{ - \, 4 \, A_H \, \Delta t \, \log(1-q_1) } \, \sin{(2 \, \pi \, q_2)}
\label{eq:diff_ disp_z}
z_d = \sqrt{ - \, 4 \, A_V \, \Delta t \, \log(1-q_3) } \, \cos{(2 \, \pi \, q_4)}
where A_H and A_V are the horizontal and vertical eddy viscosity coefficients. They can be chosen to be the same as those used in the momentum equation by the circulation model so the random motions are on the same scale of the model subgrid-scale turbulence. q_1, q_2, q_3 and q_4 are random numbers between 0 and 1. The added displacement in the horizontal plane will hence lye on a distance from the original position of
\label{eq:diff_ disp_z}
r_H = \sqrt{ x_d^2 \, + y_d^2 } \, = \, \sqrt{ 4 \, \Delta t \, A_H \, \log(1-q_1)}
The added horizontal displacement will therefore lye on a random position on a Gaussian distribution where one standard deviation is
\label{eq:diff_ disp_z}
R_H = \, \sqrt{ 4 \, \Delta t \, A_H }
The horizontal diffusion in the previous section was isotropic and the added random displacement lied on a circular disk. In this section we introduce an anistropic diffusion that is higher along the isobaths and weaker perpendicular to the isobaths. The random displacement will therefore be based on an elliptic disk instead of a circular disk as illustrated by Fig. [fig:diffusion:sub:ellipse]. This diffusion will hence take into account the fact that observations of trajectories of floats and drifters tend to follow isolines of constant planetary vorticity (f/h).
The implementation of this sub-grid parameterisation is the following. The slope, i.e. the maximum depth gradient is first calculated
\label{eq:gradient}
G = \sqrt{ \left( \frac{\Delta h}{ \Delta x}\right)^2 + \left( \frac{\Delta h}{ \Delta y}\right)^2 }
The circular disk from previous section is then stretched into an elliptic disk by
\label{eq:eliptrans}
x_d' = (c_1 G+1) x_d
\label{eq:eliptrans}
y_d' = (c_1 G+1) y_d
where c_1 is a constant. The ellipse will be a circular disk only if the bottom of the ocean is flat in the grid cells beneath the trajectory particle.
The angle from the x-axis (eastward) and the isopleths is
\label{eq:theta}
\theta = \arcsin{ \left( \frac{ \frac{\Delta h}{ \Delta x} } { G } \right) }
The ellipse is then inclined with the computed angle \theta from the x-axis (eastward) to be along the isopleth by the coordinate transform:
\label{eq:coordtrans}
x_d''= x_d' \cos(\theta)- y_d' \sin(\theta) \\
\label{eq:coordtrans}
y_d''=-x_d' \sin(\theta)+ y_d' \cos(\theta)
The Lagrangian dispersion has been shown to increase and be more realistic despite that the injected parameterised diffusion remains the same (in preparation, Döös, 2010).
The particles are taken to enter the model domain along the shoreline, and are from there transported by the motion of the water, at the same time being affected by sedimentation and resuspension processes. The horizontal motion is prescribed solely by the GCM field, while the vertical displacements are derived from the GCM field together with the settling velocity in Eq. [eq:stokes]. If a particle reaches the lower boundary of the deepest grid box in the water column, i.e. the sea floor, it will settle. Once settled it will stay at the settling position and can only leave it by resuspension. If no resuspension occurs, the particle will remain at its position until the simulation ends.
The concept underlying the sedimentation model is that suspended particulate matter is bound to follow the movements of the water. If the sinking motion of the particles in quiescent water is known, and the paths of the water can be calculated, then the particle trajectories will be a combination of these displacements.
There are two general modes of settling; viscous settling for particles smaller than 0.2 mm and inertial settling for particles larger than approximately 2 mm. Between these sizes there is a transition zone. In this study the particles modelled are silt and clay, i.e particles with a diameter much smaller than 0.2 mm, and the settling is thus viscous. In practice, the settling velocity of a particle has a basic relation to its size and shape. Since it is not possible to account for all different shapes a particle can have, the concept of equivalent size is used, viz. the size of a quartz sphere having the same settling velocity as a less spherical natural grain.
To the GCM-derived vertical velocity of the water, a settling velocity w_{s} for the particle is added. This is calculated on the basis of Stokes law, using the particle density \rho_{s}, diameter d, as well as water density and viscosity, \rho_{w} and \mu:
\label{eq:stokes}
w_{s} =\frac{\rho_{s}-\rho_{w}}{18\mu}gd^{2}.
This formula is valid when the Reynolds number is smaller than 1, corresponding to an equivalent diameter of 0.2 mm and smaller. Since the viscosity {\mu} is included in the equation, the viscous settling is temperature-dependent.
The particles are taken to enter the model domain along the shoreline, and are from there transported by the motion of the water, at the same time being affected by sedimentation and resuspension processes. The horizontal motion is prescribed solely by the GCM field, while the vertical displacements are derived from the GCM field together with the settling velocity in Eq. [eq:stokes]. If a particle reaches the lower boundary of the deepest grid box in the water column, i.e. the sea floor, it will settle. Once settled it will stay at the settling position and can only leave it by resuspension. If no resuspension occurs, the particle will remain at its position until the simulation ends.
Resuspension of a settled particle will take place if the shear stress at the bottom where the particle is located exceeds a threshold value. When this occurs, the particle will be lifted up a short distance above the bottom. There it will catch on to the water flow field again, and continue its motion in the water body.
The shear stress at the bottom is dependent on the turbulent kinetic energy. Since this is not included in the GCM-derived data set, the horizontal velocities at the lower boundary of the bottom grid box are used instead, based on the view that the water velocitie give rise to the kinetic energy that leads to the shear stress. It states the relation between the grain diameter in micrometres and the mean velocity 15 cm above the bottom in cm/s for silt and clay. For a water content of 100 % the critical velocity is 10 cm/s for the whole fraction, the relationship being valid for cohesive material of 0.1 mm and smaller.
[h]
image
[h]
image
[h]
image
[h]
image
A Lagrangian stream function can be calculated by summing over trajectories representing the desired path. A particular water mass can be isolated by following a set of trajectories between specific initial and final sections. Each trajectory, indexed by n, is associated with a volume transport T_n given by the velocity, initial area, and number of trajectories released. During transit from the initial to the final section the volume transport remains unchanged; the transport/velocity field is non-divergent, permit- ting stream-function representations. The volume transport linked to each trajectory is inversely proportional to the number of trajectories released, viz. the Lagrangian resolution (which should be sufficiently high to ensure that the stream function does not change when the number of trajectories is further increased). A non-divergent 3-D volume-transport field is obtained by recording every instance of a trajectory passing a grid-box wall. Every trajectory entering a grid-box also exits, and hence this field exactly satisfies
\label{eq:trans}
T_{i,j,k,n}^{x}-T_{i-1,j,k,n}^{x}+T_{i,j,k,n}^{y}-T_{i,j-1,k,n}^{y}+T_{i,j,k,n}^{z}-T_{i,j,k-1,n}^{z}=0
where T_{i,j,k,n}^{x}, T_{i,j,k,n}^{y} and T_{i,j,k,n}^{z}, are the trajectory-derived volume transports in the zonal (i), meridional (j), and vertical (k) directions, respectively.
By integrating vertically over the transports and over the trajectories one obtains the Lagrangian barotropic stream function \Psi_{i,j}^{LB}
\label{eq:psixy}
\Psi_{i,j}^{LB}-\Psi_{i-1,j}^{LB}={\displaystyle \sum_{k}}{\displaystyle \sum_{n}}T_{i,j,k,n}^{y}\quad or\quad\Psi_{i,j}^{LB}-\Psi_{i,j-1}^{LB}=-{\displaystyle \sum_{k}}{\displaystyle \sum_{n}}T_{i,j,k,n}^{x}
By instead integrating zonally one obtains the Lagrangian meridional overturning stream function \Psi_{j,k}^{LM}
\label{eq:psiyz}\cdot
\Psi_{j,k}^{LM}-\Psi_{j,k-1}^{LM}=-{\displaystyle \sum_{i}}{\displaystyle \sum_{n}}T_{i,j,k,n}^{y}\quad or\quad\Psi_{j,k}^{LM}-\Psi_{j-1,k}^{LM}={\displaystyle \sum_{i}}{\displaystyle \sum_{n}}T_{i,j,k,n}^{z}
Finally by integrating meridionally one obtains the Lagrangian zonal overturning stream function \Psi_{i,k}^{LZ}
\label{eq:psiyz}
\Psi_{i,k}^{LZ}-\Psi_{i,k-1}^{LZ}=-{\displaystyle \sum_{i}}{\displaystyle \sum_{n}}T_{i,j,k,n}^{y}\quad or\quad\Psi_{i,k}^{LZ}-\Psi_{i-1,k}^{LZ}={\displaystyle \sum_{i}}{\displaystyle \sum_{n}}T_{i,i,k,n}^{z}
An example of a zonal Lagrangian stream function is shown in Fig. [fig:gulf:sub:finland].
The vertical index k can be the level depth, temperature, salinity or density for ocean GCMs. For atmospheric GCMs it can be the model level, temperature, specific humidity, geopotential height or pressure.
A simulated tracer concentration can also be calculated by activating the key -Dtracer. It calculates the trajectory particle concentration in each grid at a given time. This is not a true tracer but is directly comparable with a tracer. An example of this is shown in Figure [fig:tracer] from .
The TRACMASS code is written in Fortran 95. It consists of a number of subroutines, which are preprocessed and compiled with two Makefiles. The entire unix file tree is shown in Figure [fig: unixfiletree]. The flow chart of TRACMASS is shown in Figure [fig:tracmass:sub:flowchart] and the flowchart of the subroutine loop in Figure [fig:tracmass:sub:loopflowchart].
There are two makefiles. This first one is under tracsvn/Makefile compiles TRACMASS with the chosen Fortran compiler and libraries (netcdf). It needs to be adapted to work with your computer platform. There are a few examples that work with MacOSX and different linux platforms. You will also need to set PROJECT to the the chosen "project", which is the GCM data sets you want to run TRACMASS with. There provided Makefile is set up to work with gfortran and g95 fortran and is named Makefile\_template. It needs to be adapted and renamed Makefile
Here are a few of the many possible projects that can be used and set by the very first line in this Makfile. PROJECT = orc NEMO with one of the ORCA grids. PROJECT = ifs The AGCM (IFS) from ECMWF, which is part of EC-Earth. PROJECT = tes Test basin with analytical time dependent velocity fields. PROJECT = baltix NEMO version of the Baltic and North sea.
These are the input files and subroutines that need to be adapted for each run and changed for each new project.
This input file defines the dimensions of the grid IMT = number of zonal grid points JMT = number of meridional grid points KM = number of vertical grid points NTRACMAX = maximum number of trajectories. TRACMASS will stop if this number is exceeded but will consume a lot of central memory if set to a high number. Should be able to be set to a few millions. ngcm = hours between GCM datasets iter = number of iterations between the two gcm data sets. This sets the number of calculations between the two GCM data sets explained in section [sec:timesteping], where dtmin=the time between two GCM data sets divided by iter. The trajectory solution becomes more accurate if set to a high number. 10 is however generally enough. intmax = maximum number of GCM fields
This input file sets many values of the actual run.
This subroutine reads the stored velocity, temperature and the salinity/humidity fields from the GCM. It is this subroutine that needs to be adapted when a new project is introduced.
This makefile is under tracsvn/projects/<project>/Makefile and sets the chosen pre-processing options.
-Dstat Stationary velocity fields (does not work at the moment) -Dtimestep Time changing velocity fields but using stationary solution to the differential equations with time step changing velocities -Dtimeanalyt Time changing velocity fields with time analytical solutions to the differential equations based on . This scheme is not as robust as -Dtimestep.
-Dstreamxy Calculates the Lagrangian barotropic stream function. -Dstreamv Calculates the Lagrangian vertical stream function as a function of model level -Dstreamr —————————————– " ——————————- temperature, salinity/humidity, density/pressure. The the latter two options have to be used with -Dtempsalt Calculates the temperature, salinity and density for the ocean and temperature, humidity and pressure for the atmosphere -Ddensity Calculates only the density along the trajectory. -Dtracer Stores trajectory particle positions as a simulated tracer
These options will make the code non-Lagrangian, i.e. the trajectories will not be passively advected by the velocity field as neutrally buoyant particles. -Dsediment Sediment code developed for RCO -Dtwodim Two-dimensional trajectories, which do not change depth. The vertical velocity is simply set to zero. -Dturb Adds parameterised sub-grid turbulent velocities u’ and v’ -Ddiffusion Adds a diffusive random position to the trajectory
-Drerun This option is used in order to stores the Lagrangian stream functions as a function of the end positions. The code is run a first time to find the end positions of the trajectories at chosen sections. The code is then rerun with -Drerun so that the Lagrangian stream functions will be divided into the different paths, which has been calculated in the previous trajectory run.
The following subroutines are always used by TRACMASS and consists of the pure Lagrangian trajectory calculated.
This is the program main where the code starts and ends and from which the subroutine loop is called.
The modules located in modules.f95 is where the dimensions are defined for the whole code.
This subroutine reads in the data from <project>_grid.in and <project>_run.in.
This subroutine reads the data for the initial positions of the trajectories. That is where they are "seeded".
This subroutine defines the horizontal and vertical grids as well as the pressure, density, temperature, salinity and humidity coordinates by specifying their minimum and maximum values. These need to be changed in order to fit the chosen model domain and study.
This is the central subroutine of TRACMASS. The flow chart of this subroutine is shown in Figure [fig:tracmass:sub:loopflowchart]. The subroutine consists of several loops over both time and the trajectories.
This subroutine computes the vertical velocities in the boxes of the water/air column where the trajectory particle is located by integrating Eq. [eq:cont].
This subroutine to computes the time (sp,sn) from Eq. [eq:s1], when the trajectory crosses the face of box (ia,ja,ka). Two crossings are considered for each direction: east and west for longitudinal directions, north and southward for latitudinal directions and up and down for vertical directions. The subroutine is called for each of the three directions from subroutine loop in order to determine which will be the shortest of the 6 times, which will determine which of the 6 grid walls the trajectory crosses.
The shortest of the times obtained by the subroutine cross is here used to calculate the new position of the trajectory with Eq. [eq:rs].
These subroutines are attempts to include parameterisations of the the non-resolved sub-grid scales. There are two two different options in the TRACMASS code. These subroutines are activated by using -Dturb or -Ddiffusion in the Makefile. Don’t use both of them!
It is important to realise that these scales are not included in the original velocity fields from the GCM. There is no general consensus wether it is physical to include these sub-grid parameterisations since we are adding scales that are not present in the GCM.
This subroutine computes a random horizontal turbulent velocity u', v', which are used in the subroutines vertvel, pos and cross in order to incorporate a sub-grid parameterisation of the non resolved scales. It is activated by the preprocessing option -Dturb in the Makile. Don’t use together with -Ddiffusion.
This subroutine adds a random position to the new trajectory position in order to incorporate a sub-grid parameterisation of the non resolved scales. It was originally written by Richard Levine and David Webb. It is activated by the preprocessing option -Ddiffusion in the Makile. Don’t use together with -Dturb.
The sedimentation subroutines have been written and used by Hanna Corell in her PhD thesis. They are activated by the -Dsediment in the Makefile. See the theoretical section for more details about the sedimentation model.
This subroutine calculates the settling velocity of the modelled particle size.
This subroutine checks if water velocities are large enough for resuspention of sedimentated trajectories. If large enough it puts the particle on a given level in the water column. generally chosen to be the middle of the bottom box.
This subroutine to sedimentation model for the Baltic sea. The routine calculates the a meta parameter "orb" that is multiplied with an approximation of the wave amplitude in loop, and added to the velocity in the bottom box to simmulate the water movements due to short surface waves.
The diagnostic subroutines are only to generate outputs for the trajectories and Lagrangian stream functions. They do not affect the trajectory solutions themselves.
This subroutine computes the density of the sea water from the equation of state.
This function computes the density of the sea water from the equation of state.
This subroutine computes the density of the sea water from the equation of state but with better flexibility for the pressure.
Interpolates the temperature, salinity/humidity, density/pressure from the grid box walls to the trajectory particle position for the storage of the trajectory characteristics.
Calculates pressure in dbars from depth in meters.
This subroutine writes the Lagrangian stream functions, which are computed when any of the preprocessing options -Dstreamxy, -Dstreamv or -Dstreamr are activated.
This subroutine writes the particle tracer field, which is activated with -Dtracer.
image
image
[hbt]
image
Download tracmass.zip from http://doos.misu.su.se/tracmass/ and put where you want it.
Make necessary changes in the Makefile by setting the appropriate preprocessing options described in previous chapter. Set the apropriate parameters in <project>_grid.in and <project>_run.in Compile by
or
Run TRACMASS by executing
The output data from the TRACMASS runs are stored in the directory <outDataDir>, which is set in <project>_run.in. There are basically two types od output data: trejctories and integrated quantitites such as stream functions.
The trajectories are stored in <outDataDir>/<outDataFile>_run.asc, which has to be specified in <project>_run.in The trajectories are stored in colums of ntrac,niter,x1,y1,z1,tt,t0,subvol,temp,salt,dens where ntrac is the trajectory number niter is the TRACMASS code iteration (only important for TRACMASS modellers) x1 is the zoonal position of the trajectory particle y1 is the meridional position of the trajectory particle z1 is the vertical position of the trajectory particle tt the time of the trajectory particle (in days) t0 the initial time of the trajectory particle subvol is the the "volume" transport in m3/s of the trajectory temp is the temperature of the trajectory particle salt is the salinity/specific humidity of the trajectory particle dens is the density of the trajectory particle
In addition to the trajectories there are integrated quantities stored in the code. These are activated by the preprocessing options in the Makefile. The Lagrangian stream funtions are descibed by Equation [eq:psixy] and [eq:psiyz]. They should be read in the same way as they are written in writepsi.f95. The same goes for the simulated tracer field written in writetracer.f95.
You will need to penetrate the tracmass code in case you want to implement a new setup from a GCM with a grid that do not exist in the present TRACMASS. This requires good knowlidge in the numerical discretisation of your gcm and of Fortran programing. You do not get any support from the TRACMASS team for this. You will here just get some basic guidence how to do this. 1.) Make a new <project>, which we call "mymodel" here. Copy the directory from an existing project that has a similar grid (A,B or C-grid) or has a similar format (netcdf, grib, etc.). So that or or You will now have to modify and adapt the files under tracsvn/projects/mymodel, which are readfield.f95, mymodel_grid.in and mymodel_run.in. 2.) Modify tracsvn/coord.f95 if necssary unless you specify your new grid in tracsvn/mymodel/readfield.f95 3.) You will now have to add your project into the TRACMASS code itself with the many preprocessing options depending on the project you are running. Search for a project with a similar vertical grid as yours (depth, sigma, atmospheric/ifs): or 4) Adapt the Makefile and set PROJECT = mymodel
x
natexlab#1#1url#1#1urlprefix
Blanke, B and Raynaud S, 1997. Kinematics of the Pacific Equatorial Undercurrent: a Eulerian and Lagrangian approach from GCM results. J. Phys. Oceanogr., 27, 1038-1053.
Blanke, B., M. Arhan, G. Madec, and S. Roche, 1999: Warm water paths in the equatorial Atlantic as diagnosed with a general circulation model. J. Phys. Oceanogr., 29, 2753-2768.
Blanke B., S. Speich, G. Madec and K. Döös 2001:A Global Diagnostic of Interocean Mass Transfers. J. Phys. Oceanogr. vol. 31, No. 6, 1623-1642
Döös, K. 1995. Inter-ocean exchange of water masses. J. Geophys. Res. 100 (C7), 13499- 13514.
Döös K. and A.C. Coward 1997: The Southern Ocean as the major upwelling zone of the North Atlantic Deep Water. WOCE Newsletter, No. 27, July 1997, 3-17.
Döös, K., M. Meier and R. Döscher, 2004: The Baltic Haline Conveyor Belt or The Overturning Circulation and Mixing in the Baltic. Ambio, Vol 23, No 4-5, 261-266.
Döös K., A.Engqvist, 2007: Assessment of water exchange between a discharge region and the open sea - A comparison of different methodological concepts. Estuarine, Coastal and Shelf Science. Estuarine, Coastal and Shelf Science 74 (2007) 585-597.
Döös, K., J. Nycander and A.C. Coward, 2008: Lagrangian decomposition of the Deacon Cell. J. Geophys. Res., DOI: 10.1029/ 2007JC004351
Drijfhout S., P. De Vries, K. Döös, A. Coward 2003: Impact of eddy-induced transport of the Lagrangian structure of the upper branch of the thermohaline circulation. J. Phys. Oceanogr. 33, 2141-2155.
Engqvist A., K. Döös and O. Andrejev, 2006: Modeling Water Exchange and Contaminant Transport through a Baltic Coastal Region. Ambio Vol. 35, No. 8.
Huber M., H. Brinkhuis, C. Stickley, K. Döös, A. Sluijs, J. Warnaar, S. Schellenberg, G. Williams, 2004: Eocene circulation of the Southern Ocean: was Antarctica kept warm by subtropical waters? Paleoceanography Vol. 19. PA4026.
Jönsson, B , Lundberg P. and Döös K, 2004. Baltic Sub-Basin Turnover Times Examined Using the Rossby Centre Ocean Model. Ambio, Vol 23, No 4-5, 2257-260.
Kimura S, Döös K, Coward AC, 1999: Numerical simulation to resolve the issue of downstream migration of the Japanese eel. Marine Ecology Progress Series, Vol. 186, p. 303-306.
Levine, R.C. 2005. Changes in shelf waters due to air-sea fluxes and their influence on the Arctic Ocean circulation as simulated in the OCCAM global ocean model. University of Southampton, Faculty of Engineering Science and Mathematics, School of Ocean and Earth Science, PhD Thesis, 225pp.
Marsh, R. and A. P. Megann, 2002: Tracing water masses with particle trajectories in an isopycnic-coordinate model of the global ocean. Ocean Modelling, 4, 27-53.
Mesinger, F. and Arakawa, A., 1976. Numerical methods used in atmospheric models. GARP publications series 17 I 64 pp. .
Nycander J., K Döös and A. C. Coward 2002: Chaotic and regular trajectories in the Antarctic Circumpolar. Tellus - Series A, Vol. 54, issue 1, p.99-106.
Simmons, A.J. and D.M. Burridge, 1981. An Energy and Angular-Momentum Conserving Vertical Finite-Difference Scheme and Hybrid Vertical Coordinates. Monthly Weather Review. 109, 4, 758-766.
Speich S., B. Blanke, P. de Vries , K. Döös, S. Drijfhout , A. Ganachaud and R. Marsh 2002: Tasman leakage: a new route in the global ocean conveyor belt. Geophysical Res. Letters, VOL. 29, NO. 10.
Vries, P. de and K. Döös, 2001. Calculating Lagrangian trajectories using time-dependent velocity fields. J. Atmos. Oceanic Technology. Vol. 18, No. 6, 1092-1101.