c based on sedas subroutine mbmag (jcl 3/30/2005)
      	dimension q(108,17)
	real mb
      	logical begin
      	data begin/.true./
	
	mb = 0.

1       print *, 'Compute mb, given distance (deg), depth (km),',
     *    ' amplitude (microns) and period (seconds) of P arrival.'
        print *, 'What is the distance, depth, amplitude & period?'
        read(5, *) delta, depth, amp, t

      	if(begin) then
          begin = .false.
          open(unit=1, status='old',file='qtab.txt') 
	  do i=1,17
	    read(1,20)(q(j,i),j=1,108)
20 	    format(4x,32f4.1)
	  enddo
	  close(1)
      	endif
c
	if(t .lt. .1 .or. t.gt.9. .or. 
     *    amp .le. 0. .or.
     *    (delta.lt.5. .or. delta.gt.109.).or.
     *    depth.gt.800.) go to 11

	if(.true.) goto 5

c older interpolation code to find qval
      	if(depth.lt.100.) then
          k=int(depth/25.)+2
          s1=.04*(depth-25.*(k-2))
      	else
          k=min0(int(depth/50.)+4,17)
          s1=.02*(depth-50.*(k-4))
      	endif
c
	j=max0(min0(int(delta),108),2)
       	s2=delta-j
c
c	print *,'Mbmag:  depth k s1 delta j s2 =',depth,k,s1,delta,j,s2
	q1=q(j-1,k-1)+s1*(q(j-1,k)-q(j-1,k-1))
      	qval=q1+s2*((q(j,k-1)+s1*(q(j,k)-q(j,k-1)))-q1)

c newer interpolation code to find qval
5      	do k=1,17
      	  if(k.gt.5) go to 6
      	  dep=(k-1)*25
      	  if(depth.ge.dep) go to 7
      	  s1=(depth+25.0-dep)*0.04
      	  go to 8
6 	  dep=(k-3)*50
      	  if(depth.gt.dep) go to 7
      	  s1=(depth+50.0-dep)*0.02
      	  go to 8
7     	enddo

      	k=17
      	s1=(depth+50.0-700.)*0.02
8 	do j=1,108
      	  del=j+1
      	  if(delta.ge.del) go to 9
      	  s2=delta+1.0-del
      	  go to 10
9     	enddo
      	j=108
      	s2=0.0
10 	q1=q(j-1,k-1)+s1*(q(j-1,k)-q(j-1,k-1))
      	q2=q(j,k-1)  +s1*(q(j,k)  -q(j,k-1))
      	qval=q1+s2*(q2-q1)

      	mb=alog10(amp/t)+qval
	print *, 'qval = ', qval
c      	amp = 10.**(mb - qval +3)
      	write(6, '(a, f5.1)') ' mb via NEIC formulation = ', mb

	mb = alog10(amp/t) + 0.01*delta + 5.9
	write(6, '(a, f5.1)') ' mb via Braile equation = ', mb
	goto 1
11 	amp1 = 0.0
     	print *, 'Distance must be between 5 and 109 degrees.'
	print *, 'Period must be between 0.1 and 9.0 seconds.'
	print *, 'Amplitude must be greater than zero.'
      	stop
      	end