#include "DAGH.h"
#include "t2.h"


#define DIM 3
#define MaxLev 1
/* this routines will just initialize the grid variables and then
  invoke the mg solver, and determine the error after each cycle
*/

void main(INTEGER argc, CHARACTER *argv[])  
{

   MPI_Init( &argc, &argv );

   DOUBLE start = MPI_Wtime();
   INTEGER nx ,ny ,nz,nplot;
   DOUBLE xmin,ymin,zmin,xmax,ymax,zmax,tmin,tmax;
   DOUBLE mass,tol;

   Array(2)<DOUBLE> bhr(1,1,2,3);
   Array(2)<DOUBLE> bhp(1,1,2,3);
   Array(2)<DOUBLE> bhs(1,1,2,3);

//  Characteristic functions

   INTEGER const nlmax = 6*(33*33+33*33+33*33) ;
   INTEGER nl;
//   DOUBLE d[3*nlmax];
//   INTEGER m[3*nlmax];
   INTEGER xi[nlmax],yj[nlmax],zk[nlmax];
   DOUBLE  fx[nlmax],fy[nlmax],fz[nlmax];


   INTEGER i,j,k;
   DOUBLE dx,dy,dz;
   INTEGER it,nt,itsave;
   DOUBLE dt,dt2,dtdx,time;
   DOUBLE dtsave,timesave;
   DOUBLE bb[2*DIM];
   DOUBLE bha[2];

// Debugging variables
	DOUBLE n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13,n14,n15,n16,n17,n18;

  // Set up the Grid Structure
  INTEGER shape[DIM]; 
   cout << "Begin TCODE run." << "\n" <<flush;
   INTEGER hdfflag,asciiflag;

   f_input(&nx,&ny,&nz,&nt,&nplot,&xmin,&xmax,&ymin,&ymax,&zmin,&zmax,
           &hdfflag,&asciiflag,bha,FDA(bhr),FDA(bhp),FDA(bhs),&tol);
  shape[0] = nx-1;
  bb[0] = xmin;
  bb[1] = xmax;
  shape[1] = ny-1;
  bb[2] = ymin;
  bb[3] = ymax;
  shape[2] = nz-1;
  bb[4] = zmin;
  bb[5] = zmax;


   dx = ( xmax - xmin ) / (DOUBLE) ( nx - 1 );
   dy = ( ymax - ymin ) / (DOUBLE) ( ny - 1 );
   dz = ( zmax - zmin ) / (DOUBLE) ( nz - 1 );

   time = (DOUBLE) 0.0;
   it = 0; 

   DOUBLE amp, width;
   amp = 0.5;
   width = 0.0001;


  GridHierarchy GH(DIM,NON_CELL_CENTERED,MaxLev);
  SetBaseGrid (GH,bb,shape);
  ComposeHierarchy(GH);

   INTEGER t_sten = 0; INTEGER s_sten = 2; INTEGER t_sten2 = 1;

          /* Total # of variables is now =78 */
 
   GridFunction(DIM)<DOUBLE> g11("g11",t_sten2,s_sten,GH);
   SetTimeAlias(g11,1,-1);
   GridFunction(DIM)<DOUBLE> g12("g12",t_sten2,s_sten,GH);
   SetTimeAlias(g12,1,-1);
   GridFunction(DIM)<DOUBLE> g13("g13",t_sten2,s_sten,GH);
   SetTimeAlias(g13,1,-1);
   GridFunction(DIM)<DOUBLE> g22("g22",t_sten2,s_sten,GH);
   SetTimeAlias(g22,1,-1);
   GridFunction(DIM)<DOUBLE> g23("g23",t_sten2,s_sten,GH);
   SetTimeAlias(g23,1,-1);
   GridFunction(DIM)<DOUBLE> g33("g33",t_sten2,s_sten,GH);
   SetTimeAlias(g33,1,-1);

   GridFunction(DIM)<DOUBLE> gup11("gup11",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> gup12("gup12",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> gup13("gup13",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> gup22("gup22",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> gup23("gup23",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> gup33("gup33",t_sten,s_sten,GH);

   GridFunction(DIM)<DOUBLE> k11("k11",t_sten2,s_sten,GH);
   SetTimeAlias(k11,1,-1);
   GridFunction(DIM)<DOUBLE> k12("k12",t_sten2,s_sten,GH);  
   SetTimeAlias(k12,1,-1);
   GridFunction(DIM)<DOUBLE> k13("k13",t_sten2,s_sten,GH);  
   SetTimeAlias(k13,1,-1);
   GridFunction(DIM)<DOUBLE> k21("k21",t_sten2,s_sten,GH);  
   SetTimeAlias(k21,1,-1);
   GridFunction(DIM)<DOUBLE> k22("k22",t_sten2,s_sten,GH);  
   SetTimeAlias(k22,1,-1);
   GridFunction(DIM)<DOUBLE> k23("k23",t_sten2,s_sten,GH);  
   SetTimeAlias(k23,1,-1);
   GridFunction(DIM)<DOUBLE> k31("k31",t_sten2,s_sten,GH);  
   SetTimeAlias(k31,1,-1);
   GridFunction(DIM)<DOUBLE> k32("k32",t_sten2,s_sten,GH);
   SetTimeAlias(k32,1,-1);
   GridFunction(DIM)<DOUBLE> k33("k33",t_sten2,s_sten,GH);  
   SetTimeAlias(k33,1,-1);
 
   GridFunction(DIM)<DOUBLE> c111("c111",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c112("c112",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c113("c113",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c211("c211",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c212("c212",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c213("c213",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c222("c222",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c223("c223",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c233("c233",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c311("c311",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c312("c312",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c313("c313",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c322("c322",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c323("c323",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c333("c333",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c122("c122",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c123("c123",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> c133("c133",t_sten,s_sten,GH);
 
   GridFunction(DIM)<DOUBLE> r11("r11",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> r12("r12",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> r13("r13",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> r22("r22",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> r23("r23",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> r33("r33",t_sten,s_sten,GH);

   GridFunction(DIM)<DOUBLE> alpha("alpha",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> beta1("beta1",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> beta2("beta2",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> beta3("beta3",t_sten,s_sten,GH);

   GridFunction(DIM)<DOUBLE> mask("mask",t_sten2,s_sten,GH);

   GridFunction(DIM)<DOUBLE> t11("t11",t_sten,s_sten,GH,NO_COMM);
   GridFunction(DIM)<DOUBLE> t12("t12",t_sten,s_sten,GH,NO_COMM);
   GridFunction(DIM)<DOUBLE> t13("t13",t_sten,s_sten,GH,NO_COMM);
   GridFunction(DIM)<DOUBLE> t21("t21",t_sten,s_sten,GH,NO_COMM);
   GridFunction(DIM)<DOUBLE> t22("t22",t_sten,s_sten,GH,NO_COMM);
   GridFunction(DIM)<DOUBLE> t23("t23",t_sten,s_sten,GH,NO_COMM);
   GridFunction(DIM)<DOUBLE> t31("t31",t_sten,s_sten,GH,NO_COMM);
   GridFunction(DIM)<DOUBLE> t32("t32",t_sten,s_sten,GH,NO_COMM);
   GridFunction(DIM)<DOUBLE> t33("t33",t_sten,s_sten,GH,NO_COMM);
   GridFunction(DIM)<DOUBLE> temp1("temp1",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> temp2("temp2",t_sten,s_sten,GH);
   GridFunction(DIM)<DOUBLE> temp3("temp3",t_sten,s_sten,GH);

/* Note that these arrays should only be one dimensional arrays */

   GridFunction(DIM)<INTEGER> m("m",t_sten,s_sten,GH,NO_COMM);
   GridFunction(DIM)<DOUBLE> d("d",t_sten,s_sten,GH,NO_COMM);

   DOUBLE hamnorm,mom1norm,mom2norm,mom3norm;

   INTEGER me = MY_PROC(GH);
   INTEGER num = NUM_PROC(GH);


  INTEGER t = 0, l = 0;

  time = (DOUBLE) 0.0;
  

/*******                      INITIALIZATION                          ******/
   if (me == 0) cout << "...Initialize data" << "\n" <<flush;

/*
   forall(g11,t,l,c)
      f_initdata(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(k11(t,l,c)),FDA(k12(t,l,c)),
                FDA(k13(t,l,c)),FDA(k21(t,l,c)),FDA(k22(t,l,c)),
                FDA(k23(t,l,c)),FDA(k31(t,l,c)),FDA(k32(t,l,c)),
                FDA(k33(t,l,c)),bha,FBA(bhr),FDA(bhr),FDA(temp1(t,l,c)),
                bb,&dx,&dy,&dz);
   end_forall
*/
   forall(g11,t,l,c)
            f_linearwave(FBA(g11(t,l,c)),
            FDA(g11(t,l,c)), FDA(g12(t,l,c)), FDA(g13(t,l,c)),
            FDA(g22(t,l,c)), FDA(g23(t,l,c)), FDA(g33(t,l,c)),
            FDA(k11(t,l,c)), FDA(k12(t,l,c)), FDA(k13(t,l,c)), 
            FDA(k21(t,l,c)), FDA(k22(t,l,c)), FDA(k23(t,l,c)), 
            FDA(k31(t,l,c)), FDA(k32(t,l,c)), FDA(k33(t,l,c)),
            FDA(t11(t,l,c)), FDA(t12(t,l,c)), FDA(t13(t,l,c)), 
            FDA(t21(t,l,c)), FDA(t22(t,l,c)), FDA(t23(t,l,c)), 
            FDA(t31(t,l,c)),bb, &dx, &dy, &dz, &amp, &width);
   end_forall

   if (me == 0) cout << "...Generate initial horizon mask" << "\n" <<flush;
   forall(g11,t,l,c)
      f_horizon(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(k11(t,l,c)),FDA(k12(t,l,c)),
                FDA(k13(t,l,c)),FDA(k21(t,l,c)),FDA(k22(t,l,c)),
                FDA(k23(t,l,c)),FDA(k31(t,l,c)),FDA(k32(t,l,c)),
                FDA(k33(t,l,c)),bha,FBA(bhr),FDA(bhr),FDA(mask(t,l,c)),
                bb,&dx,&dy,&dz);

      mask(-1,l,c) = mask(0,l,c);
   end_forall

   forall(g11,t,l,c)
      f_remask(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(k11(t,l,c)),FDA(k12(t,l,c)),
                FDA(k13(t,l,c)),FDA(k21(t,l,c)),FDA(k22(t,l,c)),
                FDA(k23(t,l,c)),FDA(k31(t,l,c)),FDA(k32(t,l,c)),
                FDA(k33(t,l,c)),FDA(mask(t,l,c)),FDA(mask(-1,l,c)));
   end_forall

   if (me == 0) cout << "...Compute excise mask and boundary points" <<
                         "\n" <<flush;
   Sync(mask,t,l);
   forall(mask,t,l,c)
      f_flags(FBA(mask(t,l,c)),&nlmax,xi,yj,zk,fx,fy,fz,
              FDA(mask(t,l,c)),FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)),&nl);
   end_forall

   if (me == 0) cout << "...Compute gauge" << "\n" <<flush;

   forall(mask,t,l,c)
      f_setgauge(FBA(mask(t,l,c)),FDA(mask(t,l,c)),bha,FBA(bhr),FDA(bhr),
                 FDA(alpha(t,l,c)),FDA(beta1(t,l,c)),FDA(beta2(t,l,c)),
                 FDA(beta3(t,l,c)),bb,&dx,&dy,&dz);
   end_forall

   if (me == 0) cout << "...Compute connections and Ricci" << "\n" <<flush;

   forall(g11,t,l,c)
      f_inverse(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                FDA(gup33(t,l,c)),FDA(temp1(t,l,c)),FDA(mask(t,l,c)));
   end_forall

   Sync(g11,t,l);
   Sync(g12,t,l);
   Sync(g13,t,l);
   Sync(g22,t,l);
   Sync(g23,t,l);
   Sync(g33,t,l);
   Sync(gup11,t,l);
   Sync(gup12,t,l);
   Sync(gup13,t,l);
   Sync(gup22,t,l);
   Sync(gup23,t,l);
   Sync(gup33,t,l);
	n1 = Norm2(g11,t,l);
	n2 = Norm2(g12,t,l);
	n3 = Norm2(g13,t,l);
	n4 = Norm2(g22,t,l);
	n5 = Norm2(g23,t,l);
	n6 = Norm2(g33,t,l);

	if (me == 0) {
	cout << "Norms of 3-metric" << "\n" << flush ;
	cout << n1 << " " << n2 << " " << n3 << " " << n4 << " " << n5 << " " << n6
        << "\n" << flush;
	}

   forall(g11,t,l,c)
      f_connect(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                FDA(gup33(t,l,c)),
                FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)), 
                FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)), 
                FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)), 
                FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)), 
                FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                &dx,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));

   end_forall

   Sync(c111,t,l);
   Sync(c112,t,l);
   Sync(c113,t,l);
   Sync(c122,t,l);
   Sync(c123,t,l);
   Sync(c133,t,l);
   Sync(c211,t,l);
   Sync(c212,t,l);
   Sync(c213,t,l);
   Sync(c222,t,l);
   Sync(c223,t,l);
   Sync(c233,t,l);
   Sync(c311,t,l);
   Sync(c312,t,l);
   Sync(c313,t,l);
   Sync(c322,t,l);
   Sync(c323,t,l);
   Sync(c333,t,l);


	n1 = Norm2(c111,t,l);
	n2 = Norm2(c112,t,l);
	n3 = Norm2(c113,t,l);
	n4 = Norm2(c122,t,l);
	n5 = Norm2(c123,t,l);
	n6 = Norm2(c133,t,l);
	n7 = Norm2(c211,t,l);
	n8 = Norm2(c212,t,l);
	n9 = Norm2(c213,t,l);
	n10 = Norm2(c222,t,l);
	n11 = Norm2(c223,t,l);
	n12 = Norm2(c233,t,l);
	n13 = Norm2(c311,t,l);
	n14 = Norm2(c312,t,l);
	n15 = Norm2(c313,t,l);
	n16 = Norm2(c322,t,l);
	n17 = Norm2(c323,t,l);
	n18 = Norm2(c333,t,l);

	if (me == 0) {
	cout << "Norms of connections \n" << flush ;
	cout << n1 << " " << n2 << " " << n3 << " " << n4 << " " << n5 << " " << n6
        << "\n" << flush;
	cout << n7 << " " << n8 << " " << n9 << " " << n10 << " " << n11 << " " << n12
        << "\n" << flush;
	cout << n13 << " " << n14 << " " << n15 << " " << n16 << " " << n17 << " " << n18
        << "\n" << flush;
	}

   forall(g11,t,l,c)

      f_riccidd(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                FDA(gup33(t,l,c)),
                FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)), 
                FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)), 
                FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)), 
                FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)), 
                FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                &dx,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
   end_forall
   Sync(r11,t,l);
   Sync(r12,t,l);
   Sync(r13,t,l);
   Sync(r22,t,l);
   Sync(r23,t,l);
   Sync(r33,t,l);
	n1 = Norm2(r11,t,l);
	n2 = Norm2(r12,t,l);
	n3 = Norm2(r13,t,l);
	n4 = Norm2(r22,t,l);
	n5 = Norm2(r23,t,l);
	n6 = Norm2(r33,t,l);
	if (me == 0) {
	cout << "Ricci tensor component norms \n" << flush ;
	cout << n1 << " " << n2 << " " << n3 << " " << n4 << " " << n5 << " " << n6
        << "\n" << flush;
	}

   f_timestep(&dx,&dy,&dz,&dt);
   dtsave = dt;
   timesave = time;
   itsave = it;
   dt2= -dt/8.0;
   time = time + dt2;

   forall(g11,t,l,c)
      g11(-1,l,c) = g11(0,l,c);
      g12(-1,l,c) = g12(0,l,c);
      g13(-1,l,c) = g13(0,l,c);
      g22(-1,l,c) = g22(0,l,c);
      g23(-1,l,c) = g23(0,l,c);
      g33(-1,l,c) = g33(0,l,c);
   end_forall

   if (me == 0)
      cout << "... Take Euler step backwards to generate n+1 level"
           << "\n" <<flush;

   dt = (DOUBLE)0.5 * dt2;

   Sync(beta1,t,l);
   Sync(beta2,t,l);
   Sync(beta3,t,l);
   forall(g11,t,l,c)

      f_evolve_g(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                FDA(gup33(t,l,c)),
                FDA(g11(-1,l,c)),FDA(g12(-1,l,c)),FDA(g13(-1,l,c)),
                FDA(g22(-1,l,c)),FDA(g23(-1,l,c)), FDA(g33(-1,l,c)),
                FDA(k11(t,l,c)),FDA(k12(t,l,c)),FDA(k13(t,l,c)),
                FDA(k21(t,l,c)),FDA(k22(t,l,c)),FDA(k23(t,l,c)),
                FDA(k31(t,l,c)),FDA(k32(t,l,c)), FDA(k33(t,l,c)),
                FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                FDA(beta1(t,l,c)),FDA(beta2(t,l,c)),FDA(beta3(t,l,c)),
                FDA(alpha(t,l,c)),
                FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                FDA(temp1(t,l,c)),FDA(temp2(t,l,c)),FDA(temp3(t,l,c)),
                FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                &dx,&dt,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
   end_forall


   Sync(alpha,t,l);
   forall(g11,t,l,c)
      f_evolve_ka(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                FDA(gup33(t,l,c)),
                FDA(k11(-1,l,c)),FDA(k12(-1,l,c)),FDA(k13(-1,l,c)),
                FDA(k21(-1,l,c)),FDA(k22(-1,l,c)),FDA(k23(-1,l,c)),
                FDA(k31(-1,l,c)),FDA(k32(-1,l,c)),FDA(k33(-1,l,c)),
                FDA(k11(t,l,c)),FDA(k12(t,l,c)),FDA(k13(t,l,c)),
                FDA(k21(t,l,c)),FDA(k22(t,l,c)),FDA(k23(t,l,c)),
                FDA(k31(t,l,c)),FDA(k32(t,l,c)),FDA(k33(t,l,c)),
                FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                FDA(beta1(t,l,c)),FDA(beta2(t,l,c)),FDA(beta3(t,l,c)),
                FDA(alpha(t,l,c)),
                FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                FDA(temp1(t,l,c)),FDA(temp2(t,l,c)),FDA(temp3(t,l,c)),
                FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                &dx,&dt,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
   end_forall

   Sync(temp1,t,l);
   Sync(temp2,t,l);
   Sync(temp3,t,l);
   Sync(k11,t,l);
   Sync(k12,t,l);
   Sync(k13,t,l);
   Sync(k21,t,l);
   Sync(k22,t,l);
   Sync(k23,t,l);
   Sync(k31,t,l);
   Sync(k32,t,l);
   Sync(k33,t,l);
   forall(g11,t,l,c)
      f_evolve_kb(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                FDA(gup33(t,l,c)),
                FDA(k11(-1,l,c)),FDA(k12(-1,l,c)),FDA(k13(-1,l,c)),
                FDA(k21(-1,l,c)),FDA(k22(-1,l,c)),FDA(k23(-1,l,c)),
                FDA(k31(-1,l,c)),FDA(k32(-1,l,c)),FDA(k33(-1,l,c)),
                FDA(k11(t,l,c)),FDA(k12(t,l,c)),FDA(k13(t,l,c)),
                FDA(k21(t,l,c)),FDA(k22(t,l,c)),FDA(k23(t,l,c)),
                FDA(k31(t,l,c)),FDA(k32(t,l,c)),FDA(k33(t,l,c)),
                FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                FDA(beta1(t,l,c)),FDA(beta2(t,l,c)),FDA(beta3(t,l,c)),
                FDA(alpha(t,l,c)),
                FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                FDA(temp1(t,l,c)),FDA(temp2(t,l,c)),FDA(temp3(t,l,c)),
                FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                &dx,&dt,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
   end_forall

   forall(g11,t,l,c)
      temp1(t,l,c) = g11(t,l,c);
      g11(t,l,c)   = g11(-1,l,c);
      g11(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = g12(t,l,c);
      g12(t,l,c)   = g12(-1,l,c);
      g12(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = g13(t,l,c);
      g13(t,l,c)   = g13(-1,l,c);
      g13(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = g22(t,l,c);
      g22(t,l,c)   = g22(-1,l,c);
      g22(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = g23(t,l,c);
      g23(t,l,c)   = g23(-1,l,c);
      g23(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = g33(t,l,c);
      g33(t,l,c)   = g33(-1,l,c);
      g33(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = k11(t,l,c);
      k11(t,l,c)   = k11(-1,l,c);
      k11(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = k12(t,l,c);
      k12(t,l,c)   = k12(-1,l,c);
      k12(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = k13(t,l,c);
      k13(t,l,c)   = k13(-1,l,c);
      k13(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = k21(t,l,c);
      k21(t,l,c)   = k21(-1,l,c);
      k21(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = k22(t,l,c);
      k22(t,l,c)   = k22(-1,l,c);
      k22(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = k23(t,l,c);
      k23(t,l,c)   = k23(-1,l,c);
      k23(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = k31(t,l,c);
      k31(t,l,c)   = k31(-1,l,c);
      k31(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = k32(t,l,c);
      k32(t,l,c)   = k32(-1,l,c);
      k32(-1,l,c)  = temp1(t,l,c);

      temp1(t,l,c) = k33(t,l,c);
      k33(t,l,c)   = k33(-1,l,c);
      k33(-1,l,c)  = temp1(t,l,c);

   end_forall


   dt2 = 2.0 * dt;
   if (me == 0) cout << "...Complete initial backward timestep" << "\n" <<flush;

   for (it=2;it<=8;it++)
   {
   time = time + dt2;
      forall(g11,t,l,c)

         f_inverse(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                   FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                   FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                   FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                   FDA(gup33(t,l,c)),FDA(temp1(t,l,c)),FDA(mask(t,l,c)));
      end_forall

      Sync(g11,t,l);
      Sync(g12,t,l);
      Sync(g13,t,l);
      Sync(g22,t,l);
      Sync(g23,t,l);
      Sync(g33,t,l);
      forall(g11,t,l,c)
         f_connect(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                   FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                   FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                   FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                   FDA(gup33(t,l,c)),
                   FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                   FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                   FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                   FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                   FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                   FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                   FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                   FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                   FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                   FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                   FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                   &dx,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));

   
         end_forall

         Sync(c111,t,l);
         Sync(c112,t,l);
         Sync(c113,t,l);
         Sync(c122,t,l);
         Sync(c123,t,l);
         Sync(c133,t,l);
         Sync(c211,t,l);
         Sync(c212,t,l);
         Sync(c213,t,l);
         Sync(c222,t,l);
         Sync(c223,t,l);
         Sync(c233,t,l);
         Sync(c311,t,l);
         Sync(c312,t,l);
         Sync(c313,t,l);
         Sync(c322,t,l);
         Sync(c323,t,l);
         Sync(c333,t,l);
         forall(g11,t,l,c)
            f_riccidd(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                      FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                      FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                      FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                      FDA(gup33(t,l,c)),
                      FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                      FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                      FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                      FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                      FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                      FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                      FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                      FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                      FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                      FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                      FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                      &dx,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
         end_forall
   		Sync(r11,t,l);
   		Sync(r12,t,l);
   		Sync(r13,t,l);
   		Sync(r22,t,l);
   		Sync(r23,t,l);
   		Sync(r33,t,l);

         Sync(beta1,t,l);
         Sync(beta2,t,l);
         Sync(beta3,t,l);
         forall(g11,t,l,c)
            f_evolve_g(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                      FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                      FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                      FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                      FDA(gup33(t,l,c)),
                      FDA(g11(-1,l,c)),FDA(g12(-1,l,c)),
                      FDA(g13(-1,l,c)),FDA(g22(-1,l,c)),FDA(g23(-1,l,c)),
                      FDA(g33(-1,l,c)),
                      FDA(k11(t,l,c)),FDA(k12(t,l,c)),
                      FDA(k13(t,l,c)),FDA(k21(t,l,c)),FDA(k22(t,l,c)),
                      FDA(k23(t,l,c)),FDA(k31(t,l,c)),FDA(k32(t,l,c)),
                      FDA(k33(t,l,c)),
                      FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                      FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                      FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                      FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                      FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                      FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                      FDA(beta1(t,l,c)),FDA(beta2(t,l,c)),FDA(beta3(t,l,c)),
                      FDA(alpha(t,l,c)),
                      FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                      FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                      FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                      FDA(temp1(t,l,c)),FDA(temp2(t,l,c)),FDA(temp3(t,l,c)),
                      FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                      FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                      &dx,&dt2,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
   
      end_forall
      Sync(alpha,t,l);
      forall(g11,t,l,c)
         f_evolve_ka(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                   FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                   FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                   FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                   FDA(gup33(t,l,c)),
                   FDA(k11(-1,l,c)),FDA(k12(-1,l,c)),FDA(k13(-1,l,c)),
                   FDA(k21(-1,l,c)),FDA(k22(-1,l,c)),FDA(k23(-1,l,c)),
                   FDA(k31(-1,l,c)),FDA(k32(-1,l,c)),FDA(k33(-1,l,c)),
                   FDA(k11(t,l,c)),FDA(k12(t,l,c)),FDA(k13(t,l,c)),
                   FDA(k21(t,l,c)),FDA(k22(t,l,c)),FDA(k23(t,l,c)),
                   FDA(k31(t,l,c)),FDA(k32(t,l,c)),FDA(k33(t,l,c)),
                   FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                   FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                   FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                   FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                   FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                   FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                   FDA(beta1(t,l,c)),FDA(beta2(t,l,c)),FDA(beta3(t,l,c)),
                   FDA(alpha(t,l,c)),
                   FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                   FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                   FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                   FDA(temp1(t,l,c)),FDA(temp2(t,l,c)),FDA(temp3(t,l,c)),
                   FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                   FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                   &dx,&dt2,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
   
      end_forall
      Sync(temp1,t,l);
      Sync(temp2,t,l);
      Sync(temp3,t,l);
      Sync(k11,t,l);
      Sync(k12,t,l);
      Sync(k13,t,l);
      Sync(k21,t,l);
      Sync(k22,t,l);
      Sync(k23,t,l);
      Sync(k31,t,l);
      Sync(k32,t,l);
      Sync(k33,t,l);
      forall(g11,t,l,c)
         f_evolve_kb(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                   FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                   FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                   FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                   FDA(gup33(t,l,c)),
                   FDA(k11(-1,l,c)),FDA(k12(-1,l,c)),FDA(k13(-1,l,c)),
                   FDA(k21(-1,l,c)),FDA(k22(-1,l,c)),FDA(k23(-1,l,c)),
                   FDA(k31(-1,l,c)),FDA(k32(-1,l,c)),FDA(k33(-1,l,c)),
                   FDA(k11(t,l,c)),FDA(k12(t,l,c)),FDA(k13(t,l,c)),
                   FDA(k21(t,l,c)),FDA(k22(t,l,c)),FDA(k23(t,l,c)),
                   FDA(k31(t,l,c)),FDA(k32(t,l,c)),FDA(k33(t,l,c)),
                   FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                   FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                   FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                   FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                   FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                   FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                   FDA(beta1(t,l,c)),FDA(beta2(t,l,c)),FDA(beta3(t,l,c)),
                   FDA(alpha(t,l,c)),
                   FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                   FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                   FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                   FDA(temp1(t,l,c)),FDA(temp2(t,l,c)),FDA(temp3(t,l,c)),
                   FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                   FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                   &dx,&dt2,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
      end_forall


      forall(g11,t,l,c)
            temp1(t,l,c) = g11(t,l,c);
            g11(t,l,c)   = g11(-1,l,c);
            g11(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = g12(t,l,c);
            g12(t,l,c)   = g12(-1,l,c);
            g12(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = g13(t,l,c);
            g13(t,l,c)   = g13(-1,l,c);
            g13(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = g22(t,l,c);
            g22(t,l,c)   = g22(-1,l,c);
            g22(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = g23(t,l,c);
            g23(t,l,c)   = g23(-1,l,c);
            g23(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = g33(t,l,c);
            g33(t,l,c)   = g33(-1,l,c);
            g33(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = k11(t,l,c);
            k11(t,l,c)   = k11(-1,l,c);
            k11(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = k12(t,l,c);
            k12(t,l,c)   = k12(-1,l,c);
            k12(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = k13(t,l,c);
            k13(t,l,c)   = k13(-1,l,c);
            k13(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = k21(t,l,c);
            k21(t,l,c)   = k21(-1,l,c);
            k21(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = k22(t,l,c);
            k22(t,l,c)   = k22(-1,l,c);
            k22(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = k23(t,l,c);
            k23(t,l,c)   = k23(-1,l,c);
            k23(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = k31(t,l,c);
            k31(t,l,c)   = k31(-1,l,c);
            k31(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = k32(t,l,c);
            k32(t,l,c)   = k32(-1,l,c);
            k32(-1,l,c)  = temp1(t,l,c);
      
            temp1(t,l,c) = k33(t,l,c);
            k33(t,l,c)   = k33(-1,l,c);
            k33(-1,l,c)  = temp1(t,l,c);
      
         end_forall
   }
   

   time = timesave;
   dt = dtsave;

/* swap back the 8th... */

   forall(g11,t,l,c)

         temp1(t,l,c) = g11(t,l,c);
         g11(t,l,c)   = g11(-1,l,c);
         g11(-1,l,c)  = temp1(t,l,c);
  
         temp1(t,l,c) = g12(t,l,c);
         g12(t,l,c)   = g12(-1,l,c);
         g12(-1,l,c)  = temp1(t,l,c);
  
         temp1(t,l,c) = g13(t,l,c);
         g13(t,l,c)   = g13(-1,l,c);
         g13(-1,l,c)  = temp1(t,l,c);
  
         temp1(t,l,c) = g22(t,l,c);
         g22(t,l,c)   = g22(-1,l,c);
         g22(-1,l,c)  = temp1(t,l,c);
  
         temp1(t,l,c) = g23(t,l,c);
         g23(t,l,c)   = g23(-1,l,c);
         g23(-1,l,c)  = temp1(t,l,c);
  
         temp1(t,l,c) = g33(t,l,c);
         g33(t,l,c)   = g33(-1,l,c);
         g33(-1,l,c)  = temp1(t,l,c);
  
         temp1(t,l,c) = k11(t,l,c);
         k11(t,l,c)   = k11(-1,l,c);
         k11(-1,l,c)  = temp1(t,l,c);
  
         temp1(t,l,c) = k12(t,l,c);
         k12(t,l,c)   = k12(-1,l,c);
         k12(-1,l,c)  = temp1(t,l,c);
  
         temp1(t,l,c) = k13(t,l,c);
         k13(t,l,c)   = k13(-1,l,c);
         k13(-1,l,c)  = temp1(t,l,c);

         temp1(t,l,c) = k21(t,l,c);
         k21(t,l,c)   = k21(-1,l,c);
         k21(-1,l,c)  = temp1(t,l,c);

         temp1(t,l,c) = k22(t,l,c);
         k22(t,l,c)   = k22(-1,l,c);
         k22(-1,l,c)  = temp1(t,l,c);

         temp1(t,l,c) = k23(t,l,c);
         k23(t,l,c)   = k23(-1,l,c);
         k23(-1,l,c)  = temp1(t,l,c);

         temp1(t,l,c) = k31(t,l,c);
         k31(t,l,c)   = k31(-1,l,c);
         k31(-1,l,c)  = temp1(t,l,c);

         temp1(t,l,c) = k32(t,l,c);
         k32(t,l,c)   = k32(-1,l,c);
         k32(-1,l,c)  = temp1(t,l,c);

         temp1(t,l,c) = k33(t,l,c);
         k33(t,l,c)   = k33(-1,l,c);
         k33(-1,l,c)  = temp1(t,l,c);

   end_forall

/* Restore times and time levels */

   dt = dtsave;
   time = timesave;
   it = itsave;

/*
   forall(g11,t,l,c)
      f_initdata(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(k11(t,l,c)),FDA(k12(t,l,c)),
                FDA(k13(t,l,c)),FDA(k21(t,l,c)),FDA(k22(t,l,c)),
                FDA(k23(t,l,c)),FDA(k31(t,l,c)),FDA(k32(t,l,c)),
                FDA(k33(t,l,c)),bha,FBA(bhr),FDA(bhr),FDA(temp1(t,l,c)),
                bb,&dx,&dy,&dz);
   end_forall
*/
   forall(g11,t,l,c)
            f_linearwave(FBA(g11(t,l,c)),
            FDA(g11(t,l,c)), FDA(g12(t,l,c)), FDA(g13(t,l,c)), 
            FDA(g22(t,l,c)), FDA(g23(t,l,c)), FDA(g33(t,l,c)),
            FDA(k11(t,l,c)), FDA(k12(t,l,c)), FDA(k13(t,l,c)), 
            FDA(k21(t,l,c)), FDA(k22(t,l,c)), FDA(k23(t,l,c)), 
            FDA(k31(t,l,c)), FDA(k32(t,l,c)), FDA(k33(t,l,c)),
            FDA(t11(t,l,c)), FDA(t12(t,l,c)), FDA(t13(t,l,c)), 
            FDA(t21(t,l,c)), FDA(t22(t,l,c)), FDA(t23(t,l,c)), 
            FDA(t31(t,l,c)),bb, &dx, &dy, &dz, &amp, &width);
   end_forall
   forall(g11,t,l,c)
      f_inverse(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                FDA(gup33(t,l,c)),FDA(temp1(t,l,c)),FDA(mask(t,l,c)));
   end_forall
   Sync(g11,t,l);
   Sync(g12,t,l);
   Sync(g13,t,l);
   Sync(g22,t,l);
   Sync(g23,t,l);
   Sync(g33,t,l);
     forall(g11,t,l,c)
      f_connect(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                FDA(gup33(t,l,c)),
                FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                &dx,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
   end_forall

   Sync(c111,t,l);
   Sync(c112,t,l);
   Sync(c113,t,l);
   Sync(c122,t,l);
   Sync(c123,t,l);
   Sync(c133,t,l);
   Sync(c211,t,l);
   Sync(c212,t,l);
   Sync(c213,t,l);
   Sync(c222,t,l);
   Sync(c223,t,l);
   Sync(c233,t,l);
   Sync(c311,t,l);
   Sync(c312,t,l);
   Sync(c313,t,l);
   Sync(c322,t,l);
   Sync(c323,t,l);
   Sync(c333,t,l);
     forall(g11,t,l,c)
      f_riccidd(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                FDA(gup33(t,l,c)),
                FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                &dx,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
   end_forall
   Sync(r11,t,l);
   Sync(r12,t,l);
   Sync(r13,t,l);
   Sync(r22,t,l);
   Sync(r23,t,l);
   Sync(r33,t,l);

/* Begin main loop, time-centered, space-centered. */
   if (me == 0) cout << "...Begin main loop, forward in time" << "\n" <<flush;
   if (me == 0)
      cout << "Iter Time HamErr Mom1Err Mom2err Mom3err" << "\n" << flush;

   forall(mask,t,l,c)
      f_constraint(FBA(mask(t,l,c)),FDA(mask(t,l,c)),
                   FDA(k11(t,l,c)),FDA(k12(t,l,c)),FDA(k13(t,l,c)),
                   FDA(k21(t,l,c)),FDA(k22(t,l,c)),FDA(k23(t,l,c)),
                   FDA(k31(t,l,c)),FDA(k32(t,l,c)),FDA(k33(t,l,c)),
                   &hamnorm,&mom1norm,
                   &mom2norm,&mom3norm,
                   FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                   FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                   FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                   FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                   FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                   FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                   FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                   FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                   FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),FDA(gup13(t,l,c)),
                   FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),FDA(gup33(t,l,c)),
                   FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                   FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                   FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                   FDA(temp1(t,l,c)),&dx,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
   end_forall
	Sync(temp1,t,l);
   hamnorm  = Norm2(temp1 ,t,l);
   mom1norm = Norm2(t31 ,t,l);
   mom2norm = Norm2(t32 ,t,l);
   mom3norm = Norm2(t33 ,t,l);
   it = 0; 
   if (me == 0){
   cout << it << " " << time << " " << hamnorm << " " <<
         mom1norm << " "<< mom2norm << " " << mom3norm << "\n" << flush;
   }
   while (hamnorm < tol && it < nt)
   {
   it ++;
      forall(mask,t,l,c)
         f_setgauge(FBA(mask(t,l,c)),FDA(mask(t,l,c)),bha,FBA(bhr),FDA(bhr),
                    FDA(alpha(t,l,c)),FDA(beta1(t,l,c)),FDA(beta2(t,l,c)),
                    FDA(beta3(t,l,c)),bb,&dx,&dy,&dz);
      end_forall

      f_timestep(&dx,&dy,&dz,&dt);

      time+= dt;

      Sync(beta1,t,l);
      Sync(beta2,t,l);
      Sync(beta3,t,l);
      forall(g11,t,l,c)
         f_evolve_g(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                   FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                   FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                   FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                   FDA(gup33(t,l,c)),
                   FDA(g11(-1,l,c)),FDA(g12(-1,l,c)),
                   FDA(g13(-1,l,c)),FDA(g22(-1,l,c)),FDA(g23(-1,l,c)),
                   FDA(g33(-1,l,c)),
                   FDA(k11(t,l,c)),FDA(k12(t,l,c)),
                   FDA(k13(t,l,c)),FDA(k21(t,l,c)),FDA(k22(t,l,c)),
                   FDA(k23(t,l,c)),FDA(k31(t,l,c)),FDA(k32(t,l,c)),
                   FDA(k33(t,l,c)),
                   FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                   FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                   FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                   FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                   FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                   FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                   FDA(beta1(t,l,c)),FDA(beta2(t,l,c)),FDA(beta3(t,l,c)),
                   FDA(alpha(t,l,c)),
                   FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                   FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                   FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                   FDA(temp1(t,l,c)),FDA(temp2(t,l,c)),FDA(temp3(t,l,c)),
                   FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                   FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                   &dx,&dt,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
      end_forall

      forall(g11,t,l,c)
         f_evolve_ka(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                     FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                     FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                     FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                     FDA(gup33(t,l,c)),
                     FDA(k11(-1,l,c)),FDA(k12(-1,l,c)),FDA(k13(-1,l,c)),
                     FDA(k21(-1,l,c)),FDA(k22(-1,l,c)),FDA(k23(-1,l,c)),
                     FDA(k31(-1,l,c)),FDA(k32(-1,l,c)),FDA(k33(-1,l,c)),
                     FDA(k11(t,l,c)),FDA(k12(t,l,c)),FDA(k13(t,l,c)),
                     FDA(k21(t,l,c)),FDA(k22(t,l,c)),FDA(k23(t,l,c)),
                     FDA(k31(t,l,c)),FDA(k32(t,l,c)),FDA(k33(t,l,c)),
                     FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                     FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                     FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                     FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                     FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                     FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                     FDA(beta1(t,l,c)),FDA(beta2(t,l,c)),FDA(beta3(t,l,c)),
                     FDA(alpha(t,l,c)),
                     FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                     FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                     FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                     FDA(temp1(t,l,c)),FDA(temp2(t,l,c)),FDA(temp3(t,l,c)),
                     FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                     FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                     &dx,&dt,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
      end_forall

   Sync(temp1,t,l);
   Sync(temp2,t,l);
   Sync(temp3,t,l);
   Sync(k11,t,l);
   Sync(k12,t,l);
   Sync(k13,t,l);
   Sync(k21,t,l);
   Sync(k22,t,l);
   Sync(k23,t,l);
   Sync(k31,t,l);
   Sync(k32,t,l);
   Sync(k33,t,l);
   forall(g11,t,l,c)
      f_evolve_kb(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                FDA(gup33(t,l,c)),
                FDA(k11(-1,l,c)),FDA(k12(-1,l,c)),FDA(k13(-1,l,c)),
                FDA(k21(-1,l,c)),FDA(k22(-1,l,c)),FDA(k23(-1,l,c)),
                FDA(k31(-1,l,c)),FDA(k32(-1,l,c)),FDA(k33(-1,l,c)),
                FDA(k11(t,l,c)),FDA(k12(t,l,c)),FDA(k13(t,l,c)),
                FDA(k21(t,l,c)),FDA(k22(t,l,c)),FDA(k23(t,l,c)),
                FDA(k31(t,l,c)),FDA(k32(t,l,c)),FDA(k33(t,l,c)),
                FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                FDA(beta1(t,l,c)),FDA(beta2(t,l,c)),FDA(beta3(t,l,c)),
                FDA(alpha(t,l,c)),
                FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                FDA(temp1(t,l,c)),FDA(temp2(t,l,c)),FDA(temp3(t,l,c)),
                FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                &dx,&dt,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
   end_forall

      forall(g11,t,l,c)
         temp1(t,l,c) = g11(t,l,c);
         g11(t,l,c)   = g11(-1,l,c);
         g11(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = g12(t,l,c);
         g12(t,l,c)   = g12(-1,l,c);
         g12(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = g13(t,l,c);
         g13(t,l,c)   = g13(-1,l,c);
         g13(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = g22(t,l,c);
         g22(t,l,c)   = g22(-1,l,c);
         g22(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = g23(t,l,c);
         g23(t,l,c)   = g23(-1,l,c);
         g23(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = g33(t,l,c);
         g33(t,l,c)   = g33(-1,l,c);
         g33(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = k11(t,l,c);
         k11(t,l,c)   = k11(-1,l,c);
         k11(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = k12(t,l,c);
         k12(t,l,c)   = k12(-1,l,c);
         k12(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = k13(t,l,c);
         k13(t,l,c)   = k13(-1,l,c);
         k13(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = k21(t,l,c);
         k21(t,l,c)   = k21(-1,l,c);
         k21(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = k22(t,l,c);
         k22(t,l,c)   = k22(-1,l,c);
         k22(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = k23(t,l,c);
         k23(t,l,c)   = k23(-1,l,c);
         k23(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = k31(t,l,c);
         k31(t,l,c)   = k31(-1,l,c);
         k31(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = k32(t,l,c);
         k32(t,l,c)   = k32(-1,l,c);
         k32(-1,l,c)  = temp1(t,l,c);
   
         temp1(t,l,c) = k33(t,l,c);
         k33(t,l,c)   = k33(-1,l,c);
         k33(-1,l,c)  = temp1(t,l,c);
   
         mask(-1,l,c) = mask(t,l,c);
      end_forall
      forall(g11,t,l,c)
   
         f_horizon(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                FDA(g33(t,l,c)),FDA(k11(t,l,c)),FDA(k12(t,l,c)),
                FDA(k13(t,l,c)),FDA(k21(t,l,c)),FDA(k22(t,l,c)),
                FDA(k23(t,l,c)),FDA(k31(t,l,c)),FDA(k32(t,l,c)),
                FDA(k33(t,l,c)),bha,FBA(bhr),FDA(bhr),FDA(mask(t,l,c)),
                bb,&dx,&dy,&dz);

         f_remask(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                   FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                   FDA(g33(t,l,c)),FDA(k11(t,l,c)),FDA(k12(t,l,c)),
                   FDA(k13(t,l,c)),FDA(k21(t,l,c)),FDA(k22(t,l,c)),
                   FDA(k23(t,l,c)),FDA(k31(t,l,c)),FDA(k32(t,l,c)),
                   FDA(k33(t,l,c)),FDA(mask(t,l,c)),FDA(mask(-1,l,c)));

         f_remask(FBA(g11(-1,l,c)),FDA(g11(-1,l,c)),FDA(g12(-1,l,c)),
                   FDA(g13(-1,l,c)),FDA(g22(-1,l,c)),FDA(g23(-1,l,c)),
                   FDA(g33(-1,l,c)),FDA(k11(-1,l,c)),FDA(k12(-1,l,c)),
                   FDA(k13(-1,l,c)),FDA(k21(-1,l,c)),FDA(k22(-1,l,c)),
                   FDA(k23(-1,l,c)),FDA(k31(-1,l,c)),FDA(k32(-1,l,c)),
                   FDA(k33(-1,l,c)),FDA(mask(t,l,c)),FDA(mask(-1,l,c)));

      end_forall
      Sync(mask,t,l);
      forall(mask,t,l,c)
         f_flags(FBA(mask(t,l,c)),&nlmax,xi,yj,zk,fx,fy,fz,
              FDA(mask(t,l,c)),FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)),&nl);
      end_forall
      forall(g11,t,l,c)
         f_inverse(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                   FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                   FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                   FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                   FDA(gup33(t,l,c)),FDA(temp1(t,l,c)),FDA(mask(t,l,c)));
      end_forall
   
      Sync(g11,t,l);
      Sync(g12,t,l);
      Sync(g13,t,l);
      Sync(g22,t,l);
      Sync(g23,t,l);
      Sync(g33,t,l);
        forall(g11,t,l,c)
         f_connect(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                   FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                   FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                   FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                   FDA(gup33(t,l,c)),
                   FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                   FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                   FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                   FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                   FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                   FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                   FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                   FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                   FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                   FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                   FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                   &dx,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));

      end_forall
      Sync(c111,t,l);
      Sync(c112,t,l);
      Sync(c113,t,l);
      Sync(c122,t,l);
      Sync(c123,t,l);
      Sync(c133,t,l);
      Sync(c211,t,l);
      Sync(c212,t,l);
      Sync(c213,t,l);
      Sync(c222,t,l);
      Sync(c223,t,l);
      Sync(c233,t,l);
      Sync(c311,t,l);
      Sync(c312,t,l);
      Sync(c313,t,l);
      Sync(c322,t,l);
      Sync(c323,t,l);
      Sync(c333,t,l);
        forall(g11,t,l,c)
         f_riccidd(FBA(g11(t,l,c)),FDA(g11(t,l,c)),FDA(g12(t,l,c)),
                   FDA(g13(t,l,c)),FDA(g22(t,l,c)),FDA(g23(t,l,c)),
                   FDA(g33(t,l,c)),FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),
                   FDA(gup13(t,l,c)),FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),
                   FDA(gup33(t,l,c)),
                   FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                   FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                   FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                   FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                   FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                   FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                   FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                   FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                   FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                   FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                   FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                   &dx,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
   
      end_forall
		Sync(r11,t,l);
		Sync(r12,t,l);
		Sync(r13,t,l);
		Sync(r22,t,l);
		Sync(r23,t,l);
		Sync(r33,t,l);

        forall(g11,t,l,c)
         f_constraint(FBA(mask(t,l,c)),FDA(mask(t,l,c)),
                      FDA(k11(t,l,c)),FDA(k12(t,l,c)),FDA(k13(t,l,c)),
                      FDA(k21(t,l,c)),FDA(k22(t,l,c)),FDA(k23(t,l,c)),
                      FDA(k31(t,l,c)),FDA(k32(t,l,c)),FDA(k33(t,l,c)),
                      &hamnorm,&mom1norm,
                      &mom2norm,&mom3norm,
                      FDA(c111(t,l,c)),FDA(c112(t,l,c)),FDA(c113(t,l,c)),
                      FDA(c122(t,l,c)),FDA(c123(t,l,c)),FDA(c133(t,l,c)),
                      FDA(c211(t,l,c)),FDA(c212(t,l,c)),FDA(c213(t,l,c)),
                      FDA(c222(t,l,c)),FDA(c223(t,l,c)),FDA(c233(t,l,c)),
                      FDA(c311(t,l,c)),FDA(c312(t,l,c)),FDA(c313(t,l,c)),
                      FDA(c322(t,l,c)),FDA(c323(t,l,c)),FDA(c333(t,l,c)),
                      FDA(r11(t,l,c)),FDA(r12(t,l,c)),FDA(r13(t,l,c)),
                      FDA(r22(t,l,c)),FDA(r23(t,l,c)),FDA(r33(t,l,c)),
                      FDA(gup11(t,l,c)),FDA(gup12(t,l,c)),FDA(gup13(t,l,c)),
                      FDA(gup22(t,l,c)),FDA(gup23(t,l,c)),FDA(gup33(t,l,c)),
                      FDA(t11(t,l,c)),FDA(t12(t,l,c)),FDA(t13(t,l,c)),
                      FDA(t21(t,l,c)),FDA(t22(t,l,c)),FDA(t23(t,l,c)),
                      FDA(t31(t,l,c)),FDA(t32(t,l,c)),FDA(t33(t,l,c)),
                      FDA(temp1(t,l,c)),&dx,&nl,FBA(m(t,l,c)),FDA(m(t,l,c)),FDA(d(t,l,c)));
      end_forall
   
      hamnorm = Norm2(temp1 ,t,l);
      mom1norm = Norm2(t31 ,t,l);
      mom2norm = Norm2(t32 ,t,l);
      mom3norm = Norm2(t33 ,t,l);
      if (me == 0){
         cout << it << " " << time << " " << hamnorm << " " <<
          mom1norm << " "<< mom2norm << " " << mom3norm << "\n" << flush;
      }


   }
   DOUBLE stop = MPI_Wtime();
   cout << "(" << MY_PROC(GH) << ") Exec Time = " << stop-start << endl;
      
  MPI_Finalize(); 
 }