anchor.diff:            gauss(4,l)=fourth*tetvol
           end do
         end if
  +c
  +      infetype(1)=ngauss
  +      infetype(2)=nen
  +      infetype(3)=nec
  +      infetype(4)=nee
  +c
         do l=1,ngauss
           do i=1,nen
             rr=r(i)*gauss(1,l)
  
positive.code:             drr=eight*r(i)*max(ss,one)*max(tt,one)*max(uu,one)
              dss=eight*s(i)*max(rr,one)*max(tt,one)*max(uu,one)
              dtt=eight*t(i)*max(rr,one)*max(ss,one)*max(uu,one)
              duu=-eight*u(i)*max(rr,one)*max(ss,one)*max(tt,one)
              sh(1,i,l)=drr-duu
              sh(2,i,l)=dss-duu
              sh(3,i,l)=dtt-duu
  
negative.code: c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  c
  c
        subroutine pquadtet(sh,shj,gauss,infetype,intord)
  c
  c... Subroutine to compute shape functions in natural coordinates,
  c    integration points, and weights for a quadratic tetrahedron.
  c
        include "implicit.inc"
  c
  c...  parameter definitions
  c
        include "ndimens.inc"
        include "nshape.inc"
        include "nconsts.inc"
        include "rconsts.inc"
  c
  c...  subroutine arguments
  c
        integer intord
        integer infetype(4)
        double precision sh(nsd+1,nenmax,ngaussmax)
        double precision shj(nsd+1,nenmax,ngaussmax)
        double precision gauss(nsd+1,ngaussmax)
  c
  c...  local constants
  c
        double precision r(10),s(10),t(10),u(10)
        data r/ 1d0, 0d0, 0d0, 0d0,0.5d0,0.0d0,0.5d0,0.5d0,0.0d0,0.0d0/
  