program f90_buddhabrot_gmt
    implicit none
    integer,parameter::nmax=10000
    integer,parameter::nrmax=1000000
    real(8)::randx,randy
    integer::ih,jw
    real(8)::x0,y0,rr
    real(8)::xmin,xmax,ymin,ymax
    real(8)::a,b,x1,x2,y1,y2,z
    integer::irand,ii,kx,ky
    integer,allocatable::tempa(:,:),kgrid(:,:)
    integer::i,j
    integer::col,colmax
    character(len=50)::dummy
    character(len=20)::fout

    integer::clock                !for random_seed()
    integer::nrand                ! for random_seed()
    integer,allocatable::seed(:) ! for random_seed()

    call random_seed(size=nrand)
    allocate(seed(1:nrand))
    call system_clock(count=clock)
    seed=clock
    call random_seed(put=seed)

    call getarg(1,dummy) ; read(dummy,*) jw
    call getarg(2,dummy) ; read(dummy,*) ih
    call getarg(3,dummy) ; read(dummy,*) x0
    call getarg(4,dummy) ; read(dummy,*) y0
    call getarg(5,dummy) ; read(dummy,*) rr
    call getarg(6,fout)

    write(6,'(a)') '* Job was started.'

    allocate(tempa(1:jw+1,1:ih+1),kgrid(1:jw+1,1:ih+1))
    do j=1,jw+1
        do i=1,ih+1
            kgrid(j,i)=0
        end do
    end do
    xmin=x0-rr
    xmax=x0+rr
    ymin=y0-rr
    ymax=y0+rr

    do irand=1,nrmax
        write(6,*) irand
        call random_number(randx)
        call random_number(randy)
        a=xmin+randx*(xmax-xmin)
        b=ymin+randy*(ymax-ymin)
        x1=0.0D0
        y1=0.0D0
        z=0.0D0
        do j=1,jw+1
            do i=1,ih+1
                tempa(j,i)=0
            end do
        end do
        do ii=1,nmax
            x2=x1*x1-y1*y1+a
            y2=2.0D0*x1*y1+b
            z=x2*x2+y2*y2
            if(4.0D0<z)exit
            kx=idnint((x2-xmin)*dble(jw+1)/(xmax-xmin))+1
            ky=idnint((y2-ymin)*dble(ih+1)/(ymax-ymin))+1
            tempa(kx,ky)=tempa(kx,ky)+1
            x1=x2
            y1=y2
        end do
        if(z<=4.0D0)then
            do j=1,jw+1
                do i=1,ih+1
                    tempa(j,i)=0
                end do
            end do
        end if
        do j=1,jw+1
            do i=1,ih+1
                kgrid(j,i)=kgrid(j,i)+tempa(j,i)
            end do
        end do
    end do

    colmax=0
    open(12,file=fout,status='replace')
    do i=1,ih+1
        write(12,'(a)') '>'
        do j=1,jw+1
            col=0
            if(1<=kgrid(j,i))col=kgrid(j,i)
            write(12,'(3(i10))') i-1,j-1,col
            if(colmax<col)colmax=col
        end do
    end do
    write(12,'()')
    close(12)

    write(6,*) 'colmax=',colmax
    write(6,'(a)') '** Calculation was completed.'

stop

end program f90_buddhabrot_gmt
