program quake_per_day

** Note: We use arrays here. It is much faster than example 12

* Count the number of quakes per day in an area, and write the results against the decimal year
*	The region is defined by flat1, flat2, flon1, flon2
*
*       3/22/2010. Steve Gao
*
* The input file is:
*       /home/sgao/compgeop/quakes.dat 
*  Note: This is an ANSS catalog and the first two lines in the catalog must be removed
*
* The output file is count.out
**

	parameter(max=1000000) !! this number needs to be larger than the actual number of quakes
	parameter(flat1=-90.0, flat2=90.0)  !! minimum and maximum latitudes of area of interest
	parameter(flon1=-180.0, flon2=180.0)!! minimum and maximum longitude of area of interest
	parameter(iyear1=2000) !! the first year to be counted

	integer iy(max), mo(max), id(max)
	real  flat(max), flon(max)

*1980/01/01 16:41:45.70 -31.2110   58.7170  10.00  5.00
*123456789012345678901234567890123456789012345678901234
*         1         2         3         4         5

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

** First, read the data into arrays:

	open(2,file='/home/sgao/compgeop/quakes.dat')

	do i=1, max 

	read(2,21,end=222) iy(i), mo(i), id(i), ih, im, 
     &            flat(i), flon(i), dep, fmag
21	format(i4, 1x, i2, 1x, i2, 1x, i2, 1x, i2, 7x,
     &            f8.4, 1x, f9.4, f7.2, f6.2)

	enddo
222	close(2)

	npts=i-1

	print*,'The total number of events in the catalog is ', npts

** The data were read into arrays

	do k = iyear1, 2018  !! from iyear1 to 2018

	print*,'Doing year ', k

	leep =0
	if(mod(k,4).eq.0) leep =1

	do j = 1, 365+leep

	nn=0

	do i=1,npts

** The following statement skips the quakes outside the region:

	if(flat(i).lt.flat1.or.flat(i).gt.flat2.or.
     &	   flon(i).lt.flon1.or.flon(i).gt.flon2) goto 221

** The following statement skips the quake if it did not happen in year k:

	if(iy(i).ne.k) goto 221

** This call is to the subroutine converts the date to day of the year

	call date2day(mo(i), id(i), jday, iy(i), leep)

	if(jday.eq.j) nn=nn+1

221	continue
	enddo

        yr=k*1.0 +  (1.0*j-0.5)/(365.0+leep)
	write(1,11) yr, nn
11	format(f12.7, 1x, i8)

	enddo !! j ends
	enddo !! k ends

	close(1)

	stop
	end


**************************************************************************
        subroutine date2day(m,iday,jday,iy, l)

        l=0
        if(mod(iy,4).eq.0) l=1    !! deal with leep year (e.g., 2008 has 366 days)

        if (m.eq.01) jstart=0
        if (m.eq.02) jstart=31
        if (m.eq.03) jstart=59+l
        if (m.eq.04) jstart=90+l
        if (m.eq.05) jstart=120+l
        if (m.eq.06) jstart=151+l
        if (m.eq.07) jstart=181+l
        if (m.eq.08) jstart=212+l
        if (m.eq.09) jstart=243+l
        if (m.eq.10) jstart=273+l
        if (m.eq.11) jstart=304+l
        if (m.eq.12) jstart=334+l

        jday=jstart+iday

        end