anchor.diff:              ! 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               
  
positive.code:          InitInData_BD%RootDisp     = 0.0_ReKi                                              ! {:} - - "Initial root displacement"
           InitInData_BD%RootOri      = InitInData_BD%GlbRot                                  ! {:}{:} - - "Initial root orientation"
           InitInData_BD%RootVel      = 0.0_ReKi                                              ! {:} - - "Initial root velocities and angular veolcities"                  
           
           CALL BD_Init( InitInData_BD, BD%Input(1,k), BD%p(k),  BD%x(k,STATE_CURR), BD%xd(k,STATE_CURR), BD%z(k,STATE_CURR), &
                              BD%OtherSt(k), BD%y(k), dt_BD, InitOutData_BD(k), ErrStat2, ErrMsg2 )
  
negative.code:                
                 
                 DO j=1,p_FAST%SizeJac_ED_SD_HD_BD(4+nb) !call BD_CalcOutput            
                    i = i + 1 ! i = j + sum(p_FAST%SizeJac_ED_SD_HD_BD(2:3+nb))
                    ! perturb u_BD:
                    CALL BD_CopyInput(  u_BD(nb), u_BD_perturb, MESH_UPDATECOPY, ErrStat2, ErrMsg2 )
                       CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName  )
  