#include "grid/multigrid.h"
#include "spherical.h"
#include "layered/hydro.h"
#include "layered/implicit.h"
#include "equilibrium_tide.h"

const double Omega = 7.292205e-5;
#define F0() (2.*Omega*sin(y*pi/180.))

const double Cb = 2e-3;
#define K0() (point.l > 0 && h[0,0,-1] > 10. ? 0. : h[] < dry ? HUGE : Cb*norm(u)/h[])
#include "layered/coriolis.h"

#define NL 1

#include "terrain.h" // for bathymetry
#include "layered/perfs.h"
const int aspect = 2; // the aspect ratio of the domain, here 2x1

const double hour = 3600., day = 86400., month = 30.*day, year = 365.*day;

const double BetaSAL = 0.1121;

double tspinup = 8*day;

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

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

  if (npe() > 1) {
    dimensions (ny = sqrt(npe()/aspect));
    size (360);
  }
  else
    size (360);
  periodic (right);

  Radius = 6371220.;

  origin (0, - 180./aspect);

  G = 9.8*(1. - BetaSAL);

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

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

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

  nl = NL;
  TOLERANCE = 1e-3;  
  theta_H = 0.55;

  run();

  system ("ffmpeg -y -i eta.mp4 eta-%04d.png;"
	  "ffmpeg -y -i tide.mp4 tide-%04d.png;"
	  "for f in eta-*.png; do"
	  "  f1=`echo $f | sed 's/eta//'`;"
	  "  montage -tile 1x2 -geometry 1024x456+0+5 tide$f1 eta$f1 montage$f1;"
	  "done;"
	  "ffmpeg -y -i montage-%04d.png -crf 18 -pix_fmt yuv420p -movflags +faststart montage.mp4;"
	  "rm -f eta-*.png tide-*.png montage-*.png;");
}

scalar zbs[], tide[];

event init (i = 0)
{

  equilibrium_tide_constituents();

  if (restore ("restart"))
    event ("metric");
  else {

    terrain (zbs, "~/terrain/etopo2", NULL);

    equilibrium_tide (tide, 0);
    foreach() {
      zb[] = y < -80 || y > 80 || zbs[] > - 10 ? 100 : zbs[];
      h[] = max(tide[]/(1. - BetaSAL) - zb[], 0.);
    }
  }

  u.t[top] = dirichlet(0);
  u.t[bottom] = dirichlet(0);
  zb[top] = 1000;
  zb[bottom] = 1000;
  h[top] = 0;
  h[bottom] = 0;
}

event acceleration (i++)
{
  equilibrium_tide (tide, t);
  foreach_face() {
    double ax = G/(1. - BetaSAL)*gmetric(0)*(tide[-1] - tide[])/Delta;
    foreach_layer()
      ha.x[] -= hf.x[]*ax;
  }
}

#include "harmonic.h"

scalar etaE[];

event harmonic (t = tspinup; i++)
{
  harmonic_decomposition (eta, t, (double[]){
      Etide.Q1.omega,
      Etide.O1.omega,
      Etide.K1.omega,
      Etide.N2.omega,
      Etide.M2.omega,
      Etide.S2.omega,
      Etide.K2.omega,
      0
    }, etaE);
}

double nu_H = 10; // m^2/s

event viscous_term (i++)
{
  if (nu_H > 0.) {
    vector d2u[];
    foreach_layer() {
      double dry = 1.;
      foreach()
	foreach_dimension()
	d2u.x[] = 2.*(sq(fm.x[1])/(cm[1] + cm[])*u.x[1]*(h[1] > dry) +
		      sq(fm.x[])/(cm[-1] + cm[])*u.x[-1]*(h[-1] > dry) +
		      sq(fm.y[0,1])/(cm[0,1] + cm[])*u.x[0,1]*(h[0,1] > dry) +
		      sq(fm.y[0,-1])/(cm[0,-1] + cm[])*u.x[0,-1]*(h[0,-1] > dry))
	/(sq(Delta)*cm[]);
      foreach()
	foreach_dimension() {
	double n = 2.*(sq(fm.x[1])/(cm[1] + cm[])*(1. + (h[1] <= dry)) +
		       sq(fm.x[])/(cm[-1] + cm[])*(1. + (h[-1] <= dry)) +
		       sq(fm.y[0,1])/(cm[0,1] + cm[])*(1. + (h[0,1] <= dry)) +
		       sq(fm.y[0,-1])/(cm[0,-1] + cm[])*(1. + (h[0,-1] <= dry)))
	  /(sq(Delta)*cm[]);
	u.x[] = (u.x[] + dt*nu_H*d2u.x[])/(1. + dt*nu_H*n);
      }
    }
  }
}

event outputs (t += hour)
{
  double ke = 0., keb = 0., vol = 0., volb = 0.;
  scalar etad[], nu[];
  foreach(reduction(+:ke) reduction(+:vol) reduction(+:keb) reduction(+:volb)) {
    point.l = 0;
    keb += dv()*h[]*(sq(u.x[]) + sq(u.y[]));
    volb += dv()*h[];
    foreach_layer() {
      ke += dv()*h[]*(sq(u.x[]) + sq(u.y[]));
      vol += dv()*h[];
    }
    point.l = nl - 1;
    etad[] = h[] > dry ? eta[] : 0.;
    nu[] = h[] > dry ? norm(u) : 0.;
  }

  if (i == 0) {
    fprintf (stderr, "t ke/vol keb/vol dt "
	     "mgH.i mgH.nrelax etad.stddev nu.stddev");
    for (int l = 0; l < nl; l++)
      fprintf (stderr, " d%s%d.sum/dt", h.name, l);
    fputc ('\n', stderr);
  }
  fprintf (stderr, "%g %g %g %g %d %d %g %g", t/day, ke/vol/2., keb/volb/2., dt,
	   mgH.i, mgH.nrelax,
	   statsf (etad).stddev, statsf(nu).stddev);

  static double s0[NL], t0 = 0.;
  foreach_layer() {
    double s = statsf(h).sum;
    if (i == 0)
      fprintf (stderr, " 0");
    else
      fprintf (stderr, " %g", (s - s0[_layer])/(t - t0));
    s0[_layer] = s;
  }
  fputc ('\n', stderr);
  t0 = t;
}

event movies (t = 20*day; t <= 24*day; t += 600)
{
  scalar etad[], m[];
  foreach() {
    point.l = nl - 1;
    etad[] = h[] > dry ? eta[] : 0.;
    m[] = etad[] - zbs[];
  }

  output_ppm (etad, mask = m, file = "eta.mp4", n = clamp(N,1024,2048),
	      min = -2, max = 2,
	      box = {{X0,-80},{X0+L0,80}},
	      map = jet);
  output_ppm (tide, mask = m, file = "tide.mp4", n = clamp(N,1024,2048),
	      min = -0.4, max = 0.4,
	      box = {{X0,-80},{X0+L0,80}},
	      map = jet);
}

Gauge gauges[] = {
  // file      longitude   latitude    description
  {"dahouet", -2.25 + 360., 49.25, "Dahouët, Bretagne"},
  {NULL}
};

event gauges1 (t = tspinup; i++) {
  scalar Q1a = eta.harmonic.A[0], Q1b = eta.harmonic.B[0];
  scalar O1a = eta.harmonic.A[1], O1b = eta.harmonic.B[1];
  scalar K1a = eta.harmonic.A[2], K1b = eta.harmonic.B[2];
  scalar N2a = eta.harmonic.A[3], N2b = eta.harmonic.B[3];
  scalar M2a = eta.harmonic.A[4], M2b = eta.harmonic.B[4];
  scalar S2a = eta.harmonic.A[5], S2b = eta.harmonic.B[5];
  scalar K2a = eta.harmonic.A[6], K2b = eta.harmonic.B[6];
  scalar Z = eta.harmonic.Z;
  output_gauges (gauges, {eta, Z, etaE,
			  Q1a, Q1b, O1a, O1b, K1a, K1b, N2a, N2b, M2a, M2b, S2a, S2b, K2a, K2b});
}

event harmonic_outputs (t = tspinup; t += hour)
{
  if (eta.harmonic.invertible) {
    scalar K1a = eta.harmonic.A[2], K1b = eta.harmonic.B[2];
    scalar M2a = eta.harmonic.A[4], M2b = eta.harmonic.B[4];
    scalar S2a = eta.harmonic.A[5], S2b = eta.harmonic.B[5];
    scalar K1[], M2[], S2[], m[];
    foreach() {
      point.l = nl - 1;
      K1[] = pow(sq(K1a[]) + sq(K1b[]), 1./4.);
      M2[] = pow(sq(M2a[]) + sq(M2b[]), 1./4.);
      S2[] = pow(sq(S2a[]) + sq(S2b[]), 1./4.);
      m[] = (h[] > dry ? eta[] : 0.) - zbs[];
    }
    output_ppm (M2, mask = m, file = "M2.mp4", n = clamp(N,1024,2048),
		linear = false,
		min = 0, max = 1.5,
		box = {{X0,-80},{X0+L0,80}},
		map = jet);
    output_ppm (S2, mask = m, file = "S2.mp4", n = clamp(N,1024,2048),
		linear = false,
		min = 0, max = 1.5,
		box = {{X0,-80},{X0+L0,80}},
		map = jet);
    output_ppm (K1, mask = m, file = "K1.mp4", n = clamp(N,1024,2048),
		linear = false,
		min = 0, max = 1.5,
		box = {{X0,-80},{X0+L0,80}},
		map = jet);
  }
}

event final (t = 3*month)
{
  output_field (all, box = {{X0,-80},{X0+L0,80}});
  dump();  
}

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