anchor.diff:     INTEGER(IntKi),               INTENT(IN   ):: node_total ! Total number of nodes
      INTEGER(IntKi),               INTENT(IN   ):: ngp ! Number of Gauss points
      REAL(ReKi),                   INTENT(  OUT):: Acc(:)
  -
  -   ! Local variables
  -
  -   REAL(ReKi):: MassM(dof_total,dof_total)
  -   REAL(ReKi):: RHS(dof_total)
  -   REAL(ReKi):: F_PointLoad(dof_total)
  -   REAL(ReKi):: sol_temp(dof_total)
  -   REAL(ReKi):: d
  -   INTEGER(IntKi):: indx(dof_total)
  -   INTEGER(IntKi):: j
  -   INTEGER(IntKi):: k
  -   INTEGER(IntKi):: temp_id
  -   REAL(ReKi):: temp6(6)
  -
  -   RHS(:)     = 0.0D0
  -   MassM(:,:) = 0.0D0
  +   INTEGER(IntKi),               INTENT(  OUT):: ErrStat       ! Error status of the operation
  +   CHARACTER(*),                 INTENT(  OUT):: ErrMsg        ! Error message if ErrStat /= ErrID_None
  +
  +   REAL(ReKi),                     ALLOCATABLE:: MassM(:,:)
  +   REAL(ReKi),                     ALLOCATABLE:: RHS(:)
  +   REAL(ReKi),                     ALLOCATABLE:: F_PointLoad(:)
  +   INTEGER(IntKi),                 ALLOCATABLE:: indx(:)
  +   INTEGER(IntKi)                             :: j
  +   INTEGER(IntKi)                             :: k
  +   INTEGER(IntKi)                             :: temp_id
  +   REAL(ReKi)                                 :: temp6(6)
  +   INTEGER(IntKi)                             :: ErrStat2                     ! Temporary Error status
  +   CHARACTER(ErrMsgLen)                       :: ErrMsg2                      ! Temporary Error message
  +   CHARACTER(*), PARAMETER                    :: RoutineName = 'BD_SolutionForceAcc'
  +
  +   ErrStat = ErrID_None
  +   ErrMsg  = ""
  +
  +   CALL AllocAry(MassM,dof_total,dof_total,'Mass Matrix',ErrStat2,ErrMsg2)
  +      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +   CALL AllocAry(RHS,dof_total,'Right-hand-side vector',ErrStat2,ErrMsg2)
  +      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +   CALL AllocAry(F_PointLoad,dof_total,'F_PointLoad vector',ErrStat2,ErrMsg2)
  +      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +   CALL AllocAry(indx,dof_total,'Index vector for LU',ErrStat2,ErrMsg2)
  +      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +   if (ErrStat >= AbortErrLev) then
  +       call Cleanup()
  +       return
  +   end if
   
      CALL BD_GenerateDynamicElementAcc(uuN0,uuN,vvN,Stif0,Mass0,gravity,u,&
                                        damp_flag,beta,&
                                        elem_total,node_elem,dof_total,dof_node,ngp,&
  -                                     RHS,MassM)
  +                                     RHS,MassM,ErrStat2,ErrMsg2)
  +      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      DO j=1,node_total
          temp_id = (j-1)*dof_node
          F_PointLoad(temp_id+1:temp_id+3) = u%PointLoad%Force(1:3,j)
  
positive.code:    INTEGER(IntKi)                                :: j
     INTEGER(IntKi)                                :: k
     INTEGER(IntKi)                                :: piter
     INTEGER(IntKi)                                :: niter
     INTEGER(IntKi)          :: ErrStat2                     ! Temporary Error status
     CHARACTER(ErrMsgLen)    :: ErrMsg2                      ! Temporary Error message
     CHARACTER(*), PARAMETER :: RoutineName = 'BD_Static'
  
     ErrStat = ErrID_None
     ErrMsg  = ""
  
negative.code:    INTEGER(IntKi),     INTENT(IN   ):: node_total ! Total number of nodes
     INTEGER(IntKi),     INTENT(IN   ):: ngp ! Number of Gauss points
     REAL(ReKi),         INTENT(  OUT):: Force(:)
     INTEGER(IntKi),     INTENT(  OUT):: ErrStat       ! Error status of the operation
     CHARACTER(*),       INTENT(  OUT):: ErrMsg        ! Error message if ErrStat /= ErrID_None
  
     INTEGER(IntKi)                   :: ErrStat2                     ! Temporary Error status
     CHARACTER(ErrMsgLen)             :: ErrMsg2                      ! Temporary Error message
     CHARACTER(*), PARAMETER          :: RoutineName = 'BD_DynamicSolutionForce'
  
     ErrStat = ErrID_None
     ErrMsg  = ""
  
  
     CALL BD_GenerateDynamicElementForce(uuN0,uuN,vvN,aaN,     &
                                         Stif0,Mass0,gravity,u,&
                                         damp_flag,beta,&
                                         elem_total,node_elem,dof_node,ngp,Force,&
                                         ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  
     if (ErrStat >= AbortErrLev) RETURN
  
  
     END SUBROUTINE BD_DynamicSolutionForce
  
  