program quake_time

**
* Convert the earthquake times (2010/03/22 21:00:00.00 ) into decimal year (e.g., 2010.2215753)
*	for earthquakes in a region 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 cat.out
**

	parameter(flat1=30.0, flat2=40.0)   !! minimum and maximum latitudes of area of interest
	parameter(flon1=-120.0, flon2=-110.0)!! minimum and maximum longitude of area of interest


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

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

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

	do i=1,10000000
	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


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

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

        yr=iy*1.0 +  (1.0*jday-1.0)/(365.0+leep) + ih/24.0/(365+leep)
     &		     +  im/60.0/24.0/(365+leep)

	write(1,11) yr, flat, flon, fmag, dep
11	format(f12.7, 1x, f11.4, 1x, f9.4, f7.2, 1x, f7.2)


221	continue
	enddo
222	close(2)

	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