anchor.diff:                  ! distance from downstream hub location to wake center
                   deltay = turbineYw(turbI) - (turbineYw(turb) + wake_offset)
                   
  -                ! magnitude term of gaussian at x0
  -                deltav0m = wind_speed * (                                         &
  -                           (1.0_dp - sqrt(1.0_dp - Ct(turb) *                     &
  -                           cos(yaw(turb)) / (8.0_dp * sigmay * sigmaz /           &
  -                                                        (rotorDiameter(turb) ** 2)))))
  -                                                        
  -                ! linearized gaussian magnitude term for near wake
  -                deltav = (deltav0m/x0) * x *  &
  -                    exp(-0.5_dp * ((deltay) / sigmay) ** 2) *                   &
  -                    exp(-0.5_dp * ((turbineZ(turbI) - turbineZ(turb)) / sigmaz) ** 2)
  -                
  +                ! vertical distance from downstream hub location to wake center
  +                deltaz = turbineZ(turbI) - turbineZ(turb)
  +                
  +                ! horizontal spread at far wake onset
  +                call sigmay_func(ky, 0.0_dp, rotorDiameter(turb), yaw(turb), sigmay0)
  +                
  +                ! vertical spread at far wake onset
  +                call sigmaz_func(kz, 0.0_dp, rotorDiameter(turb), sigmaz0)
  +                
  +                call deltav_near_wake_lin_func(deltay, deltaz, wake_offset, wind_speed, &
  +                                 & Ct(turb), yaw(turb), sigmay, sigmaz, & 
  +                                 & rotorDiameter(turb), x, x0, sigmay0, sigmaz0, deltav)
  +                                 
                   ! linear wake superposition (additive)
                   wtVelocity(turbI) = wtVelocity(turbI) - deltav
   
  
positive.code: 
      intrinsic cos, sqrt
      
  	! determine the initial wake angle at the onset of far wake
  	theta_c_0 = 0.3_dp * yaw * (1.0_dp - sqrt(1.0_dp - Ct * cos(yaw))) / cos(yaw)
  
  end subroutine theta_c_0_func
  
  ! calculates the horizontal spread of the wake at a given distance from the onset of far 
  ! wake condition
  
negative.code:     ! in
      real(dp), intent(in) :: deltay, deltaz, wake_offset, wind_speed, Ct, yaw, sigmay
      real(dp), intent(in) :: sigmaz, rotor_diameter
      
      ! out
      real(dp), intent(out) :: deltav
  
      intrinsic cos, sqrt, exp
      
      ! velocity difference in the wake
      deltav = wind_speed * (                                                 &
          (1.0_dp - sqrt(1.0_dp - Ct *                                                     &
                         cos(yaw) / (8.0_dp * sigmay * sigmaz / (rotor_diameter ** 2)))) * &
          exp(-0.5_dp * ((deltay) / sigmay) ** 2) * exp(-0.5_dp * ((deltaz) / sigmaz) ** 2)&
      )
      
  end subroutine deltav_func
  
  
  ! calculates the velocity difference between hub velocity and free stream for a given wake
  ! for use in the near wake region only
  subroutine deltav_near_wake_lin_func(deltay, deltaz, wake_offset, wind_speed, Ct, yaw, &
                                   & sigmay, sigmaz, rotor_diameter, x, x0, &
                                   & sigmay0, sigmaz0, deltav) 
                         
      implicit none
  
      ! define precision to be the standard for a double precision ! on local system
      integer, parameter :: dp = kind(0.d0)
  
      ! in
      real(dp), intent(in) :: deltay, deltaz, wake_offset, wind_speed, Ct, yaw, sigmay
      real(dp), intent(in) :: sigmaz, rotor_diameter, x, x0, sigmay0, sigmaz0
      
      ! local
      real(dp) :: deltav0m, deltavs
      
      ! out
      real(dp), intent(out) :: deltav
  
      intrinsic cos, sqrt, exp
      
      deltavs = 0.9*wind_speed
  
      ! magnitude term of gaussian at x0
      deltav0m = wind_speed * (                                         &
                  (1.0_dp - sqrt(1.0_dp - Ct *                     &
                  cos(yaw) / (8.0_dp * sigmay0 * sigmaz0 /           &
                                              (rotor_diameter ** 2)))))
                                              
      ! linearized gaussian magnitude term for near wake
      deltav = (((deltav0m - deltavs)/x0) * x + deltavs) *  &
          exp(-0.5_dp * (deltay / sigmay) ** 2) *                   &
          exp(-0.5_dp * (deltaz / sigmaz) ** 2)
      
  end subroutine deltav_near_wake_lin_func
  
  
  !        Generated by TAPENADE     (INRIA, Ecuador team)
  