anchor.diff:             temp_id = (j-1) * p%dof_node
              TmpPos(1:3) = p%GlbPos(1:3) + MATMUL(p%GlbRot,p%uuN0(temp_id+1:temp_id+3,i))
              temp_CRV(:) = MATMUL(p%GlbRot,p%uuN0(temp_id+4:temp_id+6,i))
  -           CALL BD_CrvMatrixR(temp_CRV,TmpDCM)
  +           CALL BD_CrvMatrixR(temp_CRV,TmpDCM,ErrStat2,ErrMsg2)
              TmpDCM(:,:) = TRANSPOSE(TmpDCM)
              temp_id = (i-1)*(p%node_elem-1)+j
              CALL MeshPositionNode ( Mesh    = u%PointLoad  &
  
positive.code:    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)
         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)
         mmm  = 0.0D0
         mEta = 0.0D0
         rho  = 0.0D0
  
negative.code:    ENDDO
     rot_temp = 0.0D0
     CALL BD_CrvCompose(rot_temp,rot0_temp,rotu_temp,0,ErrStat2,ErrMsg2)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     DO i=1,3
         uu0(i+3) = rot_temp(i)
     ENDDO
  