On the Computational Efficiency of Ocean Models: A High Performance Equation of State Dov Kruger1, Tal Ezer2, Anne M. Pence1 and Alan F. Blumberg1 1Davidson Laboratory, Civil, Environmental and Ocean Engineering Stevens Institute of Technology, Hoboken, NJ 07030 2Program in Atmospheric and Oceanic Sciences, Princeton University P.O.Box CN710, Sayre Hall, Princeton, NJ 08544-0710 Corresponding author email: [email protected] To be submitted to: Journal of Atmospheric and Oceanic Technology December 20, 20042 Abstract The numerical representation of the equation of state (EOS) in ocean models is reevaluated in terms of its computational efficiency. Ways to speedup the EOS as well as basic codes in ocean models are discussed and demonstrated. It is shown that even a simple rearrangement of a calculation order can significantly speedup the codes currently used in several community terrain-following ocean models. Using a new EOS with particular fit to data that takes into account the computational cost, can speedup the code even further. Various EOS codes were tested on different computer platforms, demonstrating the potential of speeding up codes by as much as 2-3 times on a supercomputer and up to 6-16 times on laptops and personal computers. 1. Introduction The equation of state, used to relate density of sea water to salinity, temperature and pressure, has concerned oceanographers for a long time (e.g., Fofonoff, 1956). While earlier efforts have used simple linear fits to data, a more accurate but complex equation of state (EOS) became the standard (UNESCO, 1981). However, the formulation was very computationally intensive, so ocean modelers have attempted to reduce the cost with improved formulations that may trade off accuracy for speed. For example, early z-level models, whose vertical grid follows constant depth levels, have used a table of coefficients for each level (Bryan and Cox, 1972). Simplified approach of this nature can not be used in non-aligned vertical grids such as in terrain-following (“s” or “sigma” coordinates) models (e.g., the Estuarine, Coastal and Ocean Model, ECOM, Blumberg et al., 1993; the Princeton Ocean Model, POM, Blumberg and Mellor, 1987; the Regional Ocean Modeling System, ROMS, Shchepetkin and McWilliams, 2005), or in generalized coordinate ocean models (Ezer and Mellor, 2004). Note that early sigma models were used mostly for shallow coasts and estuaries, so the pressure dependency of the EOS was often neglected. The extension of sigma ocean model applications to deep-ocean and basin-scale problems (Ezer and Mellor, 1997; Haidvogel et al., 2000) require to include3the pressure effect and led, for example, to the EOS suggested by Mellor (1991), which simplified the pressure term in the UNESCO formulation. There are two factors limiting computational fluid dynamics models, memory and computational time. In the case of ocean modeling, the memory size is usually less of a problem (except in very high resolution or global models), but time is a limiting factor (e.g., calculations of a real-time daily forecast should take much less than a day to be useful). Thus, a faster computation means the ability to have a larger domain, a finer grid, or the ability to perform more experiments. There are three basic approaches for speeding up codes. The first approach is by using parallel codes (e.g., the ROMS code) that take advantage of multiprocessor supercomputers. The second approach is by using numerical schemes that allow longer time steps between calculations (Ezer et al., 2002; Shchepetkin and McWilliams, 2005). The third approach, discussed here, is by optimizing the basic code structure so more calculations can be done per computational time, taking into account the fact that some operations are faster than others. This approach is synergistic with the other approaches. Parallel computing is more efficient when each parallelized unit is faster, and in particular shared memory architectures are much more efficient when the algorithms are efficient of memory bandwidth and are CPU-limited. While large multiprocessor supercomputers are widely used in ocean modeling, at the same time, the speed and memory available on personal computers, PCs (workstations, laptops, etc.), have increased dramatically in recent years, allowing users to run complex ocean models on local PCs. However, surprisingly little attention has been given to improvement in the performance of the basic code itself when using a single processor. It was discovered, that common ocean models may not be optimized for PCs, which was the motivation behind this study. This deficiency in the computational efficiency occur for example in the formulations of the EOS used in three community terrain-following ocean models: ROMS (Haidvogel et al., 2000; Shchepetkin and McWilliams, 2005) which uses the Jackett and McDougall (1995) formulation, POM (Blumberg and Mellor, 1987) which uses the Mellor (1991) formulation, and ECOM (Blumberg et al., 1993) which uses a version of Fofonoff (1956) without the pressure terms. Over 3000 users world wide use these community models and can benefit from improved codes. The above formulations are also compared with a new polynomial fit. The EOS in the above models accounts for ~10% of the total calculations, so optimizing4the EOS alone would speed the entire model by less than 10%. The most computationally demanding part in POM for example is the turbulence scheme (Mellor and Yamada, 1982) which takes ~20% of the entire model calculations, but preliminary tests (not shown) indicate the potential for speeding up the turbulence scheme and other parts of ocean model codes using similar concepts as demonstrated here with the EOS example. 2. The formulation of the equation of state The original UNESCO expression, without the pressure dependent terms can be written in the form, 2412/3233323142532422322215154143132121100)()(),(ScScccScccccccccccS++++++++++++++=θθθθθθθθθθθθρ (1) where ρ, S and θ are the density, salinity and potential temperature, respectively. This formulation was based on fitting the expression to a large sample of oceanic observations in order to find the optimal coefficients cij. Note however, that the original formula used the in-situ temperature, T, while ocean models need the potential temperature, θ, where ρ(S,T) = ρ(S,θ,p=0), and p is pressure. The EOS in the standard POM, is