anchor.diff:          OVz = 0.0_dp
       end if
   
  -    if (OVdYd < (OVr+OVRR)) then ! if the rotor overlaps the wake zone
  +    if (OVdYd < (OVr+OVRR)) then ! if the rotor overlaps the wake
   
           if (OVL < OVRR .and. (OVdYd-OVL) < OVr) then
               wake_overlap = OVRR*OVRR*dacos(OVL/OVRR) + OVr*OVr*dacos((OVdYd-OVL)/OVr) - OVdYd*OVz
  
positive.code:     real(dp), parameter :: pi = 3.141592653589793_dp, tol = 0.000001_dp
      real(dp) :: OVdYd, OVr, OVRR, OVL, OVz
      
      intrinsic dacos
      
      OVdYd = wake_center_y-turbine_y        ! distance between wake center and rotor center
      OVr = rotor_diameter/2.0_dp            ! rotor diameter
      OVRR = wake_diameter/2.0_dp            ! wake diameter
      OVdYd = abs(OVdYd)
      if (OVdYd >= 0.0_dp + tol) then
          ! calculate the distance from the wake center to the vertical line between
          ! the two circle intersection points
          OVL = (-OVr*OVr+OVRR*OVRR+OVdYd*OVdYd)/(2.0_dp*OVdYd)
      else
          OVL = 0.0_dp
  
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
  