subroutine st_dml(elatx,elonx,edepx,amla,semla,nmla,amlb,semlb,
     1                  nmlb,fnamam,fnamml,epi_id)
c
c
************************************************************************
*
*  by R. A. Uhrhammer
*    
*  University of California Seismographic Stations
*  Berkeley, Calif. 94720
*
*
*  This program calculates the local magnitude (ML) 
*     from data in file ampli.file generated by the Seistool.
*
*  Modified 21 Nov 94 to acquire file names from config.h file
*
*  Modified 19 July 95 to acquire station ML adjustment file names from
*     config.h file and to fix bug in ML se determination
*
*  Modified 28 Mar 97 to include 2 & 3 labeled horizontal components
*     in calculation of ML
*
************************************************************************
c
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      character*80 fnamam,fnamml
      character*9 epi_id
c
      character*4 stn(1000)
      character*3 scomp(1000)
      real*8 slat(1000),slon(1000),samp(1000), mla(1000), mlb(1000)
      real*8 delta(1000),steaz(1000),dmla(1000),resa(1000),resb(1000)
      real*8 dmlb(1000)
      integer*4 f_a(1000), f_b(1000)
c
      open (19,file=fnamml(1:lnblnk(fnamml)),status='unknown')
c
       elatd=elatx
       elond=elonx
c       write(*,*) 'st_dml - elat & elon = ',elatx,elonx
c
c ...... read station amplitude data
c
      call sadrd(fnamam,fnamml,epi_id,stn,slat,slon,scomp,samp,nstn)
c
c ...... compute distance and local magnitude
c
	call ml_calc(elatx,elonx,edepx,f_a,amla,semla,nmla,resa,
     1    f_b,amlb,semlb,nmlb,resb,mla,mlb,dmla,dmlb,delta,steaz,
     2    fnamam,fnamml,epi_id)
c
c ...... write *.ml file
c
      write(19,'(a)') fnamml(1:lnblnk(fnamml))
      do 7 i=1,nstn
      if((f_a(i).gt.0).and.(f_b(i).gt.0)) then
        write(19,'(a4,1x,a3,1x,f7.2,1x,f6.2,1x,1pe10.4,1x,0pf6.3,1x,
     1            f6.4,1x,f6.3,1x,f6.3,1x,f6.3)')
     3            stn(i),scomp(i),delta(i),steaz(i),samp(i),
     4            dlog10(samp(i)),floga0(delta(i)),dmla(i),mla(i),
     5            resa(i)
      elseif(f_a(i).gt.0) then
c        write(19,'(a4,1x,a3,1x,f7.2,1x,f6.2,1x,1pe10.4,1x,0pf6.3,1x,
c     1            f6.4,1x,f6.3,1x,f6.3,1x,f6.3)')
c     2            stn(i),scomp(i),delta(i),steaz(i),samp(i),
c     3            dlog10(samp(i)),floga0(delta(i)),dmla(i),mla(i),
c     4            resa(i)
      elseif(f_b(i).gt.0) then
c        write(19,'(a4,1x,a3,1x,f7.2,1x,f6.2,1x,1pe10.4,1x,0pf6.3,1x,
c     1            f6.4,22x,f6.3,1x,f6.3,1x,f6.3)')
c     3            stn(i),scomp(i),delta(i),steaz(i),samp(i),
c     4            dlog10(samp(i)),floga0(delta(i)),
c     5            dmlb(i),mlb(i),resb(i)
      endif
    7 continue
      write(19,'(a,f4.2,a,f5.3,a,i2,a)') 
     1      'ML(brk ) = ',amlb,' +/- ',semlb,' (',nmlb,')'
c      write(19,'(a,f4.2,a,f5.3,a,i2,a)') 
c     1      'ML(bdsn) = ',amla,' +/- ',semla,' (',nmla,')'
c      
      close (19)
c
      return
      end

	subroutine ml_calc(elatx,elonx,edepx,f_a,amla,semla,nmla,resa,
     1    f_b,amlb,semlb,nmlb,resb,mla,mlb,dmla,dmlb,delta,steaz,
     2    fnamam,fnamml,epi_id)
c
c ...... added to calculate bdsn & brk ML using station ML adjustments
c
	implicit none
c
	include "config.h"
c
	character*80 fnamam,fnamml
	character*9 epi_id
c
	integer*4 f_a(1000), f_b(1000), nmla, nmlb, i
	character*4 stn(1000), stn_t
	character*3 scomp(1000)
	real*8 ml_t, amla, amlb, rdistx
	real*8 slat(1000), slon(1000), samp(1000), mla(1000), mlb(1000)
	real*8 delta(1000), steaz(1000), dmla(1000), dmlb(1000)
        real*8 resa(1000),resb(1000), varmla, varmlb 
	real*8 elatx, elonx, edepx, semla, semlb, elatd, elond, floga0
	integer*4 nstn, i, lnblnk
c
	elatd=elatx
	elond=elonx
c
c ...... read station amplitude data
c
	call sadrd(fnamam,fnamml,epi_id,stn,slat,slon,scomp,samp,
     1             nstn)
c
	do 50 i=1,nstn
	  f_a(i)=0
	  f_b(i)=0
   50   continue
	if(nstn.lt.1) go to 20
c
	open(22,file=mladj_bdsn(1:lnblnk(mladj_bdsn)),status='old')
c
	amla=0.d0
	nmla=0
	do 4 i=1,nstn
c
c ...... compute distance and local magnitude
c
	call dist(elatd,elond,slat(i),slon(i),delta(i),steaz(i))
	rdistx=edepx*edepx+delta(i)*delta(i)
	if(rdistx.lt.0.1d0) then
	  rdistx=0.1d0
	else
	  rdistx=dsqrt(rdistx)
	endif
c
c ...... use attenuation relation tabulated from Nordquist Diagram 
c        ('42 BSSA)
c
	mla(i) = dlog10(samp(i)) + floga0(delta(i))
c
c ...... add BDSN station ML adjustments
c
	rewind(22)
    1   read(22,'(a4,1x,f6.3)',end=2) stn_t, ml_t
	if(stn(i).eq.stn_t) then
	  nmla=nmla+1
	  f_a(i)=1
	  dmla(i)=ml_t
	  mla(i)=mla(i)+dmla(i)
	  amla=amla+mla(i)
c	  write(*,*) 'bdsn - ',stn(i),' ',nmla,' ',f_a(i),' ',
c     1      dmla(i),' ',mla(i)
	else
	  go to  1
	endif
	go to 3
    2   continue
	nmla=nmla+1
	f_a(i)=1
	dmla(i)=0.d0
	mla(i)=mla(i)+dmla(i)
	amla=amla+mla(i)
    3   continue
c
    4   continue
	close(22)
	amla=amla/dble(nmla)
c
	semla=0.d0
	if(nstn.gt.2) then
	  do 5 i=1,nstn
	    if(f_a(i).gt.0) then
	      resa(i)=mla(i)-amla
	      varmla=varmla+resa(i)*resa(i)
	    endif
    5     continue
	  semla=dsqrt(varmla/dble(nmla*(nmla-1)))
	endif
c	write(*,*) 'amla, semla & nmla = ',amla, semla,' ',nmla
c
	open(23,file=mladj_brk (1:lnblnk(mladj_brk )),status='old')
c
	amlb=0.d0
	nmlb=0
	do 13 i=1,nstn
c
c ...... compute distance and local magnitude
c
	call dist(elatd,elond,slat(i),slon(i),delta(i),steaz(i))
	rdistx=edepx*edepx+delta(i)*delta(i)
	if(rdistx.lt.0.1d0) then
	  rdistx=0.1d0
	else
	  rdistx=dsqrt(rdistx)
	endif
c
c ...... use attenuation relation tabulated from Nordquist Diagram 
c        ('42 BSSA)
c
	mlb(i) = dlog10(samp(i)) + floga0(delta(i))
c
c ...... add BRK station ML adjustments
c
	rewind(23)
   11   read(23,'(a4,1x,f6.3)',end=12) stn_t, ml_t
	if(stn(i).eq.stn_t) then
	  nmlb=nmlb+1
	  f_b(i)=1
	  mlb(i) = dlog10(samp(i)) + floga0(delta(i))
	  dmlb(i)=ml_t
	  mlb(i)=mlb(i)+dmlb(i)
	  amlb=amlb+mlb(i)
c	  write(*,*) 'brk - ',stn(i),' ',nmlb,' ',f_b(i),' ',
c     1      dmlb(i),' ',mlb(i)
	else
	  go to 11
	endif
   12   continue
c
   13   continue
	close(23)
	amlb=amlb/dble(nmlb)
c
	varmlb=0.d0
	if(nstn.gt.2) then
	  do 14 i=1,nstn
	    if(f_b(i).gt.0) then
	      resb(i)=mlb(i)-amlb
	      varmlb=varmlb+resb(i)*resb(i)
	    endif
   14     continue
	  semlb=dsqrt(varmlb/dble(nmlb*(nmlb-1)))
	endif
c	write(*,*) 'amlb, semlb & nmlb = ',amlb, semlb,' ',nmlb
c
   20   continue
c
	return
	end

      subroutine sadrd(fnamam,fnamml,epi_id,stn,slat,slon,
     1                 scomp,samp,nstn)
c
c ...... read in station amplitude data
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      character*80 fnamam,fnamml,aline
      character*9 pattern,epi_id
      character*4 stn(1000)
      character*3 scomp(1000)
      real*8 slat(1000),slon(1000),samp(1000)
c      character*2 dummy
c
      pattern=epi_id(1:9)
      lml=lnblnk(fnamml)
      open(16,file=fnamam(1:lnblnk(fnamam)),status='old')
      nstn=0
    2 continue
      read(16,'(a)',end=3) aline
      nstn=nstn+1
      if(aline(1:1).eq.'#') then
        nstn=nstn-1
        go to 2
      endif
      call aldec(aline,stn(nstn),scomp(nstn),samp(nstn),ans)
      if(ans.lt.0.) then
        nstn=nstn-1
        go to 2
      endif
c	write(*,*) 'nstn, stn, scomp, samp = ',
c     1    nstn, " ", stn(nstn), scomp(nstn), samp(nstn)
c
c ...... include 2 & 3 as well as N & E orientations (28 Mar 97)
c
c      if((scomp(nstn)(3:3).ne.'N').and.
c     1   (scomp(nstn)(3:3).ne.'E')) then
c        nstn=nstn-1
c        go to 2
c      endif
c
	if((scomp(nstn)(3:3).ne.'N').and.
     1     (scomp(nstn)(3:3).ne.'E').and.
     2     (scomp(nstn)(3:3).ne.'2').and.
     3     (scomp(nstn)(3:3).ne.'3'))  then
	  nstn=nstn-1
	  go to 2
	endif
c
      call slfind(stn(nstn),slat(nstn),slon(nstn),ans)
      if(ans.lt.0.) then
        write(*,*) 'Warning - ',stn(nstn),' not in list'
        nstn=nstn-1
      endif
c      write(*,*) 'sadrd - nstn, stn, slat, slon, scomp, samp = ',
c     1      nstn,stn(nstn),slat(nstn),slon(nstn),scomp(nstn),
c     2      samp(nstn)
      go to 2
c
    3 continue
	call sd_chk(nstn,stn,scomp,samp,slat,slon)
      close(16)
      return
    4 continue
      close(16)
      write(*,*) epi_id,' not found in ampli.file'
      stop
      end

	subroutine sd_chk(nstn,stn,scomp,samp,slat,slon)
c
c ...... check for invalid stations and multiple amplitude entries
c
	implicit none
	integer*4 nstn, i, j, n, ns, ans, indx(100)
	real*8 samp(*), slat(*), slon(*)
	character*4 stn(*)
	character*3 v_stn(8), t_stn(100)
	character*3 scomp(*), t_scomp(100)
	real*8 t_samp(100), t_slat(100), t_slon(100)
c
	data v_stn /'BRK','BKS','MHC','MIN','ARC','PAC','SFB','VIN'/
	data ns /8/
c
	if(nstn.lt.1) return
	if(nstn.eq.1) then
	  ans=0
	  do 1 i=1,ns
	    if(stn(1)(1:3).eq.v_stn(i)) ans=1
    1     continue
	  if(ans.lt.1) nstn=0
	  return
	endif
	n=0
	do 3 i=1,nstn
	  do 2 j=1,ns
	    if(stn(i)(1:3).eq.v_stn(j)) then
	      n=n+1
	      t_stn(n)=stn(i)(1:3)
	      t_scomp(n)=scomp(i)
	      t_samp(n)=samp(i)
	      t_slat(n)=slat(i)
	      t_slon(n)=slon(i)
	      indx(n)=1
	    endif
    2     continue
    3   continue
	do 5 i=1,n-1
	  if(indx(i).ne.0) then
	    do 4 j=i+1,n
	      if(indx(j).ne.0) then
	        if(t_stn(i).eq.t_stn(j)) then
	          if(t_scomp(i).eq.t_scomp(j)) then
	            if(t_samp(i).gt.t_samp(j)) then
	              indx(j)=0
	            else
	              indx(i)=0
	            endif
	          endif
	        endif
	      endif
    4       continue
	  endif
    5   continue
	j=0
	do 6 i=1,n
	  if(indx(i).eq.1) then
	    j=j+1
	    stn(j)=t_stn(i)//' '
	    scomp(j)=t_scomp(i)
	    samp(j)=t_samp(i)
	    slat(j)=t_slat(i)
	    slon(j)=t_slon(i)
	  endif
    6   continue
	nstn=j
c
	return
	end

      subroutine fnxtr1(fninpt,epi_id)
c
c ...... extract date and time from pickfile name
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      character*80 fninpt
      character*9 epi_id
      character*1 achf
c
c      write(*,*) 'fnxtr1 - fninpt = ',fninpt(1:lnblnk(fninpt))
      do 1 i=1,9
      epi_id(i:i)=' '
    1 continue
      j=10
      lfa=lnblnk(fninpt)
      i=lfa+1
    2 i=i-1
      if(i.lt.1) then
        write(*,*) 'ILLEGAL file name - ',fninpt(1:lfa)
        stop
      endif
      if(j.eq.1) go to 3
      achf=fninpt(i:i)
      ichn=ichar(achf)
      if((ichn.ge.48).and.(ichn.le.57)) then
        j=j-1
        epi_id(j:j)=achf
      endif
      go to 2
    3 continue
c      write(*,*) 'fnxtr1 - epi_id = ',epi_id(1:lnblnk(epi_id))
c
      return
      end

      subroutine aldec(aline,stn,scomp,samp,ans)
c
c ...... decode line from seistool ampli.file
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
	integer*4 npar, lpar(20)
      character*80 aline
      character*4 stn
      character*3 scomp
	character*20 par(20)
      ans=-1.
	call parser(aline,npar,lpar,par)
	read(par(1)(1:3),'(a)') stn
	read(par(1)(4:6),'(a)') scomp
	read(par(3)(1:lpar(3)),'(f)') samp
      ans=+1.
c
c      write(*,*) 'scomp = ',scomp
c      write(*,*) 'samp = ',samp
c      write(*,*) 'ans = ',ans
c
      return
      end

	subroutine parser(aline,npar,lpar,par)
c
c ...... parse aline
c
	implicit none
	integer*4 i, j, lnb, lnblnk, itb(20), ite(20)
	integer*4 npar, lpar(20)
	character*80 aline
	character*20 par(20)
	character*1 tab, blank, quotes
c
	data tab,blank,quotes /Z'09',' ','"'/
c
	do 10 i=1,20
	  par(i)='                    '
	  lpar(i)=0
   10   continue
	lnb=lnblnk(aline)
	npar=0
	j=0
    1   j=j+1
	if(j.gt.lnb) go to 4
	if((aline(j:j).eq.tab).or.(aline(j:j).eq.blank)) go to 1
	npar=npar+1
	itb(npar)=j
	if(aline(j:j).eq.quotes) then
    2     j=j+1
	  if(aline(j:j).ne.quotes) go to 2
	  ite(npar)=j
	  go to 1
	endif
    3   j=j+1
	if(j.gt.lnb) then
	  ite(npar)=j-1
	  go to 4
	endif
	if((aline(j:j).ne.tab).and.(aline(j:j).ne.blank)) go to 3
	ite(npar)=j-1
	go to 1
    4   continue
	do 5 i=1,npar
	  read(aline(itb(i):ite(i)),'(a)') par(i)
	  lpar(i)=lnblnk(par(i))
    5   continue
c
	end

      subroutine slfind(stn,slat,slon,ans)
c
c ...... extract station location from instr_coord file
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
c
      include "config.h"
      character*4 stn,astn
c
      ans=-1.d0
      open(17,file=instr_coord(1:lnblnk(instr_coord)),status='old')
    1 read(17,'(a4,1x,2f9.4)',end=2) astn,slat,slon
      if(astn.ne.stn) go to 1
      slat=(datand(0.993277d0*dtand(slat)))
      ans=+1.d0
    2 continue
      close(17)
c
      return
      end

      real*8 function floga0(delta)
c
************************************************************************
*   DEFINE LOG(Ao) FUNCTION
*   Based on linear interpolation of Nordquist Nomograph (BSSA, 1948)
************************************************************************
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 d(19),a(19)
      data a /1.4d0,1.4d0,1.5d0,1.7d0,2.3d0,2.7d0,2.85d0,
     1        3.0d0,3.65d0,4.2d0,4.6d0,5.0d0,5.3d0,5.7d0,
     2        6.0d0,6.2d0,6.3d0,6.4d0,6.5d0/
      data d /0.d0,5.d0,10.d0,20.d0,35.d0,55.d0,73.d0,96.d0,
     1        220.d0,334.d0,445.d0,600.d0,750.d0,1010.d0,
     2        1340.d0,1600.d0,1750.d0,1960.d0,2230.d0/
      floga0=0.d0
      if ((delta.LT.0.d0).OR.(delta.GT.2224.d0)) return
      j=0
      do 1 i=1,19
      j=j+1
      if (delta.LT.d(i)) goto 2
    1 continue
      return
    2 j0=j-1
      floga0=a(j0)+(a(j)-a(j0))*(delta-d(j0))/(d(j)-d(j0))
      return
      end

      subroutine dist(elat,elon,slat,slon,delta,steaz)
c
c ...... compute distance & azimuth from station to epicenter
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
c
      dy = 111.195d0*(elat-slat)
      dx = 111.195d0*(elon-slon)*dcosd((elat+slat)/2.d0)
      delta=dx*dx+dy*dy
      nnn=59
      if(delta.gt.0.d0) delta=dsqrtf(nnn,delta)
      if((dx.eq.0.d0).and.(dy.eq.0.d0)) then
        steaz=0.d0
      else
        steaz=datan2d(dx,dy)
      endif
      if(steaz.lt.0.d0) steaz=steaz+360.d0
c
      return
      end

	subroutine bdsn_ml(stn,delta,wa_amp,bml,ml_logAo,ml_adj,ierr)
c
c ...... determine ML using BDSN dML & -logAo
c
c	input:
c		stn    - BDSN station code (a4)
c		delta  - distance (km)
c		wa_amp - WA maximum trace amplitude (mm)
c	output:
c		bml      - BDSN estimate of local magnitude
c		ml_logAo - -logAo
c		ml_adj   - ML station adjustment
c		ierr     - error flag
c			     ierr = 0 - returned if no error
c			       	    1 - no station adjustment used
c			            2 - distance greater than 1200 km
c
	implicit none
	integer*4 i, ii, nd, ierr, nbs
	real*4 delta, wa_amp, bml, dladd
	real*4 dml, ml_adj, ml_logAo
	real*4 logAo
	character*4 stn, bdsn_stn
	common /bdsnml/ logAo(241), dml(13), bdsn_stn(13)
c
	data nd, nbs / 241, 13/
c
	ierr=0
	ii=int(delta/5.)+1
	if(ii.gt.nd) then
	  ierr=2
	  return
	endif
	ml_adj=0.
	do 1 i=1,nbs
	  if(stn.eq.bdsn_stn(i)) then
	    ml_adj=dml(i)
	    go to 2
	  endif
    1   continue
	ierr=1
    2   continue
	dladd=(logAo(ii+1)-logAo(ii))/5.
	ml_logAo = logAo(ii) + dladd*(delta-real(5*(ii-1)))
	bml = ml_logAo + alog10(wa_amp) + ml_adj
c
	return
	end

	block data ml_pars
c
	real*4 logAo, dml
	character*4 bdsn_stn
	common /bdsnml/ logAo(241), dml(13), bdsn_stn(13)
c
	data bdsn_stn /'ARC ', 'BKS ', 'CMB ', 'HOPS', 'JRSC',
     1                 'MHC ', 'MIN ', 'ORV ', 'PKD1', 'SAO ',
     2                 'STAN', 'WDC ', 'YBH '/
	data dml /+0.209, -0.035, +0.244, +0.325, +0.146,
     1            +0.128, -0.107, +0.416, -0.219, +0.305,
     2            -0.225, +0.451, +0.519/
c
	data (logAo(i),i=1,101) / 1.489,
     1 1.489,1.588,1.685,1.782,1.976,2.168,2.359,2.448,2.537,2.625,
     2 2.713,2.744,2.776,2.809,2.841,2.870,2.901,2.934,2.969,3.000,
     3 3.031,3.064,3.098,3.132,3.167,3.202,3.237,3.272,3.307,3.341,
     4 3.374,3.407,3.439,3.470,3.500,3.530,3.560,3.589,3.617,3.645,
     5 3.672,3.699,3.725,3.751,3.775,3.798,3.821,3.844,3.866,3.889,
     6 3.911,3.933,3.955,3.976,3.998,4.020,4.041,4.063,4.085,4.107,
     7 4.129,4.151,4.173,4.195,4.218,4.240,4.261,4.278,4.294,4.311,
     8 4.328,4.344,4.361,4.378,4.395,4.412,4.429,4.446,4.463,4.480,
     9 4.497,4.515,4.532,4.549,4.567,4.584,4.602,4.619,4.637,4.649,
     1 4.662,4.674,4.687,4.699,4.712,4.725,4.737,4.750,4.762,4.775/
	data (logAo(i),i=102,201) /
     1 4.788,4.800,4.813,4.826,4.839,4.851,4.864,4.877,4.889,4.902,
     2 4.915,4.927,4.940,4.952,4.965,4.978,4.990,5.003,5.015,5.028,
     3 5.038,5.047,5.057,5.067,5.076,5.086,5.096,5.105,5.115,5.125,
     4 5.134,5.144,5.153,5.163,5.173,5.182,5.192,5.201,5.211,5.220,
     5 5.230,5.240,5.249,5.259,5.268,5.278,5.287,5.297,5.306,5.316,
     6 5.323,5.330,5.337,5.345,5.352,5.359,5.366,5.373,5.380,5.388,
     7 5.395,5.402,5.409,5.416,5.423,5.431,5.438,5.445,5.452,5.459,
     8 5.466,5.473,5.480,5.487,5.495,5.502,5.509,5.516,5.523,5.530,
     9 5.537,5.544,5.551,5.558,5.565,5.572,5.579,5.587,5.594,5.601,
     1 5.608,5.615,5.622,5.629,5.636,5.643,5.650,5.657,5.664,5.671/
	data (logAo(i),i=202,241) /
     7 5.678,5.685,5.689,5.693,5.697,5.701,5.704,5.708,5.712,5.716,
     8 5.720,5.724,5.728,5.731,5.735,5.739,5.743,5.747,5.751,5.755,
     9 5.758,5.762,5.766,5.770,5.774,5.778,5.781,5.785,5.789,5.793,
     1 5.797,5.801,5.804,5.808,5.812,5.816,5.820,5.824,5.827,5.831/
c
	end