anchor.diff:  c
   c...  initialize arrays
   c
  -      call iclear(nen,netypes)
  -      call iclear(ngauss,netypes)
  -      call clear(gauss,ngssize*netypes)
  -      call clear(sh,nshsize*netypes)
  +      call ifill(infetype,izero,ifour*netypes)
  +      call fill(gauss,zero,ngssize*netypes)
  +      call fill(sh,zero,nshsize*netypes)
  +      call fill(shj,zero,nshsize*netypes)
   c
   c... First type:  linear hex
   c
  -      call plinhex(sh(1,1,1,1),shj(1,1,1,1),gauss(1,1,1),
  -     & ngauss(1),nen(1),nenmax,ngaussmax,intord)
  +      call plinhex(sh(1,1,1,1),shj(1,1,1,1),gauss(1,1,1),infetype(1,1),
  +     & intord)
   c
   c...  Types 2-27:  linear hex + infinite boundaries
   c     Use mod3 arithmetic to cover all cases
   c
         do n=2,27
  -        nen(n)=nen(1)
  -        ngauss(n)=ngauss(1)
  +        do i=1,4
  +          infetype(i,n)=infetype(i,1)
  +        end do
           call dcopy(ngssize,gauss(1,1,1),ione,gauss(1,1,n),ione)
           call dcopy(nshsize,sh(1,1,1,1),ione,sh(1,1,1,n),ione)
           io(3)=(n-1)/inine
           io(2)=(n-1-io(3)*inine)/ithree
           io(1)=mod(n-1,ithree)
  -        do l=1,ngauss(n)
  -          do i=1,nen(n)
  +        do l=1,infetype(1,n)
  +          do i=1,infetype(2,n)
               call infellh(gauss(1,l,n),rlh(i),rr,drr,io(1))
               call infellh(gauss(2,l,n),slh(i),ss,dss,io(2))
               call infellh(gauss(3,l,n),tlh(i),tt,dtt,io(3))
  
positive.code: c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  c
  c
        subroutine preshape(sh,shj,gauss,ngauss,nen,nec,nee,ndof,nsd,
       & nenmax,ngaussmax,netypes,intord)
  c
  c...program to compute the shape functions in natural coordinates.
  c   The program also assigns integration points and weights.
  
negative.code: c
  c...  Type 28:  linear hex with one set of collapsed nodes (wrick)
  c
        infetype(1,28)=infetype(1,1)
        infetype(2,28)=iseven
        infetype(3,28)=nsd*infetype(2,28)
        infetype(4,28)=ndof*infetype(2,28)
        call dcopy(ngssize,gauss(1,1,1),ione,gauss(1,1,28),ione)
        call dcopy(nshsize,sh(1,1,1,1),ione,sh(1,1,1,28),ione)
        do l=1,infetype(1,28)
          do i=1,4
            sh(i,7,l,28)=sh(i,7,l,28)+sh(i,8,l,28)
            sh(i,8,l,28)=zero
  