** Computing seismic b-value using the maximum likelihood method of Aki
* By Steve Gao. All right reserved.

	parameter(max=100000)
	real*4 fmca(max), fba(max)
	integer mm(max)
	character cmd*100, ncat*50

	write(6,60)
60	format('Enter the full name of the catalog: ',$)
	read(5,51) ncat
51	format(a50)

	write(6,61)
61	format('Enter min and max latitude of your study area: ',$)
	read(5,*) flatmin, flatmax

	write(6,62)
62	format('Enter min and max longitude of your study area: ',$)
	read(5,*) flonmin, flonmax

	write(6,63)
63	format('Enter block size in deg.: ',$)
	read(5,*) dblock

	write(6,64)
64	format('Enter moving distance in deg.: ',$)
	read(5,*) dmove


** First of all, choose quakes in the study area and write them to cat.cut

	open(4,file='cat.cut')

	open(3,file=ncat)

	do i=1,10000000
	read (3,31, end=5651) evlat, evlon, evmag
31      format(21x, f10.4, f10.4, 9x, f4.2)

	if(evlat.lt.flatmin.or.evlat.gt.flatmax.or.
     &     evlon.lt.flonmin.or.evlon.gt.flonmax) goto 5600

	write(4,31) evlat, evlon, evmag

5600	continue

	enddo
5651	close(3)

	close(4)

** Start the iteration for blocks

	nlat=(flatmax-flatmin-dblock)/dmove+1
	nlon=(flonmax-flonmin-dblock)/dmove+1

	open(8,file='output.dat')

	do kk=1, nlon

	flon1=(kk-1)*dmove + flonmin
	flon2=flon1+dblock

	print*,'Doing longitude range ', flon1, ' ', flon2

	do k=1, nlat

	flat1=(k-1)*dmove + flatmin
	flat2=flat1+dblock


	cmd='/bin/rm tmp.out'
	call system(cmd)	

	floge=0.43429

	open(1,file='tmp.out')

	do fmc=2.0,5,0.05

	nn=0
	fbar=0

*2000.0382080 2000 014   25.5910  101.1640  33.00  5.50
*12345678901234567890123456789012345678901234567890123456
*         1         2         3         4         5
*1960/01/05 18:01        34.0148 -117.7107   6.22  3.03        


	open(3,file='cat.cut')

	do i=1,10000000
	read (3,31, end=5656) evlat, evlon, evmag

	if(evlat.lt.flat1.or.evlat.gt.flat2.or.
     &     evlon.lt.flon1.or.evlon.gt.flon2) goto 5655

	if(evmag.ge.fmc) then
		nn=nn+1
		fbar=fbar+evmag
	endif

5655	continue

	enddo !!i ends
5656    close(3)

	if(nn.ge.5.and.fbar.gt.fmc) then
	         fbar=fbar/nn
		 b=floge/(fbar-fmc)
  		write(1,202) fmc,b, nn
202     	format(f5.2, 1x, f10.4, 1x, i7)
	endif

	enddo !! fmc ends here

	close(1)


*****All loops closed, now calculate the difference of the B values between steps

c	print*,'nn is ', nn

c	if(nn.lt.1) goto 1214

	open(1, file='tmp.out')

	do i=1,max
	read(1,202,end=1212) fmca(i), fba(i), mm(i)
	enddo
1212	close(1)

	npts=i-1

	if(npts.le.2) goto 1215

	do i=1,npts
	dif=fba(i+1)-fba(i)
	if(dif.le.0.00) then
		real_b=fba(i)
		real_fmc=fmca(i)
		mm_real=mm(i)
		goto 1213
	endif
	enddo

1213    continue


c	print*,'The b value is ', real_b,'; ', real_fmc


	if(mm_real.ge.5) then
	write(8, 81) 0.5*(flat1+flat2), 0.5*(flon1+flon2), 
     &		real_b, real_fmc, mm_real
81	format(4(f10.4,1x), i8)
	write(6, 81) 0.5*(flat1+flat2), 0.5*(flon1+flon2), 
     &		real_b, real_fmc, mm_real
	endif


1214 	continue

	write(6, 71) 0.5*(flat1+flat2), 0.5*(flon1+flon2)
71	format('Done with ', 2(f10.4,1x))

1215	continue

	enddo !! lat loop
	enddo !! lon loop

	close(8)

	stop
	end