anchor.diff:                      ! calculate wake spreading parameter at each turbine if desired
                       if (calc_k_star .eqv. .true.) then
                           call k_star_func(TIturbs(turb), k_star(turb))
  -                        ky_local = k_star(turb)
  -                        kz_local = k_star(turb)
  +                        ky_local(turb) = k_star(turb)
  +                        kz_local(turb) = k_star(turb)
                       end if
   
                       ! determine the initial wake angle at the onset of far wake
                       call theta_c_0_func(yaw(turb), Ct(turb), theta_c_0)
           
                       ! horizontal spread
  -                    call sigmay_func(ky_local, deltax0, rotorDiameter(turb), yaw(turb), sigmay)
  +                    call sigmay_func(ky_local(turb), deltax0, rotorDiameter(turb), yaw(turb), sigmay)
                       !print *, "sigmay ", sigmay
                       ! vertical spread
  -                    call sigmaz_func(kz_local, deltax0, rotorDiameter(turb), sigmaz)
  +                    call sigmaz_func(kz_local(turb), deltax0, rotorDiameter(turb), sigmaz)
                       !print *, "sigmaz ", sigmaz
                       ! horizontal cross-wind wake displacement from hub
                       call wake_offset_func(rotorDiameter(turb), theta_c_0, x0, yaw(turb), &
  -                                         & ky_local, kz_local, Ct(turb), sigmay, sigmaz, wake_offset)
  +                                         & ky_local(turb), kz_local(turb), Ct(turb), sigmay, sigmaz, wake_offset)
                       !print *, "wake_offset ", wake_offset       
                                 
                       ! cross wind distance from downstream point location to wake center
  
positive.code:                 !print *, "here 1"
                  
                  ! find the final point where the original model is undefined
                  call discontinuity_point_func(x0, rotorDiameter(turb), ky, kz, yaw, Ct, & 
                                               & discontinuity_point)
                                               
                  ! far wake region
  
negative.code:         ! 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_dp))
          !print *, "TI_added = ", TI_added
          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
  