program quake_per_day

** Note: This is not an optimal code -- it takes too much CPU time. See example 13 for a better one
*
* 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(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=1990) !! the first year to be counted

*1234567890123456789012345678901234567890123456789012345
*1962/01/05 00:23:32.10 -15.5000 -177.7000  24.00  6.50

	open(1,file='count.out')
	
	do k = iyear1, 2011  !! from iyear1 (which is 1990) to 2011

	print*,'Doing year ', k

	leep =0
	if(mod(k,4).eq.0) leep =1  !! Note: mod is an intrinsic function

	do j = 1, 365+leep

	nn=0

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

	do i=1,10000000 !! this number cannot be smaller than the actual number of quakes

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

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

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

** Skip the quake if it is not in year k
	
	if(iy.ne.k) goto 221

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

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

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

        yr=k*1.0 +  (1.0*j-0.5)/(365.0+leep)

221	continue
	enddo
222	close(2)

	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