anchor.diff:              ! downstream distance from far wake onset to downstream turbine
               deltax0 = x - x0
               
  -            if (x > 0.0_dp) then
  +            if (x > (0.0_dp + tol)) then
               
                   ! horizontal spread
                   call sigmay_func(ky, deltax0, rotorDiameter(turb), yaw(turb), sigmay)
  
positive.code:                                          & ky_local, kz_local, Ct(turb), sigmay, sigmaz, wake_offset)
                      !print *, "wake_offset ", wake_offset       
                                
                      ! cross wind distance from downstream hub location to wake center
                      deltay = LocalRotorPointY + turbineYw(turbI) - (turbineYw(turb) + wake_offset)
              
                      ! cross wind distance from hub height to height of point of interest
                      deltaz = LocalRotorPointZ + turbineZ(turbI) - turbineZ(turb)
                  
                      ! far wake region
                      if (x >= x0) then
      
  
negative.code:             
              ! find velocity at point p due to the wake of turbine turb
              point_velocity = wind_speed - deficit_sum
              
              !print *, "point velocity, deficit_sum, turbI, p: ", point_velocity, deficit_sum, turbI, p    
          
              ! put sample point height in global reference frame
              point_z = LocalRotorPointZ + turbineZ(turbI)
          
              !print *, "point_z, turbI, p: ", point_z, turbI, p    
              ! adjust sample point velocity for shear
              call wind_shear_func(point_z, point_velocity, z_ref, z_0, shear_exp, point_velocity_with_shear)
              !print *, "v, vs, x, turb, turbI, p: ", point_velocity, point_velocity_with_shear, x, turb, turbI, p
              ! add sample point velocity to turbine velocity to be averaged later
              wtVelocity(turbI) = wtVelocity(turbI) + point_velocity_with_shear
          
  