REAL*8 FUNCTION Mandel(x,y,niter)
  IMPLICIT NONE
  ! Arguments
  REAL*8, intent(in)  :: x
  REAL*8, intent(in)  :: y
  INTEGER, intent(in) :: niter
  ! Local variables
  COMPLEX*16 :: z, z0
  INTEGER    :: i

  z0=dcmplx(x,y)
  z=0
  DO i=1,niter
     z = z**2 + z0
     IF (abs(z)>2.) THEN
        Mandel = 1./(i+1.)
        return
     ENDIF
  ENDDO
  Mandel = 0
  RETURN
END FUNCTION Mandel

SUBROUTINE Cmp_Mandel(z, x, y, niter, nx, ny)
  IMPLICIT NONE
  ! Arguments
  REAL*8, intent(out) :: z(nx,ny)
  REAL*8, intent(in)  :: x(nx)
  REAL*8, intent(in)  :: y(ny)
  INTEGER, intent(in) :: niter
  INTEGER, intent(in) :: nx, ny
  !f2py integer intent(hide), depend(x)  :: nx = shape(x,0)
  !f2py integer intent(hide), depend(y)  :: ny = shape(y,0)
  !f2py integer optional,intent(in) :: niter = 30
   ! Local variables
  REAL*8 :: Mandel
  INTEGER :: i, j
  DO i=1,nx
     DO j=1,ny
        z(i,j) = Mandel(x(i), y(j), niter)
     ENDDO
  ENDDO
END SUBROUTINE Cmp_Mandel

