anchor.diff:    REAL(dp), DIMENSION(nbdirs), INTENT(OUT) :: wake_overlapd
   ! local
     REAL(dp), PARAMETER :: pi=3.141592653589793_dp, tol=0.000001_dp
  -  REAL(dp) :: ovdyd, ovr, ovrr, ovl, ovz
  -  REAL(dp), DIMENSION(nbdirs) :: ovdydd, ovrd, ovrrd, ovld, ovzd
  +  REAL(dp) :: ovdyd, ovr, ovrr, ovl, ovz, ovz2
  +  REAL(dp), DIMENSION(nbdirs) :: ovdydd, ovrd, ovrrd, ovld, ovzd, &
  +& ovz2d
   ! load intrinsic functions
     INTRINSIC ACOS, SQRT
     REAL(dp) :: arg1
  
positive.code: ! line between the two circle intersection points
  !if (OVdYd >= 0.0_dp + tol) then ! check case to avoid division by zero
  !     print *, "OVdYd ", OVdYd
    IF (ovdyd .GT. 0.0_dp + tol) THEN
      DO nd=1,nbdirs
  ! check case to avoid division by zero
        ovld(nd) = ((ovrrd(nd)*ovrr-ovr*ovrd(nd)-ovrd(nd)*ovr+ovrr*ovrrd(&
  &       nd)+ovdydd(nd)*ovdyd+ovdyd*ovdydd(nd))*2.0_dp*ovdyd-(-(ovr*ovr)+&
  &       ovrr*ovrr+ovdyd*ovdyd)*2.0_dp*ovdydd(nd))/(2.0_dp*ovdyd)**2
      END DO
      ovl = (-(ovr*ovr)+ovrr*ovrr+ovdyd*ovdyd)/(2.0_dp*ovdyd)
    ELSE
      ovl = 0.0_dp
      DO nd=1,nbdirs
        ovld(nd) = 0.0_8
      END DO
    END IF
    DO nd=1,nbdirs
      ovzd(nd) = ovrrd(nd)*ovrr + ovrr*ovrrd(nd) - ovld(nd)*ovl - ovl*ovld&
  &     (nd)
    END DO
    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 .GT. 0.0_dp + tol) THEN
      DO nd=1,nbdirs
        IF (ovz .EQ. 0.0) THEN
          ovzd(nd) = 0.0_8
        ELSE
          ovzd(nd) = ovzd(nd)/(2.0*SQRT(ovz))
        END IF
      END DO
      ovz = SQRT(ovz)
    ELSE
      ovz = 0.0_dp
      DO nd=1,nbdirs
        ovzd(nd) = 0.0_8
      END DO
    END IF
  !print *, "OVRR, OVL, OVRR, OVr, OVdYd, OVz ", OVRR, OVL, OVRR, OVr, OVdYd, OVz
    IF (ovdyd .LT. ovr + ovrr) THEN
  ! if the rotor overlaps the wake
  !print *, "OVL: ", OVL
      IF (ovl .LT. ovrr .AND. ovdyd - ovl .LT. ovr) THEN
        arg1 = ovl/ovrr
        result1 = ACOS(arg1)
        arg2 = (ovdyd-ovl)/ovr
        result2 = ACOS(arg2)
        DO nd=1,nbdirs
  !         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
          arg1d(nd) = (ovld(nd)*ovrr-ovl*ovrrd(nd))/ovrr**2
          IF (arg1 .EQ. 1.0 .OR. arg1 .EQ. (-1.0)) THEN
            result1d(nd) = 0.0_8
          ELSE
            result1d(nd) = -(arg1d(nd)/SQRT(1.0-arg1**2))
          END IF
          arg2d(nd) = ((ovdydd(nd)-ovld(nd))*ovr-(ovdyd-ovl)*ovrd(nd))/ovr&
  &         **2
          IF (arg2 .EQ. 1.0 .OR. arg2 .EQ. (-1.0)) THEN
            result2d(nd) = 0.0_8
          ELSE
            result2d(nd) = -(arg2d(nd)/SQRT(1.0-arg2**2))
          END IF
          wake_overlapd(nd) = (ovrrd(nd)*ovrr+ovrr*ovrrd(nd))*result1 + &
  &         ovrr**2*result1d(nd) + (ovrd(nd)*ovr+ovr*ovrd(nd))*result2 + &
  &         ovr**2*result2d(nd) - ovdydd(nd)*ovz - ovdyd*ovzd(nd)
        END DO
        wake_overlap = ovrr*ovrr*result1 + ovr*ovr*result2 - ovdyd*ovz
      ELSE IF (ovrr .GT. ovr) THEN
        DO nd=1,nbdirs
          wake_overlapd(nd) = pi*(ovrd(nd)*ovr+ovr*ovrd(nd))
        END DO
        wake_overlap = pi*ovr*ovr
  !print *, "wake ovl: ", wake_overlap
      ELSE
        DO nd=1,nbdirs
          wake_overlapd(nd) = pi*(ovrrd(nd)*ovrr+ovrr*ovrrd(nd))
        END DO
        wake_overlap = pi*ovrr*ovrr
      END IF
    ELSE
      wake_overlap = 0.0_dp
      DO nd=1,nbdirs
        wake_overlapd(nd) = 0.0_8
  
negative.code:   END IF
  !     print *, "wake overlap in func: ", wake_overlap/(pi*OVr**2)
  !     print *, "wake overlap in func: ", wake_overlap/(pi*OVRR**2)
    IF (wake_overlap/(pi*ovr*ovr) .GT. 1.0_dp + tol .OR. wake_overlap/(pi*&
  &     ovrr*ovrr) .GT. 1.0_dp + tol) THEN
      PRINT*, 'wake overlap in func: ', wake_overlap/(pi*ovr*ovr)
      PRINT*, 'wake overlap in func: ', wake_overlap/(pi*ovrr*ovrr)
      STOP
    END IF
  END SUBROUTINE OVERLAP_AREA_FUNC_DV
  
  