anchor.diff:      real(dp) :: OVdYd, OVr, OVRR, OVL, OVz
       
       ! load intrinsic functions
  -    intrinsic dacos, sqrt
  +    intrinsic acos, sqrt
       
       ! distance between wake center and rotor center
       OVdYd = sqrt((wake_center_y-turbine_y)**2_dp + (wake_center_z - turbine_z)**2_dp)
  
positive.code:     
  end subroutine discontinuity_point_func
  
  
negative.code:     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) 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
  !         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
  