anchor.diff: +! Implementation of the Bastankhah and Porte Agel gaussian-shaped wind turbine wake 
  +! model (2016) with various farm modeling (TI and wake combination) methods included
  +! Created by Jared J. Thomas, 2017.
  +! FLight Optimization and Wind Laboratory (FLOW Lab)
  +! Brigham Young University
  +
   ! implementation of the Bastankhah and Porte Agel (BPA) wake model for analysis
   subroutine porteagel_analyze(nTurbines, nRotorPoints, turbineXw, sorted_x_idx, turbineYw, turbineZ, &
                                rotorDiameter, Ct, wind_speed, &
  
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:     
      ! TODO add other TI calculation methods
          
      ! wake combination method error 
      else
          print *, "Invalid added TI calculation method. Must be one of [1,2,3]."
          stop 1
      end if                       
      
  end subroutine added_ti_func
  
  ! compute wake spread parameter based on local turbulence intensity
  subroutine k_star_func(TI_ust, k_star_ust)
                                   
      implicit none
  