anchor.diff:      OVz = OVRR*OVRR-OVL*OVL
   
       ! Finish calculating the distance from the intersection line to the outer edge of the wake
  -    if (OVz > 0.0_dp + tol) then
  +    !if (OVz > 0.0_dp + tol) then
  +    if (OVz > 0.0_dp) then
           OVz = sqrt(OVz)
       else
           OVz = 0.0_dp
       end if
  +    
  +    !print *, "OVRR, OVL, OVRR, OVr, OVdYd, OVz ", OVRR, OVL, OVRR, OVr, OVdYd, OVz
  +    
  +    
   
       if (OVdYd < (OVr+OVRR)) then ! if the rotor overlaps the wake
           !print *, "OVL: ", OVL
  -        !if (OVL < OVRR .and. (OVdYd-OVL) < OVr) then
  -        if (OVdYd > 0.0_dp + tol) then
  -            wake_overlap = OVRR*OVRR*dacos(OVL/OVRR) + OVr*OVr*dacos((OVdYd-OVL)/OVr) - OVdYd*OVz
  +        if (OVL < OVRR .and. (OVdYd-OVL) < OVr) then
  +!         if (OVdYd > 0.0_dp + tol) then
  +!         if ((OVdYd > 0.0_dp) .and. (OVdYd > (OVRR - OVr))) then
  +            ! print *, "acos(OVL/OVRR), acos((OVdYd-OVL)/OVr), OVRR, OVL, OVr, OVdYd, OVL/OVRR, (OVdYd-OVL)/OVr ", &
  +!     & acos(OVL/OVRR), acos((OVdYd-OVL)/OVr), OVRR, OVL, OVr, OVdYd, OVL/OVRR, (OVdYd-OVL)/OVr
  +            wake_overlap = OVRR*OVRR*acos(OVL/OVRR) + OVr*OVr*acos((OVdYd-OVL)/OVr) - OVdYd*OVz
           else if (OVRR > OVr) then
               wake_overlap = pi*OVr*OVr
               !print *, "wake ovl: ", wake_overlap
  
positive.code:     real(dp) :: OVdYd, OVr, OVRR, OVL, OVz
      
      ! load intrinsic functions
      intrinsic dacos, sqrt
      
      ! distance between wake center and rotor center
      OVdYd = sqrt((wake_center_y-turbine_y)**2_dp + (wake_center_z - turbine_z)**2_dp)
  
negative.code:         
          ! Check if this is the max and use it if it is
          if (TI_tmp > TI_dst) then
             TI_dst = TI_tmp
          end if
          
      ! Niayifar and Porte Agel 2015, 2016 with smooth max
       else if (TI_calculation_method == 3) 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
          beta = 0.5_dp*((1.0_dp + sqrt(1.0_dp - Ct_ust))/sqrt(1.0_dp - Ct_ust))
          epsilon = 0.2_dp*sqrt(beta)
          
          ! calculate wake spread for TI calcs
          sigma = k_star_ust*x + rotor_diameter_ust*epsilon
          wake_diameter = 4.0_dp*sigma
          
  !         print *, "sigma, k_star_ust, x, rotor_diameter_ust, epsilon ", sigma, k_star_ust, x, rotor_diameter_ust, epsilon
          
          ! print *, "deltay, turbine_height, rotor_diameter_dst, wake_height, wake_diameter", &
  !                 & deltay, turbine_height, rotor_diameter_dst, &
  !                             wake_height, wake_diameter
          
          ! calculate wake overlap ratio
          call overlap_area_func(deltay, turbine_height, rotor_diameter_dst, &
                              0.0_dp, wake_height, wake_diameter, &
                              wake_overlap)
                              
          ! Calculate the turbulence added to the inflow of the downstream turbine by the 
          ! wake of the upstream turbine
          TI_added = 0.73_dp*(axial_induction_ust**0.8325_dp)*(TI_ust**0.0325_dp)* & 
                      ((x/rotor_diameter_ust)**(-0.32_dp))
          
          ! Calculate the total turbulence intensity at the downstream turbine based on 
          ! current upstream turbine
          rotor_area_dst = 0.25_dp*pi*rotor_diameter_dst**2_dp
          TI_tmp = sqrt(TI**2.0_dp + (TI_added*(wake_overlap/rotor_area_dst))**2.0_dp)
          
          !print *, "TI, TI_added, wake_overlap, rotor_area_dst: ", TI, TI_added, wake_overlap, rotor_area_dst
          
          ! Check if this is the max and use it if it is
          !if (TI_tmp > TI_dst) then
          !    TI_dst = TI_tmp
          !end if
  !         print *, "before: ", TI_dst, TI_tmp
          call smooth_max(TI_dst, TI_tmp, TI_dst)
  !         print *, "after:: ", TI_dst, TI_tmp
       
      !print *, "sigma: ", sigma
      ! TODO add other TI calculation methods
          
      ! wake combination method error 
      else
          !print *, "Invalid added TI calculation method. Must be one of [0,1,2,3]."
          stop 1
      end if            
      
  