subroutine element(x,y,z,e,nds)
*                =======
*     Finite Element of ACR equation
*     ------------------------------
      real x(8),y(8),z(8),e(8,8)
      real ddg(3,4),T(3,3),vel(3),flux(4)
* ... subroutine "coeffs" (User!):
      double precision xx,yy,zz
      double precision axx,ayy,azz,axy,ayz,azx,ax,ay,az,a0
      integer look(4,8+6)

* ... Local vertex numbers of tetrahedra inside bricks:
      data look / 1,2,3,5 , 2,4,1,6 , 3,1,4,7 , 4,3,2,8 ,
     /            5,7,6,1 , 6,5,8,2 , 7,8,5,3 , 8,6,7,4 ,
     /            1,2,3,4 , 2,3,1,5 , 3,1,2,6 ,
     /            4,6,5,1 , 5,4,6,2 , 6,5,4,3 /
* ... Hexahedron:
      more=0
* ... Prism:
      if(nds.eq.6) more=8
*
* ... Clear element matrix:
      do 5 i=1,nds
      do 5 j=1,nds
5     e(i,j)=0.
*
      do 10 n=1,nds
*     =============
* ... Quantities at
*     brick corners:
      xx=x(n)
      yy=y(n)
      zz=z(n)
* ... According to User:
      call coeffs(xx,yy,zz,
     , axx,ayy,azz,axy,ayz,azx,ax,ay,az,a0)
* ... Diffusion tensor:
      T(1,1)=-axx
      T(2,2)=-ayy
      T(3,3)=-azz
      T(1,2)=-axy
      T(2,1)=-axy
      T(2,3)=-ayz
      T(3,2)=-ayz
      T(3,1)=-azx
      T(1,3)=-azx
* ... First & zero order:
      vel(1)=-ax
      vel(2)=-ay
      vel(3)=-az
      zero=-a0
*
*     Differentiation Matrix
*     ======================
* ... Vertex as origin:
      x0=x(n)
      y0=y(n)
      z0=z(n)
* ... Tetrahedron:
      do 15 k=1,3
      node=look(k+1,n+more)
      ddg(k,2)=x(node)-x0
      ddg(k,3)=y(node)-y0
      ddg(k,4)=z(node)-z0
15    continue
*
* ... Inverse matrix:
      call dinv3(ddg(1,2),det)
      weight=det/nds
* ... Must be positive!
*
* ... First column:
      do 35 i=1,3
      help=0.
      do 40 j=1,3
40    help=help-ddg(i,j+1)
35    ddg(i,1)=help
*
*     Diffusion part
*     ==============
      do 20 i=1,4
      ie=look(i,n+more)
      do 20 j=1,4
      je=look(j,n+more)
      help=0.
      do 25 k=1,3
      do 25 m=1,3
25    help=help+ddg(k,i)*T(k,m)*ddg(m,j)
20    e(ie,je)=e(ie,je)-help*weight
*
*     Convection part
*     ===============
* ... Fluxes:
      do 45 i=2,4
      help=0.
      do 50 j=1,3
50    help=help+vel(j)*ddg(j,i)
45    flux(i)=2.*help*weight
*
* ... Upwind principle:
      help=0.
      do 55 i=2,4
      flux(i)=min(flux(i),0.)
55    help=help+flux(i)
      flux(1)=-help

* ... Account for it:
      do 80 j=1,4
      je=look(j,n+more)
80    e(n,je)=e(n,je)+flux(j)
*
*     Zero'th part
*     ============
      e(n,n)=e(n,n)+zero*weight
*
10    continue
*
      return
      end