anchor.diff:      logical, intent(in) :: calc_k_star
   
       ! local (General)
  -    real(dp), dimension(nTurbines) :: yaw, TIturbs, k_star, wtVelocity
  +    real(dp), dimension(nTurbines) :: yaw, TIturbs, wtVelocity
       real(dp) :: x0, deltax0, deltay, theta_c_0, sigmay, sigmaz, wake_offset, discontinuity_point
  -    real(dp) :: x, deltav, deltaz, sigmay_dp, sigmaz_dp, deltax0_dp, deficit_sum, tol
  +    real(dp) :: x, deltav, deltaz, sigmay_dp, sigmaz_dp, deltax0_dp, deficit_sum, tol, k_star
       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
  
positive.code:         ! print *, "sum of squares = ", sum_of_squares
  !         TI_dst = sqrt(sum_of_squares)
  !         print *, "TI_dst = ", TI_dst
          TI_dst = sqrt(TI_dst**2 + (TI_added*wake_overlap/rotor_area_dst)**2)
          
      
      ! Niayifar and Porte Agel 2015, 2016
  
negative.code:                     
                      if (x > discontinuity_point) then
                      
                          !print *, x
      
                          ! velocity difference in the wake
                          call deltav_func(deltay, deltaz, Ct(turb), yaw(turb), &
  