c      main program
       parameter(nx=101,ny=101,nx2=nx*2,nxy=nx*ny,ico=7)
       parameter(nx0=nx-1,ny0=ny-1,nxy0=nx0*ny0)
       parameter(nxa=nx0-1,nya=ny0-1)
       parameter(x0=0.0,y0=0.0,z0=0.0)
       parameter(xll=4.0,yll=4.0,zll=4.0)
       dimension v(nxy0),u(nxy),p(nx2)
c      icc=1:x-y plane, 2:y-z plane, 3:x-z plane
c
       pi=3.1415926
       amp=1.0
c
       do 10 j=1,ny
       do 10 i=1,nx
       x=2.0*pi*float(i-1)/float(nx-1)
       y=2.0*pi*float(j-1)/float(ny-1)
       i1=i+nx*(j-1)
       u(i1)=amp*sin(x)*sin(y)
       u(i1)=0.5*(u(i1)+amp)
c      u(i1)=0.5*(u(i1)+1.0)
   10  continue
c
       do 12 j=1,ny0
       do 12 i=1,nx0
       x=2.0*pi*(float(i-1)+0.5)/float(nx-1)
       y=2.0*pi*(float(j-1)+0.5)/float(ny-1)
       i1=i+nx0*(j-1)
       v(i1)=amp*sin(x)*sin(y)
       v(i1)=0.5*(v(i1)+amp)
   12  continue
c
       call trans1(nx,ny,v,u)
c 
       xb=x0
       yb=y0
       xl=xll
       yl=yll
       icc=3
       zcc=0.0
c      
c
       call initvrml
c
       xb=x0
       yb=y0
       zb=z0
       xl=xll
       yl=yll
       zl=zll
       icc=1
       zcc=0.0
       zcc=yb
       call image1(nx,ny,xb,yb,xl,yl,ico,icc,zcc,v,u,p)
c
       icc=2
       zcc=0.0
       zcc=xl
       call image1(nx,ny,xb,yb,xl,yl,ico,icc,zcc,v,u,p)
c
       icc=3
       zcc=0.0
       zcc=zb
       call image1(nx,ny,xb,yb,xl,yl,ico,icc,zcc,v,u,p)
c
       stop
       end
c      subroutine image1
       subroutine image1(nx,ny,xb,yb,xl,yl,ico,icc,zcc,v,u,p)
       dimension p(1),u(1),v(1)
       ny0=ny-1
       nx2=nx*2
c
       do 100 ii=1,ny0
c
       call trans2(nx,ny,ii,u,p)
       call facebe
c
       call colorp(nx2,p,ico)
c
       call coordbe
c
c      call coordp(nx2,p)
       call coordp(ii,nx,ny,xb,yb,xl,yl,icc,zcc,p)
c
       call coorden(nx)
       call faceen
c
  100  continue
c
       return
       end
c      subroutine trans1
       subroutine trans1(nx,ny,v,u)
       dimension v(1),u(1)
c
       nx0=nx-1
       ny0=ny-1
       nxa=nx0-1
       nya=ny0-1
c
       do 10 j=1,nya
       do 10 i=1,nxa
       i1=i+1+nx*j
       j1=i+nx0*(j-1)
       u(i1)=0.25*(v(j1)+v(j1+1)+v(j1+nx0)+v(j1+nx))
   10  continue
c
       do 20 i=1,nxa
       i1=i+1
       i2=i1+nx*ny0
       j1=i+1
       j2=i+1+nx0*nya
       u(i1)=0.75*(v(j1-1)+v(j1))-0.25*(v(j1-1+nx0)+v(j1+nx0))
       u(i2)=0.75*(v(j2-1)+v(j2))-0.25*(v(j2-1-nx0)+v(j2-nx0))
   20  continue
c
       do 30 j=1,nya
       i1=1+nx*j
       i2=i1+nx0
       j1=1+nx0*j
       j2=j1+nxa
       u(i1)=0.75*(v(j1-nx0)+v(j1))-0.25*(v(j1-nx0+1)+v(j1+1))
       u(i2)=0.75*(v(j2-nx0)+v(j2))-0.25*(v(j2-nx0-1)+v(j2-1))
   30  continue
c
       i1=1
       i2=i1+nx0
       i3=1+nx*ny0
       i4=i3+nx0
       j1=1
       j2=j1+nxa
       j3=1+nx0*nya
       j4=j3+nxa
       u(i1)=0.5*(3.0*v(j1)-v(j1+nx))
       u(i2)=0.5*(3.0*v(j2)-v(j2+nxa))
       u(i3)=0.5*(3.0*v(j3)-v(j3-nxa))
       u(i4)=0.5*(3.0*v(j4)-v(j4-nx))
c
c      write(6,*) i1,i2,i3,i4
c      write(6,*) j1,j2,j3,j4
c      write(6,*) nx,nx0,nxa,ny,ny0,nya
c
       return
       end
c      subroutine trans2
       subroutine trans2(nx,ny,ii,u,p)
       dimension p(1),u(1)
       do 10 i=1,nx
       j1=i+nx*(ii-1)
       j2=j1+nx
       i1=1+2*(i-1)
       i2=i1+1
       p(i1)=u(j1)
       p(i2)=u(j2)
   10  continue
       return
       end
c      subroutine initvrml
       subroutine initvrml
       write(10,10)
   10  format("#VRML V2.0 utf8")
       return
       end
       subroutine colorp(nx,p,ico)
       dimension p(1)
       r=1.0
       g=1.0
       b=1.0
       do 100 i=1,nx
       x1=p(i)
c
       a0=1.20
       a0=1.00
       if(ico.eq.2) go to 32
       if(ico.eq.3) go to 33
       if(ico.eq.4) go to 34
       if(ico.eq.5) go to 35
       if(ico.eq.6) go to 36
       if(ico.eq.7) go to 37
       r=a0-(a0-0.5)*x1
       g=a0-(a0-0.5)*x1
       b=a0-(a0-0.5)*x1
       go to 40
   32  continue
       r=a0-a0*x1
       g=1.0
       b=a0-a0*x1
       go to 40
   33  continue
       r=1.0
       g=a0-a0*x1
       b=a0-a0*x1
       go to 40
   34  continue
       a0=2.0
       y1=x1*(a0+4.0)
       r=a0-y1
       if(y1.ge.a0) r=y1-(a0+1.0)
       g=a0+3.0-y1
       b=a0+1.0-y1
       if(y1.ge.a0+2.0) b=y1-(a0+3.0)
       go to 40
   35  continue
       a0=2.0
       y1=x1*(a0+4.0)
       b=a0-y1
       if(y1.ge.a0) b=y1-(a0+1.0)
       r=a0+3.0-y1
       g=a0+1.0-y1
       if(y1.ge.a0+2.0) g=y1-(a0+3.0)
       go to 40
   36  continue
       a0=2.0
       y1=x1*(a0+4.0)
       r=a0-y1
       if(y1.ge.a0) r=y1-(a0+1.0)
       b=a0+3.0-y1
       g=a0+1.0-y1
       if(y1.ge.a0+2.0) g=y1-(a0+3.0)
       go to 40
   37  continue
       x1=10.0*x1
       b=x1
       if(x1.ge.3.0) b=4.0-x1
       if(x1.ge.6.0) b=0.25*(x1-6.0)
       g=x1-2.0
       if(x1.ge.5.0) g=6.0-x1
       if(x1.ge.6.0) g=0.25*(x1-6.0)
       r=x1
       if(x1.ge.1.0) r=2.0-x1
       if(x1.ge.4.0) r=x1-4.0
       go to 40
c
   40  continue
       r=amax1(r,0.0)
       g=amax1(g,0.0)
       b=amax1(b,0.0)
       r=amin1(r,1.0)
       g=amin1(g,1.0)
       b=amin1(b,1.0)
c
       if(i.ne.nx) write(10,20) r,g,b
       if(i.eq.nx) write(10,30) r,g,b
   20  format(f4.2,1x,f4.2,1x,f4.2,",")
   30  format(f4.2,1x,f4.2,1x,f4.2)
  100  continue
       return
       end
       subroutine coordp(ii,nx,ny,xb,yb,xl,yl,icc,zcc,p)
       dimension p(1)
       x=1.0
       y=1.0
       z=1.0
c      xb=0.0
c      yb=0.0
c      dx=1.0
c      dy=1.0
       nxy=nx*ny
       dx=xl/float(nx-1)
       dy=yl/float(ny-1)
c      do 100 j1=1,2
       do 100 i=1,nx
       do 100 j1=1,2
       j=j1-1+ii
       ixy=i+nx*(j-1)
c
       if (icc.eq.2) go to 52
       if (icc.eq.3) go to 53
c
       x=xb+dx*float(i-1)
       y=yb+dy*float(j-1)
       z=zcc
       go to 60
   52  continue
       y=xb+dx*float(i-1)
       z=yb+dy*float(j-1)
       x=zcc
       go to 60
   53  continue
       x=xb+dx*float(i-1)
       z=yb+dy*float(j-1)
       y=zcc
   60  continue
c
       if(i.ne.nx) write(10,20) x,y,z
       if(i.eq.nx) write(10,30) x,y,z
   20  format(f4.2,1x,f4.2,1x,f4.2,",")
   30  format(f4.2,1x,f4.2,1x,f4.2)
  100  continue
       return
       end
       subroutine facebe
       write(10,21)
   21  format("Shape {")
       write(10,22)
   22  format("geometry IndexedFaceSet {")
       write(10,23)
   23  format("color Color {")
       write(10,24)
   24  format("color [")
       return
       end
       subroutine coordbe
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("coord Coordinate {")
       write(10,24)
   24  format("point [")
       return
       end
       subroutine coorden(nx)
       dimension ia(2000)
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("colorPerVertex TRUE")
       write(10,24)
   24  format("coordIndex [")
c
       do 10 i=1,nx
       i1=1+5*(i-1)
       i2=2+5*(i-1)
       i3=3+5*(i-1)
       i4=4+5*(i-1)
       i5=5+5*(i-1)
       ia(i1)=0+2*(i-1)
       ia(i2)=1+2*(i-1)
       ia(i3)=3+2*(i-1)
       ia(i4)=2+2*(i-1)
       ia(i5)=-1
   10  continue
c
       na2=5*(nx-1)-1
       write(10,102) (ia(i),i=1,na2)
  102  format(4(4i2,i3),2i2,3i3,220i3,2i3,2i4,i3,
     1        449(4i4,i3),2i4,2i5,i3,3000(4i5,i3))
c
       return
       end
       subroutine faceen
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("solid FALSE")
       write(10,23)
       write(10,23)
   23  format("}")
       return
       end