anchor.diff:     integer(IntKi),          intent(  out)  :: ErrStat                         !< Error status of the operation
      character(*),            intent(  out)  :: ErrMsg                          !< Error message if ErrStat /= ErrID_None
   
  +!   real(ReKi)                              :: x_hat_disk(3)
  +!   real(ReKi)                              :: y_hat_disk(3)
  +!   real(ReKi)                              :: z_hat_disk(3)
      integer(intKi)                          :: tIndx
      integer(intKi)                          :: k                      ! loop counter for blades
      integer(intKi)                          :: ErrStat2
  
positive.code:    ! local variables
     real(ReKi)                              :: x_hat(3)
     real(ReKi)                              :: y_hat(3)
     real(ReKi)                              :: z_hat(3)
     real(ReKi)                              :: x_hat_disk(3)
     real(ReKi)                              :: y_hat_disk(3)
     real(ReKi)                              :: z_hat_disk(3)
     real(ReKi)                              :: tmp(3)
     real(R8Ki)                              :: theta(3)
     real(R8Ki)                              :: orientation(3,3)
     real(R8Ki)                              :: orientation_nopitch(3,3)
     real(ReKi)                              :: tmp_sz, tmp_sz_y
     
     integer(intKi)                          :: j                      ! loop counter for nodes
     integer(intKi)                          :: k                      ! loop counter for blades
  
negative.code:    character(*), parameter                 :: RoutineName = 'SetInputsForBEMT'
     
     
        ! Get disk average values and orientations
     call DiskAvgValues(p, u, m, x_hat_disk, y_hat_disk, z_hat_disk, Azimuth)
     call GeomWithoutSweepPitchTwist(p,u,m,thetaBladeNds,ErrStat,ErrMsg)
     if (ErrStat >= AbortErrLev) return
  
        ! Velocity in disk normal
     m%BEMT_u(indx)%Un_disk  = m%V_dot_x
       
        ! "Angular velocity of rotor" rad/s
     m%BEMT_u(indx)%omega   = dot_product( u%HubMotion%RotationVel(:,1), x_hat_disk )    
     
        ! "Angle between the vector normal to the rotor plane and the wind vector (e.g., the yaw angle in the case of no tilt)" rad 
     tmp_sz = TwoNorm( m%V_diskAvg )
     if ( EqualRealNos( tmp_sz, 0.0_ReKi ) ) then
        m%BEMT_u(indx)%chi0 = 0.0_ReKi
     else
           ! make sure we don't have numerical issues that make the ratio outside +/-1
        tmp_sz_y = min(  1.0_ReKi, m%V_dot_x / tmp_sz )
        tmp_sz_y = max( -1.0_ReKi, tmp_sz_y )
        
        m%BEMT_u(indx)%chi0 = acos( tmp_sz_y )
        
     end if
     
        ! "Azimuth angle" rad
     m%bemt_u(indx)%psi = Azimuth
  
        ! theta, "Twist angle (includes all sources of twist)" rad
        ! Vx, "Local axial velocity at node" m/s
        ! Vy, "Local tangential velocity at node" m/s
     do k=1,p%NumBlades
        do j=1,p%NumBlNds         
           
           m%BEMT_u(indx)%theta(j,k) = thetaBladeNds(j,k) ! local pitch + twist (aerodyanmic + elastic) angle of the jth node in the kth blade
                             
           x_hat = m%WithoutSweepPitchTwist(1,:,j,k)
           y_hat = m%WithoutSweepPitchTwist(2,:,j,k)
           tmp   = m%DisturbedInflow(:,j,k) - u%BladeMotion(k)%TranslationVel(:,j) ! rel_V(j)_Blade(k)
           
           m%BEMT_u(indx)%Vx(j,k) = dot_product( tmp, x_hat ) ! normal component (normal to the plane, not chord) of the inflow velocity of the jth node in the kth blade
           m%BEMT_u(indx)%Vy(j,k) = dot_product( tmp, y_hat ) ! tangential component (tangential to the plane, not chord) of the inflow velocity of the jth node in the kth blade
           
        end do !j=nodes
     end do !k=blades
     
     
        ! "Radial distance from center-of-rotation to node" m
     
     do k=1,p%NumBlades
        do j=1,p%NumBlNds
           
              ! displaced position of the jth node in the kth blade relative to the hub:
           tmp =  u%BladeMotion(k)%Position(:,j) + u%BladeMotion(k)%TranslationDisp(:,j) &
                - u%HubMotion%Position(:,1)      - u%HubMotion%TranslationDisp(:,1)
           
              ! local radius (normalized distance from rotor centerline)
           tmp_sz_y = dot_product( tmp, y_hat_disk )**2
           tmp_sz   = dot_product( tmp, z_hat_disk )**2
           m%BEMT_u(indx)%rLocal(j,k) = sqrt( tmp_sz + tmp_sz_y )
           
        end do !j=nodes      
     end do !k=blades  
    
     m%BEMT_u(indx)%UserProp = u%UserProp
     
  end subroutine SetInputsForBEMT
  
  subroutine DiskAvgValues(p, u, m, x_hat_disk, y_hat_disk, z_hat_disk, Azimuth)
     type(AD_ParameterType),  intent(in   )  :: p                               !< AD parameters
     type(AD_InputType),      intent(in   )  :: u                               !< AD Inputs at Time
     type(AD_MiscVarType),    intent(inout)  :: m                               !< Misc/optimization variables
     real(ReKi),              intent(  out)  :: x_hat_disk(3)
     real(ReKi),              intent(  out)  :: y_hat_disk(3)
     real(ReKi),              intent(  out)  :: z_hat_disk(3)
     real(R8Ki),              intent(  out)  :: Azimuth(p%NumBlNds)
     real(ReKi)                              :: z_hat(3)
     real(ReKi)                              :: tmp(3)
     real(ReKi)                              :: tmp_sz, tmp_sz_y
     integer(intKi)                          :: j                      ! loop counter for nodes
     integer(intKi)                          :: k                      ! loop counter for blades
  
        ! calculate disk-averaged relative wind speed, V_DiskAvg
     m%V_diskAvg = 0.0_ReKi
     do k=1,p%NumBlades
        do j=1,p%NumBlNds
           tmp = m%DisturbedInflow(:,j,k) - u%BladeMotion(k)%TranslationVel(:,j)
           m%V_diskAvg = m%V_diskAvg + tmp
        end do
     end do
     m%V_diskAvg = m%V_diskAvg / real( p%NumBlades * p%NumBlNds, ReKi )
  
        ! orientation vectors:
     x_hat_disk = u%HubMotion%Orientation(1,:,1) !actually also x_hat_hub
  
     m%V_dot_x  = dot_product( m%V_diskAvg, x_hat_disk )
     tmp    = m%V_dot_x * x_hat_disk - m%V_diskAvg
     tmp_sz = TwoNorm(tmp)
     if ( EqualRealNos( tmp_sz, 0.0_ReKi ) ) then
  