Computing Surface-Wave (MS) Magnitude

The following FORTRAN code is based on a subroutine used at the USGS NEIC that was provided by Ray Buland.  A text copy is also located here:  msmag.for  and the executable version is msmag.exeBut even better, Bob McClure has converted the DOS program to a Visual Basic Windows program.  Just download and unzip MagCalc.zip and run the program MagCalc.exe.

The graph on this page, msmag.html, indicates how the table values compare with those computed from the formula 1.66 alog10(delta) for the distance range 130 to 180 degrees.  The equation does not account for the increasing amplitude due to convergence of the surface waves near the antipode.

c      subroutine msmag(t,amp,msflg,ms,depth,delta) 
c
c   Msmag returns an MS station magnitude, ms, for a reading with
c   period t (s), amplitude amp (um), and flag msflg from an event at
c   depth depth (km), and distance delta (deg).  Note that the bsdist
c   table (for distances over 160 deg) was taken for QQUAKE.
c
c     character*1 msflg

      real ms
      dimension bsdist(50)
c TABLE CONTAINING Bs VALUES TO COMPUTE MS MAG AT DISTANCES
c OVER 130 DEGREES.
      data bsdist/  3.49,3.5,3.5,3.5,3.51,3.51,3.51,3.51,3.52,3.52,
     1 3.52,3.52,3.52,3.53,3.53,3.53,3.53,3.53,3.54,3.54,3.54,3.54,
     2 3.54,3.54,3.54,3.54,3.54,3.54,3.54,3.54,3.54,3.53,3.53,3.53,
     3 3.52,3.52,3.52,3.52,3.51,3.51,3.51,3.50,3.49,3.47,3.44,3.41,
     4 3.39,3.34,3.29,3.19/
c
        do i  = 1, 50
          print *, bsdist(i)
        enddo

1       print *, ' Input the amplitude, period, and distance'
        read *, amp, t, delta
        if(amp.eq.0) then
          print *, 'amplitude may not be zero.'
          stop
        else if(delta .lt. 20. .or. delta .gt. 180.) then
          print *, 'distance must be between 20 and 180 degrees,',
     *      ' inclusive.'
          stop
        else if(t.le.9..or.t.gt.99.) then
          print *,
     *      'period must be greater than 9 and ',
     *      'less than or equal 99 seconds.'
          stop
        else if(delta.le.130.) then
          ms=alog10(amp/t)+1.66*alog10(delta)+3.3
        else
          ms=alog10(amp/t)+bsdist(int(delta-129.5))+3.3
        endif

        print *, ' Ms = ', ms
        goto 1
c     if((t.lt.18..or.t.gt.22..or.delta.gt.160.).and.msflg.ne.'R')
c     1 msflg='r'
c      return
c
 10   ms=0.
      msflg='n'
      return
      end

Back