anchor.diff:       &            um1(npro),                   um2(npro),  
        &            um3(npro)
   c
  -c.... hardcoded
  -        real*8    loc(500),                dyn_org
  -        integer   countWrong
  -c.... hardcoded
           dimension lnode(27),               sgn(npro,nshl),
        &            shape(npro,nshl),        shdrv(npro,nsd,nshl)
   c
  
positive.code:             HFlux(nsrflist(iBCB(i,2))) =
       &      HFlux(nsrflist(iBCB(i,2))) + heat(i)
            enddo
  c... END DEBUGGING
  c
  c          if (countWrong .gt. 0) write(*,*) 'countWrong = ', countWrong
  c... OLD lines
  c          Force(1,nsrflist(iBCB(:,2))) = Force(1,nsrflist(iBCB(:,2))) + sum(tau1n)
  c          Force(2,nsrflist(iBCB(:,2))) = Force(2,nsrflist(iBCB(:,2))) + sum(tau2n)
  c          Force(3,nsrflist(iBCB(:,2))) = Force(3,nsrflist(iBCB(:,2))) + sum(tau3n)
  c          PresFor(1,nsrflist(iBCB(:,2))) =
  c     &    PresFor(1,nsrflist(iBCB(:,2))) + sum(pres1n)
  c          PresFor(2,nsrflist(iBCB(:,2))) =
  c     &    PresFor(2,nsrflist(iBCB(:,2))) + sum(pres2n)
  c          PresFor(3,nsrflist(iBCB(:,2))) =
  c     &    PresFor(3,nsrflist(iBCB(:,2))) + sum(pres3n)
  c          HFlux(nsrflist(iBCB(:,2)))   = HFlux(nsrflist(iBCB(:,2)))   + sum(heat)
  c... End OLD lines
          endif
  c
  c.... end of integration loop
  
negative.code: c
  c.... compute the forces on the body
  c
          do i = 1, npro
            if (nsrflist(iBCB(i,2)).ge.1) then
              tau1n(i) = ( pres(i) * bnorm(i,1) - tau1n(i) ) * WdetJb(i)
              tau2n(i) = ( pres(i) * bnorm(i,2) - tau2n(i) ) * WdetJb(i)
              tau3n(i) = ( pres(i) * bnorm(i,3) - tau3n(i) ) * WdetJb(i)
  