anchor.diff:  
      DO i=1,p%ngp*p%elem_total+2
          TmpPos(1:3) = p%GlbPos(1:3) + MATMUL(p%GlbRot,temp_L2(1:3,i))
  +       temp_CRV(:) = MATMUL(p%GlbRot,temp_L2(4:6,i))
  +       CALL BD_CrvMatrixR(temp_CRV,TmpDCM)
  +       TmpDCM(:,:) = TRANSPOSE(TmpDCM)
          CALL MeshPositionNode ( Mesh    = u%DistrLoad  &
                                 ,INode   = i            &
                                 ,Pos     = TmpPos       &
                                 ,ErrStat = ErrStat2      &
                                 ,ErrMess = ErrMsg2      &
  -                              , Orient = ??? ) !bjj: add orient=DCM ! Orientation (direction cosine matrix) of node; identity by default
  +                              , Orient = TmpDCM ) !bjj: add orient=DCM ! Orientation (direction cosine matrix) of node; identity by default
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
          
      ENDDO
  
positive.code:        DO j=1,p%node_elem
             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_id = (i-1)*(p%node_elem-1)+j
             CALL MeshPositionNode ( Mesh    = u%PointLoad  &
                                    ,INode   = temp_id      &
                                    ,Pos     = TmpPos       &
                                    ,ErrStat = ErrStat2      &
                                    ,ErrMess = ErrMsg2       &
                                    , Orient = ??? ) !bjj: add orient=DCM ! Orientation (direction cosine matrix) of node; identity by default
              CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         ENDDO
     ENDDO
     DO i=1,p%elem_total
         DO j=1,p%node_elem
             temp_id = (j-1)*p%dof_node
  !           TmpPos(1:3) = p%uuN0(temp_id+1:temp_id+3,i)
             TmpPos(1:3) = p%GlbPos(1:3) + MATMUL(p%GlbRot,p%uuN0(temp_id+1:temp_id+3,i))
             temp_id = (i-1)*p%node_elem+j
             CALL MeshPositionNode ( Mesh    = y%BldMotion  &
                                    ,INode   = temp_id      &
                                    ,Pos     = TmpPos       &
                                    ,ErrStat = ErrStat2      &
                                    ,ErrMess = ErrMsg2       &
                                    , Orient = ??? ) !bjj: add orient=DCM ! Orientation (direction cosine matrix) of node; identity by default
              CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
             
         ENDDO
  
negative.code:    DEALLOCATE(temp_w)
     ! temp_L2: the DistrLoad mesh node location
     ! temp_L2: physical coordinates of Gauss points and two end points
     CALL AllocAry(temp_L2,6,p%ngp*p%elem_total+2,'temp_L2',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     ! Temporary Gauss point intrinsic coordinates array
     CALL AllocAry(temp_GL,p%ngp,'temp_GL',ErrStat2,ErrMsg2)
  