! Calculate the area of the Mandelbrot set by simple point counting.
! Because it is symmetric about the X-axis, we save time by doing
! that only for the upper half.  It offsets the points to correct for
! the double counting of the points on the axis.

PROGRAM Mandelbrot
    IMPLICIT NONE
    INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12),    &
        maxiters = 30000
    INTEGER :: i, j
    REAL(KIND=dp) :: step, total
    CHARACTER(LEN=20) :: arg

! Read in the parameters: the number of iterations and the step.

    IF (COMMAND_ARGUMENT_COUNT() /= 1) THEN
        WRITE (*,'("Wrong number of arguments")')
        STOP 1
    END IF
    CALL GET_COMMAND_ARGUMENT(1,arg)
    READ (arg,'(F20.0)') step
    IF (step < 10.0_dp*EPSILON(step) .OR. step > 0.1_dp .OR.    &
            1.0_dp/step >= REAL(HUGE(i),KIND=dp)) THEN
        PRINT *, 'Invalid value of STEP'
        STOP 1
    END IF

! Now scan over the area with corners (-2,0) and (2,2), which includes
! all of the Mandelbrot set.

    CALL Times('')
    total = 0.0_dp
!$OMP PARALLEL DO SCHEDULE(STATIC), DEFAULT(NONE), SHARED(step), &
    !$OMP PRIVATE(i,j), REDUCTION(+:total)
    DO i = INT(-2.0_dp/step), INT(2.0_dp/step)
        DO j = 0, INT(2.0_dp/step)
            total = total+kernel(CMPLX(i*step,(j+0.5_dp)*step,KIND=dp))
        END DO
    END DO
!$OMP END PARALLEL DO
    CALL Times('Mandelbrot')
    PRINT '(A,F6.3)', 'Area of Mandelbrot set is about',    &
        2.0_dp*total*step**2

CONTAINS

    SUBROUTINE Times (which)
! If which is not empty, print the times since the previous call.
        CHARACTER(LEN=*), INTENT(IN) :: which
        REAL(KIND=dp), SAVE :: last_wall = 0.0_dp, last_cpu = 0.0_dp
        REAL(KIND=dp) :: wall, cpu
        INTEGER :: m, n

        wall = last_wall
        cpu = last_cpu
        CALL SYSTEM_CLOCK(m,n)
        last_wall = m/REAL(n,KIND=dp)
        CALL CPU_TIME(last_cpu)
        IF (LEN(which) > 0) THEN
            wall = last_wall-wall
            cpu = last_cpu-cpu
            WRITE (*,'(A," time = ",F0.2," seconds, CPU = ",F0.2," seconds")') &
                which,wall,cpu
        END IF
    END SUBROUTINE Times

    PURE FUNCTION Kernel (value)
! See if this point is within the Mandelbrot set.
        IMPLICIT NONE
        INTEGER :: Kernel
        COMPLEX(KIND=dp), INTENT(IN) :: value
        COMPLEX(KIND=dp) :: work
        INTEGER :: n

        work = value
        DO n = 1, maxiters
            work = work**2 + value
            IF (ABS(REAL(work)) > 2.0_dp .OR. ABS(AIMAG(work)) > 2.0_dp) EXIT
        END DO
        Kernel = 0
        IF (ABS(WORK) <= 2.0) Kernel = 1
    END FUNCTION Kernel

END PROGRAM Mandelbrot
