anchor.diff:                  ! far wake region
                   if (x > discontinuity_point) then
                   
  -                    !print *, "entering second logic"
  -
                       ! velocity difference in the wake
                       call deltav_func(deltay, deltaz, Ct(turb), yaw(turb),  &
  -                                    & sigmay, sigmaz, rotorDiameter(turb), opt_exp_fac, &
  -                                    & deltav)
  -                    !print *, "here 2"
  -                    !print *, 'deficit sum fw: ', deficit_sum
  +                                    & sigmay, sigmaz, rotorDiameter(turb), & 
  +                                    & wake_model_version, kz_local, x, & 
  +                                    & opt_exp_fac, deltav)
  +                                    
                   ! near wake region (linearized)
                   else
  -                    !print *, "entering else"
                       
                       ! determine distance from discontinuity point to far wake onset
                       deltax0_dp = discontinuity_point - x0
                   
  -                    !print *, "entering horizontal spread"
                       ! horizontal spread at far wake onset
                       call sigmay_func(ky_local(turb), deltax0_dp, rotorDiameter(turb), yaw(turb), &
                                       & sigmay_dp)
  -                    !print *, "entering vertical spread"
  +                                    
                       ! vertical spread at far wake onset
                       call sigmaz_func(kz_local(turb), deltax0_dp, rotorDiameter(turb), sigmaz_dp)
  -                    !print *, "entering near wake linearized"
  +                    
                       ! velocity deficit in the nearwake (linear model)
                       call deltav_near_wake_lin_func(deltay, deltaz,        &
                                        & Ct(turb), yaw(turb), sigmay_dp, sigmaz_dp,           & 
  -                                     & rotorDiameter(turb), x, x0, sigmay_dp, sigmaz_dp,    & 
  +                                     & rotorDiameter(turb), x, discontinuity_point, & 
  +                                     & sigmay_dp, sigmaz_dp,    & 
  +                                     & wake_model_version, kz_local(turb), x0, &
                                        & opt_exp_fac, deltav)
  -                    !print *, "exited near wake linearize"
  -                                     
  -                    !print *, "here 3"
  +
                   end if
               
                   ! combine deficits according to selected method wake combination method
                   call wake_combination_func(wind_speed, wtVelocity(turb), deltav,         &
                                              wake_combination_method, deficit_sum) 
  -                !print *, 'deficit sum nw: ', deficit_sum
  -            
  +                                           
               end if
               
  -            !print *, "end logic statements"
  -            
  -            ! if (x < 0.0_dp ) then
  -!                 !print *, "deltav: ", deltav
  -!                 !print *, "deficit sum at ", turb, ": ", deficit_sum 
  -!                 !print *, "deltax: ", x
  -!                 !print *, "deltav: ", deltav
  -!             end if
  -            
           end do
  -        
  -        !print *, "end inner do loop"
           
           ! find velocity at point p due to the wake of turbine turb
           point_velocity = wind_speed - deficit_sum
  
positive.code:                         call added_ti_func(TI, Ct(turb), x, ky_local(turb), rotorDiameter(turb), & 
                                             & rotorDiameter(turbI), deltay, turbineZ(turb), &
                                             & turbineZ(turbI), TIturbs(turb), &
                                             & TI_calculation_method, &
                                             & TIturbs(turbI))
                      end if
                  
  
negative.code:     
          ! final velocity calculation for turbine turbI (average equally across all points)
          wtVelocity(turbI) = wtVelocity(turbI)/nRotorPoints
          
          if (use_ct_curve) then
              call hermite_interpolation(nCtPoints, ct_curve_wind_speed, ct_curve_ct, wtVelocity(turbI), Ct(turbI))
          end if
         
      end do
  
      do, loc=1, nSamples
           
          ! set combined velocity deficit to zero for loc
          deficit_sum = 0.0_dp
  
          do, turb=1, nTurbines ! at turbineX-locations
  
              ! set velocity deficit at loc due to turb to zero
              deltav = 0.0_dp
  
              ! downstream distance between turbines
              x = velX(loc) - turbineXw(turb)
              
              if (x > 0.0_dp) then
              
                  ! determine the onset location of far wake
                  call x0_func(rotorDiameter(turb), yaw(turb), Ct(turb), alpha, TIturbs(turb), &
                              & beta, x0)
                  
                  ! downstream distance from far wake onset to downstream turbine
                  deltax0 = x - x0
          
                  ! 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(turb), deltax0, rotorDiameter(turb), yaw(turb), sigmay)
              
                  ! vertical spread
                  call sigmaz_func(kz_local(turb), deltax0, rotorDiameter(turb), sigmaz)
              
                  ! horizontal cross-wind wake displacement from hub
                  call wake_offset_func(rotorDiameter(turb), theta_c_0, x0, yaw(turb), &
  