anchor.diff:  c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   c
   c
  -      subroutine preshape(sh,shj,gauss,ngauss,nen,nec,nee,ndof,nsd,
  -     & nenmax,ngaussmax,netypes,intord)
  +      subroutine preshape(sh,shj,gauss,infetype,intord)
   c
   c...program to compute the shape functions in natural coordinates.
   c   The program also assigns integration points and weights.
  
positive.code: c   The arguments are:
  c     Input:
  c
  c       ndof                               = number of degrees of
  c                                            freedom per node
  c       nsd                                = number of spatial
  c                                            dimensions
  c       nenmax                             = maximum number of nodes
  c                                            per element
  c       ngaussmax                          = maximum number of gauss
  c                                            points per element
  c       netypes                            = number of element types
  c                                            (currently equal to 62)
  c       intord                             = integration order
  c                                            1 = full
  c                                            2 = reduced
  
negative.code: c
  c...  initialize arrays
  c
        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),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
          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,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))
  