Define U0 2
Define V0 3
Define W0 1

# Latitude (so that rotation is only along z)
Define Phi (M_PI/4.)

# Coriolis frequency
Define Omega 7.292e-5

# Analytical solution
Define A0 U0
Define B0 (sin(Phi)*V0-cos(Phi)*W0)
Define C0 (cos(Phi)*(cos(Phi)*V0+sin(Phi)*W0))
Define D0 (sin(Phi)*(cos(Phi)*V0+sin(Phi)*W0))
Define Usol0 (A0*cos(2*Omega*t)+B0*sin(2*Omega*t))
Define Vsol0 (-A0*sin(Phi)*sin(2*Omega*t) + B0*sin(Phi)*cos(2*Omega*t) + C0)
Define Wsol0 (A0*cos(Phi)*sin(2*Omega*t) - B0*cos(Phi)*cos(2*Omega*t) + D0)

1 3 GfsSimulation GfsBox GfsGEdge {} {
  Time { end = 30000  }

  Refine 2
  
  # The reference length is set to 60 m
  PhysicalParams { L = 60 }

  # The domain is 60 x 30 x 15 m with an anisotropic mesh
  MetricStretch {} { sy = 0.5 sz = 0.25 }
  
  # Initialisation of the velocity field
  Init {} {
      U = U0
      V = V0
      W = W0
  }
  
  # Viscosity (or eddy-viscosity)
  SourceViscosity 10

  # Coriolis forces
  SourceCoriolis (2.*Omega) 0. { x = 0. y = 1. z = 1. }

  # Output of the integral of each component of the velocity field over the domain every
  # 1000 time steps
  OutputScalarStats { step = 3000 } { awk '{print $3, $7}' > u.dat } { v = U }
  OutputScalarStats { step = 3000 } { awk '{print $3, $7}' > v.dat } { v = V }
  OutputScalarStats { step = 3000 } { awk '{print $3, $7}' > w.dat } { v = W }

  # Output of some kind of error measurement
  OutputScalarStats { istep = 10 } { awk '{print $3,$7}' > error.dat } {
      v = sqrt((U - Usol0)*(U - Usol0) + (V - Vsol0)*(V - Vsol0) + (W - Wsol0)*(W - Wsol0))
  }

  # Output of the final simulation file
  OutputSimulation { start = end } end.gfs
}
GfsBox {}
# Setup of periodic boundary conditions in the 3 directions of space
1 1 top
1 1 right
1 1 front