anchor.diff:          ! Calculate the turbulence added to the inflow of the downstream turbine by the 
           ! wake of the upstream turbine
           TI_added = 0.73_dp*(axial_induction_ust**0.8325)*(TI_ust**0.0325)* & 
  -                    ((x/rotor_diameter_ust)**-0.32)
  +                    ((x/rotor_diameter_ust)**(-0.32_dp))
           !print *, "TI_added = ", TI_added
  -        rotor_area_dst = 0.25_dp*pi*rotor_diameter_dst**2
  +        rotor_area_dst = 0.25_dp*pi*rotor_diameter_dst**2_dp
           ! Calculate the total turbulence intensity at the downstream turbine
           !sum_of_squares = TI_dst**2 + (TI_added*wake_overlap)**2
           ! print *, "sum of squares = ", sum_of_squares
  
positive.code:     
      end do
      
      
      !print *, wtVelocity
  
      !! make sure turbine inflow velocity is non-negative
  
negative.code:     ! local (General)
      real(dp), dimension(nTurbines) :: yaw, TIturbs, k_star
      real(dp) :: x0, deltax0, deltay, theta_c_0, sigmay, sigmaz, wake_offset
      real(dp) :: x, deltav, deltaz, sigmay_dp, sigmaz_dp, deltax0_dp, deficit_sum
      real(dp) :: ky_local, kz_local, tol, discontinuity_point
      real(dp) :: LocalRotorPointY, LocalRotorPointZ, point_velocity, point_z, point_velocity_with_shear
      Integer :: u, d, turb, turbI, p
  