program MAIN
*             ====
* This code was developed & is CopyLefted by:
*
* Han de Bruijn; Production&Services "A little bit of Physics would be  (===)
* TUD Computing Centre; P.O. Box 354  NO Idleness in Mathematics"(HdB) @-O^O-@
* 2600 AJ  Delft; The Netherlands -- http://pubwww.tudelft.nl/~rcpshdb  #/_\#
* E-mail: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 278 3787. Tel: 2751.   ###
*
      parameter (NODES=1500,LINES=4500)
      real xx(NODES),yy(NODES),zz(NODES)
      real xd(3),yd(3),zd(3)
      integer nr(3)
      logical done,triangle
*
* ... How many contour levels:
5     write(*,'(a,$)') 'Number of contour levels: '
      read(*,'(i3)',err=5) MOST
* ... Coordinates:
      open(1,file='coor2d.dat',status='old')
      do 10 k=1,NODES
      read(1,*,end=11) xx(k),yy(k),zz(k)
      NN=k
   10 continue
      stop 'NODES too small'
   11 close(1)
* ... Connectivity:
      open(1,file='topol2d.dat')
*
* ... PostScript plotfile:
      open(7,file='contour.ps')
      call plots
      do 20 k=1,LINES
      read(1,*,end=21) (nr(j),j=1,3)
      do 25 i=1,3
      ii=nr(i)+1
      xd(i)=xx(ii)
      yd(i)=yy(ii)
25    continue
* ... Mesh plotting:
      call plot(xd(3),yd(3),'u')
      do 35 i=1,3
      call plot(xd(i),yd(i),'d')
35    continue
20    continue
      stop 'LINES too small'
* ... Second page:
21    write(7,'(a)') 
     ) 'gsave stroke showpage grestore'
*
* ... Connectivity again:
      rewind(1)
      do 40 k=1,LINES
      read(1,*,end=41) (nr(i),i=1,3)
      do 45 i=1,3
      ii=nr(i)+1
      xd(i)=xx(ii)
      yd(i)=yy(ii)
      zd(i)=zz(ii)
45    continue
*
* ... Element bounds:
      dmin=1.
      dmax=0.
      do 60 m=1,3
      dmin=amin1(dmin,zd(m))
      dmax=amax1(dmax,zd(m))
60    continue
      call bounds(dmin,MOST,n1,nop)
      call bounds(dmax,MOST,nop,n2)
* ... Loop through levels:
      do 50 n=n1,n2
      zero=float(n)/float(MOST)
* ... Very basic contouring:
      done=triangle(xd,yd,zd,zero,xb,yb,xe,ye)
      if(.not.done) goto 50
      call plot(xb,yb,'u')
      call plot(xe,ye,'d')
50    continue
*
40    continue
      stop 'LINES too small'
* ... Ending:
41    write(7,'(a)') 'stroke showpage'
      close(1)
      close(7)
      stop 'OK'
      end