c      cube main program
c      parameter(nx=2,ny=2,nz=2)
c      parameter(nx=4,ny=4,nz=4)
c      parameter(nx=10,ny=10,nz=10)
c      parameter(nx=20,ny=20,nz=20)
       parameter(nx=30,ny=30,nz=30)
c      parameter(nx=100,ny=100,nz=100)
c      parameter(nx=200,ny=200,nz=200)
       parameter(n1=nx,n2=n1*ny,n3=n2*nz)
       parameter(x0=0.0,y0=0.0,z0=0.0)
       parameter(xl=1.0,yl=1.0,zl=1.0)
c      parameter(xl=9.9,yl=9.9,zl=9.9)
       parameter(nol=12,nop=3*nol,ico=7)
       dimension u(n3),p(8),pp(nop),nf(nol,6)
       dimension ic(10),aa(20)
c
c
       call initvrml
c
c
c
       ua=0.5
c      ua=0.8
c
       ic(1)=ico
       ic(2)=nx
       ic(3)=ny
       ic(4)=nz
       aa(1)=ua
       aa(2)=x1
       aa(3)=x2
       aa(4)=y1
       aa(5)=y2
       aa(6)=z1
       aa(7)=z2
       aa(8)=xl
       aa(9)=yl
       aa(10)=zl
       aa(11)=x0
       aa(12)=y0
       aa(13)=z0
c
c
       u(1)=0.01
       u(2)=0.61
       u(3)=0.451
       u(4)=0.911
       u(5)=0.221
       u(6)=0.101
       u(7)=0.901
       u(8)=0.31
c
c      amp1=2.0**6 
       amp1=4.0**6
       dx=xl/float(nx-1)
       dy=yl/float(ny-1)
       dz=zl/float(nz-1)
       do 10 k=1,nz
       do 10 j=1,ny
       do 10 i=1,nx
       z=dz*float(k-1)
       y=dy*float(j-1)
       x=dx*float(i-1)
       i1=i+nx*(j-1)+n2*(k-1)
c      u(i1)=amp1*x*y*z*(1.0-x)*(1.0-y)*(1.0-z)
       ux=abs(x*(0.5-x)*(x-0.5)*(1.0-x))
       uy=abs(y*(0.5-y)*(y-0.5)*(1.0-y))
       uz=abs(z*(0.5-z)*(z-0.5)*(1.0-z))
       u(i1)=amp1*ux*uy*uz

       ux=x*(0.5-x)
       if(x.gt.0.5) then
       ux=(x-0.5)*(1.0-x)
       end if
       uy=y*(0.5-y)
       if(y.gt.0.5) then
       uy=(y-0.5)*(1.0-y)
       end if
       uz=z*(0.5-z)
       if(z.gt.0.5) then
       uz=(z-0.5)*(1.0-z)
       end if
       u(i1)=amp1*ux*uy*uz
   10  continue
c
c
c
       call csuba(ic,aa,u)
c
c
       stop
       end