      subroutine pmacr (ld,s,ie,d,id,x,ix,f,t,jd,f0,b,
     &                  dr,ndf,ndm,nen1,nst,prt)

c     Purpose: Command language instruction driver

c     Inputs:
c        ld(*)     - Element equation number assembly vector
c        s(*)      - Element vector/matrix storage
c        ie(*)     - Material set assembly data
c        d(*)      - Material set parameters
c        id(*)     - List of equation numbers for each dof
c        x(*)      - Nodal coordinates for mesh
c        ix(*)     - Element nodal connection lists
c        f(*)      - Nodal specified force/displacmeents
c        t(*)      - Nodal temperatures
c        jd(*)     - Pointer to end of row/columns for profiles
c        f0(*)     - Base solution
c        b(*)      - Nodal solution array
c        dr(*)     - Residual/reactions
c        ndf       - Number dof/node
c        ndm       - Spatial dimension of FEM mesh
c        nen1      - DImension of ix array
c        nst       - Dimension of element arrays
c        prt       - Flag, output results if true

c     Outputs:
c        none

      implicit  none

      logical   prt,lxst
      integer   ndf,ndm,nen1,nst
      integer   i,iz, j, k, ll, nm1,nm2,nlp,nneq

      character wd(24)*4,tary*10,lwd*3
      integer   ld(*),ie(*),id(*),ix(*),jd(*)
      real*8    ct(3,100),x(*),f(*),f0(*),t(*)
      real*8    d(*),s(*),b(*),dr(*)

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

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

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

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

      integer         ior,iow
      common /iofile/ ior,iow

      integer        l,lv,lvs   ,lve   ,jct
      common /ldata/ l,lv,lvs(9),lve(9),jct(100)

      character*4    lct
      common /ldatb/ lct(100)

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

      real*8         engy,rnmax,tol,shift
      common /rdata/ engy,rnmax,tol,shift

      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

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

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

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

      data wd/'stre','utan','tang','form','mass','reac','chec','disp',
     1        'solv','mesh','rest',
     2        'tol ','dt  ','loop','next','prop','data','time','beta',
     2        'newf','subs','eigv','iden',
     p        'plot'/

c     nmi = no. macro commands in 'pmacri'; nlp = loop number

      data nm1,nm2/11,12/,iz/0/
      nlp = nm1 + 3

c     set initial values of parameters

      call pinitc(engy,rnmax,shift,tol,dt,prop,ttim,npld)
      nw1  = nm1
      nw2  = nm2 + nw1
      nneq = ndf*numnp
      call psetm (nt,nneq*2,'d',fl(1))
      do 50 j = 0,nneq*2-1
        dm(nt+j) = 0.0d0
50    continue

c     input the macro commands

100   call pmacio (jct,lct,ct,wd,nw2+1,nlp,ll)
      if(ll.le.0) go to 300

c     execute macro instruction program

      lv = 0
      l = 1
200   j = jct(l)
      i = l - 1
      call pctime (tary)
      if(l.ne.1.and.l.ne.ll) then
       write(iow,2001) i,wd(j),lct(l),(ct(k,l),k=1,3),tary
       if(ior.lt.0) write(*,2001) i,wd(j),lct(l),(ct(k,l),k=1,3),tary
      endif
      if(j.le.nw1) call pmacr1(id,ie,ix,ld,d,s,x,f,f0,t,jd,b,dr,
     1                         lct,ct,ndf,ndm,nen1,nst,nneq,prt,j)
      if(j.ge.nw1+1.and.j.le.nw2)
     1             call pmacr2(id,x,f,f0,jd,b,dr,lct,ct,ndf,nneq,j-nw1)

c     plot macro call

      if(j.eq.nw2+1) then
        call pplotf(x,ix,b,lct(l),ct(1,l),ndf,ndm,nen1)
      endif
      l = l + 1
      if(l.le.ll) go to 200
      if (ior.lt.0) go to 100
300   call pctime(tary)
      write(iow,2000) tary
      if(ior.lt.0) write(*,2000) tary
      if(.not.fl(4)) close(4,status='delete')

c     save restart information

      if(ll.lt.-1.or.fl(7)) close(3,status='delete')
      if(ll.lt.-1.or.fl(7)) return
      inquire(file=tfile(6),exist=lxst)
      lwd = 'new'
      if(lxst) lwd = 'old'
      open (7,file=tfile(6),form='unformatted',status=lwd)
      rewind 7
      write(7) numnp,numel,nummat,ndm,ndf,nhi,nhf,nrec
      write(7) ttim,(b(i),i=1,3*nneq)
      write(7) (dm(i),i=nt,nt+nneq)
      if(fl(9)) write(7) (dm(i),i=nv,nv+4*neq)
      if(nrec.gt.0) then
c        do 400 j = 1,nrec
c          call phstio(3,j,dm(nhi),nhf-nhi+1,1,tfile(2),itrec(2))
c          call phstio(7,j,dm(nhi),nhf-nhi+1,22,tfile(6),iz)
c400     continue
c       call pdefil(tfile,2,2)
      endif
      close(7)

c     formats

2000  format(' *End of macro execution*',15x,'time=',a10)
2001  format(' *Macro ',i3,' *',2(a4,1x),
     1   'V1=',g10.3,' V2=',g10.3,' V3=',g10.3/40x,'time=',a10)
      end

      subroutine pinitc(engy,rnmx,shift,tol,dt,prop,ttim,npld)

c     Purpose: Set initial values of solution parameters

c     Outputs:
c        engy       - Solution step energy         ( 0.0d0)
c        rnmx       - First iterate residual norm  ( 0.0d0)
c        whift      - Shift to eigenproblem        ( 0.0d0)
c        tol        - Solution tolerance           (1.d-12)
c        dt         - Solution time step size      ( 0.0d0)
c        prop       - Proportional loading         ( 1.0d0)
c        ttim       - Solution time                ( 0.0d0)
c        npld       - Number of proportional loads (  0   )

      implicit   none
      integer    npld, i
      real*8     dt,ttim
      real*8     engy,rnmx,shift,tol,prop

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

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

c     set initial values of parameters

      nop   = 3
      theta = 1.0
      npld  = 0
      engy  = 0.0
      rnmx  = 0.0
      shift = 0.0
      tol   = 1.e-12
      dt    = 0.0
      prop  = 1.0
      ttim  = 0.0
      do 50 i = 1,7
      fl(i+4) = .false.
50    fl(i)   = .true.
      end

      subroutine pmacio (jct,lct,ct,wd,nwd,nlp,ll)

c     Purpose: Command language instruction input subprogram

c     Inputs:
c        wd(*)      - Command names
c        nwd        - Number of commands
c        nlp        - Number of "loop" command in list

c     Outputs
c        jct(*)     - Number of command to execute
c        lct(*)     - Command option name
c        ct(3,*)    - Command parameters
c        ll         - Number of commands to execute

      implicit   none

      logical    pcomp
      integer    nwd,nlp,ll
      integer    i, j, l

      integer    jct(*)
      real*8     ct(3,*)
      character  wd(nwd)*4,lct(*)*4,clab1*4,clab2*4,tary*10,yyy*80

      integer         ior,iow
      common /iofile/ ior,iow

c     initiate the read of macro statements

      if(ior.gt.0) call prthed(iow)
      if(ior.gt.0) write(iow,2001) 

c     read macro cards

      ll = 1
      jct(1) = nlp
      ct(1,1) = 1.0
100   if(ior.lt.0) then
       call pctime(tary)
       write(*,2002) tary,ll
      endif
      ll = ll + 1
      call pintio(yyy,15)
      read(yyy,1000,err=401) clab1,clab2,(ct(i,ll),i=1,3)
      if(ior.lt.0.and.pcomp(clab1,'help')) then
        call phelp(wd,nwd,'MACRO',1)
        ll = ll - 1
        go to 100
      endif
      if(ior.gt.0.and.pcomp(clab1,'end ')) go to 150
      if(ior.lt.0) then
        if(pcomp(clab1,'exit')) ll = -1
        if(pcomp(clab1,'q   ').or.pcomp(clab1,'quit')) ll = -2
        if(ll.lt.0) return
      endif

c     set execution flag

      lct(ll) = clab2
      do 110 j = 1,nwd
110   if(pcomp(clab1,wd(j))) go to 130
      call perror('PMACIO',yyy)
      ll = ll - 1
      go to 100
130   jct(ll) = j
      if(ior.gt.0) then
        write(iow,2000) clab1,clab2,(ct(l,ll),l=1,3)
        go to 100
      endif
      ll = ll + 1
150   jct(ll)= nlp+1

c     set loop markers

      j = 0
      do 200 l = 2,ll-1
        if(jct(l).eq.nlp)   j = j + 1
        if(j.gt.8) go to 400
        if(jct(l).eq.nlp+1) j = j - 1
        if(j.lt.0) go to 400
200   continue
      if(j.ne.0.and.ior.gt.0) go to 400
      if(j.ne.0.and.ior.lt.0) ll = ll - 1
      if(j.ne.0) go to 100
      do 230 l = 1,ll-1
        if(jct(l).ne.nlp) go to 230
        j = 1
        do 210 i = l+1,ll
          if(jct(i).eq.nlp) j = j + 1
          if(jct(i).eq.nlp+1) j = j - 1
          if(j.eq.0) go to 220
210     continue
        go to 400
220     ct(2,i) = l
        ct(2,l) = i
230   continue
      return

c     error messages

400   write(iow,3000)
      if(ior.gt.0) stop
      if(ior.lt.0) write(*,3000)
      go to 100
401   call perror('PMACIO',yyy)
      go to 100
1000  format(2(a4,11x),3f15.0)
2000  format(7x,a4,1x,a4,1x,3g12.5)
2001  format('  M a c r o   I n s t r u c t i o n s'//
     1'  Macro Statement  Variable 1  Variable 2  Variable 3')
2002  format(' Input MACRO: "exit" = stop with restart',
     1   '; "quit" = quit.'/3x,'Time = ',a10,' Macro',i3,'> ',$)
3000  format(' **ERROR**  Wrong loop/next order, or > 8 loops')
      end

      subroutine pmacr1(id,ie,ix,ld,d,s,x,f,f0,t,jd,b,dr,
     &                  lct,ct,ndf,ndm,nen1,nst,nneq,prt,j)

c     Purpose: Command language instruction executions

c     Inputs:
c        id(*)     - List of equation numbers for each dof
c        ie(*)     - Material set assembly data
c        ix(*)     - Element nodal connection lists
c        ld(*)     - Element equation number assembly vector
c        d(*)      - Material set parameters
c        s(*)      - Element vector/matrix storage
c        x(*)      - Nodal coordinates for mesh
c        f(*)      - Nodal specified force/displacmeents
c        f0(*)     - Base solution
c        t(*)      - Nodal temperatures
c        jd(*)     - Pointer to end of row/columns for profiles
c        b(*)      - Nodal solution array
c        dr(*)     - Residual/reactions
c        lct       - Command option to execute
c        ct(3,*)   - Command parameters
c        ndf       - Number dof/node
c        ndm       - Spatial dimension of FEM mesh
c        nen1      - DImension of ix array
c        nst       - Dimension of element arrays
c        nneq      - Size of id, f, f0, and b arrays
c        prt       - Flag, output results if true
c        j         - Command number to execute

c     Outputs:
c        various   - Depends on command number

      implicit  none

      integer   maxa
      parameter(maxa=32000)

      logical   fa,tr,pcomp,prt,sflg,tflg
      integer   ndf,ndm,nen1,nst,nneq,i,j
      integer   ietyp,ml1,na,nal,n1,n2,n3,n4
      integer   nnpo,nnlo,nnmo,ndmo,ndfo,nhio,nhfo,nrco
      real*8    zero,shft,rnorm,enold

      character lct(*)*4,tt*10
      integer   id(*),ie(*),ix(*),jd(*),ld(*),iz
      real*8    ct(3,*),x(*),f(*),f0(*),t(*)
      real*8    b(*),d(*),dr(*),s(*),dot

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

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

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

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

      integer         maxf
      common /frdata/ maxf

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

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

      integer         ior,iow
      common /iofile/ ior,iow

      integer        l,lv,lvs   ,lve   ,jct
      common /ldata/ l,lv,lvs(9),lve(9),jct(100)

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

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

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

      real*8         engy,rnmax,tol,shift
      common /rdata/ engy,rnmax,tol,shift

      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

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

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

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

      data zero/0.0/,fa,tr/.false.,.true./,iz/0/

c     transfer to correct process

      n1 = 1
      n3 = 1
      go to (1,2,2,3,4,5,6,5,2,9,10), j

c     print stress values

1     n1 = ct(1,l)
      n2 = ct(2,l)
      n3 = ct(3,l)
      n3 = max(n3,1)
      n4 = numnp + 1
      if (pcomp(lct(l),'node')) then
        n1 = max(1,min(numnp,n1))
        n2 = max(n1,min(numnp,n2))
        enerr = 0.0
        if(.not.fl(11)) then
          call pconsd(aa,8*numnp,0.0d0)
          call formfe(b,dr,fa,fa,fa,fa,8,1,numel,1)
          call pltstr(aa,aa(n4),numnp)
        endif
        call prtstr(aa,aa(n4),numnp,n1,n2,n3)
        fl(11) = tr
      elseif (pcomp(lct(l),'erro')) then
        n1 = max(n1,1)
        n2 = 8*numnp
        do 105 i = 1,n2
          aa(n2+i) = 0.0
105     continue
        enerr = 0.0
        do 110 i = 1,n1
          call pconsd(aa,n2,0.0d0)
          call formfe(b,dr,fa,fa,fa,fa,8,1,numel,1)
          call pltstr(aa,aa(n4),numnp)
          call addvec(aa(n2),aa(n4),n2-numnp)
110     continue
        fl(11) = tr
        eerror = 0.0
        efem   = 0.0
        ebar   = 0.05*sqrt(enerr/numel)
        ietyp  = 1
        call formfe(b,dr,fa,fa,fa,fa,7,1,numel,1)
        call prterr
      else
        if(pcomp(lct(l),'all ')) then
          n2 = numel
        else
          n1 = max(1,min(numel,n1))
          n2 = max(n1,min(numel,n2))
        endif
        call formfe(b,dr,fa,fa,fa,fa,4,n1,n2,n3)
      endif
      return

c     solution step

2     shft = c1
      sflg = fl(9)
      if(j.eq.9) then
        if(.not.fl(8)) return
        fl(7) = .false.
        tflg  = .false.

c     form tangent stiffness

      else

c      unsymmetric tangent (must use profile solver)

        if(j.eq.2) fl(6) = .true.

c      symmetric tangent

        if(j.eq.3) fl(6) = .false.
        if(ct(1,l).ne.zero) then
          fl(8) = tr
          fl(7) = fa
          call raxpb(f, f0, prop, nneq, dm(nt+nneq))
          call pload(id,dm(nt),dr,nneq,dm(nl),dm(nw))
        endif
        shift= 0.
        tflg = tr
        if(.not.fl(9).and.ct(2,l).ne.zero) then
          if(fl(2)) then
            sflg = tr
            shift= ct(2,l)
            shft = -shift
            if(ior.lt.0) write(*,2006) shift
            write(iow,2006) shift
          else
            if(ior.lt.0) write(*,2007)
            write(iow,2007)
            if(ior.gt.0) stop
            return
          endif
        endif
      endif

c     call the solve routine to assemble and solve the tangent matrix

      na = maxf + 1
      nal= (maxf*(maxf+1))/2 + na
      call psolve(b,aa(na),aa,dr,aa(nal),dm(nl),s,ld,jd,im(n11c),nst,1,
     1            tflg,fl(8),sflg,shft,ipd,rnorm,engy,4)
      call pctime(tt)
      if(fl(8)) then
        fl(8) = fa
        if(tflg) write(iow,2001) rnorm,tt
        if(tflg .and. ior.lt.0) write(*,2001) rnorm,tt
        if (rnmax.eq.0.0d0) then
          rnmax = abs(engy)
          if(ct(3,l).le.0.0) ct(3,l) = 0.6
          enold = rnmax*0.9999
        endif
        write(iow,2004) rnmax,engy,tol
        if(ior.lt.0) write(*,2004) rnmax,engy,tol
        if(abs(engy).le.tol*rnmax) then
          ct(1,lve(lv)) = ct(1,lvs(lv))
          l = lve(lv) - 1
        elseif(pcomp(lct(l),'line')) then

c     line search

          if(abs(engy).gt.ct(3,l)*enold) then
           ml1 = 1 + nneq
           call serchl(dm(nt),id,engy,aa,b,dr,ct(3,l),aa(ml1),nneq)
          endif
        endif
        call update(id,dm(nt),b,dm(nv),dr,nneq,fl(9),2)
      else
        write(iow,2002) tt
        if(ior.lt.0) write(*,2002) tt
      endif
      return

c     form out of balance force for time step/iteration

3     if(fl(8)) return
      call raxpb(f, f0, prop, nneq, dm(nt+nneq))
      call pload(id,dm(nt),dr,nneq,dm(nl),dm(nw))
      call formfe(b,dr,fa,tr,fa,fa,6,1,numel,1)
      if(pcomp(lct(l),'acce').and.ttim.eq.0.0 .and. nop.eq.3) then
        call inaccl(id,dr,dm(nl),dm(nw),nneq)
      endif
      rnorm = sqrt(dot(dr,dr,neq))
      write(iow,2003) rnorm
      if(ior.lt.0) write(*,2003) rnorm
      fl(8) = tr
      return

c     form a lumped mass approximation

4     if(fl(5)) call psetm(nl,neq,'d',fl(5))
      call pconsd(dm(nl),neq,0.0d0)
      fl(2) = tr
      call formfe(b,dm(nl),fa,tr,fa,fa,5,1,numel,1)
      return

c     compute reactions and print

5     if(pcomp(lct(l),'all ')) then
        n2 = numnp
      else
        n1 = ct(1,l)
        n2 = ct(2,l)
        n3 = ct(3,l)
        n1 = max(1,min(numnp,n1))
        n2 = max(n1,min(numnp,n2))
        n3 = max(1,n3)
      endif
      if(j.eq.6) then
        call pconsd(dr,nneq,0.0d0)
        call formfe(b,dr,fa,tr,fa,tr,6,1,numel,1)
        call prtrea(dr,ndf,numnp,n1,n2,n3)
      else
        if(nop.le.2) then
          call prtdis(x,dm(nw+nneq),ttim,prop,ndm,ndf,n1,n2,n3)
        else
          call prtdis(x,b,ttim,prop,ndm,ndf,n1,n2,n3)
        endif
      endif
      return

c     check mesh for input errors

6     call formfe(b,dr,fa,fa,fa,fa,2,1,numel,1)
      return

c     modify mesh data (cannot change profile of stiffness/mass)

9     i = -1
      call pmesh(ld,ie,d,id,x,ix,f,t,ndf,ndm,nen1,i,prt)
      if (i.gt.0) go to 400
      return

c     restart previously run problem

10    open (7,file=tfile(5),form='unformatted',status='old')
      read(7) nnpo,nnlo,nnmo,ndmo,ndfo,nhio,nhfo,nrco
      if((nnpo.eq.numnp).and.(nnlo.eq.numel).and.(nnmo.eq.nummat)
     1  .and.(ndmo.eq.ndm).and.(ndfo.eq.ndf).and.(nrco.eq.nrec)
     2  .and.(nhfo-nhio.eq.nhf-nhi) ) then
         read(7) ttim,(b(i),i=1,3*nneq)
         read(7) (dm(i),i=nt,nt+nneq)
         if(fl(9)) read(7) (dm(i),i=nv,nv+4*neq)
         if(nrec.gt.0) then
           do 101 j = 1,nrec
             call phstio(7,j,dm(nhi),nhf-nhi+1,11,tfile(5),iz)
             call phstio(3,j,dm(nhi),nhf-nhi+1,2,tfile(2),itrec(2))
101        continue
         endif
         close(7)
      else
         if(ior.gt.0) write(iow,3001)
         if(ior.lt.0) write(  *,3001)
      endif
      return

c     error diagnostics

400   write(iow,3002)
      if(ior.gt.0) stop
      write(  *,3002)

c     formats

2001  format('   | R(i) | =',1pe15.7,12x,'time=',a10)
2002  format(40x,'time=',a10)
2003  format('   | R(i) | = ',1pe15.7)
2004  format('   Energy Convergence Test'/'   E(1)=',1pe21.14,
     1  ', E(i)=',1pe21.14,', Tol.=',1pe11.4)
2006  format('   Shift of',1pe11.4,' applied with mass')
2007  format('   Shift requested but no mass matrix exists.')
3001  format(' **ERROR** File incompatible with current problem.')
3002  format(' **ERROR** Attempt to change profile during mesh.')
      end

      subroutine pmacr2(id,x,f,f0,jd,b,dr,lct,ct,ndf,nneq,j)

c     Purpose: Command language instruction executions

c     Inputs:
c        id(*)     - List of equation numbers for each dof
c        x(*)      - Nodal coordinates for mesh
c        f(*)      - Nodal specified force/displacmeents
c        f0(*)     - Base solution
c        jd(*)     - Pointer to end of row/columns for profiles
c        b(*)      - Nodal solution array
c        dr(*)     - Residual/reactions
c        lct       - Command option to execute
c        ct(3,*)   - Command parameters
c        ndf       - Number dof/node
c        nneq      - Size of id, f, f0, and b arrays
c        j         - Command number to execute

c     Outputs:
c        various   - Depends on command number

      implicit  none

      logical   pcomp,sfl
      integer   i
      integer   md,mf,mg,mh,mp,mq,mv,mdp,mdt, n,n1,n2,ndm,nlm
      real*8    propld

      integer   ndf,nneq,j
      integer   id(*),jd(*)
      real*8    x(*),f0(*),f(*),ct(3,*),xtl
      real*8    b(*),dr(*)
      character lct(*)*4,ctl(2)*4,yyy*80

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

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

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

      integer         ior,iow
      common /iofile/ ior,iow

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

      integer        l,lv,lvs   ,lve   ,jct
      common /ldata/ l,lv,lvs(9),lve(9),jct(100)

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

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

      integer        maxm,ne
      common /psize/ maxm,ne

      real*8         engy,rnmax,tol,shift
      common /rdata/ engy,rnmax,tol,shift

      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

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

      real*8          uu
      common /udata/  uu(4000)

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

c     transfer to correct process

      go to (1,2,3,4,5,6,7,8,9,10,11,12), j

c     set solution tolerance

1     tol = ct(1,l)
      return

c     set time increment

2     dt = ct(1,l)
      if(fl(9)) call setci(ior)
      return

c     set loop start indicators

3     lv = lv + 1
      lvs(lv) = l
      lve(lv) = ct(2,l)
      ct(1,lve(lv)) = 1.
      return

c     loop terminator control

4     n = ct(2,l)
      ct(1,l) = ct(1,l) + 1.0
      if(ct(1,l).gt.ct(1,n)) lv = lv - 1
      if(ct(1,l).le.ct(1,n)) l = n
      return

c     input proportional load table

5     npld = ct(1,l)
      prop = propld (ttim,npld)
      return

c     data command

6     if(ior.lt.0) write(*,3000) lct(l)
      call pintio(yyy,10)
      read(yyy,1000,err=61) (ctl(i),i=1,2),xtl
      if(.not.pcomp(lct(l),ctl(1))) go to 402
      if(pcomp(ctl(1),'tol ')) tol = xtl
      if(pcomp(ctl(1),'dt  ')) dt  = xtl
      return
61    call perror('PMACR2',yyy)
      go to 6

c     increment time - initialize force / solution vectors

7     ttim = ttim + dt
      do 71 i = 0,nneq-1
        dm(nt+i) = dm(nt+nneq+i)
71    continue
      if(npld.gt.0) prop = propld(ttim,0)
      write(iow,2002) ttim,prop
      if(ior.lt.0) write(  *,2002) ttim,prop
      engy = 0.0
      rnmax = 0.0

c     update history on the disk

      if(.not.hfl) then
        hout = .true.
        call formfe(b,dr,.false.,.false.,.false.,.false.,6,1,numel,1)
        hout = .false.
      endif

c     update dynamic vectors for time step

      if(fl(9)) then
        call setci(ior)
        call update(id,dm(nt),b,dm(nv),dr,nneq,fl(9),1)
      endif

c     zero displacement increment for next time step

      call pconsd(b(nneq+1),nneq+nneq,0.0d0)
      fl(10) = .true.
      return

c     input integration parameters and initialize vectors

8     call param(lct(l),ct(1,l))
      if(fl(9)) return
      call psetm(nv,nneq*4,'d',fl(9))
      call pconsd(dm(nv),nneq*4,0.0d0)
      nw = nv + nneq

c     set initial condition for transient solution

      n1 = nw + nneq - 1
      do 81 i = 1,nneq
        dm(n1+i) = b(i)
81    continue
      fl(9) = .true.
      return

c     update the current force vector f0

9     call raxpb(f, f0, prop, nneq, f0)
      return

c     subspace eigencomputations

10    mf = ct(1,l)
      n2 = ct(2,l)
      mf = min(neq,max(1,mf))
      mq = min(mf+mf,mf+8,neq)
      if(n2.gt.0) mq = min(mf+n2,neq)
      if(fl(2)) call numass(dm(nl),neq,mq)
      if(mq.lt.mf) write(iow,2001) mq
      if(mq.lt.mf.and.ior.lt.0) write(  *,2001) mq
      mf = min(mf,mq)
      if(mf.le.0) return
      sfl = pcomp(lct(l),'prin')

c     establish addresses for eigen solutions

      md =  1
      mv = md + mq
      mg = mv + mq*neq
      mh = mg + mq*(mq+1)/2
      mdp= mh + mq*(mq+1)/2
      mdt= mdp+ mq
      mp = mdt+ mq
      nlm= max(mp+mq*mq,neq*mq+neq)
      if(nlm.gt.maxm/ipd) then
        if(ior.lt.0) then
          write(  *,3001) nlm,maxm/ipd
          return
        else
          write(iow,3001) nlm,maxm/ipd
          stop
        endif
      endif
      nlm = ne/ipd + 1
      open(7,file=tfile(3),form='unformatted',status='new')
      write(7) (dm(i),i=1,nlm)
      close (7)
      nlm = nl - 1
      n1  = neq*ipd - ipd - ipr
      do 101 i = 1,neq
        dm(i)    = dm(i+nlm)
        im(n1+i) = jd(i)
101   continue
      nlm = neq + (neq + ipd -1)/ipd + 1
      call subsp(dm(1),uu(mv),uu(mg),uu(mh),uu(md),uu(mdp),
     1   uu(mdt),uu(mp),dm(nlm),mf,mq,neq,shift,tol,sfl,25)

c     restore the solution to continue macro executions

      nlm = ne/ipd + 1
      open(7,file=tfile(3),form='unformatted',status='old')
      rewind 7
      read(7) (dm(i),i=1,nlm)
      close(7,status='delete')
      return

c     print eigenvectors

11    if(mf.gt.0) then
        n1   = ct(1,l)
        n1   = max(1,min(mq,n1))
        n2   = mv + neq*(n1-1) -1

c     expand and move eigenvector for prints

        call pconsd(dm,nneq,0.0d0)
        do 111 i = 1,nneq
          n = id(i)
          if(n.ne.0) dm(i) = uu(n2+n)
111     continue
        call prtdis(x,dm,ttim,uu(md+n1-1),ndm,ndf,1,numnp,1)
      else
        write(iow,3002)
        if(ior.lt.0) write(*,3002)
      endif
      return

c     set identity matrix for mass

12    if(fl(5)) call psetm(nl,neq,'d',fl(5))
      fl(2) = .true.
      call pconsd(dm(nl),neq,1.d0)
      return

c     error diagnostics

402   if(ior.gt.0) write(iow,3003)
      if(ior.gt.0)  stop
      write(*,3003)
      go to 6

c     formats

1000  format(a4,6x,a4,6x,f15.0)
2001  format(' **WARNING** Subspace set to',i4,' number of mass terms.')
2002  format('   Computing solution for time',1pe12.5/
     1       '   Proportional load value is ',1pe12.5)
3000  format(' Input ',a4,' macro >',$)
3001  format(' **ERROR** Subspace memory too large'/
     1 5x,'Need =',i7,' : Available =',i7)
3002  format(' **ERROR** Need eigensolution.')
3003  format(' **ERROR** Macro label mismatch on data command.')
      end

      subroutine formfe(u,b,af,bf,cf,df,is,ne1,ne2,ne3)

c     Purpose: Form finite element arrays as required

c     Inputs:
c        u(*)      - Current solution state
c        af        - Flag, assemble matrix if true
c        bf        - Flag, assemble vector if true
c        cf        - Flag, matrix unsymmetric if true
c        df        - Flag, assemble vector in full form if true
c        is        - Switch to control actions to be taken
c        ne1       - First element to process
c        ne2       - Last  element to process
c        ne3       - Increment to ne1

c     Outputs:
c        b(*)      - Element vector

      implicit  none
      integer   is,ne1,ne2,ne3

      logical   af,bf,cf,df
      real*8    u(*),b(*)

      integer        nn,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13
      common /mdata/ nn,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13

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

      integer        ndf,ndm,nen1,nst
      common /sdata/ ndf,ndm,nen1,nst

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

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

c     form appropriate f.e. array

      afl = af
      bfl = bf
      cfl = cf
      dfl = df
      isw = is
      nn1 = ne1
      nn2 = ne2
      nn3 = ne3
      call pform(dm(nn),dm(n0),dm(n1),m(n2),dm(n3),dm(n4),m(n5),dm(n6),
     &         m(n7),dm(n8),m(n9),dm(n11),m(n11c),u,b,ndf,ndm,nen1,nst)
      end

      subroutine inaccl(id,dr,xm,a,nneq)

c     Purpose:  Set initial acceleration
c        id(*)     - Equation numbers for each dof
c        dr(*)     - Residual for initial state
c        xm(*)     - Diagonal mass matrix
c        nneq      - Number of terms in id

c     Outputs:
c        a(*)      - Initial accelerations

      implicit  none

      integer   nneq, j, n

      integer   id(nneq)
      real*8    dr(*),xm(*),a(nneq)

      do 100 n = 1,nneq
        j = id(n)
        if(j.gt.0) a(n) = dr(j)/xm(j)
100   continue

      end

      subroutine phstio(iu,irec,hh,nh,isw,tfile,itrec)

c     Purpose: I/O control for history terms on disk

c     Inputs:
c        iu         - Logical unit number for I/O
c        irec       - Record number
c        hh(*)      - History terms to write (isw = 2 or 22 only)
c        nh         - Number of history terms
c        isw        - Switch: =1 or 11 for read; = 2 or 22 for write
c        tfile      - File name for history terms
c        itrec      - Record length of data

c     Outputs:
c        hh(*)      - History terms read (isw = 1 or 11 only)

      implicit  none

      character tfile*12
      integer   iu,irec, nh,isw, itrec

      real*8    hh(nh)

c     direct access read/write

      write(*,*) 'RECL =',tfile,' ',itrec
      if(isw.lt.10) open(iu,file=tfile,access='direct',recl=itrec,
     &                   status='unknown')

      if(isw.eq.1) read (iu,rec=irec) hh
      if(isw.eq.2) write(iu,rec=irec) hh
      if(isw.lt.10) close(iu)

c     sequential access read/write

      if(isw.eq.11) read (iu) hh
      if(isw.eq.22) write(iu) hh

      end

      subroutine serchl(fn,id,g0,rsd,u,d,stol,t,nneq)

c     linear line search for nonlinear problems

      implicit  none

      integer   nneq, j, linmax
      real*8    s,stol,step

      integer   id(*)
      real*8    fn(nneq,2)
      real*8    rsd(*),u(*),d(*),t(*),gamma,g,g0,ga,gb,sa,sb

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

      integer         ior,iow
      common /iofile/ ior,iow

c     compute step size for line search in direction d

      linmax = 10
      sb     = 0.0
      sa     = 1.0
      s      = 1.0
      g      = gamma(fn,id,u,rsd,d,t,s,neq,nneq)

c     find bracket on zero

      if(g*g0.gt.0.0d0) then
        write(iow,3000)
        if(ior.lt.0) write(*,3000)
      else
        j    = 0
        gb   = g0
        ga   = g
        sb   = 0.0d0
        sa   = 1.0d0
10      j    = j + 1
        step = sa - ga*(sa-sb)/(ga-gb)
        g    = gamma(fn,id,u,rsd,d,t,step,neq,nneq)
        gb   = 0.5d0*gb
        if (g*ga.lt.0.0d0) then
          sb = sa
          gb = ga
        endif
        sa = step
        ga = g
        write(iow,3001) j,step,g
        if(ior.lt.0) write(  *,3001) j,step,g
        if (j.lt.linmax) then
          if(abs(g)     .gt. stol*abs(g0)      ) go to 10
          if(abs(sb-sa) .gt. stol*0.5d0*(sa+sb)) go to 10
        endif
        do 20 j = 1, neq
          d(j) = step*d(j)
20      continue
      endif

3000  format(4x,'No line search, end points both positive.')
3001  format(4x,'Iter =',i2,' Step Size =',e12.5,' Energy =',e12.5)

      end

      double precision function gamma(fn,id,u,dr,du,t,s,nqe,nneq)

c     Purpose:  Compute energy for line search step

c     Inputs:
c        fn(*)     - Force vector for active dof's
c        id(*)     - Equation numbers for each dof
c        u(*)      - Currenet solution
c        dr(*)     - Residual
c        du(*)     - Current increment to solution
c        t(*)      - Working space
c        nqe       - Number equations
c        nneq      - Total dof in problem

c     Outputs:
c        gamma     - Value of energy for step

      implicit  none

      integer   j, n, nqe,nneq, nneq2
      real*8    s, ctem

      logical   fa,tr
      integer   id(*)
      real*8    fn(nneq,2)
      real*8    u(*),dr(*),du(*),t(*),db,dot

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

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

      integer        nn,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13
      common /mdata/ nn,n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13

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

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

      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

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

      data fa,tr/.false.,.true./


c     get a search displacement

      nneq2 = nneq + nneq
      ctem  = theta
      theta = s*theta
      do 100 n = 1,nneq
        j = id(n)
        if(j.gt.0) then
          db         = s*du(j)
          t(n)       = u(n)      + db
          t(n+nneq ) = u(n+nneq) + db
          t(n+nneq2) = db
        else
          db = theta*fn(n,2) + (1.0-theta)*fn(n,1)
          t(n+nneq2) = db - u(n)
          t(n+nneq ) = u(n+nneq) - u(n) + db
          t(n)       = db
        endif
100   continue

c     compute a residual

      call pload(id,fn,dr,nneq,dm(nl),dm(nw))
      call formfe(t,dr,fa,tr,fa,fa,6,1,numel,1)
      theta = ctem

c     update the residual for lumped mass inertial effects

      if(fl(9)) then
        ctem = s*c1
        do 110 n = 1,nneq
          j = id(n)
          if(j.gt.0) then
            dr(j) = dr(j) - dm(nl+j-1)*(dm(nw+n-1) + ctem*du(j))
          endif
110     continue
      endif

c     compute the value of gamma

      gamma = dot (du,dr,nqe)

      end
