anchor.diff:     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)
  +   REAL(ReKi)                   :: temp_aaa(6)
      INTEGER(IntKi)               :: igp
      INTEGER(IntKi)               :: i
      INTEGER(IntKi)               :: j
  
positive.code:    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)
  
negative.code:    CHARACTER(ErrMsgLen)         :: ErrMsg2                      ! Temporary Error message
     CHARACTER(*), PARAMETER      :: RoutineName = 'BD_ElememntMatrixForce'
  
     ErrStat = ErrID_None
     ErrMsg  = ""
     elf(:)  = 0.0D0
  
     CALL AllocAry(gp,ngp,'Gauss piont array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry(gw,ngp,'Gauss piont weight function array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry(hhx,node_elem,'Shape function array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry(hpx,node_elem,'Derivative of shape function array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry(GLL_temp,node_elem,'Gauss-Lobatto-Legendre (GLL) point array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry(w_temp,node_elem,'GLL weight function array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL AllocAry(temp_Naaa,dof_node*node_elem,'Temporary elemental acceleration array',ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     if (ErrStat >= AbortErrLev) then
        call Cleanup()
        return
     end if
     temp_Naaa(:)  = 0.0D0
  
     CALL BD_GenerateGLL(ngp,GLL_temp,w_temp,ErrStat2,ErrMsg2)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
     CALL BD_GaussPointWeight(ngp,gp,gw,ErrStat2,ErrMsg2)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
  
     DO igp=1,ngp
         gpr=gp(igp)
         CALL BD_ComputeJacobian(gpr,Nuu0,node_elem,dof_node,gp,GLL_temp,ngp,igp,hhx,hpx,Jacobian,ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         CALL BD_GaussPointDataAt0(hhx,hpx,Nuu0,Nrr0,node_elem,dof_node,uu0,E10,ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         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,ErrStat2,ErrMsg2)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         mmm  = 0.0D0
         mEta = 0.0D0
         rho  = 0.0D0
  