anchor.diff:  
       ! bastankhah and porte agel 2016 define yaw to be positive clockwise, this is reversed
       yaw = - yawDeg*pi/180.0_dp
  +    
  +    ! initialize TI of all turbines to free-stream value
  +    TIturbs = TI
   
       do, d=1, nTurbines
       
  
positive.code: !         print *, turbineXw(turbI)
          
          do, u=1, nTurbines ! at turbineX-locations
              
              ! set this iterations velocity deficit to 0
              deltav = 0.0_dp
          
              ! get index of upstream turbine
              turb = sorted_x_idx(u) + 1
          
              ! determine the onset location of far wake
              call x0_func(rotorDiameter(turb), yaw(turb), Ct(turb), alpha, I, beta, x0)
          
              ! downstream distance between turbines
              x = turbineXw(turbI) - turbineXw(turb)
                  
              ! downstream distance from far wake onset to downstream turbine
              deltax0 = x - x0
  
              ! far wake region
              if (x >= x0) then
                  
                  ! 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, deltax0, rotorDiameter(turb), yaw(turb), sigmay)
                  
                  ! vertical spread
                  call sigmaz_func(kz, deltax0, rotorDiameter(turb), sigmaz)
                  
                  ! horizontal cross-wind wake displacement from hub
                  call wake_offset_func(rotorDiameter(turb), theta_c_0, x0, yaw(turb), &
                                       & ky, kz, Ct(turb), sigmay, sigmaz, wake_offset)
                                       
                  ! cross wind distance from downstream hub location to wake center
                  deltay = turbineYw(turbI) - (turbineYw(turb) + wake_offset)
                  
                  ! cross wind distance from hub height to height of point of interest
                  deltaz = turbineZ(turbI) - turbineZ(turb)
                  
                  ! velocity difference in the wake
                  call deltav_func(deltay, deltaz, wake_offset, wind_speed, Ct(turb), & 
                                   & yaw(turb), sigmay, sigmaz, rotorDiameter(turb), deltav)
      
                  ! linear supposition
                  !deficit_sum = deficit_sum + deltav
              
                  ! combine deficits according to selected method wake combination method
                  call wake_combination_func(wind_speed, wtVelocity(turb), deltav,         &
                                             wake_combination_method, deficit_sum)  
  
              ! near wake region (linearized)
              else if (x > 0.0_dp) then
                  
                  ! 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, deltax0, rotorDiameter(turb), yaw(turb), sigmay)
                  
                  ! vertical spread
                  call sigmaz_func(kz, deltax0, rotorDiameter(turb), sigmaz)
                                                  
                  ! horizontal cross-wind wake displacement from hub
                  call wake_offset_func(rotorDiameter(turb), theta_c_0, x0, yaw(turb), &
                                       & ky, kz, Ct(turb), sigmay, sigmaz, wake_offset)
                  
                  ! distance from downstream hub location to wake center
                  deltay = turbineYw(turbI) - (turbineYw(turb) + wake_offset)
                  
                  ! vertical distance from downstream hub location to wake center
                  deltaz = turbineZ(turbI) - turbineZ(turb)
                  
                  ! horizontal spread at far wake onset
                  call sigmay_func(ky, 0.0_dp, rotorDiameter(turb), yaw(turb), sigmay0)
                  
                  ! vertical spread at far wake onset
                  call sigmaz_func(kz, 0.0_dp, rotorDiameter(turb), sigmaz0)
                  
                  ! velocity deficit in the nearwake (linear model)
                  call deltav_near_wake_lin_func(deltay, deltaz, wake_offset, wind_speed, &
                                   & Ct(turb), yaw(turb), sigmay, sigmaz, & 
                                   & rotorDiameter(turb), x, x0, sigmay0, sigmaz0, deltav)
                  ! linear supposition
                  !deficit_sum = deficit_sum + deltav
              
                  ! combine deficits according to selected method wake combination method
                  call wake_combination_func(wind_speed, wtVelocity(turb), deltav,         &
                                             wake_combination_method, deficit_sum)                
              end if
              
              
              
          end do
          
          ! final velocity calculation for turbine turbI
  
negative.code:     ! in
      integer, intent(in) :: nTurbines
      integer, intent(in) :: wake_combination_method, TI_calculation_method
      logical, intent(in) :: calc_k_star
      real(dp), dimension(nTurbines), intent(in) :: turbineXw, turbineYw, turbineZ
      integer, dimension(nTurbines), intent(in) :: sorted_x_idx
      real(dp), dimension(nTurbines), intent(in) :: rotorDiameter, yawDeg
      real(dp), dimension(nTurbines), intent(in) :: Ct
      real(dp), intent(in) :: ky, kz, alpha, beta, TI, wind_speed
  
      ! 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, deltav0m, deltaz, sigmay0, sigmaz0, deficit_sum
      Integer :: u, d, turb, turbI
  