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