anchor.diff:      ! model out
       real(dp), dimension(nSamples), intent(out) :: wsArray
   
  -    intrinsic cos, atan, max, sqrt, log
  +    intrinsic sin, cos, atan, max, sqrt, log
   
       ! bastankhah and porte agel 2016 define yaw to be positive clockwise, this is reversed
       yaw = - yawDeg*pi/180.0_dp
  +    
  +    ! set tolerance for location checks
  +    tol = 0.1_dp
   
       ! initialize location velocities to free stream
       wsArray = wind_speed
  +    
  +    ! initialize wind turbine velocities to 0.0
  +    wtVelocity = 0.0_dp
       
       ! initialize TI of all turbines to free-stream values
       TIturbs = TI
  
positive.code:     ! local (General)
      real(dp), dimension(nTurbines) :: yaw, TIturbs, k_star, wtVelocity
      real(dp) :: x0, deltax0, deltay, theta_c_0, sigmay, sigmaz, wake_offset
      real(dp) :: x, deltav, deltaz, sigmay0, sigmaz0, deficit_sum
      real(dp) :: LocalRotorPointY, LocalRotorPointZ, point_velocity, point_z, point_velocity_with_shear
      real(dp), dimension(nTurbines) :: ky_local, kz_local
      Integer :: u, d, turb, turbI, loc, p
  
negative.code:             ! scale rotor sample point coordinate by rotor diameter (in rotor hub ref. frame)
              LocalRotorPointY = RotorPointsY(p)*0.5_dp*rotorDiameter(turbI)
              LocalRotorPointZ = RotorPointsZ(p)*0.5_dp*rotorDiameter(turbI)
  !             print *, "rotorDiameter after local rotor points", rotorDiameter
  !             print *, "local rotor points Y,Z: ", LocalRotorPointY, LocalRotorPointZ
          
              do, u=1, nTurbines ! at turbineX-locations
              
  