anchor.diff:  REAL(ReKi) FUNCTION GetPerturb(x)
      REAL(ReKi), INTENT(IN) :: x
         
  -   !GetPerturb = sqrt( EPSILON(x)) * max( abs(x, 1._ReKi)  
  +   !GetPerturb = sqrt( EPSILON(x)) * max( abs(x), 1._ReKi)  
   !      GetPerturb = 1.0e6
      GetPerturb = 1.0
         
  
positive.code:                   ! calculate outputs with perturbed inputs:
                    CALL BD_CalcOutput( this_time, u_BD_perturb, p_BD(nb), x_BD(nb), xd_BD(nb), z_BD(nb), OtherSt_BD(nb), y_BD_perturb(nb), ErrStat2, ErrMsg2 )
                       CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName  )
                    
                    CALL U_ED_SD_HD_BD_Residual(y_ED, y_SD, y_HD, y_BD_perturb, u_perturb, Fn_U_perturb) ! get this perturbation  
                 
                    IF ( ErrStat >= AbortErrLev ) THEN
  
negative.code:             ! Get BeamDyn's contribution: (note if p_FAST%CompElast /= Module_BD, SizeJac_ED_SD_HD_BD(5) = 0)
              !............................... 
              DO nb=1,p_FAST%nBeams
                 CALL BD_CopyOutput(y_BD(nb),y_BD_perturb(nb),MESH_UPDATECOPY,ErrStat2,ErrMsg2)
                    CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName  )
              END DO
                          
              DO nb=1,p_FAST%nBeams
                 
                 ! make sure we perturb the outputs from only the current blade (overwrite the previous perturbation with y_BD values)
                 if (nb > 1) then               
  