anchor.diff:     CALL BD_DestroyInput(u_tmp, ErrStat2, ErrMsg2 )
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   
  -   END SUBROUTINE BD_CalcIC
  -
  -   SUBROUTINE BD_CalcForceAcc( u, p, x, OtherState )
  -   !
  -   ! Routine for computing derivatives of continuous states.
  -   !........................................................................................................................
  -
  +END SUBROUTINE BD_CalcIC
  +
  +SUBROUTINE BD_CalcForceAcc( u, p, x, OtherState, ErrStat, ErrMsg)
  +!-------------------------------------------------------------------
  +! Routine for computing derivatives of continuous states.
  +! This subroutine calculates the reaction forces/moments at the root
  +! and the accelerations at all other FE points along the beam
  +!-------------------------------------------------------------------
      TYPE(BD_InputType),           INTENT(IN   ):: u           ! Inputs at t
      TYPE(BD_ParameterType),       INTENT(IN   ):: p           ! Parameters
      TYPE(BD_ContinuousStateType), INTENT(IN   ):: x           ! Continuous states at t
      TYPE(BD_OtherStateType),      INTENT(INOUT):: OtherState  ! Other/optimization states
  -
  -   ! local variables
  -   REAL(ReKi)                                 :: Acc(p%dof_total)
  +   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_CalcForceAcc'
  +
  +   ErrStat = ErrID_None
  +   ErrMsg  = ""
   
      CALL BD_SolutionForceAcc(p%uuN0,x%q,x%dqdt,p%Stif0_GL,p%Mass0_GL,p%gravity,u,&
                               p%damp_flag,p%beta,&
                               p%node_elem,p%dof_node,p%elem_total,p%dof_total,p%node_total,p%ngp,&
  -                            Acc)
  -   OtherState%Acc(:) = Acc(:)
  -
  -   END SUBROUTINE BD_CalcForceAcc
  -
  -   SUBROUTINE BD_SolutionForceAcc(uuN0,uuN,vvN,Stif0,Mass0,gravity,u,                    &
  -                                  damp_flag,beta,                                        &
  -                                  node_elem,dof_node,elem_total,dof_total,node_total,ngp,&
  -                                  Acc)
  -   !***************************************************************************************
  -   ! This subroutine calls other subroutines to apply the force, build the beam element
  -   ! stiffness and mass matrices, build nodal force vector.  The output of this subroutine
  -   ! is the second time derivative of state "q".
  -   !***************************************************************************************
  +                            OtherState%Acc,ErrStat2,ErrMsg2)
  +      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  +
  +   if (ErrStat >= AbortErrLev) RETURN
  +
  +END SUBROUTINE BD_CalcForceAcc
  +
  +SUBROUTINE BD_SolutionForceAcc(uuN0,uuN,vvN,Stif0,Mass0,gravity,u,                    &
  +                               damp_flag,beta,                                        &
  +                               node_elem,dof_node,elem_total,dof_total,node_total,ngp,&
  +                               Acc,ErrStat,ErrMsg)
  +!***************************************************************************************
  +! This subroutine calls other subroutines to apply the force, build the beam element
  +! stiffness and mass matrices, build nodal force vector.  The output of this subroutine
  +! is the second time derivative of state "q".
  +!***************************************************************************************
      REAL(ReKi),                   INTENT(IN   ):: uuN0(:,:) ! Initial position vector
      REAL(ReKi),                   INTENT(IN   ):: Stif0(:,:,:) ! Element stiffness matrix
      REAL(ReKi),                   INTENT(IN   ):: Mass0(:,:,:) ! Element stiffness matrix
  
positive.code:    ENDDO
     OtherState%Xcc(:) = OtherState%Acc(:)
  
     END SUBROUTINE BD_IniAcc
  
  
  END MODULE BeamDyn
  
negative.code:               call Cleanup()
                return
            end if
  !DO j=1,dof_total
  !WRITE(QiMass,6000) MassM(:,j)
  !WRITE(QiStif,6000) StifK(:,j)
  !ENDDO
  !6000 FORMAT (18ES21.12)
  !CLOSE (QiMass)
  !CLOSE (QiStif)
  !STOP
         StifK = MassM + coef(7) * DampG + coef(8) * StifK
         DO j=1,node_total
             temp_id = (j-1)*dof_node
  