anchor.diff:  ! 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   ):: 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
  +   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
  +   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
  +   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
  
positive.code:    CHARACTER(ErrMsgLen)         :: ErrMsg2                      ! Temporary Error message
     CHARACTER(*), PARAMETER      :: RoutineName = 'BD_ElememntMatrixForce'
  
     ErrStat = ErrID_None
     ErrMsg  = ""
     elf(:)  = 0.0_BDKi
     
     fact = .FALSE.
  
  
negative.code:              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
         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)
  