anchor.diff:  c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   c
   c
  -      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.
  -c   All shape functions are stored in a single array.
  +      subroutine preshape(sh,shj,gauss,intord,ietype,nen,ngauss,ierr,
  +     & errstrng)
  +c
  +c...program to compute information for a given element type.
   c
   c   The arguments are:
   c     Input:
   c
  -c       intord                             = integration order
  -c                                            1 = full
  -c                                            2 = reduced
  -c                                            3 = Bbar
  +c       intord                  = integration order
  +c                                 1 = full
  +c                                 2 = reduced
  +c                                 3 = Bbar
  +c       ietype                  = element type
  +c       nen                     = number of element nodes
  +c       ngauss                  = number of element gauss points
   c
   c     Output:
  -c       sh(nsd+1,nenmax,
  -c          ngaussmax,netypes)   = shape functions and their derivatives
  -c       shj(nsd+1,nenmax,
  -c           ngaussmax,netypes)  = shape functions and derivatives used
  +c       sh(nsd+1,nen,ngauss)    = shape functions and their derivatives
  +c       shj(nsd+1,nen,ngauss)   = shape functions and derivatives used
   c                                 for calculating Jacobian.  Unless the
   c                                 element is an infinite element, this
   c                                 is identical to sh.
  -c       gauss(nsd+1,ngaussmax,netypes)     = Gauss point coordinates and
  -c                                            weights
  -c       infetype(4,netypes)                = element type info array
  -c         infetype(1,i) = number of gauss points for element type
  -c         infetype(2,i) = number of nodes for element type
  -c         infetype(3,i) = number of element coordinate directions for
  -c                         element type
  -c         infetype(4,i) = number of element degrees of freedom for
  -c                         element type
  +c       gauss(nsd+1,ngauss)     = Gauss point coordinates and weights
   c
         include "implicit.inc"
   c
  
positive.code: c
  c...  local variables
  c
        integer i,l,n,nshsize,ngssize
        double precision rr,ss,tt,drr,dss,dtt
        integer io(3)
  c
  
negative.code: c
  c...  subroutine arguments
  c
        integer intord,ietype,nen,ngauss,ierr
        double precision sh(nsd+1,nen,ngauss),shj(nsd+1,nen,ngauss)
        double precision gauss(nsd+1,ngauss)
        character*(*) errstrng
  c
  c...  local constants
  c
  