anchor.diff:               uu0(temp_id1+1:temp_id1+6),E10(temp_id2+1:temp_id2+3),node_elem,dof_node,&
                uuu,uup,E1,RR0,kapa,Stif,cet,ErrStat2,ErrMsg2)
             CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  -       mmm  = 0.0_BDKi
  -       mEta = 0.0_BDKi
  -       rho  = 0.0_BDKi
  +       mmm          = 0.0_BDKi
  +       rho          = 0.0_BDKi
          mmm          = EMass0_GL(1,1,igp)
  +       mEta(1)      = 0.0_BDKi
          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)
  
positive.code: ! This subroutine computes Global mass matrix and force vector to 
  ! calculate the forces along the beam
  !----------------------------------------------------------------------------------------
     REAL(BDKi),         INTENT(IN   ):: uuN(:) ! Displacement of Mass 1: m
     REAL(BDKi),         INTENT(IN   ):: vvN(:) ! Velocity of Mass 1: m/s
     REAL(BDKi),         INTENT(IN   ):: aaN(:) ! Velocity of Mass 1: m/s
     REAL(BDKi),         INTENT(IN   ):: Stif0(:,:,:) ! Element stiffness matrix
     REAL(BDKi),         INTENT(IN   ):: Mass0(:,:,:) ! Element stiffness matrix
     TYPE(BD_InputType), INTENT(IN   ):: u           ! Inputs at t
     INTEGER(IntKi),     INTENT(IN   ):: damp_flag ! Number of Gauss points
     REAL(BDKi),         INTENT(IN   ):: beta(:)
     INTEGER(IntKi),     INTENT(IN   ):: elem_total ! Total number of elements
     INTEGER(IntKi),     INTENT(IN   ):: node_elem ! Node per element
     INTEGER(IntKi),     INTENT(IN   ):: dof_node ! Degrees of freedom per node
     INTEGER(IntKi),     INTENT(IN   ):: ngp ! Number of Gauss points
     INTEGER(IntKi),     INTENT(IN   ):: quadrature ! Number of Gauss points
     REAL(BDKi),        INTENT(IN   ):: gw(:)
     REAL(BDKi),        INTENT(IN   ):: hhx(:,:)
     REAL(BDKi),        INTENT(IN   ):: hpx(:,:)
     REAL(BDKi),        INTENT(IN   ):: Jaco(:,:)
     REAL(BDKi),        INTENT(IN   ):: uu0(:,:)
     REAL(BDKi),        INTENT(IN   ):: E10(:,:)
     REAL(BDKi),         INTENT(  OUT):: RHS(:) ! Right hand side of the equation Ax=B
     INTEGER(IntKi),     INTENT(  OUT):: ErrStat       ! Error status of the operation
     CHARACTER(*),       INTENT(  OUT):: ErrMsg        ! Error message if ErrStat /= ErrID_None
  
     REAL(BDKi),           ALLOCATABLE:: Nuuu(:) ! Nodal displacement of Mass 1 for each element
     REAL(BDKi),           ALLOCATABLE:: Nrrr(:) ! Nodal rotation parameters for displacement of Mass 1
     REAL(BDKi),           ALLOCATABLE:: Nvvv(:) ! Nodal velocity of Mass 1: m/s for each element
     REAL(BDKi),           ALLOCATABLE:: Naaa(:) ! Nodal velocity of Mass 1: m/s for each element
     REAL(BDKi),           ALLOCATABLE:: EStif0_GL(:,:,:) ! Nodal material properties for each element
     REAL(BDKi),           ALLOCATABLE:: EMass0_GL(:,:,:) ! Nodal material properties for each element
     REAL(BDKi),           ALLOCATABLE:: DistrLoad_GL(:,:) ! Nodal material properties for each element
     REAL(BDKi),           ALLOCATABLE:: elf(:) ! Total element force (Fc, Fd, Fb)
     INTEGER(IntKi)                   :: dof_elem ! Degree of freedom per node
     INTEGER(IntKi)                   :: rot_elem ! Rotational degrees of freedom
     INTEGER(IntKi)                   :: nelem ! number of elements
     INTEGER(IntKi)                   :: j ! Index counter
     INTEGER(IntKi)                   :: temp_id ! Index counter
     INTEGER(IntKi)                   :: ErrStat2                     ! Temporary Error status
     CHARACTER(ErrMsgLen)             :: ErrMsg2                      ! Temporary Error message
     CHARACTER(*), PARAMETER          :: RoutineName = 'BD_GenerateDynamicElementForce'
  
     ErrStat    = ErrID_None
  
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
  