anchor.diff:  !----------------------------------------------------------------------------------------
   ! This subroutine computes gyroscopic forces
   !----------------------------------------------------------------------------------------
  -   REAL(BDKi),    INTENT(IN   ):: mEta(:) ! m\Eta resolved in inertia frame at Gauss point
  -   REAL(BDKi),    INTENT(IN   ):: rho(:,:) ! Tensor of inertia resolved in inertia frame at Gauss point
  -   REAL(BDKi),    INTENT(IN   ):: vvv(:) ! Velocities at Gauss point (including linear and angular velocities)
  -   REAL(BDKi),    INTENT(  OUT):: Fb(:) ! Gyroscopic forces
  +   REAL(BDKi),    INTENT(IN   ):: mEta(:)       ! m\Eta resolved in inertia frame at Gauss point
  +   REAL(BDKi),    INTENT(IN   ):: rho(:,:)      ! Tensor of inertia resolved in inertia frame at Gauss point
  +   REAL(BDKi),    INTENT(IN   ):: vvv(:)        ! Velocities at Gauss point (including linear and angular velocities)
  +   REAL(BDKi),    INTENT(  OUT):: Fb(:)         ! Gyroscopic forces
      INTEGER(IntKi),INTENT(  OUT):: ErrStat       ! Error status of the operation
      CHARACTER(*),  INTENT(  OUT):: ErrMsg        ! Error message if ErrStat /= ErrID_None
   
  
positive.code:        CALL BD_ElasticForce(E1,RR0,kapa,Stif,cet,fact,Fc,Fd,Oe,Pe,Qe,ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         mmm          = EMass0_GL(1,1,igp)
         mEta(2)      =-EMass0_GL(1,6,igp)
         mEta(3)      = EMass0_GL(1,5,igp)
         rho(1:3,1:3) = EMass0_GL(4:6,4:6,igp)
  
negative.code:       do l = 1,np+1
           
         if ((abs(spts(j)-1.).LE.eps).AND.(l.EQ.np+1)) then
           dPhis(l,j) = REAL((np+1)*np, BDKi)/4.0_BDKi
         elseif ((abs(spts(j)+1.).LE.eps).AND.(l.EQ.1)) then
           dPhis(l,j) = -REAL((np+1)*np, BDKi)/4.0_BDKi
         elseif (abs(spts(j)-npts(l)).LE.eps) then
           dPhis(l,j) = 0.0_BDKi
         else
           dPhis(l,j) = 0.0_BDKi
           den = 1.0_BDKi          
           do i = 1,np+1
             if (i.NE.l) then
               den = den*(npts(l)-npts(i))
             endif
             dnum = 1.0_BDKi
             do k = 1,np+1
               if ((k.NE.l).AND.(k.NE.i).AND.(i.NE.l)) then
                 dnum = dnum*(spts(j)-npts(k))
               elseif (i.EQ.l) then
                 dnum = 0.0_BDKi
               endif
             enddo
             dPhis(l,j) = dPhis(l,j) + dnum
  