src/examples/gulf-stream.c

The Gulf Stream

We seek a “minimal setup” able to reproduce a realistic oceanic circulation in the North Atlantic, with the Gulf Stream an obvious dominant feature. The setup is as close as possible to that used by Hurlburt & Hogan, 2000 but uses the layered solver described in Popinet, 2020.

See Hurlburt & Hogan, 2000 for details but the main characteristics of the setup are:

• 5 isopycnal layers (see Table 2 of H&H, 2000 for the reference densities, thicknesses etc.)
• Wind stress in the top layer given by the monthly climatology of Hellerman & Rosenstein, 1983
• “Compressed bathymetry” as in H&H, 2000
• Quadratic bottom friction (Cb = 2 x 10^-3), Laplacian horizontal viscosity (10 m²/s)
• Atlantic Meridional Overturning Circulation (AMOC) driven by fluxes at the northern and southern boundaries (see Table 2 of H&H, 2000 for the values of fluxes)

Setup

We use the generic ocean model with 5 isopycnal layers. The documentation of the generic model needs to be consulted for details.

#define NL 5
#include "ocean.h"

Isopycnal layers

The setup uses five isopycnal layers with an initial thickness of 250 meters for the top 4 layers (the dh array below). The corresponding relative density differences are given by the drho array and correspond to those in Table 2 of Hurlburt & Hogan, 2000.

double * dh   = (double [NL]){ 0., 250, 250, 250, 250 };
double * drho = (double [NL]){ 2.13/rho0, 1.72/rho0, 1.41/rho0, 1.01/rho0, 0. };

Diapycnal entrainment

See entrainment.h for explanations. The mininum and maximum layer thicknesses are given by hmin and hmax (in meters here). The average entrainement velocity is set to 0.07 cm/s (see Table 2 of H&H, 2000) and the coefficient of additional interfacial friction associated with entrainment is zero.

double * hmin = (double [NL]){ 40, 40, 40, 40, 50 };
double * hmax = (double [NL]){ HUGE, HUGE, HUGE, HUGE, HUGE };
double omr = 0.07e-2, Cm = 0.;

Boundary fluxes

The Atlantic Meridional Overturning Circulation (AMOC) is driven primarily by heat (and salinity) fluxes at the surface of the ocean, which are not included in this simplified model. Reasonably realistic North Atlantic circulations can be obtained without them but much more realistic results are obtained when including the simplified representation proposed by Hurlburt & Hogan, 2000. See bflux.h for details.

#if 1
#include "bflux.h"
#endif

main()

The simulation can be run in OpenMP/serial/GPU or with MPI. With MPI the number of processes must be of the form 2^(i+1) i.e.: 8, 32, 128, 512, 2048 etc. (see Tips for explanations). This can be done using the Makefile with:

CC='mpicc -D_MPI=8' CFLAGS=-disable-dimensions make gulf-stream.tst

it can also be run on GPUs using e.g.

OMP_NUM_THREADS=8 CFLAGS=-DSHOW make gulf-stream.gpu.tst

Command-line parameters can be passed to the code to change the spatial resolution and/or the timestep.

int main (int argc, char * argv[])
{

The domain is rectangular with a 2x1 aspect ratio.

  dimensions (2, 1);

On a spherical grid, the sizes are given in degrees (of longitude).

  size (84);

The Earth radius, here in meters, sets the length unit.

  Radius = 6371220 [1];

The longitude and latitude of the lower-left corner.

  origin (-98, 9);

The acceleration of gravity sets the time unit (seconds here).

  G = 9.8;

The default resolution is 512 (longitude) x 256 (latitude) i.e. 42/256 \approx 1/6 of a degree.

  if (argc > 1)
    N = atoi(argv[1]);
  else
    N = 512;

The default timestep is 600 seconds. Note that using a larger timestep (at low resolutions) can significantly affect the structure of the boundary current.

  DT = 600 [0,1];
  if (argc > 2)
    DT = atof(argv[2]);

We can change the spinup time using the third command-line parameter.

  if (argc > 3)
    tspinup = atof(argv[3]);

The number of layers is set to NL (five).

  nl = NL;

The “implicitness parameter” of the implicit barotropic solver is set to 0.55 rather than the default 0.5 (“Crank-Nicholson”). This causes numerical damping of the barotropic mode and is necessary to prevent “basin resonances” at low resolution (N = 512). At higher resolutions, values closer to 0.5 (e.g. 0.51) seem to work fine. At all resolutions the sensitivity of the results to this parameter is low (resonances excepted).

  theta_H = 0.55;

  run();
}

Run times

The simulation can be run on parallel machines using something like:

../qcc -source -D_MPI=1 gulf-stream.c
scp _gulf-stream.c navier.lmm.jussieu.fr:gulf-stream/
mpicc -Wall -std=c99 -D_XOPEN_SOURCE=700 -O2 _gulf-stream.c -o gulf-stream -L$HOME/lib -lkdt -lm

On 128 cores of the “navier” cluster at d’Alembert, runtimes are of the order of 89 simulated years per day (ypd) for a resolution of 512 x 256, 28 ypd for 1024 x 512 and 6 ypd for 2048 x 1024.

The simulation also runs fine on GPUs using e.g.

OMP_NUM_THREADS=16 make gulf-stream.gpu.tst

On an RTX 4090 runtimes are approx. 48 ypd for a resolution of 2048 x 1024.

Spinup takes at least 5 years and statistics (for e.g. the standard deviations of SSH) at least 10 years.

Results

All the reference generated files (log, perfs, *.mp4 etc.) are accessible in the following folders (you will need to specify the exact access path e.g. ../4096/eta.mp4):

• gulf-stream/512/
• gulf-stream/1024/
• gulf-stream/2048/
• gulf-stream/4096/

Todo

• Postprocessing (as in the sandbox).
• Boundary fluxes cause large velocities which may unnecessarily restrict the timestep: smoother conditions may then improve runtimes.
• Real (“uncompressed”) bathymetry runs but the Gulf Stream stays attached: this may be caused by a tuning of isopycnal layers, fluxes etc. which is specific to the compressed bathymetry.

References