anchor.diff:                          deltax0_dp = discontinuity_point - x0
                   
                           ! horizontal spread at discontinuity point
  -                        call sigmay_func(ky_local, deltax0_dp, rotorDiameter(turb), yaw(turb), sigmay_dp)
  +                        call sigmay_func(ky_local(turb), deltax0_dp, rotorDiameter(turb), yaw(turb), sigmay_dp)
                   
                           ! vertical spread at discontinuity point
  -                        call sigmaz_func(kz_local, deltax0_dp, rotorDiameter(turb), sigmaz_dp)
  +                        call sigmaz_func(kz_local(turb), deltax0_dp, rotorDiameter(turb), sigmaz_dp)
                   
                           ! velocity deficit in the nearwake (linear model)
                           call deltav_near_wake_lin_func(deltay, deltaz, &
  
positive.code:     
      end do
      
      
      !print *, wtVelocity
  
      !! make sure turbine inflow velocity is non-negative
  
negative.code:                 !print *, "here 1"
                  
                  ! find the final point where the original model is undefined
                  call discontinuity_point_func(x0, rotorDiameter(turb), ky_local(turb), kz_local(turb), yaw(turb), Ct(turb), & 
                                               & discontinuity_point)
                                               
                  ! far wake region
  