anchor.diff:  
   END SUBROUTINE BD_AssembleRHS
   
  -SUBROUTINE ludcmp(a,n,indx,d)
  -!***************************************************************************************
  -! This subroutine performs LU Decomposition
  -!***************************************************************************************
  -
  -INTEGER(IntKi),INTENT(IN):: n ! DOF total - 6
  -INTEGER(IntKi),INTENT(OUT):: indx(:)
  -REAL(ReKi),INTENT(INOUT):: a(:,:) ! Mass matrix
  -REAL(ReKi),INTENT(OUT):: d ! Determinate
  -
  -REAL(ReKi),PARAMETER:: tolf = 1.0D-20
  -
  -INTEGER(IntKi):: i,imax,j,k
  -REAL(ReKi):: aamax,dum,summ,vv(n)
  -
  -   d=1.d0
  -do i=1,n
  -    aamax=0.
  -    do j=1,n
  -        if (abs(a(i,j)) .gt. aamax) aamax=abs(a(i,j))
  -    enddo
  -    if (aamax==0.) then
  -        write(*,*) "singular matrix"
  -        stop
  -    endif
  -    vv(i)=1./aamax
  -enddo
  -
  -do j=1,n
  -   do i=1,j-1
  -      summ=a(i,j)
  -        do k=1,i-1
  -         summ=summ-a(i,k)*a(k,j)
  -        enddo
  -        a(i,j)=summ
  -    enddo
  -    aamax=0.
  -    do i=j,n
  -        summ=a(i,j)
  -        do k=1,j-1
  -         summ=summ-a(i,k)*a(k,j)
  -      enddo
  -        a(i,j)=summ
  -        dum=vv(i)*abs(summ)
  -        if (dum.ge.aamax) then
  -            imax=i
  -            aamax=dum
  -        endif
  -    enddo
  -    if (j.ne.imax) then
  -        do k=1,n
  -            dum=a(imax,k)
  -            a(imax,k)=a(j,k)
  -            a(j,k)=dum
  -        enddo
  -        d=-d
  -        vv(imax)=vv(j)
  -    endif
  -    indx(j)=imax
  -    if (a(j,j).eq.0.) a(j,j)=tolf
  -    if (j.ne.n) then
  -        dum=1./a(j,j)
  -        do i=j+1,n
  -            a(i,j)=a(i,j)*dum
  -        enddo
  -    endif
  -enddo
  -
  -END SUBROUTINE ludcmp
  -
  -   SUBROUTINE lubksb(a,n,indx,b,ui)
  -!***************************************************************************************
  -! This subroutine solves Ax=b by backward substitution.
  -!***************************************************************************************
  -    INTEGER(IntKi),INTENT(IN):: indx(:)
  -    INTEGER(IntKi),INTENT(IN):: n ! DOF total - 6
  -    REAL(ReKi),INTENT(IN):: a(:,:) ! Mass matrix
  -    REAL(ReKi),INTENT(INOUT):: b(:) ! Unknowns in Ax=b
  -    REAL(ReKi),INTENT(OUT):: ui(:) ! Unknowns in Ax=b
  -
  -    INTEGER(IntKi):: i,ii,j,ll
  -    REAL(ReKi):: summ
  -
  -    ii = 0
  -    DO i = 1 , n
  -      ll = indx(i)
  -        summ = b(ll)
  -        b(ll) = b(i)
  -        IF (ii .ne. 0) THEN
  -         DO j = ii , i-1
  -               summ = summ - a(i,j) * b(j)
  -         ENDDO
  -        ELSE IF (summ .ne. 0.0) THEN
  -         ii = i
  -        ENDIF
  -        b(i) = summ
  -   ENDDO
  -
  -    DO i = n,1,-1
  -      summ = b(i)
  -        DO j = i+1 , n
  -         summ = summ - a(i,j) * b(j)
  -        ENDDO
  -        b(i) = summ / a(i,i)
  -    ENDDO
  -
  -    do i=1,n
  -      ui(i)=b(i)
  -    enddo
  -
  -   END SUBROUTINE lubksb
  -
   SUBROUTINE BD_UpdateDynamicGA2(ainc,uf,vf,af,xf,coef,node_total,dof_node,ErrStat,ErrMsg)
   !---------------------------------------------------------------------------
   ! This subroutine updates the 1) displacements/rotations(uf)
  
positive.code:        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,ErrStat2,ErrMsg2)
         CALL BD_ElasticForce(E1,RR0,kapa,Stif,cet,Fc,Fd,Oe,Pe,Qe,ErrStat2,ErrMsg2)
         IF(damp_flag .EQ. 1) THEN
             CALL BD_DissipativeForce(beta,Stif,vvv,vvp,E1,Fc,Fd,Sd,Od,Pd,Qd,&
                                      betaC,Gd,Xd,Yd,ErrStat2,ErrMsg2)
         ENDIF
         CALL BD_GyroForce(mEta,rho,uuu,vvv,Fb)
  
         DO i=1,node_elem
             DO j=1,dof_node
  
negative.code:        ENDDO
     ENDDO
     temp33 = 0.0D0
     temp33 = MATMUL(BD_Tilde(ome),rho)
     DO i=1,3
         DO j=1,3
  