program f90_julia_gmt3d
    implicit none
    integer,parameter::nmax=10000
    integer::ih,jw
    real(8)::x0,y0,rr
    real(8)::a,b
    real(8)::xx,yy
    integer::i,j
    real(8)::col
    character(len=50)::dummy
    character(len=50)::fnameF
    character(len=20)::fout='_out.txt'

    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,dummy) ; read(dummy,*) a
    call getarg(7,dummy) ; read(dummy,*) b
    call getarg(8,fnameF)

    write(6,'(a)') '* Job for "'//trim(adjustl(fnameF))//'" was started.'

    open(12,file=fout,status='replace')
    do i=0,ih
        write(12,'(a)') '>'
        yy=y0+rr-2.0D0*rr/dble(ih)*dble(i)
        do j=0,jw
            xx=x0-rr+2.0D0*rr/dble(jw)*dble(j)
            col=conv(xx,yy,a,b,nmax)
            write(12,'(3(e15.7))') xx,yy,col
        end do
        write(12,'(3(e15.7))') x0+rr,yy,0.0D0
    end do
    write(12,'()')
    close(12)

    write(6,'(a)') '** Calculation was completed.'
    call GMT3DSH(x0,y0,rr,fnameF,fout,jw,ih)
    call system('bash _gmt3d.sh')
    write(6,'(a)') '*** All works were completed.'

stop

contains

    real(8) function conv(xx,yy,a,b,nmax)
        real(8),intent(in)::xx,yy
        real(8),intent(in)::a,b
        integer,intent(in)::nmax
        real(8)::x1,y1,x2,y2,z
        integer::i

        conv=log10(dble(nmax))+1.0D0
        x1=xx
        y1=yy
        z=0.0d0
        do i=1,nmax
            x2=x1*x1-y1*y1+a
            y2=2.0D0*x1*y1+b
            z=x2*x2+y2*y2
            if(z>4.0D0)exit
            x1=x2
            y1=y2
        end do
        if(z>4.0D0)conv=log10(dble(i))
    end function conv

    subroutine GMT3DSH(x0,y0,rr,fnameF,fout,jw,ih)
        real(8),intent(in)::x0,y0,rr
        character(len=50),intent(in)::fnameF
        character(len=20),intent(in)::fout
        integer,intent(in)::jw,ih
        character(len=50)::dummy1,dummy2,dummy3,dummy4
        character(len=50)::fnameW='_gmt3d.sh'
        character(len=200)::strdat

        write(dummy1,'(e15.7)') x0-rr
        write(dummy2,'(e15.7)') x0+rr
        write(dummy3,'(e15.7)') y0-rr
        write(dummy4,'(e15.7)') y0+rr
        strdat=trim(adjustl(dummy1))//'/'//trim(adjustl(dummy2))
        strdat=trim(adjustl(strdat))//'/'//trim(adjustl(dummy3))//'/'//trim(adjustl(dummy4))

        open(12,file=fnameW,status='replace')
        write(12,'(a)') 'range='//trim(adjustl(strdat))//'/0/7'
        write(12,'(a)') 'scale=10/10'
        write(12,'(a)') 'inp_dat='//trim(adjustl(fout))
        write(12,'(a)') 'fig=_gmt.eps'
        write(12,'(a)') 'gmt psxyz $inp_dat -R$range -JX$scale -JZ5 -G255 -W0.25 -p150/30 -P > $fig'
        write(12,'(a)') 'convert -trim -density 300 $fig '//trim(adjustl(fnameF))
        close(12)
    end subroutine GMT3DSH
end program f90_julia_gmt3d