anchor.diff:          ! calculate BPA spread parameters Bastankhah and Porte Agel 2014
           beta = 0.5_dp*((1.0_dp + sqrt(1.0_dp - Ct_ust))/sqrt(1.0_dp - Ct_ust))
           epsilon = 0.2_dp*sqrt(beta)
  -        print *, "epsilon = ", epsilon
  +        !print *, "epsilon = ", epsilon
           ! calculate wake spread for TI calcs
           sigma = k_star_ust*x + rotor_diameter_ust*epsilon
           wake_diameter = 4.0_dp*sigma
  -        print *, "sigma = ", sigma
  +        !print *, "sigma = ", sigma
           ! calculate wake overlap ratio
           call overlap_area_func(deltay, turbine_height, rotor_diameter_dst, &
                               0.0_dp, wake_height, wake_diameter, &
                               wake_overlap)
  -        print *, "wake_overlap = ", wake_overlap   
  +        !print *, "wake_overlap = ", wake_overlap   
           ! 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)
  -        print *, "TI_added = ", TI_added
  +        !print *, "TI_added = ", TI_added
           rotor_area_dst = 0.25_dp*pi*rotor_diameter_dst**2
           ! Calculate the total turbulence intensity at the downstream turbine
           !sum_of_squares = TI_dst**2 + (TI_added*wake_overlap)**2
  
positive.code:             
                      ! horizontal spread
                      call sigmay_func(ky_local, deltax0, rotorDiameter(turb), yaw(turb), sigmay)
                      print *, "sigmay ", sigmay
                      ! vertical spread
                      call sigmaz_func(kz_local, 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)
                      print *, "wake_offset ", wake_offset                 
                      ! cross wind distance from downstream hub location to wake center
                      deltay = turbineYw(turbI) - (turbineYw(turb) + wake_offset)
                  
  
negative.code:     
      ! intrinsic functions
      intrinsic sqrt
      !print *, "TI_dst in: ", TI_dst
      ! Niayifar and Porte Agel 2015, 2016 (adjusted by Annoni and Thomas for SOWFA match 
      ! and optimization)
      if (TI_calculation_method == 1) then
  