anchor.diff:     DO piter=1,p%niter
   
         ! compute the finite differenced stiffness matrix
  -      IF ( p%tngt_stf_fd .or. p%tngt_stf_comp ) CALL BD_FD_Stat( x, gravity, p, m, ErrStat2, ErrMsg2 )
  +      IF ( p%tngt_stf_fd .or. p%tngt_stf_comp ) CALL BD_FD_Stat( x, gravity, p, m )
   
         CALL BD_QuadraturePointData( p,x,m )         ! Calculate QP values uuu, uup, RR0, kappa, E1
         CALL BD_GenerateStaticElement(gravity, p, m) ! Calculate RHS and analytical tangent stiffness matrix
  
positive.code: 
  !-----------------------------------------------------------------------------------------------------------------------------------
  !> This subroutine computes the finite differenced tangent stiffness matrix
  SUBROUTINE BD_FD_Stat( x, gravity, p, m, ErrStat, ErrMsg )
  
      ! Function arguments
      TYPE(BD_ContinuousStateType),    INTENT(INOUT) :: x            !< Continuous states at t on input at t + dt on output
      REAL(BDKi),                      INTENT(IN   ) :: gravity(:)   !< not the same as p%gravity (used for ramp of loads and gravity)
      TYPE(BD_ParameterType),          INTENT(IN   ) :: p            !< Parameters
      TYPE(BD_MiscVarType),            INTENT(INOUT) :: m            !< misc/optimization variables
      INTEGER(IntKi),                  INTENT(  OUT)  :: ErrStat     !< Error status of the operation
      CHARACTER(*),                    INTENT(  OUT)  :: ErrMsg      !< Error message if ErrStat /= ErrID_None
  
      ! local variables
      INTEGER(IntKi)                                 :: i
      INTEGER(IntKi)                                 :: idx_dof
      INTEGER(IntKi)                                 :: ErrStat2 ! Temporary Error status
      REAL(BDKi), allocatable                        :: RHS_m(:,:), RHS_p(:,:)
      CHARACTER(ErrMsgLen)                           :: ErrMsg2  ! Temporary Error message
      CHARACTER(*), PARAMETER                        :: RoutineName = 'BD_FD_Stat'
  
      ErrStat = ErrID_None
      ErrMsg  = ""
  
      ! zero out the local matrices.
      m%RHS_m    = 0.0_BDKi
  
negative.code: 
        fact = MOD(piter-1,p%n_fact) .EQ. 0  ! when true, we factor the jacobian matrix
  
        IF ( (p%tngt_stf_fd .OR. p%tngt_stf_comp) .AND. fact ) CALL BD_FD_GA2( x, OtherState, p, m )
  
           ! Apply accelerations using F=ma ?  Is that what this routine does?
           ! Calculate Quadrature point values needed
  