anchor.diff:     DO igp=1,ngp
          gpr=gp(igp)
   
  -       CALL BD_ComputeJacobian(gpr,Nuu0,node_elem,dof_node,gp,GLL_temp,ngp,igp,hhx,hpx,Jacobian)
  -       CALL BD_GaussPointDataAt0(hhx,hpx,Nuu0,Nrr0,node_elem,dof_node,uu0,E10)
  +       CALL BD_ComputeJacobian(gpr,Nuu0,node_elem,dof_node,gp,GLL_temp,ngp,igp,hhx,hpx,Jacobian,ErrStat2,ErrMsg2)
  +       CALL BD_GaussPointDataAt0(hhx,hpx,Nuu0,Nrr0,node_elem,dof_node,uu0,E10,ErrStat2,ErrMsg2)
          Stif(:,:) = 0.0D0
          Stif(1:6,1:6) = EStif0_GL(1:6,1:6,igp)
  -       CALL BD_GaussPointData(hhx,hpx,Nuuu,Nrrr,uu0,E10,node_elem,dof_node,uuu,uup,E1,RR0,kapa,Stif,cet)
  +       CALL BD_GaussPointData(hhx,hpx,Nuuu,Nrrr,uu0,E10,node_elem,dof_node,&
  +                              uuu,uup,E1,RR0,kapa,Stif,cet,ErrStat2,ErrMsg2)
          mmm  = 0.0D0
          mEta = 0.0D0
          rho  = 0.0D0
  
positive.code:    CHARACTER(ErrMsgLen)        :: ErrMsg2                      ! Temporary Error message
     CHARACTER(*), PARAMETER     :: RoutineName = 'BD_GaussPointDataAt0'
  
     uu0 = 0.0D0
     E10 = 0.0D0
     DO inode=1,node_elem
         hhi = hhx(inode)
         hpi = hpx(inode)
  
negative.code:    DO nelem=1,elem_total
         Nuu0(:) = uuN0(:,nelem)
         CALL BD_ElemNodalDisp(uuN,node_elem,dof_node,nelem,Nuuu,ErrStat2,ErrMsg2)
         CALL BD_NodalRelRot(Nuu0,node_elem,dof_node,Nrr0,ErrStat2,ErrMsg2)
         CALL BD_NodalRelRot(Nuuu,node_elem,dof_node,Nrrr,ErrStat2,ErrMsg2)
         CALL BD_ElemNodalDisp(vvN,node_elem,dof_node,nelem,Nvvv,ErrStat2,ErrMsg2)
         CALL BD_ElemNodalDisp(aaN,node_elem,dof_node,nelem,Naaa,ErrStat2,ErrMsg2)
         temp_id = (nelem-1)*ngp
  