subroutine tetraedr(x,y,z,g,zero)
*                ========
*     Isosurfaces in a tetrahedron
*     ----------------------------
* ... (x,y,z) coordinates and function values
*     at tetrahedron vertices:
      real x(4),y(4),z(4),g(4),f(4)
* ... Permuted & calculated coordinates / values:
      real cv(3,4),fv(4),ca(3,4)
* ... Permutations & actions:
      integer np(4,24),must(81),move(81)
      logical first
* ... Save the following:
      data np,must,move / 96*0 , 81*0 , 81*0 /
      data first / .true. /
*
* ... Initialize actions & permutations:
      if(first) then
      call setup4(np,must,move)
      first=.false.
      endif

* ... Case encoding (base 3):
      more=0
      do 20 k=1,4
      f(k)=g(k)-zero
      if(f(k).eq.0.) more=more*3
      if(f(k).lt.0.) more=more*3+1
      if(f(k).gt.0.) more=more*3+2
20    continue
      label=must(more+1)
*
* ... Permutation of case:
      keus=move(more+1)
      do 30 i=1,4
      ik=np(i,keus)
      cv(1,i)=x(ik)
      cv(2,i)=y(ik)
      cv(3,i)=z(ik)
30    fv(i)=f(ik)
*
cHdB  write(13,*) fv
cHdB  write(13,*) label
*
      goto (1,2,3,4,5,6,7,8,9),label
*
* ... Nothing else to do:
1     continue
* ... ====
3     continue
* ... ==<<
5     continue
* ... =<<<
7     continue
* ... <<<<
      return
*
* ... Coplanar with boundary:
2     continue
* ... ===<
*     Polygon:
      do 22 j=1,3
22    write(7,*) (cv(k,j),k=1,3)
      write(7,*) '-'
      return
*
* ... Through two vertices:
4     continue
* ... ==<>
      fac=fv(3)/(fv(3)-fv(4))
      do 41 k=1,3
      ca(k,1)=cv(k,1)
      ca(k,2)=cv(k,2)
41    ca(k,3)=cv(k,3)+fac*(cv(k,4)-cv(k,3))
*     Polygon:
      do 44 j=1,3
44    write(7,*) (ca(k,j),k=1,3)
      write(7,*) '-'
      return
*
* ... Through one vertex:
6     continue
* ... =<<>
      do 63 k=1,3
63    ca(k,1)=cv(k,1)
      do 61 i=2,3
      fac=fv(i)/(fv(i)-fv(4))
      do 62 k=1,3
62    ca(k,i)=cv(k,i)+fac*(cv(k,4)-cv(k,i))
61    continue
*     Polygon:
      do 66 j=1,3
66    write(7,*) (ca(k,j),k=1,3)
      write(7,*) '-'
      return
*
* ... Triangular surface:
8     continue
* ... <<<>
      do 81 i=1,3
      fac=fv(4)/(fv(4)-fv(i))
      do 82 k=1,3
82    ca(k,i)=cv(k,4)+fac*(cv(k,i)-cv(k,4))
81    continue
*     Polygon:
      do 88 j=1,3
88    write(7,*) (ca(k,j),k=1,3)
      write(7,*) '-'
      return
*
* ... Quadrilateral surface:
9     continue
* ... <<>>
      do 91 i=1,2
      fac=fv(3)/(fv(3)-fv(i))
      do 92 k=1,3
92    ca(k,i)=cv(k,3)+fac*(cv(k,i)-cv(k,3))
91    continue
      do 93 i=1,2
      fac=fv(4)/(fv(4)-fv(i))
      do 94 k=1,3
94    ca(k,5-i)=cv(k,4)+fac*(cv(k,i)-cv(k,4))
93    continue
*     Polygon:
      do 99 j=1,4
99    write(7,*) (ca(k,j),k=1,3)
      write(7,*) '-'
      return
*
      end