      subroutine acheck(x,y,n0,nl)

c     Purpose: Data parser

c     Inputs:
c        x        - Character input data
c        n0       - Field widths for parse
c        nl       - Number of characters to check

c     Outputs:
c        y        - Parsed array

      implicit none

      integer  n0,nl, i,ii,il, k

      character*1 x(nl),y(nl)

      do 100 ii = nl,1,-1
         if(x(ii).ne.' ') go to 110
100   continue

110   do 150 i = 1,nl
         y(i) = ' '
150   continue

      k = 0
      il= 0
      do 200 i = 1,ii
      if(x(i).eq.',') then
         k  = k + n0
         if(k.gt.nl-n0) go to 210
         il = k - i
      else
         y(i+il) = x(i)
      endif
200   continue
      k  = k + n0

210   call just(y,k,n0)
      do 220 i = n0,nl,n0
        if(y(i).eq.' ') y(i) = '0'
220   continue

      end

      subroutine ckisop(ix,xl,shp,ndm)

c     Purpose: Check isoparametric elements for errors

c     Inputs:
c        ix(*)     - Node numbers on element
c        xl(ndm,*) - Nodal coordinates for element
c        shp(3,*)  - Shape function storage
c        ndm       - Spatial dimension of mesh

c     Outputs:
c        none

      implicit  none

      integer   ndm
      integer   i,ineg, l

      integer   xn(9),yn(9),ic(18),ix(*)
      real*8    xl(ndm,*)
      real*8    shp(3,*),ss,tt,xsj

      real*8          dm
      integer            n,ma,mct,iel,nel
      common /eldata/ dm,n,ma,mct,iel,nel

      integer         ior,iow
      common /iofile/ ior,iow

      data xn/-1,1,1,-1,0,1,0,1,0/, yn/-1,-1,1,1,-1,0,1,0,0/

c     check the element for input errors

      ineg = 0
      do 100 l = 1,nel
        if(xl(1,l).eq. -999.0 .and. ix(l).ne.0) then
          ic(ineg+1) = l
          ic(ineg+2) = abs(ix(l))
          ineg = ineg + 2
        endif
100   continue
      if(ineg.gt.0) then
        write(iow,2000) n,(ic(i),i=1,ineg)
        if(ior.lt.0) write(*,2000) n,(ic(i),i=1,ineg)
      else
        do 110 l = 1,nel
          ss = xn(l)
          tt = yn(l)
          call  shape (ss,tt,xl,shp,xsj,ndm,nel,ix,.false.)
          if(xsj.le.0.0d0) then
            ic(ineg+1) = l
            ic(ineg+2) = abs(ix(l))
            ineg = ineg + 2
          endif
  110   continue
        if(ineg.gt.0) then
          write(iow,2001) n,(ic(i),i=1,ineg)
          if(ior.lt.0) write(*,2001) n,(ic(i),i=1,ineg)
        endif
      endif
2000  format(' >Element',i4,' coordinates not input for nodes:'/
     1      ('                Local =',i3,' Global =',i4))
2001  format(' >Element',i4,' has negative jacobian at nodes:'/
     1      ('                Local =',i3,' Global =',i4))

      end

      double precision function dot (a,b,n)

c     Purpose: Compute dot product of two vectors

c     Inputs:
c        a(*)     - Vector 1
c        b(*)     - Vector 2
c        n        - Length of vectors

c     Outputs:
c        dot      - Dot product value

      implicit none

      integer  n, k
      real*8   a(*),b(*)

c     dot product function 

      dot = 0.0d0
      do 10 k=1,n
        dot = dot + a(k)*b(k)
   10 continue

      end

      subroutine elmlib(d,u,x,ix,t,s,p,i,j,k,iel,isw)

c     Purpose: Element library driver

c     Inputs:
c        d(*)      - Material set parameters
c        u(*)      - Element nodal solution parameters
c        x(*)      - Element nodal coordinates
c        ix(*)     - Element node numbers
c        t(*)      - Element nodal temperatures
c        i         - Number dof/node
c        j         - Spatial dimension of mesh
c        k         - Size of element arrays
c        iel       - Element type number
c        isw       - Switch for action to be taken by element routine

c     Outputs:
c        s(k,*)    - Element array
c        p(*)      - Element vector

      implicit  none

      integer   i, j, k, iel, isw
      real*8    d(*),p(*),s(*),u(*)
      real*8    x(*),t(*)
      integer   ix(*)

      integer         ior,iow
      common /iofile/ ior,iow

      if(isw.ge.3.and.isw.le.6) then
        call pconsd(p,k,0.0d0)
        call pconsd(s,k*k,0.0d0)
      endif
      if(iel.eq.1) then
        call elmt01(d,u,x,ix,t,s,p,i,j,k,isw)
      elseif(iel.eq.2) then
        call elmt02(d,u,x,ix,t,s,p,i,j,k,isw)
      elseif(iel.eq.3) then
        call elmt03(d,u,x,ix,t,s,p,i,j,k,isw)
      elseif(iel.eq.4) then
        call elmt04(d,u,x,ix,t,s,p,i,j,k,isw)
      elseif(iel.eq.5) then
        call elmt05(d,u,x,ix,t,s,p,i,j,k,isw)
      elseif(iel.eq.6) then
        call elmt06(d,u,x,ix,t,s,p,i,j,k,isw)
      else
        write(iow,4000) iel
        if(ior.lt.0) write(*,4000) iel
        stop
      endif

4000  format('  **ERROR** Element type number',i3,' found.')

      end

      subroutine pconsd(v,nn,cc)

c     Purpose: Set real*8 array to constant value

c     Inputs:
c        nn        - Length of array
c        cc        - Value to assign

c     Outputs:
c        v(*)      - Array set to 'cc'

      implicit  none

      integer   n,nn
      real*8    v(nn),cc

      do 100 n = 1,nn
        v(n) = cc
100   continue

      end

      subroutine pconsi(iv,nn,ic)

c     Purpose: Set integer array to constant value

c     Inputs:
c        nn        - Length of array
c        ic        - Value to assign

c     Outputs:
c        iv(*)     - Array set to 'ic'

      implicit  none

      integer   n,nn,ic
      integer   iv(nn)

      do 100 n = 1,nn
        iv(n) = ic
100   continue

      end

      subroutine pconsr(v,nn,cr)

c     Purpose: Set real*4 array to constant value

c     Inputs:
c        nn        - Length of array
c        cr        - Value to assign

c     Outputs:
c        v(*)      - Array set to 'cr'

      implicit none

      integer  n,nn
      real     v(nn),cr
      do 100 n = 1,nn
        v(n) = cr
100   continue

      end

      subroutine pdisk(disk,files)

c     Purpose: Set disk name character string

c     Inputs:
c        disk        - Disk label for files (e.g., c for DOS hard)

c     Output:
c        files(*)    - File name with disk label added

      implicit    none

      integer     i
      character*1 disk,files(*)

      i = 1
      if((files(1).ne.' ') .and. (files(2).eq.':')) i = 3
      files(i) = disk

      end

      subroutine pend(subnam)

c     Purpose: End of file found

c     Inputs:
c        subnam    - Name of subroutine where EOF encountered

      implicit  none
      character subnam*6

      integer         ior,iow
      common /iofile/ ior,iow

      if(ior.gt.0) write(iow,4000) subnam
      if(ior.lt.0) write(  *,4000) subnam
      stop

4000  format(' ** ERROR in ',a6,' ** end of file encountered')

      end

      subroutine perror(subnam,yy)

c     Purpose: Read error encountered

c     Inputs:
c        subnam    - Name of subroutine where EOF encountered
c        yy        - String of data read

      implicit  none
      character subnam*6, yy*80

      integer         ior,iow
      common /iofile/ ior,iow

      if(ior.gt.0) then
        write(iow,4000) subnam,yy
        stop
      else
        write(*,4000) subnam
      endif

4000  format(' **ERROR in ',a6,'** Reinput last record:'/1x,a80)

      end

      subroutine pgauss(l,lint,r,z,w)

c     Purpose: Gauss points and weights for two dimensions

c     Inputs:
c        l         - Order of quadrature

c     Outputs:
c        lint      - Number of points
c        r(*),z(*) - Points
c        w(*)      - Weights

      implicit  none

      integer   i,l,lint
      integer   lr(9),lz(9),lw(9)
      real*8    r(*),z(*),w(*),g,h

      real*8          dm
      integer            n,ma,mct,iel,nel
      common /eldata/ dm,n,ma,mct,iel,nel

      data lr/-1,1,1,-1,0,1,0,-1,0/,lz/-1,-1,1,1,-1,0,1,0,0/
      data lw/4*25,4*40,64/

      lint = l*l

c     1x1 integration

      if(l.eq.1) then
        r(1) = 0.
        z(1) = 0.
        w(1) = 4.

c     2x2 integration

      elseif(l.eq.2) then
        g = 1.0/sqrt(3.d0)
        do 21 i = 1,4
          r(i) = g*lr(i)
          z(i) = g*lz(i)
          w(i) = 1.
21      continue

c     3x3 integration

      elseif(l.eq.3) then
        g = sqrt(0.60d0)
        h = 1.0/81.0d0
        do 31 i = 1,9
          r(i) = g*lr(i)
          z(i) = g*lz(i)
          w(i) = h*lw(i)
31      continue
      endif

      end

      subroutine phelp(wd,nwd,wrd,isw)

c     Purpose: Help file for macro command list

c     Inputs:
c        wd(*)     - List of FEAP command names
c        nwd       - Number of command names
c        wrd       - Request name for help
c        isw       - Switch: =1 for solution commands; otherwise mesh.

      implicit  none
      integer   nwd,isw
      character wd(nwd)*4,wrd*5

      integer         ior,iow
      common /iofile/ ior,iow

      if(ior.gt.0) return
      if(isw.eq.1) write(*,2000)
      if(isw.ne.1) write(*,2001) wrd
      write(*,2002) wd
      write(*,2003) wrd

2000  format(//' The following macro commands are available'//
     &         ' use of loop must terminate with a matching next'//
     &         ' multiple loop-next pairs may occur up to depth of 8')
2001  format(//' The following ',a5,'commands are available:')
2002  format(/8(3x,a4))
2003  format(/' Terminate ',a5,' execution with an "end" command.')

      end

      subroutine psetm(na,nl,typ,afl)

c     Purpose:  Allocate array storage pointer from blank common
c
c     Inputs:
c        nl        - length of array requested
c        typ       - type of array requested ('i', 'r', 'd')

c     Outputs:
c        na        - Pointer to first word in array
c        afl       - Flag, error if true

      implicit  none

      logical   afl
      integer   na,nl, np,ns
      real*8    amx
      character typ*1

      integer         iodr,iodw,ipd,ipr,ipi
      common /iofild/ iodr,iodw,ipd,ipr,ipi

      integer         ior,iow
      common /iofile/ ior,iow

      integer         maxm,ne
      common /psize/  maxm,ne

c     Set data management pointers for arrays

c     alignment for the blank common d(1),r(1),i(xx): xx = size
c     precisions: 'd' = ipd ; 'r' = ipr ; 'i' = ipi

c     N.B. ipd >= ipr >= ipi

      if( typ .eq. 'd' ) then
        np = ipd
        ns = 0
      elseif( typ .eq. 'r' ) then
        np = ipr
        ns = ipd/ipr
      elseif( typ .eq. 'i' ) then
        np = ipi
        ns = (ipd+ipr)/ipi
      endif
      na = ne
      ne = na + nl*np + mod(ipd - mod(nl*np,ipd),ipd)
      na   = (na + np - 1)/np - ns
      afl = .false.
      amx = maxm
      amx = ne/amx
      if(amx.gt.0.90) write(*,1001) ne,maxm,amx
      if(ne.le.maxm) return
      write(iow,1000) ne,maxm
      if(ior.lt.0) write(*,1000) ne,maxm
      stop

1000  format(2x,'**ERROR** Insufficient storage in blank common'/
     &      10x,'Required  =',i6/10x,'Available =',i6/)
1001  format('  **Memory warning** Used =',i6,' Avail =',i6,' % =',f6.3)

      end

      subroutine pstres(sig,p1,p2,p3)

c     Purpose: Principal stresses (2 dimensions): sig = sig-xx,sig-xy,sig-yy

c     Inputs:
c        sig(*)   - Stress state (as above)

c     Outputs:
c        p1,p2    - Principal stresses
c        p3       - Angle p1 makes with x-1 axis

      implicit none

      real*8   sig(3),p1,p2,p3,xi1,xi2,rho

      xi1 = (sig(1) + sig(3))/2.
      xi2 = (sig(1) - sig(3))/2.
      rho = sqrt(xi2*xi2 + sig(2)*sig(2))
      p1 = xi1 + rho
      p2 = xi1 - rho
      p3 = 45.0
      if(xi2.ne.0.0d0) p3 = 22.5*atan2(sig(2),xi2)/atan(1.0d0)

      end

      subroutine saxpb (a,b,x,n,c)

c     Purpose:  Scalar times vector plus vector

c     Inputs:
c        a(*)      - Vector to be multiplied by scalar
c        b(*)      - Vector to add
c        x         - Scalar value multiplying 'a'
c        n         - Length of vectors

c     Outputs:
c        c(*)      - Result of operation
      implicit none

      integer  k,n
      real*8   a(*),b(*),c(*),x

c    vector times scalar added to second vector 

      do 10 k=1,n
        c(k) = a(k)*x +b(k)
   10 continue

      end

      subroutine shap2(s,t,shp,ix,nel)

c     Purpose: Add quadratic shape functions as necessary

c     Inputs:
c       s,t     - Natural coordinate point for evaluation
c       ix(*)   - Element connection list
c       nel     - Number of nodes on element

c     Outputs:
c       shp(3,*) - Quadratic shape functions generated

      implicit  none

      integer   i, j, k, l, nel
      real*8    shp(3,1),s,t,s2,t2

      integer   ix(*)

      s2 = (1.-s*s)/2.
      t2 = (1.-t*t)/2.
      do 120 i=5,9
        do 100 j = 1,3
          shp(j,i) = 0.0d0
100     continue
120   continue

c     midside nodes (serendipity)

      if(ix(5).eq.0) go to 101
      shp(1,5) = -s*(1.-t)
      shp(2,5) = -s2
      shp(3,5) = s2*(1.-t)
101   if(nel.lt.6) go to 108
      if(ix(6).eq.0) go to 102
      shp(1,6) = t2
      shp(2,6) = -t*(1.+s)
      shp(3,6) = t2*(1.+s)
102   if(nel.lt.7) go to 108
      if(ix(7).eq.0) go to 103
      shp(1,7) = -s*(1.+t)
      shp(2,7) = s2
      shp(3,7) = s2*(1.+t)
103   if(nel.lt.8) go to 108
      if(ix(8).eq.0) go to 104
      shp(1,8) = -t2
      shp(2,8) = -t*(1.-s)
      shp(3,8) = t2*(1.-s)

c     interior node (lagrangian)

104   if(nel.lt.9) go to 108
      if(ix(9).eq.0) go to 108
      shp(1,9) = -4.*s*t2
      shp(2,9) = -4.*t*s2
      shp(3,9) = 4.*s2*t2

c     correct edge nodes for interior node (lagrangian)

      do 107 j= 1,3
        do 105 i = 1,4
          shp(j,i) = shp(j,i) - 0.25*shp(j,9)
105     continue
        do 106 i = 5,8
          if(ix(i).ne.0) shp(j,i) = shp(j,i) - .5*shp(j,9)
106     continue
107   continue

c     correct corner nodes for presense of midside nodes

108   do 110 i = 1,4
        k = mod(i+2,4) + 5
        l = i + 4
        do 109 j = 1,3
          shp(j,i) = shp(j,i) - 0.5*(shp(j,k)+shp(j,l))
109     continue
110   continue

      end

      subroutine shape(ss,tt,xl,shp,xsj,ndm,nel,ix,flg)

c     Purpose: Shape function routine for two dimensional elements

c     Inputs:
c       ss,tt     - Natural coordinates for evaluation
c       xl(ndm,*) - Nodal coordinates for element
c       ndm       - Spatial dimension of mesh
c       nel       - Number of nodes on element
c       ix(*)     - Element nodal connections
c       flg       - Flag, compute globagl derivatives if true

c     Outputs:
c       shp(3,*)  - Shape functions and derivatives
c       xsj       - Jacobian determinant at point

      implicit none

      logical  flg
      integer  i, j, k, ndm,nel
      real*8   xl(ndm,1),s(4),t(4)
      real*8   shp(3,*),xs(2,2),sx(2,2),ss,tt,xsj,tp
      integer  ix(*)

      data s/-.5d0,.5d0,.5d0,-.5d0/, t/-.5d0,-.5d0,.5d0,.5d0/

c     form 4-node quadrilateral shape functions

      do 100 i = 1,4
        shp(3,i) = (0.5+s(i)*ss)*(0.5+t(i)*tt)
        shp(1,i) = s(i)*(0.5+t(i)*tt)
        shp(2,i) = t(i)*(0.5+s(i)*ss)
100   continue
      if(nel.eq.3) then

c     form triangle by adding third and fourth together

        do 110 i = 1,3
          shp(i,3) = shp(i,3)+shp(i,4)
110     continue
      endif

c     add quadratic terms if necessary

      if(nel.gt.4) call shap2(ss,tt,shp,ix,nel)

c     construct jacobian and its inverse

      do 125 i = 1,2
      do 125 j = 1,2
        xs(i,j) = 0.0
        do 120 k = 1,nel
          xs(i,j) = xs(i,j) + xl(i,k)*shp(j,k)
120     continue
125   continue
      xsj = xs(1,1)*xs(2,2)-xs(1,2)*xs(2,1)
      if(flg) return
      if(xsj.le.0.0d0) xsj = 1.0
      sx(1,1) = xs(2,2)/xsj
      sx(2,2) = xs(1,1)/xsj
      sx(1,2) =-xs(1,2)/xsj
      sx(2,1) =-xs(2,1)/xsj

c     form global derivatives

      do 130 i = 1,nel
        tp        = shp(1,i)*sx(1,1)+shp(2,i)*sx(2,1)
        shp(2,i)  = shp(1,i)*sx(1,2)+shp(2,i)*sx(2,2)
        shp(1,i) = tp
130   continue

      end
