c      subroutime auba2.f
c      subroutine initvrml
       subroutine initvrml
c      initialization of VRML
c      write(10,*) "#VRML V2.0 utf8"
       write(10,10)
   10  format("#VRML V2.0 utf8")
       return
       end
       subroutine point3d(nx,ny,nz,x0,y0,z0,xl,yl,zl,ico,u)
c
c      nx: number in x-direction
c      ny: number in y-direction
c      nz: number in z-direction
c      x0: x-value at the left and bottom point
c      y0: y-value at the left and bottom point
c      z0: z-value at the left and bottom point
c      xl: length in x-direction
c      yl: length in y-direction
c      zl: length in z-direction
c      ico: color code
c      u(nx*ny): input data [0.0,1.0]
c
       dimension u(1)
c
       hx=xl/float(nx-1)
       hy=yl/float(ny-1)
       hz=zl/float(nz-1)
c
       call colorb
       call colorpb(ico,nx,ny,nz,u)
       call colore
c
       call pointb
       call pointp(nx,ny,nz,x0,y0,z0,hx,hy,hz)
       call pointe
c
       return
       end
       subroutine colorb
c
       write(10,*) " Shape {"
       write(10,*) " geometry PointSet {"
       write(10,*) " color Color {"
       write(10,*) " color ["
c
       return
       end
       subroutine colore
       write(10,*) " ]"
       write(10,*) " }"
c
       return
       end
       subroutine pointb
       write(10,*) " coord Coordinate {"
       write(10,*) " point ["
       return
       end
       subroutine pointe
       write(10,*) " ]"
       write(10,*) " }"
       write(10,*) " }"
       write(10,*) " }"
c
       return
       end
       subroutine colorpb(ico,nx,ny,nz,u)
       dimension u(1)
c
       n1=nx
       n2=n1*ny
       n3=n2*nz
       last=nx*ny*nz
c
       do 10 k=1,nz
       do 10 j=1,ny
       do 10 i=1,nx
       x=hx*float(i-1)
       y=hy*float(j-1)
       z=hz*float(k-1)
       ii=i+n1*(j-1)+n2*(k-1)
c
       x1=u(ii)
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(ii.lt.last) then
       write(10,20) r,g,b
   20  format(1h ,f3.1,1x,f3.1,1x,f3.1,1h,)
       else
       write(10,22) r,g,b
   22  format(1h ,f3.1,1x,f3.1,1x,f3.1)
       end if
   10  continue
c
       return
       end
       subroutine pointp(nx,ny,nz,x0,y0,z0,hx,hy,hz)
c
       n1=nx
       n2=n1*ny
       n3=n2*nz
       last=nx*ny*nz
c
       do 10 k=1,nz
       z=z0+hz*float(k-1)
       do 10 j=1,ny
       y=y0+hy*float(j-1)
       do 10 i=1,nx
       x=x0+hx*float(i-1)
       ii=i+n1*(j-1)+n2*(k-1)
       if(ii.lt.last) then
       write(10,20) x,z,y
   20  format(1h ,f5.2,1x,f5.2,1x,f5.2,1h,)
       else
       write(10,22) x,z,y
   22  format(1h ,f5.2,1x,f5.2,1x,f5.2)
       end if
   10  continue
       return
       end
c      subroutine line3d(lasl,r,g,b,po)
       subroutine line3d(lasl,r,g,b,po)
c
c      lasl: number of points on line
c      r: red color   [0.0,1.0]
c      g: green color [0.0,1.0]
c      b: blue color  [0.0,1.0]
c      po(3*lasl): x,y,z-values in position
c
       dimension po(1)
c
       call linebe(r,g,b)
c
       do 40 ii=1,lasl
       lin=ii
       li1=3*(lin-1)+1
       li2=3*(lin-1)+2
       li3=3*lin
       x=po(li1)
       y=po(li2)
       z=po(li3)
       call linep(x,y,z,ii,lasl)
   40  continue
c
       call lineen(lasl)
       return
       end
c      subroutine symbvb
       subroutine symbvb
       write(10,10)
   10  format("Viewpoint { position 10 4 8 }")
       return
       end
       subroutine symblv(x,y,z,h,chr,n)
c
c      x: x-position of start point of character
c      y: y-position of start point of character
c      z: z-position of start point of character
c      h: height of character
c      chr: character
c      n: number of character
c
       character chr*80,ica*80,ich(80)*1
       equivalence (ica,ich)
       ica=chr
       write(10,10)
c  10  format("Viewpoint { position 10 4 8 }"/"Transform {")
   10  format("Transform {")
       write(10,12) x,y,z
   12  format("translation ",3(1x,f5.1))
       write(10,14)
   14  format("children ["/"Shape {"/"geometry Text {")  
c      write(10,16) (ich(i),i=1,n)
c  16  format("string ",'"',80a1,'"')
       write(10,*) "string ",'"',(ich(i),i=1,n),'"'
       write(10,18)
   18  format("fontStyle FontStyle {")
       write(10,20) h
   20  format("size ",f5.1)      
       write(10,22)
   22  format('family "SERIF"'/'style "ROMAN"')
       write(10,24)
   24  format("}"/"}"/"}"/"]"/"}")
       return
       end
       subroutine symblvc(x,y,z,h,r,g,b,chr,n)
c
c      x: x-position of start point of character
c      y: y-position of start point of character
c      z: z-position of start point of character
c      h: height of character
c      r: red   [0.0,1.0]
c      g: green [0.0,1.0]
c      b: blue  [0.0,1.0]
c      chr: character
c      n: number of character
c
       character chr*80,ica*80,ich(80)*1
       equivalence (ica,ich)
       ica=chr
       write(10,10)
c  10  format("Viewpoint { position 10 4 8 }"/"Transform {")
   10  format("Transform {")
       write(10,12) x,y,z
   12  format("translation ",3(1x,f5.1))
       write(10,14) r,g,b
   14  format("children ["/"Shape {"/"appearance Appearance {"
     1  /"material Material {"/"diffuseColor ",3(1x,f3.1)/
     2  "}"/"}"/"geometry Text {")
c      write(10,16) (ich(i),i=1,n)
c  16  format("string ",'"',80a1,'"')
       write(10,*) "string ",'"',(ich(i),i=1,n),'"'
       write(10,18)
   18  format("fontStyle FontStyle {")
       write(10,20) h
   20  format("size ",f5.1)
       write(10,22)
   22  format('family "SERIF"'/'style "ROMAN"')
       write(10,24)
   24  format("}"/"}"/"}"/"]"/"}")
       return
       end
c      subroutine image1
       subroutine image1(nx,ny,xb,yb,xl,yl,ico,icc,zcc,u,p)
c
c      nx: number in x-direction
c      nx: number in x-direction
c      x0: x-value at the left and bottom point
c      y0: y-value at the left and bottom point
c      xl: length in x-direction
c      yl: length in y-direction
c      ico: color code
c      icc: selection of plane, i
c         : 1 for xy plane, s2 for xz plane, 3 for yz plane
c      zcc: position of selected plane
c      u(nx*ny): input data
c      p: work array
c
       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 image2
       subroutine image2(nx,ny,xb,yb,xl,yl,ico,icc,zcc,v,u,p)
c
c      nx: number in x-direction
c      nx: number in x-direction
c      x0: x-value at the left and bottom point
c      y0: y-value at the left and bottom point
c      xl: length in x-direction
c      yl: length in y-direction
c      ico: color code
c      icc: selection of plane, i
c         : 1 for xy plane, s2 for xz plane, 3 for yz plane
c      zcc: position of selected plane
c      v: (nx-1)*(ny-1) additional input data, unnecessary
c      u(nx*ny): input data
c      p: work array
c
       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
       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
       subroutine linep(x,y,z,ii,lasl)
c
       if(ii.lt.lasl) then
       write(10,10) x,z,y
c  10  format(f5.1,1x,f5.1,1x,f5.1,1h,)
c  10  format(f5.2,1x,f5.2,1x,f5.2,1h,)
   10  format(f6.2,1x,f6.2,1x,f6.2,1h,)
       else
       write(10,20) x,z,y
c  20  format(f5.1,1x,f5.1,1x,f5.1)
c  20  format(f5.2,1x,f5.2,1x,f5.2)
   20  format(f6.2,1x,f6.2,1x,f6.2)
       end if
       return
       end
       subroutine linep2(x,y,z,ii,lasl)
c
       if(ii.lt.lasl) then
       write(10,10) x,z,y
   10  format(f5.1,1x,f5.1,1x,f5.1,1h,)
       else
       write(10,20) x,z,y
   20  format(f5.1,1x,f5.1,1x,f5.1)
       end if
       return
       end
       subroutine linebe(r,g,b)
c
       write(10,21)
   21  format("Shape {")
       write(10,22)
   22  format("geometry IndexedLineSet {")
       write(10,23)
   23  format("color Color {")
       write(10,20) r,g,b
   20  format("color [ ",f3.1,1x,f3.1,1x,f3.1," ]")
       write(10,25)
   25  format("}")
       write(10,26)
   26  format("coord Coordinate {")
       write(10,27)
   27  format("point [")
       return
       end
       subroutine lineen(lasl)
       integer ic(3000)
c
       do 10 i=1,3000
       ic(i)=i-1
   10  continue
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("colorIndex [ 0 ]")
       write(10,24)
   24  format("coordIndex [")
c
       write(10,20) (ic(i),i=1,lasl)
   20  format(10i2,90i3,900i4/(500i5))
c
       write(10,25)
   25  format("]")
       write(10,26)
   26  format("colorPerVertex FALSE")
       write(10,27)
   27  format("}")
       write(10,28)
   28  format("}")
       return
       end
       subroutine rect1(icc,zcc,r,g,b,x1,y1,x2,y2)
       lasl=5
       call linebe(r,g,b)
c
       if (icc.eq.2) go to 52
       if (icc.eq.3) go to 53
c
       z=zcc
       call linep(x1,y1,z,1,lasl)
       call linep(x2,y1,z,2,lasl)
       call linep(x2,y2,z,3,lasl)
       call linep(x1,y2,z,4,lasl)
       call linep(x1,y1,z,5,lasl)
       go to 60
c
   52  continue
       x=zcc
       call linep(x,x1,y1,1,lasl)
       call linep(x,x2,y1,2,lasl)
       call linep(x,x2,y2,3,lasl)
       call linep(x,x1,y2,4,lasl)
       call linep(x,x1,y1,5,lasl)
       go to 60
c
   53  continue
       y=zcc
       call linep(x1,y,y1,1,lasl)
       call linep(x2,y,y1,2,lasl)
       call linep(x2,y,y2,3,lasl)
       call linep(x1,y,y2,4,lasl)
       call linep(x1,y,y1,5,lasl)
   60  continue
       call lineen(lasl)
c
       return
       end       
c      zcsub1.f subroutine program
       subroutine zcoordp(nx,p)
       dimension p(1)
c
       do 100 i=1,nx
       i1=1+3*(i-1)
       i2=i1+1
       i3=i2+1
       x=p(i1)
       y=p(i2)
       z=p(i3)
       if(i.ne.nx) write(10,20) x,y,z
       if(i.eq.nx) write(10,30) x,y,z
c  20  format(f3.1,1x,f3.1,1x,f3.1,",")
c  30  format(f3.1,1x,f3.1,1x,f3.1)
   20  format(f5.3,1x,f5.3,1x,f5.3,",")
   30  format(f5.3,1x,f5.3,1x,f5.3)
  100  continue
       return
       end
       subroutine zcoorden(nx,ib)
       dimension ib(1)
       write(10,21)
   21  format("]")
       write(10,22)
   22  format("}")
       write(10,23)
   23  format("colorPerVertex TRUE")
       write(10,24)
   24  format("coordIndex [")
c
c      do 10 i=1,nx
c      ib(i)=i-1
c  10  continue
c
       write(10,102) (ib(i),i=1,nx)
  102  format(10i5)
c
       return
       end
c
c      zsub12.f  cube subroutine program
       subroutine csub2(ic,ib,aa,p,pt,ct)
       parameter(nol=12,nop=3*nol)
c      dimension pp(nop),p(8),nf(nol,6),ii(6),jj(6)
       dimension pp(nop),p(8),nf(nol,6),ii(10),jj(10)
       dimension ic(1),ib(1),aa(1),pt(1),ct(1)
c
c
       ico=ic(1)
       ua=aa(1)
       x1=aa(2)
       x2=aa(3)
       y1=aa(4)
       y2=aa(5)
       z1=aa(6)
       z2=aa(7)
c
       do 10 j1=1,6
       do 10 i1=1,nol
       nf(i1,j1)=0
   10  continue
c
       do 12 i1=1,nop
       pp(i1)=-1.0
   12  continue
c
       u1=(ua-p(1))*(ua-p(2))
       if(u1.lt.0.0) then
       nf(1,1)=1
       nf(1,3)=1
       x=(x1*(p(2)-ua)+x2*(ua-p(1)))/(p(2)-p(1))
       y=y1
       z=z1
       i1=1
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(2))*(ua-p(3))
       if(u1.lt.0.0) then
       nf(2,1)=1
       nf(2,4)=1
       x=x2
       y=(y1*(p(3)-ua)+y2*(ua-p(2)))/(p(3)-p(2))
       z=z1
       i1=1+3*1
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
ccc
c      do 242 j1=1,12
c      i1=1+3*(j1-1)
c      i2=3*j1
c      write(6,22) j1,(pp(ia),ia=i1,i2)
c 242  continue
c
       u1=(ua-p(3))*(ua-p(4))
       if(u1.lt.0.0) then
       nf(3,1)=1
       nf(3,5)=1
       x=(x2*(p(4)-ua)+x1*(ua-p(3)))/(p(4)-p(3))
       y=y2
       z=z1
       i1=1+3*2
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(4))*(ua-p(1))
       if(u1.lt.0.0) then
       nf(4,1)=1
       nf(4,6)=1
       x=x1
       y=(y2*(p(1)-ua)+y1*(ua-p(4)))/(p(1)-p(4))
       z=z1
       i1=1+3*3
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c
       u1=(ua-p(5))*(ua-p(6))
       if(u1.lt.0.0) then
       nf(5,2)=1
       nf(5,3)=1
       x=(x1*(p(6)-ua)+x2*(ua-p(5)))/(p(6)-p(5))
       y=y1
       z=z2
       i1=1+3*4
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(6))*(ua-p(7))
       if(u1.lt.0.0) then
       nf(6,2)=1
       nf(6,4)=1
       x=x2
       y=(y1*(p(7)-ua)+y2*(ua-p(6)))/(p(7)-p(6))
       z=z2
       i1=1+3*5
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(7))*(ua-p(8))
       if(u1.lt.0.0) then
       nf(7,2)=1
       nf(7,5)=1
       x=(x2*(p(8)-ua)+x1*(ua-p(7)))/(p(8)-p(7))
       y=y2
       z=z2
       i1=1+3*6
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(8))*(ua-p(5))
       if(u1.lt.0.0) then
       nf(8,2)=1
       nf(8,6)=1
       x=x1
       y=(y2*(p(5)-ua)+y1*(ua-p(8)))/(p(5)-p(8))
       z=z2
       i1=1+3*7
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c
       u1=(ua-p(1))*(ua-p(5))
       if(u1.lt.0.0) then
       nf(9,3)=1
       nf(9,6)=1
       x=x1
       y=y1
       z=(z1*(p(5)-ua)+z2*(ua-p(1)))/(p(5)-p(1))
       i1=1+3*8
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(2))*(ua-p(6))
       if(u1.lt.0.0) then
       nf(10,3)=1
       nf(10,4)=1
       x=x2
       y=y1
       z=(z1*(p(6)-ua)+z2*(ua-p(2)))/(p(6)-p(2))
       i1=1+3*9
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(3))*(ua-p(7))
       if(u1.lt.0.0) then
       nf(11,4)=1
       nf(11,5)=1
       x=x2
       y=y2
       z=(z1*(p(7)-ua)+z2*(ua-p(3)))/(p(7)-p(3))
       i1=1+3*10
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
       u1=(ua-p(4))*(ua-p(8))
       if(u1.lt.0.0) then
       nf(12,5)=1
       nf(12,6)=1
       x=x1
       y=y2
       z=(z1*(p(8)-ua)+z2*(ua-p(4)))/(p(8)-p(4))
       i1=1+3*11
       i2=i1+1
       i3=i2+1
       pp(i1)=x
       pp(i2)=y
       pp(i3)=z
       end if
c
c
c      write(6,20) (p(i1),i1=1,8)
   20  format(1h ,4f10.4)
       do 24 j1=1,12
       i1=1+3*(j1-1)
       i2=3*j1
c      write(6,22) j1,(pp(ia),ia=i1,i2)
   22  format(1h ,i5,3f10.4)
   24  continue
       do 28 j1=1,12
c      write(6,26) j1,(nf(j1,ia),ia=1,6)
   26  format(1h ,i5,6i5)
   28  continue
c
c
       call csub1(ic,ib,p,pp,nf,pt,ct)
c
c 100  continue
       return
       end
c
c       
c
c      cube subroutine program
       subroutine csub1(ic,ib,p,pp,nf,pt,ct)
       parameter(nol=12,nop=3*nol)
c      dimension p(8),pp(nop),nf(nol,6),ii(6),jj(6),ic(1)
       dimension p(8),pp(nop),nf(nol,6),ii(10),jj(10),ic(1)
       dimension ib(1),pa(nop),pt(1),ct(1)
c
       ico=ic(1)
       iin=0
       do 10 i=1,4
       if(nf(i,1).eq.1) then
       iin=iin+1
       ii(iin)=i
       jj(iin)=1
       end if
   10  continue
c
       do 12 i=5,8
       if(nf(i,2).eq.1) then
       iin=iin+1
       ii(iin)=i
       jj(iin)=2
       end if
   12  continue
cc
cc
       if(iin.eq.0) go to 100
       if(iin.eq.2.and.ii(2).le.4) go to 200
       if(iin.eq.2.and.ii(2).ge.5) go to 300
       if(iin.eq.4) go to 400
       go to 990
c
  100  continue
       do 110 i=9,12
       j=i-6
       if(nf(i,j).eq.1) then
       iin=iin+1
       ii(iin)=i
       jj(iin)=j
       end if
  110  continue
       go to 1000
c
c
  200  continue
       i1=ii(1)
       i2=ii(2)
       j1=jj(1)
       j2=jj(2)
c
       do 210 j=3,6
       if(nf(i2,j).eq.1) then
       j3=j
       end if
  210  continue
       do 212 i=9,12
       if(i.eq.i2) go to 212
       if(nf(i,j3).eq.1) then
       iin=iin+1
       i3=i
       ii(iin)=i3
       jj(iin)=j3
       end if
  212  continue
c
       i4=0
       do 220 j=3,6
       if(j.eq.j3) go to 220
       if(nf(i3,j).eq.1) then
       j4=j
       end if
  220  continue
       do 222 i=9,12
       if(i.eq.i2) go to 222
       if(i.eq.i3) go to 222
       if(nf(i,j4).eq.1) then
       iin=iin+1
       i4=i
       ii(iin)=i4
       jj(iin)=j4
       end if
  222  continue
       if(i4.eq.0) go to 1000
c
       do 230 j=3,6
       if(j.eq.j3) go to 230
       if(j.eq.j4) go to 230
       if(nf(i4,j).eq.1) then
       j5=j
       end if
  230  continue
       do 232 i=9,12
       if(i.eq.i2) go to 232
       if(i.eq.i3) go to 232
       if(i.eq.i4) go to 232
       if(nf(i,j5).eq.1) then
       iin=iin+1
       i5=i
       ii(iin)=i5
       jj(iin)=j5
       end if
  232  continue
       go to 1000
c
c
  300  continue
       i1=ii(1)
       i2=ii(2)
       j1=jj(1)
       j2=jj(2)
c
       do 310 j=3,6
       if(nf(i2,j).eq.1) then
       j3=j
       end if
  310  continue
       do 312 i=9,12
       if(i.eq.i2) go to 312
       if(nf(i,j3).eq.1) then
       iin=iin+1
       i3=i
       ii(iin)=i3
       jj(iin)=j3
       end if
  312  continue
c
       i4=0
       do 320 j=3,6
       if(j.eq.j3) go to 320
       if(nf(i3,j).eq.1) then
       j4=j
       end if
  320  continue
       do 322 i=9,12
       if(i.eq.i2) go to 322
       if(i.eq.i3) go to 322
       if(nf(i,j4).eq.1) then
       iin=iin+1
       i4=i
       ii(iin)=i4
       jj(iin)=j4
       end if
  322  continue
       if(i4.eq.0) go to 1000
c
       do 330 j=3,6
       if(j.eq.j3) go to 330
       if(j.eq.j4) go to 330
       if(nf(i4,j).eq.1) then
       j5=j
       end if
  330  continue
       do 332 i=9,12
       if(i.eq.i2) go to 332
       if(i.eq.i3) go to 332
       if(i.eq.i4) go to 332
       if(nf(i,j5).eq.1) then
       iin=iin+1
       i5=i
       ii(iin)=i5
       jj(iin)=j5
       end if
  332  continue
       go to 1000
c
c
  400  continue
       iin=2
       i1=ii(1)
       i2=ii(2)
       j1=jj(1)
       j2=jj(2)
c
       do 410 j=3,6
       if(nf(i2,j).eq.1) then
       j3=j
       end if
  410  continue
       do 412 i=5,12
       if(i.eq.i2) go to 412
       if(nf(i,j3).eq.1) then
       iin=iin+1
       i3=i
       ii(iin)=i3
       jj(iin)=j3
       end if
  412  continue
c
       i4=0
       do 420 j=2,6
       if(j.eq.j3) go to 420
       if(nf(i3,j).eq.1) then
       j4=j
       end if
  420  continue
       do 422 i=5,12
       if(i.eq.i2) go to 422
       if(i.eq.i3) go to 422
       if(nf(i,j4).eq.1) then
       iin=iin+1
       i4=i
       ii(iin)=i4
       jj(iin)=j4
       end if
  422  continue
       if(i4.eq.0) go to 1000
c
       i5=0
       do 430 j=2,6
       if(j.eq.j3) go to 430
       if(j.eq.j4) go to 430
       if(nf(i4,j).eq.1) then
       j5=j
       end if
  430  continue
       do 432 i=5,12
       if(i.eq.i2) go to 432
       if(i.eq.i3) go to 432
       if(i.eq.i4) go to 432
       if(nf(i,j5).eq.1) then
       iin=iin+1
       i5=i
       ii(iin)=i5
       jj(iin)=j5
       end if
  432  continue
       if(i5.eq.0) go to 1000
c
       do 440 j=2,6
       if(j.eq.j3) go to 440
       if(j.eq.j4) go to 440
       if(j.eq.j5) go to 440
       if(nf(i5,j).eq.1) then
       j6=j
       end if
  440  continue
       do 442 i=5,12
       if(i.eq.i2) go to 442
       if(i.eq.i3) go to 442
       if(i.eq.i4) go to 442
       if(i.eq.i5) go to 442
       if(nf(i,j6).eq.1) then
       iin=iin+1
       i6=i
       ii(iin)=i6
       jj(iin)=j6
       end if
  442  continue
c
c
 1000  continue
       do 1010 j=1,iin
       ii1=ii(j)
       jj1=jj(j)
       i1=1+3*(ii1-1)
       i2=i1+1
       i3=i2+1
       j1=1+3*(j-1)
       j2=j1+1
       j3=j2+1
c
       ic(5)=ic(5)+1
       ic(6)=ic(6)+1
       k=ic(5)
       ka=ic(6)
       k1=1+3*(k-1)
       k2=k1+1
       k3=k2+1
       ct(k)=p(j)
       ib(ka)=k-1
c
       x=pp(i1)
       y=pp(i2)
       z=pp(i3)
       pt(k1)=x
       pt(k2)=y
       pt(k3)=z
c      write(6,1020) j,ii1,jj1,x,y,z
 1020  format(1h ,3i4,3f10.4)
 1010  continue
c
       if(iin.eq.0) go to 990
       ic(6)=ic(6)+1
       ka=ic(6)
       ib(ka)=-1
c
c
  990  continue
       return
       end
c      cube subroutime program
       subroutine csuba(ic,aa,u)
c      isosurface plot
c
c      ic(10):
c      ic(1)=ico,ic(2)=nx,ic(3)=ny,ic(4)=nz
c      aa(20):
c      aa(1)=ua,aa(2)=x1,aa(3)=x2,aa(4)=y1,aa(5)=y2
c      aa(6)=z1,aa(7)=z2,aa(8)=xl,aa(9)=yl,aa(10)=zl
c      aa(11)=x0,aa(12)=y0,aa(13)=z0
c      u(nx*ny*nz): input data
c
       parameter(ipt=9999,ipp=ipt*3)
       parameter(nol=12,nop=3*nol)
c      dimension pp(nop),p(8),nf(nol,6)
       dimension pp(nop),p(8),nf(nol,6)
       dimension ic(1),aa(1),u(1)
       dimension ib(ipt),ct(ipt),pt(ipp)
c
c
c      dimension ic(1),aa(1),u(1),pp(nop),nf(nol,6)
c      dimension p(8)
c
       ico=ic(1)
       nx=ic(2)
       ny=ic(3)
       nz=ic(4)
       ua=aa(1)
       x1=aa(2)
       x2=aa(3)
       y1=aa(4)
       y2=aa(5)
       z1=aa(6)
       z2=aa(7)
       xl=aa(8)
       yl=aa(9)
       zl=aa(10)
       x0=aa(11)
       y0=aa(12)
       z0=aa(13)
c
       n1=nx
       n2=n1*ny
       n3=n2*nz
       nxs=nx-1
       nys=ny-1
       nzs=nz-1
       dx=xl/float(nxs)
       dy=yl/float(nys)
       dz=zl/float(nzs)
c
c
       do 200 k=1,nzs
       z1=z0+dz*float(k-1)
       z2=z1+dz
       ic(5)=0
       ic(6)=0
c
       do 100 j=1,nys
       y1=y0+dy*float(j-1)
       y2=y1+dy
       do 100 i=1,nxs
       x1=x0+dx*float(i-1)
       x2=x1+dx
c
       i1=i+n1*(j-1)+n2*(k-1)
       i2=i1+1
       i3=i2+nx
       i4=i1+nx
       i5=i1+n2
       i6=i5+1
       i7=i6+nx
       i8=i5+nx
       p(1)=u(i1)
       p(2)=u(i2)
       p(3)=u(i3)
       p(4)=u(i4)
       p(5)=u(i5)
       p(6)=u(i6)
       p(7)=u(i7)
       p(8)=u(i8)
c
c
       aa(1)=ua
       aa(2)=x1
       aa(3)=x2
       aa(4)=y1
       aa(5)=y2
       aa(6)=z1
       aa(7)=z2
c
c
c      write(6,12) (ic(i1),i1=1,10)
c      write(6,14) (aa(i1),i1=1,20)
c      write(6,16) (p(i1),i1=1,8)
   12  format(1h ,10i5)
   14  format(1h ,5f10.5)
   16  format(1h ,4f10.5)
c
       call csub2(ic,ib,aa,p,pt,ct)
c
  100  continue
c
       iin1=ic(5)
       iin2=ic(6)
c      write(6,12) k,iin1,iin2
c
       if(iin1.eq.0) go to 200
c
       call facebe
c
       call colorp(iin1,ct,ico)
c
       call coordbe
c
       call zcoordp(iin1,pt)
c
       iin21=iin2-1
       call zcoorden(iin21,ib)
       call faceen

c
  200  continue
       return
       end