T-UGOm model developments
Spectral solver discretisations
New discretisations have been added in a still-going effort to improve currents representation. Availble paires are the following:
LGP0xLGP1, DGP1xLGP2, DNP1xLGP2, LGP2xLGP2 (using quadrature integration on Gauss points)
First discretisation given in above paires indicates curents discretisation, the second one indicates elevetion discretisation (DNP1xLGP2 means discontinuous non-conforming P1 discretisation for currents, and continuous Lagrange P2 discretisation for elevations).
Quite extensive testing have shown that the more efficient paire in FES2012 is the DNP1xLGP2. However, some LGP2xLGP2 variants were not fully explored and may require further attention.
Open MPI parallel implementation:
Latest implementation includes HIPS solver. It proved to be more scalable than the MUMPS's one.
Status: available for the time stepping simulation only (real-valued implicit systems)
Open MP parallel implementation:
As most recent computers offer multiple cores, Open MP parallel implementation is considered as really worthy. Contrary to MPI implementation, it keeps most of the code unchanged if one except parallel directives.
Spectral linear solver implementations:
UMFPACK solver tends to get unstable as the number of dof increases (1 500 000 complex unknowns and above). The SpDomestic solver (T-UGOm embarked solver) does not show such vulnerability, but is considerably slower than UMFPACK. In consequence, a new solver is being tested (PASTIX), that is naturally multi-threads optimised (contrary to UMFPACK that is multi-threads capable only through BLAS muti-threading).
Thanks to the gret work of Cyril Nguyen (POC computer engineer), PASTIX is now operational in T-UGOm. Latest benchmarks show improved numerical performances (compared to UMFPACK) and no failure were evidenced yet (29/12/2011).
Forthcoming actions: more testing is needed to optimize PASTIX performances (compilation options, BLAS library etc...). Also use of additional memory during system resolution has been detected (rising use up to 13 Go with FES2012 configuration), A full proofing should be done ASAP.
Ices shelf enhanced friction:
The increase of friction below permanent ice shelf (such as Weddell Sea ice shelf) has been reworked to allow easier generation of frictional members (ensemble data assimilation)
Implicit zero-flux boundary condition:
The zero-flux boundary condition (rigid boundaries) is easily implicitely enforced in the continuity equation. Its implementation in the momentum equation is a bit more tricky, and has been left aside for a while. In the latest T-UGOm version, it has been generalized to improve currents solutions in coastal areas.
Linear friction (still underway):
In some shallow ocean regions, the presence of fine sediments and/or mud flow can dramatically change the physics of bottom friction. In such locations (such as the Northern Amazon River estuary, see Yoann Lebars PhD report)), a linear parameterisation can be more adequate than the usual quadratic one. To comply with those peculiar regions, a new linear, regionally constrained parameterisation has been added to the T-UGOm spectral friction procedures.
Some additional adjustement are needed to avoid solution discontinuities due to the sudden change in friction (smoothed by horizontal diffusion in time-stepping simulation, but left abrupt in spectral simulation).
S2 atmospheric surface pressure forcing
The radiational contribution (i.e. the ocean response to atmospheric pressure S2 tide) represents about 15% of the total S2 ocean tide (astronomical plus radiational). Separation of the astronomical and radiational component in a S2 ocean tide observation is quite tricky, even if "smoothness" hypothesis can be used to deduce the S2 astronomic tides from priorly known semi-diurnal constituents. As assimilated data are a measure of the total S2 ocean tide, it was decided to add the S2 (surface) pressure tide to the forcing of the prior S2 ocean tide simulation, so the prior S2 solution get more consistent with the assimilated data and closer to the (forthcoming) optimal one.
S2 atmospheric tide (units in Pascals), 2010 analysis
Presently, 10 years of ECMWF DCDA fields have harmonically analysed to extract a global surface pressure S2 tide. (equations...) Preliminary tests have shown that std error is dived by a factor 2 when taking atmospheric pressure in to account in the S2 tidal forcing.
Non-linear tides determination
The spectral approach is clearly less convenient for the simulation of non-linear tides. Contrary to the astronomical components, the non-linear tides are mutually inter-dependent. M4 is the product of the interaction of M2 with itself, but also with M6. Reversely M6 is mainly due to the interaction of M2 with M4. An iterative approach is used to solve for a coherent set of of non-linear tides, so that mutual interactions are mostly taken into acccount.
Current tidal prediction (from the intermediate state of resolved tidal spectrum) is used to reconstruct fully non-linear bottom friction. It might be useful to extend this approach to the 2 other non-liear tide producing terms (transport divergence in mass conservation and advection).
Additional tests (based on M4 simulation) have shown that a linear parameterisation of non-linear tides loading potential and deformation might be more adequate than to use a prior LSA field (i.e. computed from a prior tidal solution).
Latitude-dependent gravitational acceleration g
Due to Earth rotation and geometry, the gravitational acceleration g varies with latitude. In combination, the equatorial bulge and the effects of the earths inertia mean that sea-level gravitational acceleration increases from about 9.780 m·s−2 at the Equator to about 9.832 m·s−2 at the poles, so an object will weigh about 0.5% more at the poles than at the Equator.
This point was raised to us by Claire Perigaud some years ago, but was not really pushed forward as the tidal model errors were largely dominated by others terms. However, because of the present reduced error budget in FES2012/T-UGOm (about 9% of M2 amplitude), it was thought that the implementation of the latitude-dependent gravitational constant was worth a trial.
Preliminary tests were successful, with a gain of 10% in M2 simulations.
Internal wave drag above critical latitude
In a linear case, internal tides can not propagate above critical latitudes (Gill, Ocean-Atmosphere Dynamics). Non-linearities might occasionally allow for partially breaking this rule, but numerical experiments tend to show that such exceptions are ... quite exceptionals. Consequently, energy transfer between barotropic tides and baroclinic tides were assumed to be negligibles. However, reduced propagation means only reduced energy transfer, so the internal wave drag has been modified to keep a minimum energy transfer even above critical latitude.
Internal wave drag parameters
The internal tide drag parameters such as Brünt-Vaissala frequency and characteristic wavenumber (tidal excursion) have been explored to get closer to local conditions. The Brünt-Vaissala frequency must be representative of the local dominant vertical mode (usually the first one), and various attempt have been done from climatology or model simulation outputs (solving the classical Sturm-Liouville problem for vertical eigenmode).
Equivalent Brünt-Vaissala frequency for the first vertical mode (computed from WOA2005)
Generic spectral energy estimates and budget
Tidal energy examination is a powerful way to analyse and validate tidal solution. the historical LGP1xLGP1 routines have been reformed to be able to deal with all type of discretisations. Basic diagnostics are exactly computed (integrated) over the elements and then projected on LGP1 discretisation for display purposes.
Ocean partition and toponyms
A new way to prescribe regional-specific values for various parameters has been introduced in complement of polygons partitioning methods. It is based on a gridded file indicating geographical toponyms.