* This program makes a list of all the vertical component (0z) SAC files
*	under a directory, and read some info from the header of the files

** By Steve Gao (sgao@mst.edu), 11/16/2018. For Computational Geophysics

* The parameter statement defines the value of the variable max
	parameter(max=1000000)

* The following line defines two character variables. 
*  Both variables (cmd and nin) have a maximum length of 100 characters.
*  cmd is for calling subroutine "system", and nin is the name of
*	the SAC input file(s).

	character cmd*100, nin*100

* stname in the following statement is a character variable to hold 
*	the station name. It has a maximum length of 10 characters.

	character stname*10

* The following statement defines a 1-dimension real array with a 
*	maximum number of elements of max (which is defined in the
*	parameter statement above.

	real x(max)

** All the definitions and other "headers" of the program have been defined.

* Remove the existing list file (to prepare for the "ls" command in
*	the next statements).

* The "call system(cmd)" statement will execute the system command defined
*	in "cmd" directly from inside this Fortran program.

	cmd='/bin/rm sac.list'
	call system(cmd)

* Make a list of the 0z (i.e., vertical) component of the seismograms.
c Note: Please change this to your own path. Note that
c		you can use relative path name (e.g. ../../data etc.).

	cmd='ls /home/sgao/demo/04_sac/data/EQ*/*0z* >sac.list'
	call system(cmd)

* Open the output file. This will hold the station locations:

	open(2,file='stloc.out')

* Work on each of the files in sac.list:

	open(1,file='sac.list')
	do k=1,max

* Read the k-th file into variable nin:
* Note that the lines in sac.list cannot be longer than 100 characters
* Read from unit 1, with a format stated in statement 11. When the read
*	reaches the end of the file, go to statement 111 (which closes unit 1).

	read(1,11,end=111) nin

* The format for the above read statement. a100 means 100 characters:

11	format(a100)

* Prepare for reading the header values by calling rsac1 (which is a
*	subroutine in the /home/sgao/lib/sacio.a library file):
	
	call rsac1(nin, x, npts, beg, dt, max, nerr)

* In the above, nin is the input SAC file name.
*               x is the array to hold the data points of the SAC file
*               npts is the number of points in the SAC file returned by rsac1
*               beg is the beginning tme in the SAC file returned by rsac1
*               dt is the sampling interval in second in the SAC file returned by rsac1
*               max is the max value defined in the "parameter" statement
*               nerr is a "flag" for the correctness of the operation. Ignore it.

* Get the station latitude and longitude by calling getfhv.
* getfhv is a subroutine provided by the SAC package (in the sacio.a library together
*	with a bunch of other subroutines including rsac1 and wsac1 etc.).

	call getfhv('STLA', stlat, nerr)
	call getfhv('STLO', stlon, nerr)

* In the above statement, STLA and STLO are SAC header items for station latitude
*	and station longitude, respectively. You can find other SAC header items by
*	starting SAC, reading a SAC file into sac, and using the lh command.
*   stlat is the variable to hold station latitude.
*   stlon is the variable to hold station longitude.

* Get station name by calling getkhv. KSTNM is a SAC header for station name.
* "k" here stand for character. getkhv is a subroutine to get character type
*   SAC header values

	call getkhv('KSTNM', stname, nerr)

* Write the output to stloc.out

	write(2,21) stname, stlat, stlon
21	format(a10, 1x, f10.4, 1x, f10.4)

	enddo !! k loop ends here

111	close(1)  !! close unit 1

	close(2)  !! close unit 2

	stop
	end

** Note: You should copy the Makefile to your dir from this directory
*	and compile your code by typing "make".