      subroutine addvec(a,b,nn)

c     Purpose: Add two vectors

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

c     Outputs:
c        a(*)     - Sum of vector 1 and 2

      implicit none
      integer  n,nn
      real*8   a(*),b(*)

      do 100 n = 1,nn
        a(n) = a(n) + b(n)
100   continue

      end

      subroutine just(y,k,n0)

c     Purpose: Complete parser alignment of input data

c     Inputs:
c        y(*)        - Input string
c        k           - Length of input string
c        n0          - Field width for parsed data

c     Outputs:
c        y(*)        - Adjusted data to field widths of n0

      implicit    none
      integer     i,j,k,kl,l,n0,n1,n2
      character*1 y(k),yi,ze,ni,mi,pl,dt,sp

      data ze,ni,mi,pl,dt,sp/'0','9','-','+','.',' '/

      n1 = n0 - 1
      n2 = n1 - 1
      do 140 i = 1,k,n0
        do 100 j = i,i+n1
          if(y(j).ne.sp) go to 115
100     continue
        y(i+n1) = ze
115     if(y(i+n1).ne.sp) go to 140
        yi = y(i)
        if((yi.ge.ze.and.yi.le.ni).or.(yi.eq.mi).or.(yi.eq.pl)
     &                            .or.(yi.eq.dt)) then
          do 110 j = i+n2,i,-1
            if(y(j).ne.sp) go to 120
110       continue
120       kl = n1 + i - j
          do 130 l = j,i,-1
            y(l+kl) = y(l)
            y(l) = sp
130       continue
        endif
140   continue

      end

      subroutine modify(b,s,dul,nst)

c     Purpose: Modify for non-zero displacement boundary conditions

c     Inputs:
c        s(nst,*) - Element tangent matrix
c        dul(*)   - Element vector of specified displacements
c        nst      - Dimension of element vectors
c     Outputs:
c        b(*)     - Residual modified for specified displacements

      implicit none

      integer  nst,i,j
      real*8   b(*),s(nst,*),dul(*)

      do 110 i = 1,nst
        do 100 j = 1,nst
          b(i) = b(i) - s(i,j)*dul(j)
100     continue
110   continue

      end

      subroutine pangl(ix,nen,angl,angg,nrot)

c     Purpose: Set up table of rotation angles

c     Inputs:
c        ix(*)     - Element nodal connections
c        nen       - Number of nodes on element
c        angg(*)   - Angle for sloping b.c. of each node point

c     Outputs:
c        angl(*)   - Angle for each element node
c        nrot      - Number of nodes requiring rotation

      implicit  none

      integer   nen,nrot, ii,n
      integer   ix(nen)
      real*8    angl(nen),angg(*)

      nrot = 0
      do 100 n = 1,nen
        angl(n) = 0.0
        ii = abs(ix(n))
        if (ii.gt.0) then
          angl(n) = angg(ii)
          if (angg(ii).ne.0.0) nrot = nrot + 1
        endif
100   continue

      end

      subroutine pform(ul,xl,tl,ld,p,s,ie,d,id,x,ix,t,idl,u,b,
     &                 ndf,ndm,nen1,nst)

c     Purpose: Compute element arrays and assemble global arrays

c     Inputs:
c        ul(ndf,*) - Element nodal solution parameters
c        xl(ndm,*) - Element nodal coordinates
c        tl(*)     - Element nodal temperatures
c        ld(*)     - Element equation numbers for each dof
c        p(*)      - Element vector
c        s(nst,*)  - Element matrix
c        ie(9,*)   - Material set assembly data
c        d(*)      - Material set parameters
c        id(*)     - Equation numbers for each dof
c        x(ndm,*)  - Nodal coordinates for mesh
c        ix(nen1,*)- Element nodal connection lists
c        t(*)      - Nodal temperatures for mesh
c        idl(*)    - Working vector
c        u(*)      - Nodal displacements and increments
c        ndf       - Number dof/node
c        ndm       - Spatial dimension of mesh
c        nen1      - Dimension of ix array
c        nst       - Dimension of element arrays

c     Outputs:
c        b(*)      - Residual/reactions or diagonal mass

      implicit  none

      logical   efl
      integer   ndf,ndm,nen1,nst
      integer   i,ii,iid,ijj,isx,ixi, j,jj, k
      integer   ne2,ne3,ne4, nu,nt0,nt1,nrot,numnp2
      real*8    tm
      real*8    dun,un

      integer   ld(ndf,*),ie(9,*),id(ndf,*),ix(nen1,*),idl(*)
      real*8    xl(ndm,*),x(ndm,*),tl(*),t(*)
      real*8    d(18,*),p(*),s(nst,*),b(*),ul(ndf,*),u(ndf,*)

      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

      logical        fl    ,pfr
      common /fdata/ fl(11),pfr

      integer         ior,iow
      common /iofile/ ior,iow

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

      integer        nv,nw,nl
      common /ndata/ nv,nw,nl

      integer        nh1,nh2
      common /hdata/ nh1,nh2

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

      real*8         prop,ap
      integer                      iexp    ,ik    ,npld
      common /prlod/ prop,ap(6,10),iexp(10),ik(10),npld

      real*8         beta,gamm,theta
      integer                        nop,nt
      common /tbeta/ beta,gamm,theta,nop,nt

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

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

      integer        isw,nn1,nn2,nn3
      logical                        afl,bfl,cfl,dfl
      common /xdata/ isw,nn1,nn2,nn3,afl,bfl,cfl,dfl

      real*8 dm
      real*8       rm
      integer            im
      common dm(1),rm(1),im(1000)

c     Set up local arrays before calling element library

      iel = 0
      tm  = 1. - theta
      efl = .false.
      if(.not.dfl.and.isw.eq.6) efl = .true.
      if(bfl.and.isw.eq.3)      efl = .true.
      if(isw.ne.3.or.nn1.eq.1)  irec = 0
      ne2 = nen + nen
      ne3 = nen + ne2
      ne4 = nen + ne3
      nt0 = nt  - 1
      nt1 = nt0 + ndf*numnp
      numnp2 = numnp + numnp
      do 110 nu = 1,numel
       n = idl(nu)
       if( (n.ge.nn1 .and. n.le.nn2) .and. (mod(n-nn1,nn3).eq.0) ) then

c        Set up history terms

         ma  = ix(nen1,n)
         nh1 = ix(nen+1,n)
         nh2 = nh1

c        write(iow,*) 'HFL,HOUT:',hfl,hout
c        write(iow,*) 'IREC,JREC:',irec,ix(nen+2,n)
c        write(iow,*) 'NHI,NHF:',nhi,nhf,nhf-nhi+1,itrec(2)
c        write(*,*) 'HFL,HOUT:',hfl,hout

         if(.not.hfl) then
           jrec= ix(nen+2,n)
           if(jrec.ne.irec) then
             if(hout .and. irec.ne.0) then
               call phstio(3,irec,dm(nhi),nhf-nhi+1,2,tfile(2),itrec(2))
             endif
             call phstio(3,jrec,dm(nhi),nhf-nhi+1,1,tfile(2),itrec(2))
             irec = jrec
           endif
         endif
         call pconsd(ul,5*nen*ndf,0.0d0)
         call pconsd(xl,nen*ndm,0.0d0)
         call pconsd(tl,nen,0.0d0)
         call pconsd(dm(n11a),nen,0.0d0)
         call pconsi(ld,nst,0)
         un = 0.0
         dun= 0.0
         call pangl(ix(1,n),nen,dm(n11a),dm(n11b),nrot)
         do 108 i = 1,nen
           ixi= ix(i,n)
           ii = abs(ixi)
           if(ii.gt.0) then
             iid = ii*ndf - ndf
             nel = i
             tl(i) = t(ii)
             do 106 j = 1,ndm
               xl(j,i) = x(j,ii)
106          continue
             do 107 j = 1,ndf
               jj = ie(j,ma)
               if(jj.le.0) go to 107
               ijj     = iid + jj
               k       = id(jj,ii)
               ul(j,i) = u(jj,ii)
               ul(j,i+nen) = u(jj,ii+numnp)
               ul(j,i+ne2) = u(jj,ii+numnp2)
               if(fl(9)) ul(j,i+ne3) = dm(nv+ijj)
               if(k.le.0) then
                 ul(j,i+ne4) = theta*dm(nt1+ijj)+tm*dm(nt0+ijj)-u(jj,ii)
                 dun = max(dun,abs(ul(j,i+ne4)))
               endif
               un = max(un,abs(ul(j,i)))
               if(dfl) then
                 ld(j,i) = ijj
               else
                 ld(j,i) = sign(k,ixi)
               endif
107          continue
           endif
108      continue

c        Form element array

         dq  = prop
         if(ie(7,ma).ne.iel) mct = 0
         iel = ie(7,ma)
         isx = isw
         if(efl .and. dun.gt.0.0000001d0*un .and. .not.afl) isx = 3
         if(nrot.gt.0)
     &   call ptrans(ia(1,iel),dm(n11a),ul,p,s,nel,nen,ndf,nst,1)
         call elmlib(d(1,ma),ul,xl,ix(1,n),tl,s,p,ndf,ndm,nst,iel,isx)
         if(nrot.gt.0)
     &   call ptrans(ia(1,iel),dm(n11a),ul,p,s,nel,nen,ndf,nst,2)

c        Modify for non-zero displacement boundary conditions

         if(efl) call modify(p,s,ul(1,ne4+1),nst)

c        Assemble a vector if needed

         if(bfl) then
           do 109 i = 1,nst
             j = abs(ld(i,1))
             if(j.ne.0) b(j) = b(j) + p(i)
109        continue
         endif
       endif
110   continue

c     Put last history state on disk

      if(hout) call phstio(3,jrec,dm(nhi),nhf-nhi+1,2,tfile(2),itrec(2))

      end

      subroutine pload(id,fn,b,nn,xm,ac)

c     Purpose: Form load vector in compact form

c     Inputs:
c
c        id(*)     - Equation numbers for dof
c        fn(nn,2)  - Forces for active dofs at t_n and t_n+1
c        nn        - Size of id and fn arrays (total dof's)
c        xm(*)     - Diagonal mass array
c        ac(*)     - Acceleration vector

c     Outputs:
c        b(*)      - Residual vector

      implicit  none

      integer   j, n,nn
      integer   id(*)
      real*8    fn(nn,2)
      real*8    b(*),xm(*),ac(*)

      logical        fl    ,pfr
      common /fdata/ fl(11),pfr

      real*8         beta,gamm,theta
      integer                        nop,nt
      common /tbeta/ beta,gamm,theta,nop,nt

      fl(11) = .false.
      call pconsd(b,nn,0.0d0)
      do 100 n = 1,nn
        j = id(n)
        if(j.gt.0) then
          b(j) = theta*fn(n,2) + (1. - theta)*fn(n,1) + b(j)
          if(fl(9)) b(j) =  b(j) - xm(j)*ac(n)
        endif
100   continue

      end

      function propld(t,j)

c     Purpose: Proportional load table (j load records, maximum 10)

c     Inputs:
c        t      - Time
c        j      - Number of proportional loads to input

c     Outputs:
c        propld - Proportional loading for time 't' (when j = 0)

      implicit  none

      character yyy*80
      integer   i,j,l,m, nprop
      real*8    t,tmin,tmax
	  real*8    propld

      integer         ior,iow
      common /iofile/ ior,iow

      real*8         prop,a
      integer                     iexp    ,ik    ,npld
      common /prlod/ prop,a(6,10),iexp(10),ik(10),npld

      if(j.gt.0) then

c       Input table of proportional loads

        write(iow,2000)
        if(ior.lt.0) then
          write(  *,2000)
          write(  *,2003)
        endif
        do 100 i=1,j
101       call pintio(yyy,10)
          read(yyy,1000,err=102) ik(i),iexp(i),(a(m,i),m=1,6)
          go to 103

c         Error message

102       call perror('PROPLD',yyy)
          go to 101

c         Set a default ramp table if a type "0" input

103       if(ik(i).eq.0) then
            a(2,i) = 1.e+6
            a(4,i) = 1.
          endif
          write(iow,2001) i,ik(i),(a(m,i),m=1,6),iexp(i)
          if(ior.lt.0) write(*,2001) i,ik(i),(a(m,i),m=1,6),iexp(i)
100     continue
        nprop = j
      endif

c     Compute value at time t

      propld = 0.0
      do 220 i = 1,nprop
        tmin = a(1,i)
        tmax = a(2,i)
        if(t.lt.tmin.or.t.gt.tmax) go to 220
        l = max(iexp(i),1)
        propld = a(3,i)+a(4,i)*t+ a(5,i)*(sin(a(6,i)*t+tmin))**l
     &         + propld
220   continue

c     Formats

1000  format(2i10,6f10.0)
2000  format(30x,'P r o p o r t i o n a l   L o a d   T a b l e')
2001  format(/,' number    type      tmin',10x,'tmax',/i3,i10,7x,g10.4,
     & 4x,g10.4,/6x,'a(1)',10x,'a(2)',10x,'a(3)',10x,'a(4)',10x,
     & 'exp',/4x,4(g10.4,4x),i5/)
2003  format(' Input: type, exponent, tmin, tmax, a(i),i=1,4'/' >',$)

      end

      subroutine prtdis(x,b,ttim,prop,ndm,ndf,n1,n2,n3)

c     Purpose: Output nodal displacement values

c     Inputs:
c        x(ndm,*)  - Nodal coordinates for mesh
c        b(ndf,*)  - solution vector for nodes
c        ttim      - Current time in solution
c        prop      - Current proportional load value
c        ndm       - Spatial dimension of mesh
c        ndf       - Number dof/node
c        n1        - First node to output
c        n2        - Last  node to output
c        n3        - Increment to n1 for next node to output, etc.

c     Outputs:
c        none

      implicit  none

      character cd*6,di*6
      integer   ndm,ndf,n1,n2,n3
      integer   i, n, kount
      real*8    ttim,prop

      real*8    x(ndm,*),b(ndf,*)

      integer         ior,iow
      common /iofile/ ior,iow

      data cd/' coord'/,di/' displ'/

      kount = 0
      do 100 n = n1,n2,n3
        kount = kount - 1
        if(kount.le.0) then
        call prthed(iow)
        write(iow,2000) ttim,prop,(i,cd,i=1,ndm),(i,di,i=1,ndf)
        if(ior.lt.0) then
          write(*,2000) ttim,prop,(i,cd,i=1,ndm),(i,di,i=1,ndf)
        endif
        kount = 48
        endif
        if(x(1,n).ne. -999.) then
          write(iow,2001) n,(x(i,n),i=1,ndm),(b(i,n),i=1,ndf)
          if(ior.lt.0) then
            write(*,2001) n,(x(i,n),i=1,ndm),(b(i,n),i=1,ndf)
          endif
        else
          write(iow,2002) n
          if(ior.lt.0) write(*,2002) n
        endif
100   continue

2000  format('  N o d a l   D i s p l a c e m e n t s',5x,
     &  'time',e18.5/31x,'prop. ld. (eigenvalue)',e13.5//
     &  '  node',9(i7,a6))
2001  format(i6,1p9e13.6)
2002  format(i6,' not input.')

      end

      subroutine prterr()

c     Purpose:  Output error estimation values.

      implicit  none

      real*8    elrind,ecrind

      integer         ior,iow
      common /iofile/ ior,iow

      real*8          eerror,elproj,ecproj,efem,enerr,ebar
      common /errind/ eerror,elproj,ecproj,efem,enerr,ebar

c     Output error indicator values

      elrind = 0.0d0
      ecrind = 0.0d0
      if(elproj.ne.0.0d0) elrind = sqrt((efem-elproj)/elproj)
      if(ecproj.ne.0.0d0) ecrind = sqrt((efem-ecproj)/ecproj)
      if(ecproj.ne.0.0d0) eerror = sqrt(eerror/ecproj)
      write(iow,2000) efem,elproj,ecproj,elrind,ecrind,eerror
      if(ior.lt.0) then
        write(*,2000) efem,elproj,ecproj,elrind,ecrind,eerror
      endif

2000  format(/'   Finite Element Stress Measure              =',e15.8/
     &        '   L u m p e d  Projected Stress Measure      =',e15.8/
     &        '   Consistent   Projected Stress Measure      =',e15.8/
     &        '   L u m p e d  Error Indicator               =',e15.8/
     &        '   Consistent   Error Indicator               =',e15.8/
     &        '   Direct Error Indicator                     =',e15.8/)

      end

      subroutine prtrea(r,ndf,numnp,n1,n2,n3)

c     Purpose: Print nodal reactions

c     Inputs:
c        r(ndf,*)  - Nodal reactions
c        ndf       - Number dof/node
c        numnp     - Number of nodes in mesh
c        n1        - First node to output
c        n2        - Last  node to output
c        n3        - Increment to n1 for next node to output, etc.

c     Outputs:
c        none

      implicit  none

      integer   ndf,numnp, n1,n2,n3
      integer   i, k, n, kount

      real*8    rr(6),r(ndf,1),rsum(6),asum(6),psum(6)

      integer         ior,iow
      common /iofile/ ior,iow

      call pconsd(rsum,ndf,0.0d0)
      call pconsd(psum,ndf,0.0d0)
      call pconsd(asum,ndf,0.0d0)

      do 85 i = 1,numnp
        do 75 k = 1,ndf
          rsum(k) = rsum(k) - r(k,i)
          asum(k) = asum(k) + abs(r(k,i))
75      continue
85    continue

      kount = 0
      do 100 n = n1,n2,n3
        kount = kount - 1
        if(kount.le.0) then
          call prthed(iow)
          write(iow,2000) (k,k=1,ndf)
          if(ior.lt.0) write(*,2000) (k,k=1,ndf)
          kount = 50
        endif

        do 80 k = 1,ndf
          rr(k) = -r(k,n)
          psum(k) = psum(k) + rr(k)
80      continue

        write(iow,2001) n,(rr(k),k=1,ndf)
        if(ior.lt.0) write(*,2001) n,(rr(k),k=1,ndf)
100   continue

c     Print statics check

      write(iow,2002) (rsum(k),k=1,ndf)
      write(iow,2003) (psum(k),k=1,ndf)
      write(iow,2004) (asum(k),k=1,ndf)
      if(ior.lt.0) then
        write(*,2002) (rsum(k),k=1,ndf)
        write(*,2003) (psum(k),k=1,ndf)
        write(*,2004) (asum(k),k=1,ndf)
      endif

2000  format('  N o d a l    R e a c t i o n s'//6x,'node',6(i9,' dof'))
2001  format(i10,6e13.4)
2002  format(/7x,'sum',6e13.4)
2003  format( 3x,'prt sum',6e13.4)
2004  format( 3x,'abs sum',6e13.4)

      end

      subroutine prtstr(dt,ds,numnp,n1,n2,n3)

c     Purpose: Output projected nodal stress values
c
c     Inputs:
c        dt(*)       - Value of 8-th nodal value
c        st(numnp,*) - Value of 1-7  nodal values
c        numnp       - Number of nodes in mesh
c        n1          - First node to output
c        n2          - Last  node to output
c        n3          - Increment to n1 for next node to output, etc.

      implicit  none

      integer   i, kount, n, numnp, n1,n2,n3

      real*8    dt(numnp),ds(numnp,7)

      integer         ior,iow
      common /iofile/ ior,iow

      kount = 0
      do 200 n = n1,n2,n3
        kount = kount - 1
        if(kount.le.0) then
          call prthed(iow)
          write(iow,2000)
          if(ior.lt.0) write(*,2000)
          kount = 17
        endif
        write(iow,2001) n,(ds(n,i),i=1,7),dt(n)
        if(ior.lt.0) write(*,2001) n,(ds(n,i),i=1,7),dt(n)
200   continue

2000  format('   N o d al   S t r e s s e s'//' node',4x,'11-stress',4x,
     & '12-stress',4x,'22-stress',4x,'33-stress'/10x,'1-stress',5x,
     & '2-stress',4x,'max-shear',8x,'angle')
2001  format(i5,4e13.5/5x,4e13.5/1x)

      end

      subroutine ptrans(ia,angl,ul,p,s,nel,nen,ndf,nst,isw)

c     Purpose: Subroutine to make two-dimesional rotations

c     Inputs:
c        ia(2)       - Dof's to rotate
c        angl(*)     - Angle of sloping boundary condition at element nodes
c        nel         - Number of nodes on element
c        nen         - Maximum number of nodes on element
c        ndf         - Number dof/node
c        nst         - Dimension of element arrays
c        ist         - Switch to control actions taken

c     Outputs:
c        ul(ndf,*,*) - Rotated element solution vector/rates (isw = 1)
c        p(ndf,*)    - Rotated element vector (isw = 2)
c        s(nst,*)    - Rotated element matrix (isw = 2)

      implicit  none

      integer   nel,nen,ndf,nst,isw
      integer   i, j, i1, ij1,ij2
      real*8    ang,cs,sn,tm

      integer   ia(2)
      real*8    angl(*)
      real*8    ul(ndf,nen,*),p(ndf,*),s(nst,nst)

c     Recover dof to be rotated

      ij1 = ia(1)
      ij2 = ia(2)
      i1  = 0

c     Transform displacement quantities to element coordinates

      if(isw.eq.1) then
        do 110 i = 1,nel
          if(angl(i).ne.0.0) then
            ang = angl(i)*0.01745329252d0
            cs  = cos(ang)
            sn  = sin(ang)
            do 100 j = 1,4
              tm          = cs*ul(ij1,i,j) - sn*ul(ij2,i,j)
              ul(ij2,i,j) = sn*ul(ij1,i,j) + cs*ul(ij2,i,j)
              ul(ij1,i,j) = tm
100         continue
          endif
110     i1 = i1 + ndf

c     Transform element arrays to global coordinates

      else
        do 220 i = 1,nel
          if(angl(i).ne.0.0) then
            ang = angl(i)*0.01745329252d0
            cs  = cos(ang)
            sn  = sin(ang)

c           Transform load vector

            tm       = cs*p(ij1,i) + sn*p(ij2,i)
            p(ij2,i) =-sn*p(ij1,i) + cs*p(ij2,i)
            p(ij1,i) = tm

c           Postmultiply s by transformation

            do 210 j = 1,nst
              tm         = s(j,i1+ij1)*cs + s(j,i1+ij2)*sn
              s(j,i1+ij2)=-s(j,i1+ij1)*sn + s(j,i1+ij2)*cs
              s(j,i1+ij1)= tm
210         continue

c           Premultiply s by transformation

            do 215 j = 1,nst
              tm         = cs*s(i1+ij1,j) + sn*s(i1+ij2,j)
              s(i1+ij2,j)=-sn*s(i1+ij1,j) + cs*s(i1+ij2,j)
              s(i1+ij1,j)= tm
215         continue
          endif
220     i1 = i1 + ndf
      endif

      end

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

c     Purpose: Real*4 vector times scalar added to vector

c     Inputs:
c        a(*)     - Vector to be multiplied by scalar
c        b(*)     - Vector to add
c        x        - Scalar multiplier for a-vector
c        n        - Length of vectors

c     Outputs:
c        c(*)     - Sum of 'b' and 'x*a'

      implicit none
      integer  n, k
      real*8   a(*),b(*),c(*)
      real*8   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 param(ctl,ct)

c     Purpose: Set appropriate time integration parameters

c     Inputs:
c        ctl       - Character string of integration type
c        ct(3)     - Integration type parameters

c     Outputs:
c        none      - Output through common block

      implicit  none

      logical   pcomp
      character ctl*4
      real*8    ct(3)

      integer         ior,iow
      common /iofile/ ior,iow

      real*8         beta,gamm,theta
      integer                        nop,nt
      common /tbeta/ beta,gamm,theta,nop,nt

c     Set integration parameters

      beta = ct(1)
      gamm = ct(2)

c     SS11 algorithm

      if(pcomp(ctl,'ss11')) then
        nop = 1
        if(beta.le.0.0) beta = 0.50
        write(iow,2001) beta
        if(ior.lt.0) write(*,2001) beta
        theta = beta

c     SS22 algorithm

      elseif(pcomp(ctl,'ss22')) then
        nop = 2
        if(beta.le.0.0) beta = 0.50
        if(gamm.le.0.0) gamm = 0.50
        write(iow,2002) beta,gamm
        theta = gamm
        if(ior.lt.0) write(*,2002) beta,gamm

c     Newmark algorithm

      else
        nop = 3
        if(beta.le.0.0) beta = 0.25
        if(gamm.le.0.0) gamm = 0.50
        write(iow,2003) beta,gamm
        if(ior.lt.0) write(*,2003) beta,gamm
        theta = beta
      endif

2001  format(' SS-11 Method Parameter: theta = ',f9.4)
2002  format(' SS-22 Method Parameters'/
     &       ' theta-1 = ',f9.4,' ;  theta-2 = ',f9.4)
2003  format(' Newmark Method Parameters'/
     &       ' beta = ',f9.4,' ;  gamma = ',f9.4)

      end

      subroutine setci(ior)

c     Purpose: Compute integration constants 'c1' to 'c5' for current 'dt'

c     Inputs:
c        ior      - Data input logical unit number

c     Outputs:
c        none     - Output through common block

      implicit none

	  integer  ior

      real*8         beta,gamm,theta
      integer                        nop,nt
      common /tbeta/ beta,gamm,theta,nop,nt

      real*8         ttim,dt,c1,c2,c3,c4,c5,c6
      common /tdata/ ttim,dt,c1,c2,c3,c4,c5,c6

      if(dt.le.0.0) then
        write(*,2000)
        if(ior.gt.0) stop
        return
      endif

c     Compute integration constants 'c1' to 'c5' for current 'dt'

      if(nop.eq.1) then
        c1 = 1.d0/(beta*dt)
        c2 = c1
        c3 = 1.d0/beta
        c6 = beta
      elseif(nop.eq.2) then
        c5 = 2.d0/(gamm*dt)
        c1 = c5/dt
        c2 = c5*beta
        c3 = dt*beta
        c4 = 1.d0/gamm
        c6 = gamm
      else
        c1 = 1.d0/(beta*dt*dt)
        c2 = gamm/(dt*beta)
        c3 = 1.d0 - 1.d0/(beta+beta)
        c4 = 1.d0 - gamm/beta
        c5 = (1.d0 - gamm/(beta+beta))*dt
        c6 = dt*c1
      endif

2000  format(' **ERROR** Input DT as nonzero number.')

      end

      subroutine update(id,fn,u,v,du,nneq,fdyn,isw)

c     Purpose: Update displacements (and velocities and accelerations)

c     Inputs:
c       id(*)      - Equation numbers for each dof
c       fn(*)      - Forced solution values at restrained dof
c       u(*)       - Solution vector and incrmeents: isw = 1
c       du(*)      - Solution increment: isw = 2
c       nneq       - Number of total parameters in u and v
c       fdyn       - Flag, transient solution if true
c       isw        - Switch

c     Outputs:
c       u(*)       - Initialization for   t_n+1: isw = 1
c                    Solution iterates at t_n+1: isw = 2
c
c       v(*)       - Initialized rates for t_n+1: isw = 1
c                    Rate iterates at      t_n+1: isw = 2

      implicit  none

      logical   fdyn
	  integer   j,n,nneq,isw
	  real*8    ur1,ur2,dot

      integer   id(*)
      real*8    fn(nneq,2)
      real*8    u(nneq,*),v(nneq,*),du(*)

      integer         ior,iow
      common /iofile/ ior,iow

      real*8         beta,gamm,theta
      integer                        nop,nt
      common /tbeta/ beta,gamm,theta,nop,nt

      real*8         ttim,dt,c1,c2,c3,c4,c5,c6
      common /tdata/ ttim,dt,c1,c2,c3,c4,c5,c6

c     update solution vectors to begin a step

      if(isw.eq.1) then

        ur1 = sqrt(dot(v(1,1),v(1,1),nneq))
        ur2 = sqrt(dot(v(1,2),v(1,2),nneq))
        write(iow,2000) ur1,ur2
        if(ior.lt.0) write(*,2000) ur1,ur2

        do 100 n = 1,nneq

c         SS11 algorithm u=u-bar; v1=v-bar; v2=a=abar; v3=u;

          if(nop.eq.1) then
            u(n,1) = v(n,3)
            v(n,1) = 0.0d0
            v(n,2) = 0.0d0

c         SS22 algorithm u=u-bar; v1=v-bar; v2=a=abar; v3=u; v4=v

          elseif(nop.eq.2) then
            u(n,1) = v(n,3) + c3*v(n,4)
            v(n,1) = v(n,4)
            v(n,2) = 0.0d0
            v(n,3) = v(n,3) + dt*v(n,4)

c         Newmark algorithm

          else
            ur2    = - c6*v(n,1) + c3*v(n,2)
            v(n,1) =   c4*v(n,1) + c5*v(n,2)
            v(n,2) =   ur2
          endif
100     continue
      else

c       Update displacement and its increments within time step

        ur2 = 1. - theta
        do 200 n = 1,nneq
          j = id(n)
          if (j.gt.0) then

c           For active degrees-of-freedom compute values from solution

            u(n,1) = du(j) + u(n,1)
            u(n,2) = du(j) + u(n,2)
            u(n,3) = du(j)
          else

c           For fixed degrees-of-freedom compute values from forced inputs

            ur1    = theta*fn(n,2) + ur2*fn(n,1)
            u(n,3) = ur1    - u(n,1)
            u(n,2) = u(n,2) + u(n,3)
            u(n,1) = ur1
          endif
200     continue

c       For time dependent solutions update rate terms

        if(fdyn) then

          do 210 n = 1,nneq
            v(n,1) = v(n,1) + c2*u(n,3)
            v(n,2) = v(n,2) + c1*u(n,3)

c           SS11 algorithm

            if(nop.eq.1) then
              v(n,3) = v(n,3) + c3*u(n,3)

c           SS22 algorithm

            elseif(nop.eq.2) then
              v(n,3) = v(n,3) + c4*u(n,3)
              v(n,4) = v(n,4) + c5*u(n,3)
            endif
210       continue
        endif
      endif

2000  format('   N o r m s   f o r   D y n a m i c s'/
     &   10x,'Velocity:',e13.5,' Acceleration:',e13.5)

      end
