anchor.diff:      ! out  
       real(dp), intent(inout) :: k_star_ust
       
  +    ! calculate wake spread parameter from Niayifar and Porte Agel (2015, 2016)
       k_star_ust = 0.3837*TI_ust+0.003678
       
   end subroutine k_star_func
  
positive.code:     ! out
      real(dp), intent(out) :: deltav
  
      intrinsic cos, sqrt, exp
      
      !print *, "here 2.1", nPoints, pointsY(1), pointsZ(1)
      
      !deltavr = 0.9
      
      !print *, "inputs = ", nPoints, deltay, deltaz, Ct, yaw,  &
                !                    & sigmay, sigmaz, rotor_diameter_ust,  &
  !                                  & rotor_diameter_dst, x, x0,           &
  !                                  & sigmay0, sigmaz0, pointsY, pointsZ
  
      ! magnitude term of gaussian at x0
      deltav0m = (                                         &
                  (1.0_dp - sqrt(1.0_dp - Ct *                          &
                  cos(yaw) / (8.0_dp * sigmay0 * sigmaz0 /              &
                                              (rotor_diameter_ust ** 2)))))
      deltavr = deltav0m
      !print *, "deltav0m = ", deltav0m
      deltav_sum = 0.0_dp
      
      do, point=1, nPoints
          !print *, "point= ", point
          deltay_p = deltay + 0.5_dp*pointsY(point)*rotor_diameter_dst
          deltaz_p = deltaz + 0.5_dp*pointsZ(point)*rotor_diameter_dst    
          !print *, "dyp, dzp: ", deltay_p, deltaz_p
          !print *, "sum a = ", deltav_sum
          ! linearized gaussian magnitude term for near wake
          deltav_sum = deltav_sum + (((deltav0m - deltavr)/x0) * x + deltavr) *       &
              exp(-0.5_dp * (deltay_p / sigmay) ** 2) *                                 &
              exp(-0.5_dp * (deltaz_p / sigmaz) ** 2)
          !print *, "sum b = ", deltav_sum
      
      end do
      
      ! take the average of the velocity at all points
      deltav = deltav_sum/nPoints
      
      !print *, "here 2.2"
                  
  end subroutine deltav_near_wake_lin_func
  
  
negative.code:     
      ! Niayifar and Porte Agel 2015, 2016
      else if (TI_calculation_method == 2) then
      
          ! calculate axial induction based on the Ct value
          call ct_to_axial_ind_func(Ct_ust, axial_induction_ust)
          
          ! calculate BPA spread parameters Bastankhah and Porte Agel 2014
  