anchor.diff:      ! out
       real(dp), intent(out) :: deltav
   
  +    ! load intrinsic functions
       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)))))
  +    ! initialize the gaussian magnitude term at the rotor for the linear interpolation
       deltavr = deltav0m
  -    !print *, "deltav0m = ", deltav0m
  +    
  +    ! set the velocity sum across all points to zero
       deltav_sum = 0.0_dp
       
  +    ! calculate the velocity at each rotor sample point
       do, point=1, nPoints
  -        !print *, "point= ", point
  +    
  +        ! scale rotor sample point locations by rotor radius and place them in the wake 
  +        ! reference frame
           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
  +    ! compute effective hub velocity by taking the average of the velocity at all points
       deltav = deltav_sum/nPoints
  -    
  -    !print *, "here 2.2"
                   
   end subroutine deltav_near_wake_lin_func
   
  
positive.code:     ! load intrinsic functions
      intrinsic cos, sqrt, exp
      
      !print *, "here 1.1"
      
      deltav_sum = 0.0_dp
      
      do, point=1, nPoints
          deltay_p = deltay + 0.5_dp*pointsY(point)*rotor_diameter_dst
          deltaz_p = deltaz + 0.5_dp*pointsZ(point)*rotor_diameter_dst
          
          ! velocity difference in the wake
          deltav_sum = deltav_sum + (                                                              &
              (1.0_dp - sqrt(1.0_dp - Ct *                                                         &
                             cos(yaw) / (8.0_dp * sigmay * sigmaz / (rotor_diameter_ust ** 2)))) *     &
              exp(-0.5_dp * ((deltay_p) / sigmay) ** 2) * exp(-0.5_dp * ((deltaz_p) / sigmaz) ** 2)&
          )
  
      end do 
      
      ! take the average of the velocity at all points
      deltav = deltav_sum/nPoints
      
      !print *, "here 1.2"
      
  end subroutine deltav_func
  
  
  
negative.code:     ! 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
  