      subroutine elmt06(d,ul,xl,ix,tl,s,p,ndf,ndm,nst,isw)

c     Purpose: Two dimensional laplace equation with a reaction term

c     Inputs:
c        d(*)      - Material set parameters
c        ul(ndf,*) - Solution variables for element
c        xl(ndf,*) - Nodal coordinates for element
c        ix(*)     - Nodal connections for element
c        tl(*)     - Nodal temperatures for element
c        ndf       - Number dof/node
c        ndm       - Spatial dimension of FEM mesh
c        nst       - Size of element matrix
c        isw       - Switch parameter to control action taken by routine

c     Outputs:
c        s(nst,*)  - Element matrix
c        p(*)      - Element vector

      implicit real*8 (a-h,o-z)

      character*4 head,wlab(2)*12,yyy*80
c     integer*2   ix(1)
      integer     ix(1)
      real*8      xl(ndm,*),tl(*),dm,aa
      real*8      d(*),ul(ndf,*),s(nst,nst),p(nst),shp(3,9),sg(9),tg(9),
     1            wg(9),coord

      integer         maxa
      parameter      (maxa=32000)
      common /adata/  aa(maxa)
      common /bdata/  head(20)
      common /cdata/  numnp,numel,nummat,nen,neq
      common /eldata/ dm,n,ma,mct,iel,nel
      common /iofile/ ior,iow
      common /ydata / yyy
      data wlab/'  p l a n e ','axisymmetric'/
c.... compute the quadrature points
      if(isw.gt.2) then
        l = max(1,min(nel/2,3))
        call pgauss(l,lint,sg,tg,wg)
      endif
c.... input material properties
      if(isw.eq.1) then
  1     if(ior.lt.0) write(*,3000)
        call pintio(yyy,10)
        read(yyy,1000,err=110) (d(i),i=1,6),nn,kat
        if(kat.ne.2) kat=1
        nn = max(-1,min(1,nn))
        write(iow,2000) (d(i),i=1,6),nn,wlab(kat)
        if(ior.lt.0) write(*,2000) (d(i),i=1,6),nn,wlab(kat)
        d(7) = nn
        d(8) = kat
        d(9) = d(2)*d(3)
        return
110     call perror('PCELM6',yyy)
        go to 1
c.... compute conductivity matrix and residual
      elseif(mod(isw,3).eq.0) then
        nn  = d(7)
        kat = d(8)
        do 330 l=1,lint
          call shape(sg(l),tg(l),xl,shp,xsj,ndm,nel,ix,.false.)
          xsj = xsj*wg(l)
          if(kat.eq.2) xsj = xsj*coord(xl,shp,ndm,nel)
          call flo06(1.0d0,shp,ul,q1,q2,qm,uu,ndf,nel)
          if(nn.eq.0) then
            qq = dm*d(4)
            dq = 0.0
          elseif(nn.eq.1) then
            qq = dm*d(4)*exp(d(5)*uu)
            dq = d(5)
          elseif(uu.ne.0.0d0) then
            qq = dm*d(4)*exp(d(5)-d(5)*d(6)/uu)
            dq = d(5)*d(6)/uu**2
          else
            write(*,*) ' ** ELMT06 ERROR ** T = 0.0: stop'
            stop
          endif
          j1 = 1
          do 320 j=1,nel
            a1    = d(1)*shp(1,j)*xsj
            a2    = d(1)*shp(2,j)*xsj
            a3    = qq*shp(3,j)*xsj
c....  compute residual
            p(j1) = p(j1) + a1*q1 + a2*q2 - a3
            a3    = a3*dq
c....  compute tangent
            if(isw.eq.3) then
              i1  = 1
              do 310 i=1,nel
                s(i1,j1) = s(i1,j1)+a1*shp(1,i)+a2*shp(2,i)+a3*shp(3,i)
                i1 = i1 + ndf
310           continue
            endif
            j1 = j1 + ndf
320       continue
330     continue
c.... compute the flows in each element
      elseif(isw.eq.4) then
        call shape(0.0d0,0.0d0,xl,shp,xsj,ndm,nel,ix,.false.)
        call flo06(d,shp,ul,q1,q2,qm,uu,ndf,nel)
        rr  = coord(xl(1,1),shp,ndm,nel)
        zz  = coord(xl(2,1),shp,ndm,nel)
        mct = mct - 1
        if(mct.lt.0) then
          mct = 50
          write(iow,2001) head
          if(ior.lt.0) write(*,2001) head
        endif
        write(iow,2002) n,ma,rr,zz,q1,q2,qm,uu
        if(ior.lt.0) write(*,2002) n,ma,rr,zz,q1,q2,qm,uu
c.... compute heat capacity (mass) matrix
      elseif(isw.eq.5) then
        do 520 l=1,lint
          call shape(sg(l),tg(l),xl,shp,xsj,ndm,nel,ix,.true.)
          xsj = xsj*wg(l)
          if(kat.eq.2) xsj = xsj*coord(xl,shp,ndm,nel)
          j1 = 1
          do 510 j=1,nel
            p(j1) = p(j1) + d(9)*shp(3,j)*xsj
            j1 = j1 + ndf
510       continue
520     continue
c.... compute the nodal flow values
      elseif(isw.eq.8) then
        call stcn06(ix,d,xl,ul,shp,aa,aa(numnp+1),ndf,ndm,nel,
     1              numnp,sg,tg,wg,lint)
      endif
c.... formats
1000  format(6f10.0,2i10)
2000  format(3x,'Laplace Element with Reaction Loading'//
     1  4x,'Conductivity  ',e12.4/4x,'Specific Heat ',e12.4/
     x  4x,'Density       ',e12.4/4x,'Load Amplitude',e12.4/
     2  4x,'Reaction Exp. ',e12.4/4x,'Ambient Temp. ',e12.4/
     3  4x,'n - (Temp**n) ',i5/   4x,a12,' analysis')
2001  format(1x,20a4//'  L a p l a c e   E q u a t i o n'//
     1 ' elem  mat    1-coord    2-coord    1-flow     2-flow',
     2 '   max flow     U-value'/)
2002  format(2i5,0p2f11.3,1p4e11.3)
3000  format(' Input:K, c, rho, Q, r, Ta, nn, geom(1=plane,2=axisym)'
     1       /3x,'>',$)
      end
c
      double precision function coord(xl,shp,ndm,nel)
      real*8 xl(ndm,*)
      real*8 shp(3,*)
c
      coord = 0.d0
      do 100 i = 1,nel
        coord = coord + shp(3,i)*xl(1,i)
100   continue
      end
c
      subroutine flo06(d,shp,ul,q1,q2,qm,uu,ndf,nel)

      implicit real*8 (a-h,o-z)
      real*8   shp(3,*),ul(ndf,*)
c.... compute flows at current point
      q1 = 0.0d0
      q2 = 0.0d0
      uu = 0.0d0
      do 100 i = 1,nel
        q1 = q1 - d*shp(1,i)*ul(1,i)
        q2 = q2 - d*shp(2,i)*ul(1,i)
        uu = uu +   shp(3,i)*ul(1,i)
100   continue
      qm =    sqrt(q1*q1 + q2*q2)
      end
c
      subroutine stcn06(ix,d,xl,ul,shp,dt,st,ndf,ndm,nel,numnp,
     1                  sg,tg,wg,lint)

      implicit real*8 (a-h,o-z)

c....  project values to nodes
c     integer*2 ix(*)
      integer   ix(*)
      real*8    dt(numnp),st(numnp,*),xl(ndm,*)
      real*8    d(*),ul(ndf,*),sg(*),tg(*),wg(*),shp(3,*)
      do 110 l = 1,lint
        call shape(sg(l),tg(l),xl,shp,xsj,ndm,nel,ix,.false.)
        call flo06(d,shp,ul,q1,q2,qm,uu,ndf,nel)
        do 100 ii = 1,nel
          ll = abs(ix(ii))
          if(ll.gt.0) then
            xsji     = xsj*shp(3,ii)*wg(l)
            dt(ll)   = dt(ll)   + xsji
            st(ll,1) = st(ll,1) + q1*xsji
            st(ll,3) = st(ll,3) + q2*xsji
            st(ll,4) = st(ll,4) + qm*xsji
            st(ll,5) = st(ll,5) + ul(1,ii)*xsji
          endif
100     continue
110   continue
      end
