src/test/bleck.c

    Isopycnal gyres intersecting the bathymetry

    This test verifies that isopycnals intersecting the bottom topography do not lead to spurious imbalance. The setup is close, but not identical, to the “two-layer model, high-cut topography” case of section 3 of Bleck & Smith, 1990.

    The velocity/height field can be compared to Fig.3 of Bleck & Smith, 1990.

    Velocity field and surface height

    Velocity field and surface height

    At convergence the bottom layer must be stationary.

    set xlabel 'Time (years)'
set ylabel 'Kinetic energy per unit area (J/m^2)'
set logscale y
year = 86400.*365.
set grid
set key center right
plot [0:5]						\
   'log' u ($1/year):3 w l t 'top layer',		\
   'log' u ($1/year):2 w l t 'bottom layer'
    Evolution of the kinetic enery per unit area for both layers (script)

    Evolution of the kinetic enery per unit area for both layers (script)

    References

    [bleck1990]

    Rainer Bleck and Linda T Smith. A wind-driven isopycnic coordinate model of the north and equatorial Atlantic Ocean: 1. Model development and supporting experiments. Journal of Geophysical Research: Oceans, 95(C3):3273–3285, 1990. [ DOI ]

    		 
    #include "grid/multigrid.h"

#define hour 3600.
#define day (24.*hour)
#define year (365.*day)

    This flag changes the way face heights are computed for isopycnals intersecting the bathymetry (in hydro.h).

    #include "layered/hydro.h"
#include "layered/implicit.h"

    Coriolis

    We include Coriolis acceleration on a \beta-plane.

    double Beta = 2e-11;
#define F0() (0.83e-4 + Beta*y)
#define alpha_H 0.5 // fixme: unstable for alpha_H 1
#include "layered/coriolis.h"

    Vertical stratification

    There are only two (isopycnal) layers. The bottom layer has a relative density difference of 10-2.

    double * drho = (double []){ 1e-2, 0. };
#include "layered/isopycnal.h"

    Surface “wind stress”

    The “surface wind stress” is added directly as a (depth-weighted) acceleration (along the x-direction) in the top layer.

    double tau0 = 1e-4;
  
event acceleration (i++)
{
  foreach_face(x)
    ha.x[0,0,1] += hf.x[] > dry ? tau0*cos(2.*pi*y/L0) : 0.;
}

    Horizontal viscosity

    double nu_H = 1e5;

event acceleration (i++)
{
  if (nu_H > 0.) {
    vector d2u[];
    foreach_layer() {
      foreach()
	foreach_dimension()
	  d2u.x[] = (h[1]*u.x[1] + h[-1]*u.x[-1] +
		     h[0,1]*u.x[0,1] + h[0,-1]*u.x[0,-1] - 4.*h[]*u.x[])/
	sq(Delta);
      foreach()
	if (h[] > dry)
	  foreach_dimension()
	    u.x[] += dt*nu_H*d2u.x[]/max(h[],10.);
    }
  }
}

    main

    int main()
{
  G = 9.8;
  nl = 2;
  size (6000e3 [1]);
  origin (0, - L0/2.);
  DT = 3000 [0,1];
  N = 32;
  theta_H = 0.51;
  run();
}

    Initial conditions

    event init (i = 0)
{

    No-slip (and dry) conditions on all boundaries. Note that this is only relevant for “wet” boundaries. For example, given that there is a (dry) coastline at x = 4400 km (see below), the u.t[right] boundary condition below is useless.

      foreach_dimension() {
    u.t[left] = dirichlet(0);
    zb[left] = 1000.;
    h[left] = 0.;
    u.t[right] = dirichlet(0);
    zb[right] = 1000.;
    h[right] = 0.;
  }

  foreach() {

    The bathymetry.

    #if 0
    // low-cut case
    zb[] = x > 1400e3 ? - 4000. : - 1550. - 2450.*x/1400e3;
#elif 1
    // high-cut case
    zb[] = x > 2000e3 ? - 4000. : - 200. - 3800.*x/2000e3;
#else // flat bottom
    zb[] = -4000;
#endif

    A “coastline” at x = 4400 km.

        if (x > 4400e3)
      zb[] = 1000.;

    The bottom layer is between the bottom and - 1000 m, the top layer is between 0 and -1000 m.

        h[] = max (- 1000. - zb[], 0.);
    h[0,0,1] = max (- zb[] - h[], 0.);
  }
}

    Outputs

    #include "view.h"

event end (t = 5*year)
{
  view (fov = 19.7885,
	tx = -0.35, ty = 0, width = 440, height = 600);
  squares ("zb < 100 ? eta : nodata", spread = -1);
  vectors ("u1", scale = 8e5);
  save ("gyre.png");
}

event logfile (i += 300)
{
  double ke0 = 0., ke1 = 0., area = 0., rho0 = 1000;
  scalar h1 = lookup_field ("h1");
  vector u1 = lookup_vector ("u1");
  foreach (reduction (+:ke0) reduction(+:ke1) reduction(+:area)) {
    ke0 += dv()*h[]*(sq(u.x[]) + sq(u.y[]))/2.;
    ke1 += dv()*h1[]*(sq(u1.x[]) + sq(u1.y[]))/2.;
    area += dv();
  }
  fprintf (stderr, "%g %g %g %g %g %g %d %d\n", t, rho0*ke0/area, rho0*ke1/area,
	   statsf (h).sum, statsf (h1).sum, dt, mgH.i, mgH.nrelax);
}