anchor.diff:  c
   c.... compute the forces on the body
   c
  -c.... DEBUGGING
  -c          where (nsrflist(iBCB(:,2)).ge.1)
  -c            tau1n = ( pres * bnorm(:,1) - tau1n ) * WdetJb
  -c            tau2n = ( pres * bnorm(:,2) - tau2n ) * WdetJb
  -c            tau3n = ( pres * bnorm(:,3) - tau3n ) * WdetJb
  -c            pres1n =  pres * bnorm(:,1) * WdetJb
  -c            pres2n =  pres * bnorm(:,2) * WdetJb
  -c            pres3n =  pres * bnorm(:,3) * WdetJb
  -c            heat  = - heat * WdetJb
  -c          elsewhere
  -c            tau1n = zero
  -c            tau2n = zero
  -c            tau3n = zero
  -c            pres1n = zero
  -c           pres2n = zero
  -c            pres3n = zero
  -c            heat  = zero
  -c          endwhere
  -c
  -c.... HARDCODED
  -c        loc(1:300)   = 0.0
  -c        loc(301:377) = 1.6924 * 1e-6
  -c        loc(378:390) = 1.8338 * 1e-6
  -c        loc(391:402) = 1.8738 * 1e-6
  -c        loc(403:415) = 2.0562 * 1e-6
  -c        loc(416:427) = 2.3090 * 1e-6
  -c        loc(428:440) = 2.2736 * 1e-6
  -c        loc(441:452) = 1.7806 * 1e-6
  -
  -c        dyn_org = 0.0
  -c        do i = 300,lstep      
  -c          dyn_org = dyn_org + loc(i-1) 
  -c        enddo
  -c        countWrong = 0; 
  -c.... HARDCODED
  -c
           do i = 1, npro
             if (nsrflist(iBCB(i,2)).ge.1) then
  -c            if ((xlb(i,1,1) .ge. dyn_org) .and. (iBCB(i,2) .eq. 10)) then
  -c               iBCB(i,2)  = 20
  -c               countWrong = countWrong + 1 
  -c            endif
  -c            if ((xlb(i,1,1). lt. dyn_org) .and. (iBCB(i,2) .eq. 20)) then
  -c               iBCB(i,2)  = 10
  -c               countWrong = countWrong + 1
  -c            endif
               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)
  
positive.code:      &            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
  
negative.code:             HFlux(nsrflist(iBCB(i,2))) =
       &      HFlux(nsrflist(iBCB(i,2))) + heat(i)
            enddo
          endif
  c
  c.... end of integration loop
  