anchor.diff:              
                   ! get index of upstream turbine
                   turb = sorted_x_idx(u) + 1
  +                
  +                ! skip this loop if turb = turbI (turbines impact on itself)
  +                if (turb .eq. turbI) cycle
           
                   ! downstream distance between turbines, adjust for downstream turbine yaw
                   x = turbineXw(turbI) - turbineXw(turb) + LocalRotorPointY*sin(yaw(turbI))
  
positive.code:                 ! vertical distance from downstream hub location to wake center
                  deltaz = velZ(loc) - turbineZ(turb)
                  !print *, "here 1"
                  ! far wake region
                  if (x >= x0) then
  
                      ! velocity difference in the wake
                      call deltav_func(deltay, deltaz, Ct(turb), yaw(turb),  &
  
negative.code: 
                      ! near wake region (linearized)
                      else
                          
                          ! determine distance from discontinuity point to far wake onset
                          deltax0_dp = discontinuity_point - x0
                  
                          ! horizontal spread at discontinuity point
                          call sigmay_func(ky_local, deltax0_dp, rotorDiameter(turb), yaw(turb), sigmay_dp)
                  
                          ! vertical spread at discontinuity point
                          call sigmaz_func(kz_local, deltax0_dp, rotorDiameter(turb), sigmaz_dp)
                  
                          ! velocity deficit in the nearwake (linear model)
                          call deltav_near_wake_lin_func(deltay, deltaz, &
                                           & Ct(turb), yaw(turb), sigmay_dp, sigmaz_dp, & 
                                           & rotorDiameter(turb), x, discontinuity_point, sigmay_dp, sigmaz_dp, & 
                                           & deltav)
                
                      end if
  