%%% "integrate plane Couette flow from a given initial condition and save velocity fields to disk."
function couette(u, [T0 T1], dt, dT, Reynolds, label)

  ArgList args(argc, argv, purpose);

  const Real T0       = args.getreal("-T0", "--T0", 0.0, "start time");
  const Real T1       = args.getreal("-T1", "--T1", 100, "end time");
  const Real dtarg    = args.getreal("-dt", "--dt", 0.0, "timestep, 0 => variable");
  const Real dT       = args.getreal("-dT", "--dT", 1.0, "save interval");
  const string outdir = args.getpath("-o", "--outdir", "data/", "output directory");
  const string label  = args.getstr ("-l", "--label", "u", "output field prefix");
  const Real Reynolds = args.getreal("-R", "--Reynolds", 400, "Reynolds number");
  const string uname = args.getstr (1, "<flowfield>", "initial condition");
  
  args.check();
  args.save("./");
  args.save(outdir);
  mkdir(outdir);

  FlowField u(uname);
  const int Nx = u.Nx();
  const int Ny = u.Ny();
  const int Nz = u.Nz();
  const Real Lx = u.Lx();
  const Real Lz = u.Lz();
  const Real a = u.a();
  const Real b = u.b();


  const Real CFLmin = 0.40;
  const Real CFLmax = 0.60;
  const Real dtmin  = 0.001;
  const Real dtmax  = 0.10;
  const bool variable_dt = (dtarg == 0.0) ? true : false;
  const Real dtinit = (variable_dt) ? 0.03125 : dtarg;
 
  TimeStep dt(dtinit, dtmin, dtmax, dT, CFLmin, CFLmax);
  const bool inttime = (abs(dT - int(dT)) < 1e-12) ? true : false;


  // Construct Navier-Stoke Integrator
  DNSFlags flags; 
  flags.timestepping = SBDF3;
  flags.dealiasing   = DealiasXZ;
  flags.nonlinearity = Rotational;
  flags.constraint   = PressureGradient;
  flags.dPdx         = 0.0;

  // Define flow parameters
  const Real nu = 1.0/Reynolds;

  cout << "Constructing u,q, and optimizing FFTW..." << flush;
  FlowField q(Nx,Ny,Nz,1,Lx,Lz,a,b);

  // Construct base flow for plane Couette: U(y) = y
  Vector y = chebypoints(Ny,a,b);
  
  cout << "Nx,Ny,Nz == " << Nx << ',' << Ny << ',' << Nz << endl;
  ChebyTransform trans(Ny);
  ChebyCoeff Ubase(Ny,a,b,Physical);
  for (int ny=0; ny<Ny; ++ny) 
    Ubase[ny] = y[ny];

  cout << "constructing DNS..." << flush;
  DNS dns(u, Ubase, nu, dt, flags, Real(T0));
  if (variable_dt && dt.adjust(dns.CFL()))
    dns.reset_dt(dt);
  cout << "done" << endl;

  for (Real t=T0; t<=T1; t += dT) {
    
    Real l2norm = L2Norm(u);
    cout << "===============================================" << endl;
    cout << "           t == " << t << endl;
    cout << "          dt == " << dt << endl;
    cout << "         CFL == " << dns.CFL() << endl;
    cout << "   L2Norm(u) == " << l2norm << endl;
    cout << "       Ubulk == " << dns.Ubulk() << endl;
    cout << "       ubulk == " << Re(u.profile(0,0,0)).mean() << endl;
    cout << "        dPdx == " << dns.dPdx() << endl;

    dns.advance(u, q, dt.n());
    if (variable_dt && dt.adjust(dns.CFL()))
      dns.reset_dt(dt);

    u.binarySave(outdir + label + t2s(t, inttime));
  }
  
  cout << "done!" << endl;