	program hinv2cnss
c
c	programmer: david oppenheimer (2/8/99)

c	8/31/2004: 
c	  fixed bug in datenm that allowed seconds to slip through outside 0-60.  
c	  picks up SEED code from calsta2000.loc
c	  fixed f7.2 format in $amp output to f7.4
c	  added Y data source code for BK data, and occasional CI data sources
c	11/07/2007:
c	  added support for H=Hand magnitudes in sum2cnss.f
c
c	creates cnss format data (cnss-catalog-ver-1.0.3).

c	hypoinverse archive input read from standard input
c	summary output written to standard output
c	phase output written to hinv2cnss.phase
c	error output written to hinv2cnss.log
c	USGS-SEED conversion file read from calsta2000.loc

	implicit none
	integer		cunit				! logical unit number of file that contains USGS->SEED channel names
	integer		eunit				! logical unit number of error output
	integer		iunit				! logical unit number of phase file input
	integer		ounit				! logical unit number of phase output
	integer		sunit				! logical unit number of summary (catalog) output
	integer		mxphs				! maximum number of phases allowed in an event
	integer		ncode				! maximum number of entries in calsta2000.loc file
c
	parameter (eunit=3, iunit=5, ounit=4, sunit=6, cunit = 7)
	parameter (mxphs = 3000, ncode = 10000)
c
        real                    amag                    ! Eaton amplitude magnitude 
	real			ampmag			! station coda magnitude
	character*5		cal_name(ncode)		! station name
	character*2		cal_net(ncode)		! network code
	character*3		cal_usgs(ncode)		! usgs channel name
	character*3		cal_seed(ncode)		! seed channel name
	character*2		cal_loc(ncode)		! seed location code
	character*200           card(mxphs)		! phase cards
        real                    cmag                    ! Eaton coda magnitude 
	character*8		date			! date string
	character*8		date1(ncode)		! calsta2000 usage date string
	character*8		date2(ncode)		! calsta2000 usage date string
	character*17            dollar			! scratch variable to test for presence of shadow card & y2k
	real			durmag			! station duration magnitude
	logical                 erflag			! flag: true=errors found
	integer			id(3)			! numeric m/d/y date fields returned by system function
	integer                 i			! loop index
	integer                 iawt			! amplitude magnitude weight
	integer                 icwt			! coda magnitude weight
	integer                 ierr			! error return flag
	integer                 ios			! iostat variable
	integer                 iline			! line # of file
        integer                 iday			! day of origin time
        integer                 ihour			! hour of origin time
        integer                 imin			! minute of origin time
        integer                 imonth			! month of origin time
        integer                 iyear			! year of origin time
	integer                 j			! loop index
	integer                 leng			! function
        real                    lmag                    ! WA mag
	integer                 nampl			! number of WA amp readings used in Eaton amplitude magnitude
	integer                 nampx			! number of S amp readings used in local (ML) amplitude magnitude
	integer                 nc			! number of USGS->SEED conversions
	integer                 nch			! number of non-blank characters in line
	integer                 ncodv			! number of coda readings used in high-gain magnitude
	integer                 ncodz			! number of coda readings used in low-gain magnitude
	integer                 np			! number of phase cards
        real			sec			! sec of origin time
	character*200		shadow(mxphs)		! shadow card
	logical                 shdo			! flag:t(f)=do(not) read shadow card
	character*2		typmag			! CNSS magnitude type
	character*1             vrep			! flag:y(n)=do(not) replace a space in 4th column of station with a "v"
        real                    wmag                    ! Moment mag
	logical			y2k			! y2k format found for hypoinverse file
        real                    zmag                    ! Hirshorn Z coda mag

c	read calsta2000.loc for USGS->SEED conversion

	open (eunit, file = 'hinv2cnss.log', status = 'unknown')
	open (cunit, file='calsta2000.loc', iostat=ios, status='old')
	if (ios .ne. 0) then
	  write (eunit, *) 'calsta2000.loc not found'
	  stop
	endif
	nc = 1
90	read (cunit, 91, end = 100) cal_name(nc), cal_net(nc), 
     1 cal_usgs(nc), date1(nc), date2(nc), cal_seed(nc), cal_loc(nc)
91	format (a5, a2, 1x, a3, t90, 2(a8, 1x), t127, a3, 1x, a2)
	nc = nc + 1
	if (nc .gt. ncode) then
	  write (eunit, *) 'number of entries in calsta2000.loc > '
     1 , ncode
	  stop
	endif
	goto 90
100	nc = nc - 1
	close (cunit)

	np = 2
	read (iunit, '(a)') dollar
	if (dollar(9:9) .ne. ' ' .and.
     1 (dollar(17:17) .eq. ' ' .or. dollar(17:17) .eq. 's')) then
	  y2k = .false.
	elseif (dollar(12:12) .ne. ' ' .and. dollar(17:17).ne.' ')then
	  y2k = .true.
	else
	  write (eunit, *) 'y2k format is unclear'
	  stop
	endif

	read (iunit, '(a)') dollar
	if (dollar(1:1) .eq. '$') then
	  shdo = .true.
	else
	  shdo = .false.
	endif
	rewind(iunit)
	open (ounit, file = 'cnss.phase', status = 'unknown') 
	write (ounit, '(a)') '$fmt cnss-catalog-ver-1.0.3'
	write (sunit, '(a)') '$fmt cnss-catalog-ver-1.0.3'
	vrep = 'n'
	erflag = .false.
        call idate(id)
        write (date, '(i4, 2i2.2)'), id(3), id(2), id(1)
	iline = 0

c	next event; reset counter

20	np = 0
	ncodv = 0
	ncodz = 0
	nampx = 0
	nampl = 0
	ierr = 0
30	np = np + 1
	if (np .gt. 3000) then
	  write (eunit, *) 'number of phases in event exceeds ', mxphs
	  stop
	endif
        iline = iline + 1

c	read a line 

        read (iunit, '(a)', iostat = ios, end = 200) card(np)
	if (ios .ne. 0) then
	  write (eunit, 40) iline
40	  format (/, 1x, '***** error reading archive file at line', i)
	  stop
	end if

	if (shdo) then
	  iline = iline + 1
	  read (iunit, '(a)', iostat = ios) shadow(np)
	  if (ios .ne. 0) then
	    write (eunit, 41) iline
41	    format (/, 1x, '***** error shadow card at line', i)
	    stop
	  end if
	end if

c	end of event

	if (card(np)(1:4) .eq. '    ') then
          goto 45
	else

c	Reconstruct count of number of coda, amplitude readings for cnss output.  
c	Hypoinverse reports weighted number. CNSS needs unweighted count

	  if (np .gt. 1) then
	    if (.not. y2k) then
	      read (card(np), 42, iostat = ierr) icwt, durmag
42	      format (t67, i1, t79, f2.1)
	    else
	      read (card(np), 422, iostat = ierr) icwt, durmag
422	      format (t83, i1, t95, f3.2)
	    endif
	    if (ierr .ne. 0) then
	      write (eunit, 43) iline
43	      format ('***** error decoding card at line', i)
	      stop
	    end if
            if (durmag .gt. 0. .and. icwt .lt. 4) then
	      if (.not. y2k) then
	        if (card(np)(93:93) .eq. 'D') ncodv = ncodv + 1
	        if (card(np)(93:93) .eq. 'Z') ncodz = ncodz + 1
	      else
	        if (card(np)(110:110) .eq. 'D') ncodv = ncodv + 1
	        if (card(np)(110:110) .eq. 'Z') ncodz = ncodz + 1
	      endif
	    endif
	    if (.not. y2k) then
	      read (card(np), 44, iostat = ierr) iawt, ampmag
44	      format (t66, i1, t81, f2.1) 
	    else
	      read (card(np), 444, iostat = ierr) iawt, ampmag
444	      format (t82, i1, t98, f3.2) 
	    endif
	    if (ierr .ne. 0) then
	      write (eunit, 43) iline
	      stop
	    end if
            if (ampmag .gt. 0. .and. iawt .lt. 4) then
	      if (.not. y2k) then
	        if (card(np)(94:94) .eq. 'L') nampl = nampl + 1
	        if (card(np)(94:94) .eq. 'X') nampx = nampx + 1
	      else
	        if (card(np)(111:111) .eq. 'L') nampl = nampl + 1
	        if (card(np)(111:111) .eq. 'X') nampx = nampx + 1
	      endif
	    endif
	  endif
	  goto 30
	endif

c	summary card

45	ierr = 0
	write (ounit, '(a4)') '$beg'
	write (sunit, '(a4)') '$beg'
	call sum2cnss (card(1), ounit, sunit, eunit, date, cmag, amag, 
     1 zmag, lmag, wmag, iyear, imonth, iday, ihour, imin, sec, ncodv, 
     2 ncodz, nampl, nampx, y2k, ierr, typmag)

	if (ierr .eq. 1) then
	  if (shdo) then
	    i = iline - 2*np +1
	  else
	    i = iline - np +1
	  endif
	  write (eunit, 50) i
	  write (eunit, 50) i
50	  format ('Error on summary card on line ', i)
	  nch = leng(card(1))
	  write (eunit, '(a)') card(1)(1:nch)
	  write (eunit, '(a)') card(1)(1:nch)
	  stop
	endif

c	phase/shadow card

	do i = 2, np - 1
	  ierr = 0
	  call phs2cnss (card(i), ounit, eunit, shdo, shadow(i),
     2 cmag, amag, zmag, lmag, wmag, iyear, imonth, iday, ihour, 
     3 imin, sec, y2k, ierr, nc, cal_name, cal_net, cal_usgs, 
     4 cal_seed, cal_loc, typmag, date1, date2)
	  if (ierr .eq. 1) then
	    if (shdo) then
	      j = iline - 2*np + 2*i - 1
	    else
	      j = iline - np + i - 1
	    endif
	    write (eunit, '(a, i, 2a)') 'Error occurred on line ', j, 
     1 ' for eventid ', card(np)(63:72)
	    nch = leng(card(i))
	    write (eunit, '(a)') card(i)(1:nch)
	    write (eunit, '(a)') ' '
	  endif
	end do

	write (ounit, '(a4)') '$end'
	write (sunit, '(a4)') '$end'
	goto 20

c	end of file

200	if (np .gt. 1 ) write (eunit, *) 'unexpected EOF'
	close (eunit)
	close (iunit)
	close (ounit)
	close (sunit)
	stop
	end
