anchor.diff:          if(ngravflag.ne.izero) call daxpy(neq,one,bgravity,ione,
        &   bextern,ione)
           if(nconcflag.ne.izero) call daxpy(neq,one,bconcforce,ione,
  -     &   bextern,ione)
  -        if(nprestrflag.ne.izero) call daxpy(neq,alpha,bprestress,ione,
        &   bextern,ione)
           call dcopy(neq,bextern,ione,bresid,ione)
   cdebug        write(6,*) "bextern:"
  
positive.code: c
  c...  subroutine arguments
  c
        integer nextflag,ntractflag,ngravflag,nconcflag,nprestrflag,neq
        double precision bextern(nextflag*neq),btraction(ntractflag*neq)
        double precision bgravity(ngravflag*neq),bconcforce(nconcflag*neq)
        double precision bprestress(nprestrflag*neq),bintern(neq)
        double precision bresid(neq)
  c
  cdebug      integer idb
  c
  
negative.code: c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  c
  c
        subroutine bsum(bextern,btraction,bgravity,bconcforce,bintern,
       & bresid,nextflag,ntractflag,ngravflag,nconcflag,neq)
  c
  c...subroutine to sum all contributions to the external force vector,
  c   then subtract the internal force vector from the external vector
  