anchor.diff:     ENDIF
   
      END SUBROUTINE BD_MotionTensor
  + 
  +   SUBROUTINE BD_GenerateDynamicElementAcc(uuN0,uuN,vvN,Stif0,Mass0,gravity,u,         &
  +                                           damp_flag,beta,                             &
  +                                           elem_total,node_elem,dof_total,dof_node,ngp,&
  +                                           RHS,MassM)
  +   !----------------------------------------------------------------------------------------
  +   ! This subroutine computes Global mass matrix and force vector for the beam.
  +   !----------------------------------------------------------------------------------------
  +   REAL(ReKi),        INTENT(IN   ):: uuN0(:,:) ! Initial position vector
  +   REAL(ReKi),        INTENT(IN   ):: uuN(:) ! Displacement of Mass 1: m
  +   REAL(ReKi),        INTENT(IN   ):: vvN(:) ! Velocity of Mass 1: m/s
  +   REAL(ReKi),        INTENT(IN   ):: Stif0(:,:,:) ! Element stiffness matrix
  +   REAL(ReKi),        INTENT(IN   ):: Mass0(:,:,:) ! Element stiffness matrix
  +   REAL(ReKi),        INTENT(IN   ):: gravity(:) ! Velocity of Mass 1: m/s
  +   TYPE(BD_InputType),INTENT(IN   ):: u           ! Inputs at t
  +   INTEGER(IntKi),    INTENT(IN   ):: damp_flag ! Total number of elements
  +   REAL(ReKi),        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_total ! Degrees of freedom per node
  +   INTEGER(IntKi),    INTENT(IN   ):: dof_node ! Degrees of freedom per node
  +   INTEGER(IntKi),    INTENT(IN   ):: ngp ! Number of Gauss points
  +   REAL(ReKi),        INTENT(  OUT):: MassM(:,:) ! Mass matrix
  +   REAL(ReKi),        INTENT(  OUT):: RHS(:) ! Right hand side of the equation Ax=B
  +
  +   REAL(ReKi)                      :: Nuu0(dof_node*node_elem) ! Nodal initial position for each element
  +   REAL(ReKi)                      :: Nuuu(dof_node*node_elem) ! Nodal displacement of Mass 1 for each element
  +   REAL(ReKi)                      :: Nrr0(3*node_elem) ! Nodal rotation parameters for initial position
  +   REAL(ReKi)                      :: Nrrr(3*node_elem) ! Nodal rotation parameters for displacement of Mass 1
  +   REAL(ReKi)                      :: Nvvv(dof_node*node_elem) ! Nodal velocity of Mass 1: m/s for each element
  +   REAL(ReKi)                      :: EStif0_GL(6,6,node_elem-1) ! Nodal material properties for each element
  +   REAL(ReKi)                      :: EMass0_GL(6,6,node_elem-1) ! Nodal material properties for each element
  +   REAL(ReKi)                      :: DistrLoad_GL(6,node_elem-1) ! Nodal material properties for each element
  +   REAL(ReKi)                      :: elf(dof_node*node_elem) ! Total element force (Fc, Fd, Fb)
  +   REAL(ReKi)                      :: elm(dof_node*node_elem,dof_node*node_elem) ! Element mass matrix
  +   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_GenerateDynamicElementAcc'
  +
  +   dof_elem = dof_node * node_elem
  +   rot_elem = (dof_node/2) * node_elem
  +   DO nelem=1,elem_total
  +       Nuu0(:) = uuN0(:,nelem)
  +       CALL BD_ElemNodalDisp(uuN,node_elem,dof_node,nelem,Nuuu,ErrStat2,ErrMsg2)
  +       temp_id = (nelem-1)*ngp
  +       DO j=1,ngp
  +           EStif0_GL(1:6,1:6,j) = Stif0(1:6,1:6,temp_id+j)
  +           EMass0_GL(1:6,1:6,j) = Mass0(1:6,1:6,temp_id+j)
  +           DistrLoad_GL(1:3,j) = u%DistrLoad%Force(1:3,temp_id+j+1)
  +           DistrLoad_GL(4:6,j) = u%DistrLoad%Moment(1:3,temp_id+j+1)
  +       ENDDO
  +
  +       CALL BD_NodalRelRot(Nuu0,node_elem,dof_node,Nrr0)
  +       CALL BD_NodalRelRot(Nuuu,node_elem,dof_node,Nrrr)
  +       CALL BD_ElemNodalDisp(vvN,node_elem,dof_node,nelem,Nvvv,ErrStat2,ErrMsg2)
  +
  +       elf(:) = 0.0D0
  +       elm(:,:) = 0.0D0
  +       CALL BD_ElementMatrixAcc(Nuu0,Nuuu,Nrr0,Nrrr,Nvvv,&
  +                                EStif0_GL,EMass0_GL,gravity,DistrLoad_GL,&
  +                                ngp,node_elem,dof_node,damp_flag,beta,&
  +                                elf,elm)
  +
  +
  +       CALL BD_AssembleStiffK(nelem,node_elem,dof_elem,dof_node,&
  +                              elm,MassM)
  +       CALL BD_AssembleRHS(nelem,dof_elem,node_elem,dof_node,elf,RHS)
  +
  +   ENDDO
  +
  +   END SUBROUTINE BD_GenerateDynamicElementAcc
  +
  +   SUBROUTINE BD_ElementMatrixAcc(Nuu0,Nuuu,Nrr0,Nrrr,Nvvv,&
  +                                  EStif0_GL,EMass0_GL,gravity,DistrLoad_GL,&
  +                                  ngp,node_elem,dof_node,damp_flag,beta,&
  +                                  elf,elm)
  +
  +   !-------------------------------------------------------------------------------
  +   ! This subroutine total element forces and mass matrices
  +   !-------------------------------------------------------------------------------
  +
  +   REAL(ReKi),INTENT(IN   )    :: Nuu0(:) ! Nodal initial position for each element
  +   REAL(ReKi),INTENT(IN   )    :: Nuuu(:) ! Nodal displacement of Mass 1 for each element
  +   REAL(ReKi),INTENT(IN   )    :: Nrr0(:) ! Nodal rotation parameters for initial position
  +   REAL(ReKi),INTENT(IN   )    :: Nrrr(:) ! Nodal rotation parameters for displacement of Mass 1
  +   REAL(ReKi),INTENT(IN   )    :: Nvvv(:) ! Nodal velocity of Mass 1: m/s for each element
  +   REAL(ReKi),INTENT(IN   )    :: EStif0_GL(:,:,:) ! Nodal material properties for each element
  +   REAL(ReKi),INTENT(IN   )    :: EMass0_GL(:,:,:) ! Nodal material properties for each element
  +   REAL(ReKi),INTENT(IN   )    :: gravity(:) !
  +   REAL(ReKi),INTENT(IN   )    :: DistrLoad_GL(:,:) ! Nodal material properties for each element
  +   REAL(ReKi),INTENT(  OUT)    :: elf(:)  ! Total element force (Fd, Fc, Fb)
  +   REAL(ReKi),INTENT(  OUT)    :: elm(:,:) ! Total element mass matrix
  +   INTEGER(IntKi),INTENT(IN   ):: ngp ! Number of Gauss points
  +   INTEGER(IntKi),INTENT(IN   ):: node_elem ! Node per element
  +   INTEGER(IntKi),INTENT(IN   ):: dof_node ! Degrees of freedom per node
  +   INTEGER(IntKi),INTENT(IN   ):: damp_flag ! Degrees of freedom per node
  +   REAL(ReKi),INTENT(IN   )    :: beta(:)
  +
  +   REAL(ReKi)                  :: gp(ngp) ! Gauss points
  +   REAL(ReKi)                  :: gw(ngp) ! Gauss point weights
  +   REAL(ReKi)                  :: hhx(node_elem) ! Shape function
  +   REAL(ReKi)                  :: hpx(node_elem) ! Derivative of shape function
  +   REAL(ReKi)                  :: GLL_temp(node_elem) ! Temp Gauss-Lobatto-Legendre points
  +   REAL(ReKi)                  :: w_temp(node_elem) ! Temp GLL weights
  +   REAL(ReKi)                  :: uu0(6)
  +   REAL(ReKi)                  :: E10(3)
  +   REAL(ReKi)                  :: RR0(3,3)
  +   REAL(ReKi)                  :: kapa(3)
  +   REAL(ReKi)                  :: E1(3)
  +   REAL(ReKi)                  :: Stif(6,6)
  +   REAL(ReKi)                  :: cet
  +   REAL(ReKi)                  :: uuu(6)
  +   REAL(ReKi)                  :: uup(3)
  +   REAL(ReKi)                  :: Jacobian
  +   REAL(ReKi)                  :: gpr
  +   REAL(ReKi)                  :: Fc(6)
  +   REAL(ReKi)                  :: Fd(6)
  +   REAL(ReKi)                  :: Fg(6)
  +   REAL(ReKi)                  :: vvv(6)
  +   REAL(ReKi)                  :: vvp(6)
  +   REAL(ReKi)                  :: mmm
  +   REAL(ReKi)                  :: mEta(3)
  +   REAL(ReKi)                  :: rho(3,3)
  +   REAL(ReKi)                  :: Fb(6)
  +   REAL(ReKi)                  :: Mi(6,6)
  +   REAL(ReKi)                  :: Oe(6,6)
  +   REAL(ReKi)                  :: Pe(6,6)
  +   REAL(ReKi)                  :: Qe(6,6)
  +   REAL(ReKi)                  :: Sd(6,6)
  +   REAL(ReKi)                  :: Od(6,6)
  +   REAL(ReKi)                  :: Pd(6,6)
  +   REAL(ReKi)                  :: Qd(6,6)
  +   REAL(ReKi)                  :: betaC(6,6)
  +   REAL(ReKi)                  :: Gd(6,6)
  +   REAL(ReKi)                  :: Xd(6,6)
  +   REAL(ReKi)                  :: Yd(6,6)
  +   REAL(ReKi)                  :: temp_Naaa(dof_node*node_elem)
  +   REAL(ReKi)                  :: temp_aaa(dof_node)
  +
  +   INTEGER(IntKi)              :: igp
  +   INTEGER(IntKi)              :: i
  +   INTEGER(IntKi)              :: j
  +   INTEGER(IntKi)              :: m
  +   INTEGER(IntKi)              :: n
  +   INTEGER(IntKi)              :: temp_id1
  +   INTEGER(IntKi)              :: temp_id2
  +   INTEGER(IntKi)              :: ErrStat2                     ! Temporary Error status
  +   CHARACTER(ErrMsgLen)        :: ErrMsg2                      ! Temporary Error message
  +   CHARACTER(*), PARAMETER     :: RoutineName = 'BD_ElementMatrixAcc'
  +
  +
  +   elf(:) = 0.0D0
  +   elm(:,:) = 0.0D0
  +
  +   CALL BD_GenerateGLL(ngp,GLL_temp,w_temp,ErrStat2,ErrMsg2)
  +   CALL BD_GaussPointWeight(ngp,gp,gw,ErrStat2,ErrMsg2)
  +
  +   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
  +       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)
  +       CALL BD_GaussPointDataMass(hhx,hpx,Nvvv,temp_Naaa,RR0,node_elem,dof_node,vvv,temp_aaa,vvp,mmm,mEta,rho)
  +       CALL BD_MassMatrix(mmm,mEta,rho,Mi)
  +
  +       CALL BD_GyroForce(mEta,rho,uuu,vvv,Fb)
  +       CALL BD_GravityForce(mmm,mEta,gravity,Fg)
  +       CALL BD_ElasticForce(E1,RR0,kapa,Stif,cet,Fc,Fd,Oe,Pe,Qe)
  +       IF(damp_flag .NE. 0) THEN
  +           CALL BD_DissipativeForce(beta,Stif,vvv,vvp,E1,Fc,Fd,Sd,Od,Pd,Qd,betaC,Gd,Xd,Yd)
  +       ENDIF
  +
  +       Fd(:) = Fd(:) + Fb(:) - DistrLoad_GL(:,igp) - Fg(:)
  +
  +       DO i=1,node_elem
  +           DO j=1,node_elem
  +               DO m=1,dof_node
  +                   temp_id1 = (i-1)*dof_node+m
  +                   DO n=1,dof_node
  +                       temp_id2 = (j-1)*dof_node+n
  +                       elm(temp_id1,temp_id2) = elm(temp_id1,temp_id2) + hhx(i)*Mi(m,n)*hhx(j)*Jacobian*gw(igp)
  +                   ENDDO
  +               ENDDO
  +           ENDDO
  +       ENDDO
  +
  +       DO i=1,node_elem
  +           DO j=1,dof_node
  +               temp_id1 = (i-1) * dof_node+j
  +               elf(temp_id1) = elf(temp_id1) - hpx(i)*Fc(j)*Jacobian*gw(igp)
  +               elf(temp_id1) = elf(temp_id1) - hhx(i)*Fd(j)*Jacobian*gw(igp)
  +           ENDDO
  +       ENDDO
  +   ENDDO
  +
  +
  +   END SUBROUTINE BD_ElementMatrixAcc
   
      SUBROUTINE BD_MassMatrix(m00,mEta,rho,Mi)
   !----------------------------------------------------------------------------------------
  
positive.code:    REAL(ReKi),          ALLOCATABLE:: elf(:)
     REAL(ReKi),          ALLOCATABLE:: elm(:,:)
     REAL(ReKi),          ALLOCATABLE:: elg(:,:)
     REAL(ReKi)           ALLOCATABLE:: EStif0_GL(:,:,:)
     REAL(ReKi)           ALLOCATABLE:: EMass0_GL(:,:,:)
     REAL(ReKi)           ALLOCATABLE:: DistrLoad_GL(:,:)
     INTEGER(IntKi)                  :: dof_elem
     INTEGER(IntKi)                  :: rot_elem
     INTEGER(IntKi)                  :: nelem
  
negative.code:    elk(:,:) = 0.0D0
     elm(:,:) = 0.0D0
     elg(:,:) = 0.0D0
     EStif0_GL(:,:,:)  = 0.0D0
     EMass0_GL(:,:,:)  = 0.0D0
     DistrLoad_GL(:,:) = 0.0D0
  
     DO nelem=1,elem_total
         Nuu0(:) = uuN0(:,nelem)
  