      subroutine pmesh(idl,ie,d,id,x,ix,f,t,ndf,ndm,nen1,iii,prt)

c     Purpose: Data input routine for mesh description

c     Inputs:
c        idl(*)     - Working vector
c        ndf        - Number dof/node
c        ndm        - Spatial dimension of mesh
c        nen1       - Dimension of ix array
c        iii        - Flag to initialize arrays
c        prt        - Flag, output results if true

c     Outputs:
c        ie(9,*)    - Material set assembly data
c        d(*)       - Material set parameters
c        id(ndf,*)  - Boundary condition codes
c        x(ndm,*)   - Nodal coordinates of mesh
c        ix(nen1,*) - Element nodal connection lists
c        f(ndf,*)   - Nodal force/displacements
c        t(*)       - Nodal temperatures

      implicit    none

      logical     prt,error,pcomp
      integer     ndf,ndm,nen1
      integer     i,iii, list

      integer     ie(9,*),id(ndf,*),ix(nen1,*),idl(6)
      real*8      x(ndm,*),f(ndf,*),t(*)
      real*8      d(18,*)
      character*4 wd(13),cc,yyy*80,cds*12,tmp*12,fds*12,an*12

      integer        numnp,numel,nummat,nen,neq
      common /cdata/ numnp,numel,nummat,nen,neq

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

      integer         ior,iow
      common /iofile/ ior,iow

      integer        n11a,n11b,n11c,ia
      common /mdat2/ n11a,n11b,n11c,ia(2,11)

      real*8 dm
      common dm(1000)

      data wd/'coor','elem','mate','boun','forc','temp','prin','nopr',
     &        'bloc','pola','angl','ebou','end '/

      data an/'  angles    '/, cds/' coordinates'/, tmp/' temperature'/

      data fds/' force/displ'/, list/13/

c     Initialize arrays

      error = .false.
      if(iii.ge.0) then
        prt = .true.
        call pconsd(dm(n11b),numnp,0.0d0)
        call pconsd(   f,numnp*ndf,0.0d0)
        call pconsi(id,numnp*ndf,0)
        if(iii.eq.0) call pconsd( t,numnp,0.0d0)
        if(iii.eq.0) call pconsd( x,numnp*ndm,-999.0d0)
      endif
102   if(ior.lt.0) write(*,2000)
      call pintio(yyy,10)
      read(yyy,1000,err=901,end=900) cc
      if((ior.lt.0).and.pcomp(cc,'help')) then
         call phelp(wd,list,'MESH ',0)
         go to 102
      endif
20    do 30 i = 1,list
30    if(pcomp(cc,wd(i))) go to 40
      go to 102

c     Nodal coordinate data input

40    if(i.eq.1) call genvec(ndm,x,cds,prt,error,.true.)

c     Element data input

      if(i.eq.2) call pelin(idl,ix,nen1,nen,numnp,numel,error,prt)

c     Material data input

      if(i.eq.3) call pmatin(d,x,ix,idl,ie,nummat,ndm,ndf,prt)

c     Read in restraint conditions for each node

      if(i.eq.4) call pbcin(iii,idl,id,numnp,ndf,prt)

c     Force/displ data input

      if(i.eq.5) call genvec(ndf,f,fds,prt,error,.false.)

c     Temperature data input

      if(i.eq.6) call genvec(1,t,tmp,prt,error,.false.)

c     Set print flag

      if(i.eq.7.or.i.eq.8) prt = i.eq.7

c     Generate block of nodes and 4-node elements

      if(i.eq.9) then
        if(iii.lt.0) write(iow,3000)
        call blkgen(ndm,nen,nen1,x,ix,prt)
      endif

c     Convert polar/cylindrical to cartesian coordinates

      if(i.eq.10) call polar(x,ndm,prt)

c     Set boundary angles

      if(i.eq.11) call genvec(1,dm(n11b),an,prt,error,.false.)

c     Set edge boundary values

      if(i.eq.12) call pedges(iii,x,id,ndm,ndf,numnp,prt)

c     End of mesh

      if(i.eq.13) then
        if(error) stop
        return
      endif
      go to 102

c     End of file encountered

900   call pend('pmesh ')
      stop
901   call perror('PMESH ',yyy)
      go to 102
1000  format(a4)
2000  format('     Mesh 1 > ',$)
3000  format(' **WARNING** Element connections necessary to use '
     &      ,'block in macro program')

      end

      subroutine pmatin(d,x,ix,idl,ie,nummat,ndm,ndf,prt)

c     Purpose: Input material set data

c     Inputs:
c        idl(*)     - Working vector
c        nummat     - Number of material sets
c        ndm        - Spatial dimension of mesh
c        ndf        - Number dof/node
c        prt        - Flag, output results if true

c     Outputs:
c        d(*)       - Material set parameters
c        x(ndm,*)   - Nodal coordinates of mesh
c        ix(nen1,*) - Element nodal connection lists
c        ie(9,*)    - Material set assembly data

      implicit  none

      integer   maxa
      parameter(maxa=32000)

      integer   ndm,nummat,ndf
      integer   i

      logical   prt
      character yyy*80

      integer   ix(*),idl(*),ie(9,*)

      real*8    x(ndm,*)
      real*8    d(18,*)

      real*8          ad
      common /adata/  ad(maxa)

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

      integer         nh1,nh2
      common /hdata/  nh1,nh2

      integer         ior,iow
      common /iofile/ ior,iow

c     Material data input

      if(prt) then
        call prthed(iow)
        write(iow,2001)
        if(ior.lt.0) write(*,2001)
      endif
      do 306 n = 1,nummat
301     if(ior.lt.0) write(*,3000)
        call pintio(yyy,10)
        read(yyy,1000,err=311) ma,iel,(idl(i),i=1,ndf)
        if(ma.le.0) return

c       Set all zero inputs

        do 302 i = 1,ndf
302     ie(i,ma) = idl(i)
        do 303 i = 1,ndf
        if(idl(i).ne.0) go to 305
303     continue
        do 304 i = 1,ndf
304     ie(i,ma) = i
305     ie(7,ma) = iel
        mct = 0
        nh1 = 0
        if(prt) then
          write(iow,2000) ma,iel,(i,ie(i,ma),i=1,ndf)
          if(ior.lt.0) write(*,2000) ma,iel,(i,ie(i,ma),i=1,ndf)
        endif
        call elmlib(d(1,ma),ad,x,ix,dq,ad,ad,ndf,ndm,ndf,iel,1)
        if(nh1.eq.0 .and. mct.ne.0) nh1 = mct
        ie(8,ma) = nh1
        go to 306
311     call perror('PMESH ',yyy)
        go to 301
306   continue

c     Formats

1000  format(8i10)
2000  format(/5x,'Material Set',i3,' for Element Type',i2,5x,//
     &   10x,'degree of freedom assignments    local    global' /
     &   42x, 'number',4x,'number'/(36x,2i10))
2001  format('    M a t e r i a l    P r o p e r t i e s')
3000  format(' Input: matl. no., elmt type'/3x,'>',$)

      end

      subroutine blkgen(ndm,nel,nel1,x,ix,prt)

c     Purpose: Generate a block of nodes/elements

c     Inputs:
c        ndm        - Spatial dimension of mesh
c        nel        - Number of nodes/element
c        nel1       - Dimension of ix array
c        prt        - Flag, output results if true

c     Outputs:
c        x(ndm,*)   - Nodal coordinates of mesh
c        ix(nen1,*) - Element nodal connection lists

      implicit  none

      logical   prt
      character yyy*80
      integer   ndm,nel,nel1
      integer   i, j, k, l, ma, n, ne,nf,ng,ni,nm,nn,nr,ns,ntyp,nodinc
      real*8    r, s, t

      integer   ix(nel1,*),ixl(9)
      real*8    x(ndm,*),xl(3,9)
      real*8    shp(3,9),dr,ds

      integer        numnp,numel,nummat,nen,neq
      common /cdata/ numnp,numel,nummat,nen,neq

      integer         ior,iow
      common /iofile/ ior,iow

100   if(ior.lt.0) write(*,5000)
      call pintio(yyy,5)
      read(yyy,1000,err=101) nn,nr,ns,ni,ne,ma,nodinc,ntyp
      go to 102
101   call perror('BLKGEN',yyy)
      go to 100
102   nodinc = max(nodinc,0)
      nr = max(nr,1)
      ns = max(ns,1)
      ni = max(ni,1)
      ma = max(ma,1)
      if(prt) then
        call prthed(iow)
        write(iow,2000) nr,ns,ni,ne,ma,nodinc,ntyp
        if(ne.eq.0) write(iow,2005)
        write(iow,2002) (i,i=1,ndm)
        if(ior.lt.0) then
          write(*,2000) nr,ns,ni,ne,ma,nodinc,ntyp
          if(ne.eq.0) write(*,2005)
          write(*,2002) (i,i=1,ndm)
        endif
      endif
      do 10 n = 1,9
      do 10 j = 1,ndm
      xl(j,n) = 0.0
10    ixl(n) = 0
      nm = 0
      do 20 n = 1,nn
200   if(ior.lt.0) write(*,5001)
      call pintio(yyy,10)
      read(yyy,1001,err=201) l,r,s,t
      go to 202
201   call perror('BLKGEN',yyy)
      go to 200
202   if(l.eq.0) l = n
      nm = max(nm,l)
      ixl(l) = l
      xl(1,l) = r
      xl(2,l) = s
      xl(3,l) = t
20    if(prt) write(iow,2001) l,(xl(i,l),i=1,ndm)
      if(prt.and.ior.lt.0) write(*,2001) l,(xl(i,l),i=1,ndm)
      dr = 2.d0/nr
      ds = 2.d0/ns
      if (ntyp.eq.0) then
         nf = ne + nr*ns - 1
      elseif (ntyp.gt.7) then
         nf = ne + (nr*ns)/4 - 1
      else
         nf = ne + 2*nr*ns - 1
      endif
      if(nf.gt.numel.and.ne.gt.0) go to 401
      nr = nr + 1
      ns = ns + 1
      if(ndm.eq.1) ns = 1
      ng = nr*ns + ni -1
      if(ng.gt.numnp) go to 400

c     Form block

      call sblk(nr,ns,xl,ixl,shp,x,ix,dr,ds,ni,ne,ndm,
     &     nel1,nodinc,ntyp,nm,ma,prt)

c     Print lists if wanted

      if(.not.prt) return

c     Print element lists

      if(ne.le.0) return
      do 502 n = ne,nf,50
        call prthed(iow)
        write(iow,2003) (i,i=1,nel)
        if(ior.lt.0) write(*,2003) (i,i=1,nel)
        j = min(nf,n+49)
        do 501 i = n,j
          write(iow,2004) i,ma,(ix(k,i),k=1,nel)
          if(ior.lt.0) write(*,2004) i,ma,(ix(k,i),k=1,nel)
501     continue
502   continue
      return
400   write(iow,2006) ng,numnp
      if(ior.lt.0) write(*,2006) ng,numnp
      return
401   write(iow,2007) nf,numel
      if(ior.lt.0) write(*,2007) nf,numel

1000  format(8i5)
1001  format(i10,3f10.0)
2000  format('    N o d e   G e n e r a t i o n s'//
     &9x,'Number of r-increments:',i5/9x,'Number of s-increments:',i5/
     &9x,'First node number     :',i5/9x,'First element number  :',i5/
     &9x,'Element material type :',i5/9x,'Node line increment   :',i5/
     &   9x,'Block type (0-9)       :',i5/1x)
2001  format(i9,1p3e12.3)
2002  format(5x,'node',3(i6,' coord'))
2003  format('    E l e m e n t   C o n n e c t i o n s'//
     & 3x,'elmt  matl',8(i3,'-node')/(13x,8(i3,'-node')))
2004  format(2i6,8i8/(13x,8i8))
2005  format(' **WARNING** No elements are generated ')
2006  format(' **ERROR** insufficient storage for nodes'/
     &   10x,'final node =',i5,5x,'numnp =',i5)
2007  format(' **ERROR** insufficient storage for elements'/
     &   10x,'final element =',i5,5x,'numel =',i5)
5000  format(' Input: nn,nr,ns,ni,ne,ma,nodinc,ntyp'/3x,'>',$)
5001  format(' Input: node, x-1, x-2, x-3'/3x,'>',$)

      end

      subroutine genvec(ndm,x,cd,prt,err,prtz)

c     Purpose: Generate real data arrays by linear interpolation

c     Inputs:
c        ndm        - Dimension of x array
c        cd         - Character label for outputs
c        prt        - Flag, output results if true
c        prtz       - Flag, output zero results if true

c     Outputs:
c        x(ndm,*)   - Nodal coordinates of mesh
c        err        - Flag, true if errors occurred

      implicit  none

      logical   prt,err,prtz
      character cd*12,yyy*80
      integer   ndm
      integer   i, j, l,li,lg, mct, n,ng

      real*8    x(ndm,*),xl(6)

      integer        numnp,numel,nummat,nen,neq
      common /cdata/ numnp,numel,nummat,nen,neq

      integer         ior,iow
      common /iofile/ ior,iow

      mct = 0
      n = 0
      ng = 0
102   l = n
      lg = ng
100   if(ior.lt.0) write(*,5000) cd
      call pintio(yyy,10)
      read(yyy,1000,err=11) n,ng,(xl(i),i=1,ndm)
      go to 12
11    call perror('GENVEC',yyy)
      go to 100
12    if(n.gt.numnp) write(iow,3001) n,cd
      if(ior.lt.0.and.n.gt.numnp) write(*,3001) n,cd
      if(n.le.0.or.n.gt.numnp) go to 109
      do 103 i = 1,ndm
103   x(i,n) = xl(i)
      if(lg) 104,102,104
104   lg = sign(lg,n-l)
      li =(abs(n-l+lg)-1)/abs(lg)
      do 105 i = 1,ndm
105   xl(i) = (x(i,n)-x(i,l))/li
106   l = l + lg
      if((n-l)*lg.le.0) go to 102
      if(l.le.0.or.l.gt.numnp) go to 108
      do 107 i = 1,ndm
107   x(i,l) = x(i,l-lg) + xl(i)
      go to 106
108   write(iow,3000) l,cd
      if(ior.lt.0) write(*,3000) l,cd
      err = .true.
      go to 102
109   if(.not.prt) return
      do 113 j = 1,numnp
      if(prtz) go to 111
      do 110 l = 1,ndm
      if(x(l,j).ne.0.0) go to 111
110   continue
      go to 113
111   mct = mct - 1
      if(mct.gt.0) go to 112
      mct = 50
      write(iow,2000) cd,(l,cd,l=1,ndm)
      if(ior.lt.0) write(*,2000) cd,(l,cd,l=1,ndm)
112   if(x(1,j).eq.-999.0) write(iow,2001) j
      if(x(1,j).ne.-999.0) write(iow,2002) j,(x(l,j),l=1,ndm)
      if(ior.lt.0) then
        if(x(1,j).eq.-999.0) write(*,2001) j
        if(x(1,j).ne.-999.0) write(*,2002) j,(x(l,j),l=1,ndm)
      endif
113   continue

1000  format(2i10,6f10.0)
2000  format('    N o d a l: ',a12//6x,'node',9(i7,a6))
2001  format(i10,' has not been input or generated')
2002  format(i10,9f13.4)
3000  format(' **ERROR** attempt to generate node',i5,' in ',a12)
3001  format(' **ERROR** attempt to input node',i5,', terminate'
     & ,' input of nodes in ',a12)
5000  format(' Input ',a12,' values: node, inc, value(i),i=1,nval'/
     &       3x,'>',$)

      end

      subroutine pbcin(iii,idl,id,numnp,ndf,prt)

c     Purpose: Input restraint conditions for each node

c     Inputs:
c        iii        - Flag to test change in b.c. inputs
c        idl(*)     - Working vector
c        numnp      - Number nodes in mesh
c        ndf        - Number dof/node
c        prt        - Flag, output results if true

c     Outputs:
c        id(ndf,*)  - Boundary condition codes

      implicit  none

      logical   prt
      character yyy*80
      integer   iii,numnp,ndf
      integer   i, l,lg, n,ng

      integer   idl(*),id(ndf,*)

      integer         ior,iow
      common /iofile/ ior,iow

c     Read in restraint conditions for each node

      iii = 1
      n = 0
      ng = 0
400   l = n
      lg = ng
401   if(ior.lt.0) write(*,5000)
      call pintio(yyy,10)
      read(yyy,1000,err=412) n,ng,(idl(i),i=1,ndf)
      if(n.gt.0.and.n.le.numnp) then
        do 402 i = 1,ndf
          id(i,n) = idl(i)
          if(l.ne.0.and.idl(i).eq.0.and.id(i,l).lt.0) id(i,n) = -1
402     continue
        lg = sign(lg,n-l)
403     l = l + lg
        if((n-l)*lg.le.0) go to 400
        do 404 i = 1,ndf
          if(id(i,l-lg).lt.0) id(i,l) = -1
404     continue
        go to 403
      endif

c     Output nodes with nonzero codes

      if(prt) then
        call prthed(iow)
        write(iow,2000) (i,i=1,ndf)
        do 407 n = 1,numnp
          do 406 l = 1,ndf
            if(id(l,n).ne.0) then
              write(iow,2001) n,(id(i,n),i=1,ndf)
              go to 407
            endif
406       continue
407     continue
      endif
      return
412   call perror('PBCIN ',yyy)
      go to 401

c     Formats

1000  format(8i10)
2000  format('    N o d a l   B. C.'//6x,'node',8(i3,' b.c.')/1x)
2001  format(i10,8i8)
5000  format(' Input: node, inc, b.c. codes(i),i=1,ndf'/3x,'>',$)

      end

      subroutine pelin(idl,ix,nen1,nen,numnp,numel,error,prt)

c     Purpose: Element data input

c     Inputs:
c        idl(*)     - Working vector
c        nen1       - Dimension of ix array
c        nen        - Number nodes/element
c        numnp      - Number of nodes in mesh
c        numel      - Number of elements in mesh
c        prt        - Flag, output results if true

c     Outputs:
c        ix(nen1,*) - Element nodal connection lists
c        error      - Flag, true if error occurs

      implicit  none

      logical   error,prt
      character yyy*80
      integer   nen1,nen,numnp,numel
      integer   i, j, k, l,lk,lx, n,nx

      integer   idl(*),ix(nen1,*)

      integer         ior,iow
      common /iofile/ ior,iow

c     Element data input

      l = 0
      do 210 i = 1,numel,50
      if(prt) call prthed(iow)
      if(prt) write(iow,2001) (k,k=1,nen)
      if(ior.lt.0.and.prt) write(  *,2001) (k,k=1,nen)
      j = min(numel,i+49)
      do 209 n = i,j
      if(l-n) 200,202,204
100   call perror('PELIN ',yyy)
200   if(ior.lt.0) write(*,5000)
      call pintio(yyy,5)
      read(yyy,1001,err=100) l,lk,(idl(k),k=1,nen),lx
      if(l.eq.0) l = numel+1
      if(lx.eq.0) lx=1
      if(l-n) 201,202,204
201   write(iow,3001) l,n
      if(ior.lt.0) write(*,3001) l,n
      error = .true.
      go to 209
202   nx = lx
      do 203 k = 1,nen
      if(idl(k).gt.numnp.or.idl(k).lt.0) go to 208
      ix(k,l) = idl(k)
203   continue
      ix(nen1,l) = lk
      go to 206
204   ix(nen1,n) = ix(nen1,n-1)
      do 205 k = 1,nen
      ix(k,n) = ix(k,n-1) + nx
      if(ix(k,n-1).eq.0) ix(k,n) = 0
      if(ix(k,n).gt.numnp.or.ix(k,n).lt.0) go to 208
205   continue
206   if(prt .and. ior.lt.0) then
        write(  *,2002) n,ix(nen1,n),(ix(k,n),k=1,nen)
      endif
      if(prt) write(iow,2002) n,ix(nen1,n),(ix(k,n),k=1,nen)
      go to 209
208   write(iow,3002) n
      error = .true.
209   continue
210   continue

c     Formats

1001  format(16i5)
2001  format('    E l e m e n t   C o n n e c t i o n s'//
     & 3x,'elmt  matl',8(i3,'-node')/(13x,8(i3,'-node')))
2002  format(2i6,8i8/(13x,8i8))
3001  format(' **ERROR** element',i5,' appears after element',i5)
3002  format(' **ERROR** element',i5,' has illegal nodes')
5000  format(' Input: elm, mat, ix(i),i=1,nen, inc'/3x,'>',$)

      end

      subroutine polar(x,ndm,prt)

c     Purpose: Convert polar to cartesian coordinates

c     Inputs:
c        x(ndm,*)   - Polar form of nodal coordinates
c        ndm        - Spatial dimension of mesh
c        prt        - Flag, output results if true

c     Outputs:
c        x(ndm,*)   - Cartesian form of nodal coordinates

      implicit  none

      logical   prt
      character yyy*80
      integer   ndm
      integer   i,inc, mct, n,ni,ne
      real*8    x0,y0, r

      real*8    th
      real*8    x(ndm,*)

      integer        numnp,numel,nummat,nen,neq
      common /cdata/ numnp,numel,nummat,nen,neq

      integer         ior,iow
      common /iofile/ ior,iow

      if(ndm.eq.1) return
      mct = 0
      th = atan(1.0d0)/45.0d0
100   if(ior.lt.0) write(*,5000)
      call pintio(yyy,10)
      read(yyy,1000,err=101) ni,ne,inc,x0,y0
      go to 102
101   call perror('POLAR ',yyy)
      go to 100
102   if(ni.le.0) return
      if(ni.gt.numnp.or.ne.gt.numnp) go to 300
      inc = sign(max(abs(inc),1),ne-ni)
      if(ne.eq.0) ne = ni
      n = ni
200   r = x(1,n)
      x(1,n) = x0 + r*cos(x(2,n)*th)
      x(2,n) = y0 + r*sin(x(2,n)*th)
      if(mct.gt.0) go to 250
      if(prt) call prthed(iow)
      if(prt) write(iow,2000) x0,y0,(i,i=1,ndm)
      if(ior.lt.0.and.prt) write(*,2000) x0,y0,(i,i=1,ndm)
      mct = 50
250   if(prt) write(iow,2001) n,(x(i,n),i=1,ndm)
      if(ior.lt.0.and.prt) write(*,2001) n,(x(i,n),i=1,ndm)
      mct = mct - 1
      n = n + inc
      if((ne-n)*inc.ge.0) go to 200
      if(mod(ne-ni,inc).eq.0) go to 100
      ni = ne
      n = ne
      go to 200

c     Error

300   write(iow,3000) ni,ne
      if(ior.lt.0) write(*,3000) ni,ne
      stop

c     Formats

1000  format(3i10,2f10.0)
2000  format('    P o l a r   t o   C a r t e s i a n   C o o r d s.'/
     & 8x,'Center: x0 = ',e12.4,' y0 = ',e12.4/6x,'node',6(i7,'-coord'))
2001  format(i10,6f13.4)
3000  format(' **ERROR** attempt to convert nodes ni= ',i6,' - ne= ',i6)
5000  format(' Input: 1-node, 2-node, inc, x1(cent), x2(cent)'/3x,'>',$)

      end

      subroutine sblk(nr,ns,xl,ixl,shp,x,ix,dr,ds,ni,ne,ndm,nel1,
     &                nodinc,ntyp,nm,ma,prt)

c     Inputs:
c        nr         - Number increments in 'r' direction
c        ns         - Number increments in 's' direction
c        xl(ndm,*)  - Master element nodal coordinates
c        ixl(*)     - Master element nodal connections list
c        shp(3,*)   - Shape function storage
c        dr         - 'r' coordinate increment
c        ds         - 's' coordinate increment
c        ni         - Initial node number
c        ne         - Initial element number
c        ndm        - Spatial dimension of mesh
c        nel1       - Dimension of ix array
c        nodinc     - Increment to node numbers at end of 'r' line
c        ntyp       - Block type
c        nm         - Number of nodes on master element
c        ma         - Material set number for block
c        prt        - Flag, output results if true

c     Outputs:
c        x(ndm,*)   - Nodal coordinates of mesh
c        ix(nel1,*) - Element nodal connection lists

      implicit  none

      logical prt
      integer nr,ns,ni,ne,ndm,nel1,nodinc,ntyp,nm,ma
      integer i,inc, j, k, l, m,me,mct, n

      real*8  r,s,dr,ds,xsj

      real*8  xl(3,*),x(ndm,*)
      real*8  shp(3,*)

      integer   ixl(*),ix(nel1,*)

      integer        numnp,numel,nummat,nen,neq
      common /cdata/ numnp,numel,nummat,nen,neq

      integer         ior,iow
      common /iofile/ ior,iow

      n = ni
      mct = 0
      s = -1.0
      do 200 j = 1,ns
      r = -1.0
      do 100 i = 1,nr
      call shape(r,s,xl,shp,xsj,3,nm,ixl,.true.)
      do 55 l = 1,ndm
        x(l,n) = 0.0
        do 50 k = 1,9
          m = ixl(k)
          if(m.gt.0) x(l,n) = x(l,n) + shp(3,m)*xl(l,m)
50      continue
55    continue
      if(prt) then
         mct = mct + 1
         if(mod(mct,50).eq.1) then
           call prthed(iow)
           write(iow,2000) (k,k=1,ndm)
           if(ior.lt.0) write(*,2000) (k,k=1,ndm)
         endif
         write(iow,2001) n,(x(k,n),k=1,ndm)
         if(ior.lt.0) write(*,2001) n,(x(k,n),k=1,ndm)
      endif
      n = n + 1
100   r = r + dr
      n = n + nodinc
200   s = s + ds
      if(ne.le.0) return
      me = ne - 1
      n = ni
      inc = 1
      if(ntyp.ge.8) inc = 2
      do 400 j = 1,ns-1,inc
      do 300 i = 1,nr-1,inc
      n = n + 1
      me = me + 1
      ix(nel1,me) = ma
      if(ntyp.eq.0) then
         ix(1,me)    = n - 1
         ix(2,me)    = n
         if(ndm.ne.1) then
            ix(3,me)    = n + nr + nodinc
            ix(4,me)    = n + nr - 1 + nodinc
         endif
      elseif(ntyp.ge.8) then
         ix(1,me) = n-1
         ix(5,me) = n
         ix(2,me) = n+1
         ix(8,me) = nr+nodinc + n-1
         if(ntyp.gt.8) ix(9,me) = nr+nodinc + n
         ix(6,me) = nr+nodinc + n+1
         ix(4,me) = 2*(nr+nodinc) + n-1
         ix(7,me) = 2*(nr+nodinc) + n
         ix(3,me) = 2*(nr+nodinc) + n+1
         n = n+1
      endif
300   continue
400   n = n + (inc-1) * nr + nodinc + 1

2000  format('    N o d a l   C o o r d i n a t e s'//6x,'node',
     &     3(i7,' coord'))
2001  format(i10,3f13.4)

      end

      subroutine pedges(iii,x,id,ndm,ndf,numnp,prt)

c     Purpose: Set boundary edge for restraints

c     Inputs:
c        iii        - Flag to test change in b.c. inputs
c        x(ndm,*)   - Nodal coordinates of mesh
c        ndm        - Spatial dimension of mesh
c        ndf        - Number dof/node
c        numnp      - Number of nodes in mesh
c        prt        - Flag, output results if true

c     Outputs:
c        id(ndf,*)  - Boundary condition codes

      implicit  none

      logical   prt
      character yyy*80
      integer   iii,ndm,ndf,numnp
      integer   i, j, n
      real*8    x0,dx

      integer   id(ndf,numnp),idl(6)
      real*8    x(ndm,numnp)

      integer         ior,iow
      common /iofile/ ior,iow

c     Read input of boundary edge for restraints

      iii = 1
100   if(ior.lt.0) write(*,5000)
      call pintio(yyy,10)
      read(yyy,1000,err=101) i,x0,idl
      go to 102
101   call perror('PEDGES',yyy)
      go to 100
102   if(i.le.0.or.i.gt.ndm) go to 4
      dx = 1.e-04
      do 103 n = 1,numnp
        dx = max( dx , x(i,n)*1.e-04 )
103   continue
      do 200 n = 1,numnp
        if(x(i,n).ne.-999. .and. abs(x(i,n)-x0).le.dx) then
          do 210 j = 1,ndf
            id(j,n) = max(abs(id(j,n)),abs(idl(j)))
210       continue
        endif
200   continue
      go to 100
4     if(prt) then
        call prthed(iow)
        write(iow,2000) (i,i=1,ndf)
        do 410 n = 1,numnp
          do 408 i = 1,ndf
            if(id(i,n).ne.0) then
              write(iow,2001) n,(id(j,n),j=1,ndf)
              go to 410
            endif
408       continue
410     continue
      endif

1000  format(i10,f10.0,6i10)
2000  format(/'    E d g e    N o d a l    B. C.'/
     &       /6x,'node',8(i3,' b.c.')/1x)
2001  format(i10,8i8)
5000  format(' Input: i-dir, x(i), b.code(j),j=1,ndf'/3x,'>',$)

      end

      subroutine sethis(ie,ix,idl,ipd,nen,nen1,numel)

c     Purpose: Set up history addresses in ix array

c     Inputs:
c        ie(9,*)    - Material set assembly data
c        ix(nen1,*) - Element nodal connection lists
c        idl(*)     - Working vector
c        ipd        - Precision of 'real*8' variables
c        nen        - Number nodes/element
c        nen1       - Dimension of ix array
c        numel      - Number of elements in mesh

c     Outputs:
c        ix(nen1,*) - Element nodal connection lists with history
c                     requirements added.

      implicit  none

      integer   maxa
      parameter(maxa=32000)

      integer   ipd,nen,nen1,numel
      integer   ihsiz,ihfac,ihmin, n,nu,nh0,nhinc

      integer   ie(9,*),ix(nen1,*),idl(*)

      integer        nhi,nhf,ihbuff,irec,jrec,nrec
      logical                                      hfl,hout
      common /hdatb/ nhi,nhf,ihbuff,irec,jrec,nrec,hfl,hout

      integer         maxm,ne
      common /psize/  maxm,ne

      character*12    tfile
      common /temfl1/ tfile(6)

      integer         itrec   ,nw1,nw2
      common /temfl2/ itrec(4),nw1,nw2

      real*8          dm
      common          dm(1000)

      nh0 = (ne+ipd-1)/ipd
      nhi = nh0

c     Record length factors

      ihfac  = 8
      ihsiz  = maxa
      ihbuff = ihsiz + 1 - nh0
      hfl    = .true.
      hout   = .false.
      irec   = 0
      nrec   = 1

c     Determine buffer size needed for history terms

      ihmin  = 0
      do 50  nu = 1,numel
        n     = idl(nu)
        ihmin = ihmin + ie(8,ix(nen1,n))
50    continue

c     Set buffer length and record length to minimum possible

      ihbuff = min(ihbuff,ihmin)
      nhf    = nhi + ihbuff - 1

c     Integer values

      ihbuff = ihbuff*ihfac

c     Set history area

      do 100 nu = 1,numel
      n     = idl(nu)
      nhinc = ie(8,ix(nen1,n))
      if(hfl .and. nhinc.ne.0) then
        itrec(2) = ihbuff
        open(3,file=tfile(2),access='direct', status='new',
     &         form='unformatted',recl=itrec(2))
        close(3)
        hfl = .false.
        call pconsd(dm(nhi),ihbuff/ihfac,0.0d0)
      endif

c     Integer values

      if(nh0+nhinc .gt. ihsiz) then
        call phstio(3,nrec,dm(nhi),ihsiz-nhi+1,2,tfile(2),itrec(2))
        nrec = nrec + 1
        nh0 = nhi
      endif
      ix(nen+1,n) = nh0
      ix(nen+2,n) = nrec
100   nh0 = nh0 + nhinc

c     Check for errors and finish initialization

      if(nrec.gt.numel) then
        write(*,2000)
        call pdefil(tfile,2,2)
        stop
      elseif(nh0.gt.nhi) then
        call phstio(3,nrec,dm(nhi),nhf-nhi+1,2,tfile(2),itrec(2))
      else
        nrec = nrec - 1
      endif
      ne   = (nhf+1)*ipd + 1
      close(3)

2000  format(' **ERROR** Insufficient storage for history terms')

      end
