Define U0 (2.*M_PI)

6 12 GfsAdvection GfsBox GfsGEdge {} {
  PhysicalParams { L = 2.*M_PI/4. }
  MetricCubed M LEVEL

  # Use alternative implementation
  # MetricCubed1 M
  # MapFunction {
  #     x = atan2 (X, Z)*180./M_PI
  #     y = asin (CLAMP(Y,-1.,1.))*180./M_PI
  # }

  Time { end = 1 }
  Refine LEVEL
  VariableTracer T { 
      gradient = gfs_center_gradient 
      cfl = 1
  }
  Global {
      #define DTR (M_PI/180.)
      double lambda1 (double lambda, double theta, double lambdap, double thetap) {
	  /* eq. (8) */
	  return atan2 (cos (theta)*sin (lambda - lambdap),
	                cos (theta)*sin (thetap)*cos (lambda - lambdap) - 
                        cos (thetap)* sin (theta));
      }
      double theta1 (double lambda, double theta, double lambdap, double thetap) {
	  /* eq. (9) */
	  return asin (sin (theta)*sin (thetap) + cos (theta)*cos (thetap)*cos (lambda - lambdap));
      }
      double bell (double lambda, double theta, double lc, double tc) {
	  double h0 = 1.;
	  double R = 1./3.;
	  double r = acos (sin(tc)*sin(theta) + cos (tc)*cos (theta)*cos (lambda - lc));
	  return r >= R ? 0. : (h0/2.)*(1. + cos (M_PI*r/R));
      }
      double bell_moving (double lambda, double theta, double t) {
	  double lambda2 = lambda1 (lambda, theta, M_PI, M_PI/2. - ALPHA*DTR);
	  double theta2 = theta1 (lambda, theta, M_PI, M_PI/2. - ALPHA*DTR);
	  double lc = 3.*M_PI/2. - U0*t, tc = 0.;
	  return bell (lambda2, theta2, lc, tc);
      }      
  }
  Init {} { T = bell_moving(x*DTR,y*DTR,0) }
  VariableStreamFunction Psi -U0*(sin (y*DTR)*cos (DTR*ALPHA) - cos (x*DTR)*cos (y*DTR)*sin (DTR*ALPHA))
  AdaptGradient { istep = 1 } { cmax = 1e-4 maxlevel = LEVEL } T
  OutputTime { istep = 10 } stderr
#  OutputSimulation { istep = 10 } stdout
  OutputSimulation { start = end } end-LEVEL-ALPHA.gfs
  OutputErrorNorm { istep = 1 } { awk '{ print LEVEL,$3,$5,$7,$9}' > error-LEVEL-ALPHA } { v = T } {
      s = bell_moving(x*DTR,y*DTR,t)
      v = E
      relative = 1
  }
  OutputScalarSum { istep = 1 } t-LEVEL-ALPHA { v = T }
  OutputScalarSum { istep = 1 } area-LEVEL-ALPHA { v = 1 }
  EventScript { start = end } {
      ( cat isolines.gfv
        echo "Save isolines.gnu { format = Gnuplot }"
        echo "Clear"
        cat reference.gfv
        echo "Save reference.gnu { format = Gnuplot }"
	echo "Clear"
        cat zero.gfv
        echo "Save zero.gnu { format = Gnuplot }"
      ) | gfsview-batch2D end-LEVEL-ALPHA.gfs
      cat <<EOF | gnuplot
        set term postscript eps lw 2 18 color
        set output 'isolines-LEVEL-ALPHA.eps'
        set size ratio -1
        set xlabel 'Longitude'
        set ylabel 'Latitude'
        unset key
        plot [60:120][-30:30]'isolines.gnu' w l, 'reference.gnu' w l, 'zero.gnu' w l
EOF
      fixbb isolines-LEVEL-ALPHA.eps
      rm -f isolines.gnu reference.gnu zero.gnu

      # check mass conservation
      if awk '
        BEGIN { min = 1000.; max = -1000.; }{ 
          if ($5 < min) min = $5; 
          if ($5 > max) max = $5; 
        }
        END {
          if (max - min != 0.)
            exit (1);
        }' < t-LEVEL-ALPHA; then
	  exit 0
      else
	  exit $GFS_STOP
      fi
  }
}
GfsBox {}
GfsBox {}
GfsBox {}
GfsBox {}
GfsBox {}
GfsBox {}
1 2 right
2 3 top
3 4 right
4 5 top
5 6 right
6 1 top
1 3 top left
3 5 top left
5 1 top left
2 6 bottom right
4 2 bottom right
6 4 bottom right