anchor.diff:                  wsArray(loc) = wsArray(loc) - deltav
   
               end if
  +            
  +            ! make sure turbine flow-field velocity is non-negative
  +            if (wsArray(loc) .lt. 0.0_dp) then 
  +                wsArray(loc) = 0.0_dp
  +            end if
  +            
           end do
       end do
   
  
positive.code: 
      ! in
      integer, intent(in) :: nTurbines
      real(dp), dimension(nTurbines), intent(in) :: turbineXw, turbineYw, turbineZ
      real(dp), dimension(nTurbines), intent(in) :: rotorDiameter, yawDeg
      real(dp), dimension(nTurbines), intent(in) :: Ct
  
negative.code: ! implementation of the Bastankhah and Porte Agel (BPA) wake model for analysis
  subroutine porteagel_analyze(nTurbines, turbineXw, turbineYw, turbineZ, &
                               rotorDiameter, Ct, wind_speed, &
                               yawDeg, ky, kz, alpha, beta, I, wake_combination_method, &
                               TI_calculation_method, wtVelocity)
  
      ! independent variables: turbineXw turbineYw turbineZ rotorDiameter
      !                        Ct yawDeg
  