anchor.diff:          ! 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)*(TI_ust**0.0325)* & 
  -                    ((x/rotor_diameter_ust)**-0.32)
  +                    ((x/rotor_diameter_ust)**(-0.32_dp))
           !print *, "TI_added = ", TI_added
  -        rotor_area_dst = 0.25_dp*pi*rotor_diameter_dst**2
  +        rotor_area_dst = 0.25_dp*pi*rotor_diameter_dst**2_dp
           ! Calculate the total turbulence intensity at the downstream turbine
           !sum_of_squares = TI_dst**2 + (TI_added*wake_overlap)**2
           ! print *, "sum of squares = ", sum_of_squares
  
positive.code:                     !print *, "local y,z : ", LocalRotorPointY, LocalRotorPointZ
                      
                      ! find the final point where the original model is undefined
                      call discontinuity_point_func(x0, rotorDiameter(turb), ky, kz, yaw, Ct, & 
                                                   & discontinuity_point)
                                                   
                      ! far wake region
  
negative.code:     end if
  
      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
          else if (OVRR > OVr) then
              wake_overlap = pi*OVr*OVr
              !print *, "wake ovl: ", wake_overlap
          else
              wake_overlap = pi*OVRR*OVRR
          end if
      else
          wake_overlap = 0.0_dp
      end if
      
      !print *, "wake overlap in func: ", wake_overlap
                               
  end subroutine overlap_area_func
  
  