VARS=STRUCTURED ARRAY(p,u,v,w) OF REAL

#define uxav0(iy,iz) (0.75*u(0,iy,iz)+0.25*u(-1,iy,iz))
#define uxav1(iy,iz) (0.25*u(0,iy,iz)+0.75*u(-1,iy,iz))
#define vyav0(ix,iz) (0.75*v(ix,0,iz)+0.25*v(ix,-1,iz))
#define vyav1(ix,iz) (0.25*v(ix,0,iz)+0.75*v(ix,-1,iz))
#define wzav0(ix,iy) (0.75*w(ix,iy,0)+0.25*w(ix,iy,-1))
#define wzav1(ix,iy) (0.25*w(ix,iy,0)+0.75*w(ix,iy,-1))

SUBROUTINE interp(ARRAY(*,*,*) OF VARS c,f^) WITH c
  f.u(0,0,0)=~+uxav0(0,0)
  f.u(-1,0,0)=~+uxav1(0,0)
  f.u(0,-1,0)=~+0.5*[uxav0(0,0)+uxav0(-1,0)]
  f.u(-1,-1,0)=~+0.5*[uxav1(0,0)+uxav1(-1,0)]
  f.u(0,0,-1)=~+0.5*[uxav0(0,0)+uxav0(0,-1)]
  f.u(-1,0,-1)=~+0.5*[uxav1(0,0)+uxav1(0,-1)]
  f.u(0,-1,-1)=~+0.25*[uxav0(0,0)+uxav0(0,-1)+uxav0(-1,0)+uxav0(-1,-1)]
  f.u(-1,-1,-1)=~+0.25*[uxav1(0,0)+uxav1(0,-1)+uxav1(-1,0)+uxav1(-1,-1)]
  f.v(0,0,0)=~+vyav0(0,0)
  f.v(0,-1,0)=~+vyav1(0,0)
  f.v(-1,0,0)=~+0.5*[vyav0(0,0)+vyav0(-1,0)]
  f.v(-1,-1,0)=~+0.5*[vyav1(0,0)+vyav1(-1,0)]
  f.v(0,0,-1)=~+0.5*[vyav0(0,0)+vyav0(0,-1)]
  f.v(0,-1,-1)=~+0.5*[vyav1(0,0)+vyav1(0,-1)]
  f.v(-1,0,-1)=~+0.25*[vyav0(0,0)+vyav0(0,-1)+vyav0(-1,0)+vyav0(-1,-1)]
  f.v(-1,-1,-1)=~+0.25*[vyav1(0,0)+vyav1(0,-1)+vyav1(-1,0)+vyav1(-1,-1)]
  f.w(0,0,0)=~+wzav0(0,0)
  f.w(0,0,-1)=~+wzav1(0,0)
  f.w(-1,0,0)=~+0.5*[wzav0(0,0)+wzav0(-1,0)]
  f.w(-1,0,-1)=~+0.5*[wzav1(0,0)+wzav1(-1,0)]
  f.w(0,-1,0)=~+0.5*[wzav0(0,0)+wzav0(0,-1)]
  f.w(0,-1,-1)=~+0.5*[wzav1(0,0)+wzav1(0,-1)]
  f.w(-1,-1,0)=~+0.25*[wzav0(0,0)+wzav0(0,-1)+wzav0(-1,0)+wzav0(-1,-1)]
  f.w(-1,-1,-1)=~+0.25*[wzav1(0,0)+wzav1(0,-1)+wzav1(-1,0)+wzav1(-1,-1)]
  f.p(0,0,0)=~+p(0,0,0)
  f.p(-1,0,0)=~+0.5*[p(0,0,0)+p(-1,0,0)]
  f.p(0,-1,0)=~+0.5*[p(0,0,0)+p(0,-1,0)]
  f.p(-1,-1,0)=~+0.25*[p(0,0,0)+p(-1,0,0)+p(0,-1,0)+p(-1,-1,0)]
  f.p(0,0,-1)=~+0.5*[p(0,0,0)+p(0,0,-1)]
  f.p(-1,0,-1)=~+0.25*[p(0,0,0)+p(-1,0,0)+p(0,0,-1)+p(-1,0,-1)]
  f.p(0,-1,-1)=~+0.25*[p(0,0,0)+p(0,-1,0)+p(0,0,-1)+p(0,-1,-1)]
  f.p(-1,-1,-1)=~+0.125*[p(0,0,0)+p(-1,0,0)+p(0,-1,0)+p(-1,-1,0)+
                         p(0,0,-1)+p(-1,0,-1)+p(0,-1,-1)+p(-1,-1,-1)]
END interp

SUBROUTINE restrict(ARRAY(*,*,*) OF VARS f,c^) WITH f
  c.u(0,0,0)=~+0.125*[uxav0(0,0)+0.5*uxav0(-1,0)+0.5*uxav0(0,-1)+0.25*uxav0(-1,-1)]
  c.u(-1,0,0)=~+0.125*[uxav1(0,0)+0.5*uxav1(-1,0)+0.5*uxav1(0,-1)+0.25*uxav1(-1,-1)]
  c.u(0,-1,0)=~+0.125*[0.5*uxav0(-1,0)+0.25*uxav0(-1,-1)]
  c.u(-1,-1,0)=~+0.125*[0.5*uxav1(-1,0)+0.25*uxav1(-1,-1)]
  c.u(0,0,-1)=~+0.125*[0.5*uxav0(0,-1)+0.25*uxav0(-1,-1)]
  c.u(-1,0,-1)=~+0.125*[0.5*uxav1(0,-1)+0.25*uxav1(-1,-1)]
  c.u(0,-1,-1)=~+0.125*0.25*uxav0(-1,-1)
  c.u(-1,-1,-1)=~+0.125*0.25*uxav1(-1,-1)
  c.v(0,0,0)=~+0.125*[vyav0(0,0)+0.5*vyav0(-1,0)+0.5*vyav0(0,-1)+0.25*vyav0(-1,-1)]
  c.v(0,-1,0)=~+0.125*[vyav1(0,0)+0.5*vyav1(-1,0)+0.5*vyav1(0,-1)+0.25*vyav1(-1,-1)]
  c.v(-1,0,0)=~+0.125*[0.5*vyav0(-1,0)+0.25*vyav0(-1,-1)]
  c.v(-1,-1,0)=~+0.125*[0.5*vyav1(-1,0)+0.25*vyav1(-1,-1)]
  c.v(0,0,-1)=~+0.125*[0.5*vyav0(0,-1)+0.25*vyav0(-1,-1)]
  c.v(0,-1,-1)=~+0.125*[0.5*vyav1(0,-1)+0.25*vyav1(-1,-1)]
  c.v(-1,0,-1)=~+0.125*0.25*vyav0(-1,-1)
  c.v(-1,-1,-1)=~+0.125*0.25*vyav1(-1,-1)
  c.w(0,0,0)=~+0.125*[wzav0(0,0)+0.5*wzav0(-1,0)+0.5*wzav0(0,-1)+0.25*wzav0(-1,-1)]
  c.w(0,0,-1)=~+0.125*[wzav1(0,0)+0.5*wzav1(-1,0)+0.5*wzav1(0,-1)+0.25*wzav1(-1,-1)]
  c.w(-1,0,0)=~+0.125*[0.5*wzav0(-1,0)+0.25*wzav0(-1,-1)]
  c.w(-1,0,-1)=~+0.125*[0.5*wzav1(-1,0)+0.25*wzav1(-1,-1)]
  c.w(0,-1,0)=~+0.125*[0.5*wzav0(0,-1)+0.25*wzav0(-1,-1)]
  c.w(0,-1,-1)=~+0.125*[0.5*wzav1(0,-1)+0.25*wzav1(-1,-1)]
  c.w(-1,-1,0)=~+0.125*0.25*wzav0(-1,-1)
  c.w(-1,-1,-1)=~+0.125*0.25*wzav1(-1,-1)
  c.p(0,0,0)=~+0.125*[p(0,0,0)+0.5*p(-1,0,0)+0.5*p(0,-1,0)+0.25*p(-1,-1,0)+
                         0.5*p(0,0,-1)+0.25*p(-1,0,-1)+0.25*p(0,-1,-1)+0.125*p(-1,-1,-1)]
  c.p(-1,0,0)=~+0.125*[0.5*p(-1,0,0)+0.25*p(-1,-1,0)+0.25*p(-1,0,-1)+0.125*p(-1,-1,-1)]
  c.p(0,-1,0)=~+0.125*[0.5*p(0,-1,0)+0.25*p(-1,-1,0)+0.25*p(0,-1,-1)+0.125*p(-1,-1,-1)]
  c.p(-1,-1,0)=~+0.125*[0.25*p(-1,-1,0)+0.125*p(-1,-1,-1)]
  c.p(0,0,-1)=~+0.125*[0.5*p(0,0,-1)+0.25*p(-1,0,-1)+0.25*p(0,-1,-1)+0.125*p(-1,-1,-1)]
  c.p(-1,0,-1)=~+0.125*[0.25*p(-1,0,-1)+0.125*p(-1,-1,-1)]
  c.p(0,-1,-1)=~+0.125*[0.25*p(0,-1,-1)+0.125*p(-1,-1,-1)]
  c.p(-1,-1,-1)=~+0.125*0.125*p(-1,-1,-1)
END restrict