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