      SUBROUTINE HYMAG
C--CALCULATE MAGNITUDES FOR ALL STATIONS FOR HYPOINVERSE

      IMPLICIT NONE

      INCLUDE 'common.inc'

      REAL      RCAL(7), RSPA(26), FMCAR(6), XMCAR(6), S(7)
      INTEGER   IATN(7), IAEXP(7), IYEARI(7)
      CHARACTER STN*5,SNET*2,SCOMP*3, SLOC*2

C--ARRAYS FOR GETTING WEIGHTED MEDIAN MAGNITUDES
C     INTEGER*2 IMAG(MAXPHS),IMWT(MAXPHS)
C     INTEGER*2 IMAG2(MAXPHS),IMWT2(MAXPHS)
      INTEGER   IMAG(MAXPHS),IMWT(MAXPHS)
      INTEGER   IMAG2(MAXPHS),IMWT2(MAXPHS)

      REAL      A0MAG, AMPMM, CAL, FQ, G, GG
      REAL      SMAG, SMAGB, SMAGF, SMAGT
      REAL      PER, SLDIS
      REAL      TEMP, TAU, TLOG
      INTEGER   I, J, K, IM, IFQ, IOS
      INTEGER   LATYP, NDATE
      INTEGER   KTEMP, KFMWT, KXMWT
      INTEGER   MEDIAN, MSD, NMED, NMED2

C--MOTOR CONSTANTS FOR VARIOUS VELOCITY SEISMOMETERS IN V/CM/SEC.
C  INDEX (SUBSCRIPT) IS THE SEISMOMETER TYPE CODE JTYPE(J)
      DATA S/   1.0,    !1 L4C, ATTENUATION HISTORY
     2          1.0,    !2 UNUSED, ACTUALLY A SPRENGNETHER
     3          1.0,    !3 L4C, CAL FACTOR HISTORY
     4          0.448,  !4 HS1 (USED WITH NANOMETRICS)
     5          8.0,    !5 GURALP (USED WITH REFTEKS)
     6          24.0,   !6 STREKHEISEN STS-1 (USED WITH BDSN)
     7          15.0/   !7 STREKHEISEN STS-2 (USED WITH BDSN)
      SAVE S

C--RESPONSE CURVE OF THE USGS STANDARD HIGH GAIN (L4C 1 SEC.) RELATIVE TO WA
C  WOOD ANDERSON HAS MAGNIFICATION 2080 & DAMPING 0.8
C     FREQUENCY  .16   .20   .25   .32   .40   .50   .63   .79   1.00
C     LOG FREQ   -.8   -.7   -.6   -.5   -.4   -.3   -.2   -.1   0.0
      DATA RSPA /.288, .432, .561, .680, .786, .891, .983, 1.066,1.138,

C FREQUENCY 1.26   1.59   2.00   2.51   3.16   3.98   5.01   6.31   7.94
C LOG FREQ   .1     .2     .3     .4     .5     .6     .7     .8     .9
     2    1.205, 1.276, 1.355, 1.443, 1.535, 1.630, 1.726, 1.822, 1.916,

C FREQUENCY 10.0   12.6   15.9   20.0   25.1   31.6   39.8   50.1
C LOG FREQ  1.0    1.1    1.2    1.3    1.4    1.5    1.6    1.7
     3    2.007, 2.090, 2.145, 2.099, 1.878, 1.546, 1.172, .771/
      SAVE RSPA

C--INITIALIZE SUMS & VALUES
C     MFMAG=0
C     MXMAG=0
      NFMAG=0
      NXMAG=0
      FMMAD=0.
      XMMAD=0.
      XMAG=VNOMAG
      FMAG=VNOMAG
      NMED=0
C     MFMAG2=0
C     MXMAG2=0
      NFMAG2=0
      NXMAG2=0
      FMMAD2=0.
      XMMAD2=0.
      XMAG2=VNOMAG
      FMAG2=VNOMAG
      NMED2=0
      NDATE=KYEAR2*1000000 +KMONTH*10000 +KDAY*100 +KHOUR

C--OPTIONALLY SET FMAG2 FLAG FOR ALL STATIONS IF EPI IN LAT/LON BOX (MAR COMMAND)
      IF (USEMAR) THEN
        IF (((YLATMX-CLAT)*(CLAT-YLATMN)).GE.0. .AND.
     2  ((XLONMX-CLON)*(CLON-XLONMN)).GE.0.) THEN
          LINBOX=.TRUE.
          DO J=1,JSTA
            JFM2(J)=.TRUE.
          END DO
        ELSE
          LINBOX=.FALSE.
          DO J=1,JSTA
            JFM2(J)=.FALSE.
          END DO
        END IF
      END IF

C--LOOP OVER STATIONS TO UPDATE FMAG CORRECTIONS, CAL FACTORS & CALCULATE FMAGS
      DO 50 K=1,KSTA
      KFMAG(K)=MINMAG           !IS THIS STATEMENT NECESSARY?
      J=KINDX(K)

C++++++++++++++++++++++++++++++++++++++++++++++++++++++
C--IF A FILE OF FMAG CORRECTIONS WITH EXPIRATION DATES WAS READ, CHECK TO SEE IF
C  THE ONE FOR THIS STATION HAS EXPIRED.
      IF (JFEXP(J).GT.0) THEN
        IF (NDATE.GT.JFEXP(J)) THEN

C--READ FILE TO GET A NEW EXPIRATION DATE & MAG CORRECTION. OPEN
C  CORRECTION FILE & SEARCH FOR THIS STATION. END OF FILE SHOULD NEVER OCCUR.
          CALL OPENR (13,FMCFIL,'F',IOS)

2         IF (L2000) THEN
            READ (13,1202) STN,SNET,SLOC,SCOMP,
     2      (FMCAR(I),IYEARI(I),IAEXP(I),I=1,6)
1202        FORMAT (A5,1X,A2,A2,A3,1X, 6(F5.2,1X,I4,I6,1X))

          ELSE
            READ (13,1002) STN,SNET,SCOMP,
     2      (FMCAR(I),IYEARI(I),IAEXP(I),I=1,6)
C--READ 2-DIGIT YEAR, MO, DA, HOUR, AND IGNORE MINUTES IF PRESENT
1002        FORMAT (A5,1X,A2,2X,A3,1X, 6(F5.2,1X,I2,I6,3X))
            SLOC='  '
          END IF

          IF (STN(1:NSTLET) .NE. STANAM(J)(1:NSTLET) .OR.
     2    SNET(1:NETLET) .NE. JNET(J)(1:NETLET) .OR.
     3    (SLOC(1:NSLOC2) .NE. JSLOC(J)(1:NSLOC2) .AND.
     3    SLOC(1:NSLOC2) .NE. JSLOC2(J)(1:NSLOC2)) .OR.
     4    SCOMP(1:NCOMP) .NE. JCOMP3(J)(1:NCOMP)) GOTO 2

C--SEARCH THE LIST OF EXPIRATION DATES FOR THE FIRST ONE AFTER THE CURRENT DATE,
C  OR AN EXPIRATION DATE OF 0 (FMAG COR GOOD THROUGH THE FUTURE).
          DO I=1,6
            IF (IYEARI(I).LT.100 .AND. IYEARI(I).GT.0)
     2      IYEARI(I)=IYEARI(I)+ICENT
C--TRUNCATE YEAR TO 2146, THE LARGEST STORABLE IN I*4, STORE WHOLE DATE
            IF (IYEARI(I).GT.2146) IYEARI(I)=2146
            IAEXP(I)=IYEARI(I)*1000000 +IAEXP(I)

            IF (IAEXP(I).EQ.0 .OR. IAEXP(I).GT.NDATE) THEN

C--ADD 10 TO CORRECTION TO GIVE GAIN CORRECTION 0 WEIGHT
C--YOU CAN ADD BOTH 10 AND 5, BUT MUST DO THIS TEST FIRST
              IF (FMCAR(I).GT.7.45) THEN
                JFGWT(J)=0
                FMCAR(I)=FMCAR(I)-10.
              ELSE
                JFGWT(J)=1
              END IF

C--ADD 5 TO CORRECTION TO GIVE IT ZERO WEIGHT
              IF (FMCAR(I).GT.2.45) THEN
                JFWT(J)=0
                FMCAR(I)=FMCAR(I)-5.
              ELSE
                JFWT(J)=10
              END IF
              JFCOR(J)=NINT(100.*FMCAR(I))

C--WRITE MESSAGE THAT A NEW FMC WAS FOUND
              IF (LPRT .AND. KPRINT.GE.3) WRITE(15,1003)
     2          FMCAR(I),STN,SNET,SCOMP,SLOC, JFEXP(J),IAEXP(I)
1003          FORMAT (' * NEW FMAG CORRECTION',F6.2,' ASSIGNED TO ',
     2        A5,'-',A2,'-',A3,'-',A2,'.'/
     3        ' NEW START DATE = ',I10,', EXPIRATION DATE = ',I10)
              JFEXP(J)=IAEXP(I)
              GOTO 3
            END IF
          END DO
3         CLOSE (13)
        END IF
      END IF

C++++++++++++++++++++++++++++++++++++++++++++++++++++++
C--IF A FILE OF XMAG CORRECTIONS WITH EXPIRATION DATES WAS READ, CHECK TO SEE IF
C  THE ONE FOR THIS STATION HAS EXPIRED.
      IF (JXEXP(J).GT.0) THEN
        IF (NDATE.GT.JXEXP(J)) THEN

C--READ FILE TO GET A NEW EXPIRATION DATE & MAG CORRECTION. OPEN
C  CORRECTION FILE & SEARCH FOR THIS STATION. END OF FILE SHOULD NEVER OCCUR.
C  IGNORE STATION TYPE BECAUSE WE GOT IT ON FIRST READING
          CALL OPENR (13,XMCFIL,'F',IOS)

12        IF (L2000) THEN
            READ (13,1212) STN,SNET,SLOC,SCOMP,
     2      (XMCAR(I),IYEARI(I),IAEXP(I),I=1,6)
1212        FORMAT (A5,1X,A2,A2,A3, 3X, 6(F5.2,1X,I4,I6,1X))

          ELSE
            READ (13,1012) STN,SNET,SCOMP,
     2      (XMCAR(I),IYEARI(I),IAEXP(I),I=1,6)
1012        FORMAT (A5,1X,A2,2X,A3, 3X, 6(F5.2,1X,I2,I6,3X))
            SLOC='  '
          END IF

          IF (STN(1:NSTLET) .NE. STANAM(J)(1:NSTLET) .OR.
     2    SNET(1:NETLET) .NE. JNET(J)(1:NETLET) .OR.
     3    (SLOC(1:NSLOC2) .NE. JSLOC(J)(1:NSLOC2) .AND.
     3    SLOC(1:NSLOC2) .NE. JSLOC2(J)(1:NSLOC2)) .OR.
     4    SCOMP(1:NCOMP) .NE. JCOMP3(J)(1:NCOMP)) GOTO 12

C--SEARCH THE LIST OF EXPIRATION DATES FOR THE FIRST ONE AFTER THE CURRENT DATE,
C  OR AN EXPIRATION DATE OF 0 (XMAG COR GOOD THROUGH THE FUTURE).
          DO I=1,6
            IF (IYEARI(I).LT.100 .AND. IYEARI(I).GT.0)
     2      IYEARI(I)=IYEARI(I)+ICENT

C--TRUNCATE YEAR TO 2146, THE LARGEST STORABLE IN I*4, STORE WHOLE DATE
            IF (IYEARI(I).GT.2146) IYEARI(I)=2146
            IAEXP(I)=IYEARI(I)*1000000 +IAEXP(I)

            IF (IAEXP(I).EQ.0 .OR. IAEXP(I).GT.NDATE) THEN
C--ADJUST STATION WEIGHT DEPENDING ON RANGE OF CORRECTION
              IF (XMCAR(I).GT.2.45) THEN
                JXWT(J)=0
                XMCAR(I)=XMCAR(I)-5.
              ELSE
                JXWT(J)=10
              END IF
              JXCOR(J)=NINT(100.*XMCAR(I))

C--WRITE MESSAGE THAT A NEW XMC WAS FOUND
              IF (LPRT .AND. KPRINT.GE.3) WRITE(15,1013)
     2          XMCAR(I),STN,SNET,SCOMP,SLOC, JXEXP(J),IAEXP(I)
1013          FORMAT (' * NEW XMAG CORRECTION',F6.2,' ASSIGNED TO ',
     2        A5,'-',A2,'-',A3,'-',A2,'.'/
     3        ' NEW START DATE = ',I10,', EXPIRATION DATE = ',I10)
              JXEXP(J)=IAEXP(I)
              GOTO 13
            END IF
          END DO
13        CLOSE (13)
        END IF
      END IF

C++++++++++++++++++++++++++++++++++++++++++++++++++++++
C--IF USING A FILE OF CAL FACTORS WITH EXPIRATION DATES, CHECK TO SEE IF
C  CAL FOR THIS STATION HAS EXPIRED.  IF NOT USING A HISTORY FILE,
C  THE EXPIRATION DATE SHOULD ALWAYS BE 0.
      IF (JCEXP(J).GT.0) THEN
        IF (NDATE.GT.JCEXP(J)) THEN

C--FOR INSTRUMENT TYPE 1, USE ATTENUATION HISTORY FILE
          IF (JTYPE(J).EQ.1) THEN

C--READ FILE TO GET A NEW EXPIRATION DATE & ATTENUATION. OPEN ATTENUATION FILE
C  & SEARCH FOR THIS STATION. END OF FILE SHOULD NEVER OCCUR.
            CALL OPENR (13,ATNFIL,'F',IOS)

5           IF (L2000) THEN
              READ (13,1200) STN,SNET,SLOC,SCOMP,
     2        (IATN(I),IYEARI(I),IAEXP(I),I=1,7)
1200          FORMAT (A5,1X,A2,A2,A3,1X, 7(I2,1X,I4,I6,1X))

            ELSE
              READ (13,1000) STN,SNET,SCOMP,
     2        (IATN(I),IYEARI(I),IAEXP(I),I=1,7)
1000          FORMAT (A5,1X,A2,2X,A3,1X, 7(I2,1X,I2,I6,3X))
              SLOC='  '
            END IF

            IF (STN(1:NSTLET) .NE. STANAM(J)(1:NSTLET) .OR.
     2      SNET(1:NETLET) .NE. JNET(J)(1:NETLET) .OR.
     3      (SLOC(1:NSLOC2) .NE. JSLOC(J)(1:NSLOC2) .AND.
     3      SLOC(1:NSLOC2) .NE. JSLOC2(J)(1:NSLOC2)) .OR.
     4      SCOMP(1:NCOMP) .NE. JCOMP3(J)(1:NCOMP)) GOTO 5

C--SEARCH THE LIST OF EXPIRATION DATES FOR THE FIRST ONE AFTER THE CURRENT DATE,
C  OR AN EXPIRATION DATE OF 0 (CAL GOOD THROUGH THE FUTURE).
            DO I=1,7
              IF (IYEARI(I).LT.100 .AND. IYEARI(I).GT.0)
     2        IYEARI(I)=IYEARI(I)+ICENT

C--TRUNCATE YEAR TO 2146, THE LARGEST STORABLE IN I*4, STORE WHOLE DATE
              IF (IYEARI(I).GT.2146) IYEARI(I)=2146
              IAEXP(I)=IYEARI(I)*1000000 +IAEXP(I)

              IF (IAEXP(I).EQ.0 .OR. IAEXP(I).GT.NDATE) THEN
C--ATTEN MUST BE A MULTIPLE OF 6
                KTEMP=IATN(I)/6
                JCAL(J)=CALSV(KTEMP)*1000.

C--WRITE MESSAGE THAT A NEW ATTENUATION WAS FOUND
                IF (LPRT.AND.KPRINT.GE.3) WRITE (15,1001) CALSV(KTEMP),
     2          IATN(I), STN,SNET,SCOMP,SLOC, JCEXP(J), IAEXP(I)
1001            FORMAT (' * NEW CAL FACTOR',F6.3,' (ATTEN=',I2,
     2          ') ASSIGNED TO ',A5,'-',A2,'-',A3,'-',A2,'.'/
     4          ' NEW CAL FACTOR START DATE = ',
     3          I10,', EXPIRATION DATE = ',I10)
                JCEXP(J)=IAEXP(I)
                GOTO 6
              END IF
            END DO
6           CLOSE (13)
          END IF

C--FOR INSTRUMENT TYPES OTHER THAN 1, USE CAL FACTOR HISTORY FILE
          IF (JTYPE(J).NE.1) THEN

C--READ FILE TO GET A NEW EXPIRATION DATE & CAL FACTOR. OPEN ATTENUATION FILE
C  & SEARCH FOR THIS STATION. END OF FILE SHOULD NEVER OCCUR.
            CALL OPENR (13,CALFIL,'F',IOS)

7           IF (L2000) THEN
              READ (13,1207) STN,SNET,SLOC,SCOMP,
     2        (RCAL(I),IYEARI(I),IAEXP(I),I=1,7)
1207          FORMAT (A5,1X,A2,A2,A3,1X, 7(F7.2,1X,I4,I6,1X))

            ELSE
              READ (13,1007) STN,SNET,SCOMP,
     2        (RCAL(I),IYEARI(I),IAEXP(I),I=1,7)
1007          FORMAT (A5,1X,A2,2X,A3,1X, 7(F7.2,1X,I2,I6,1X))
              SLOC='  '
            END IF

            IF (STN(1:NSTLET) .NE. STANAM(J)(1:NSTLET) .OR.
     2      SNET(1:NETLET) .NE. JNET(J)(1:NETLET) .OR.
     3      (SLOC(1:NSLOC2) .NE. JSLOC(J)(1:NSLOC2) .AND.
     3      SLOC(1:NSLOC2) .NE. JSLOC2(J)(1:NSLOC2)) .OR.
     4      SCOMP(1:NCOMP) .NE. JCOMP3(J)(1:NCOMP)) GOTO 7

C--SEARCH THE LIST OF EXPIRATION DATES FOR THE FIRST ONE AFTER THE CURRENT DATE,
C  OR AN EXPIRATION DATE OF 0 (CAL GOOD THROUGH THE FUTURE).
            DO I=1,7
              IF (IYEARI(I).LT.100 .AND. IYEARI(I).GT.0)
     2        IYEARI(I)=IYEARI(I)+ICENT

C--TRUNCATE YEAR TO 2146, THE LARGEST STORABLE IN I*4, STORE WHOLE DATE
              IF (IYEARI(I).GT.2146) IYEARI(I)=2146
              IAEXP(I)=IYEARI(I)*1000000 +IAEXP(I)

              IF (IAEXP(I).EQ.0 .OR. IAEXP(I).GT.NDATE) THEN
                JCAL(J)=RCAL(I)*1000.

C--WRITE MESSAGE THAT A NEW CAL FACTOR WAS FOUND
                IF (LPRT .AND. KPRINT.GE.3) WRITE (15,1008) RCAL(I),
     2          STN,SNET,SCOMP,SLOC, JCEXP(J), IAEXP(I)
1008            FORMAT (' * NEW CAL FACTOR',F6.3,
     2          ' ASSIGNED TO ',A5,'-',A2,'-',A3,'-',A2,'.'/
     4          ' NEW CAL FACTOR START DATE = ',
     3          I10,', EXPIRATION DATE = ',I10)
                JCEXP(J)=IAEXP(I)
                GOTO 8
              END IF
            END DO
8           CLOSE (13)
          END IF

        END IF
      END IF
C++++++++++++++++++++++++++++++++++++++++++++++++++++++

C------- CALCULATE CODA MAGNITUDES -------------------------------------
C--MAKE SURE THERE IS A DURATION FOR THIS STATION
      IF (FMPK(K).LE.0.) GOTO 50
      TEMP=FMPK(K)

C--GET CALIBRATION FACTOR
C--USE KCAL IN PREFERENCE TO JCAL IF KCAL IS PRESENT
      IF (KCAL(K).EQ.0) THEN
        CAL=.001*JCAL(J)
      ELSE
        CAL=.01*KCAL(K)
      END IF

C--CONVERT CAL FACTOR TO GAIN CORRECTION (NO CAL FACTOR MEANS NO CORRECTION)
      G=0.
C--3.95 IS THE GAIN OF A "TYPICAL" 15 DB STATION WITH A CORRECTION OF 0.0
C--DO NOT APPLY GAIN CORRECTION TO STATIONS INSTRUCTED NOT TO
      IF (CAL.NE.0.) G=-ALOG10 (CAL/3.95) *JFGWT(J)

C--GG IS THE ACTUAL GAIN CORRECTION USED
C  IT WILL BE 0 IF COMPONENT IS NOT ON LIST AND COMPONENT SELECTION IS DESIRED
      IF (IDUG.LT.0) THEN               !ALL COMPONENTS
        GG=G
      ELSE IF (IDUG.EQ.0) THEN          !NO COMPONENTS
        GG=0.
      ELSE IF (IDUG.GT.0) THEN          !SOME COMPONENTS
        GG=0.
        DO I=1,IDUG
          IF (JCOMP3(J)(1:NCOMP) .EQ. CDUG(I)(1:NCOMP)) GG=G
        END DO
      END IF

C--CALCULATE 3 MAGNITUDE TYPES (SMAGF, SMAGT, SMAGB)
C--USE A BI-LINEAR LOG (F-P) SCALE WITH DIST & DEPTH CORRECTIONS
      TLOG=ALOG10(TEMP)
      IF (TEMP.LT.FMBRK) THEN
        SMAGF=FMA1 +FMB1*TLOG +FMD1*DIS(K) +FMZ1*Z1 +.01*JFCOR(J)
     2  +FMF1*TEMP
      ELSE
        SMAGF=FMA2 +FMB2*TLOG +FMD2*DIS(K) +FMZ2*Z1 +.01*JFCOR(J)
     2  +FMF2*TEMP
      END IF
      SMAGF=SMAGF +FMGN*GG

C--SECOND CODA MAGNITUDE
      IF (TEMP.LT.FMBRKB) THEN
        SMAGB=FMA1B +FMB1B*TLOG +FMD1B*DIS(K) +FMZ1B*Z1 +.01*JFCOR(J)
     2  +FMF1B*TEMP
      ELSE
        SMAGB=FMA2B +FMB2B*TLOG +FMD2B*DIS(K) +FMZ2B*Z1 +.01*JFCOR(J)
     2  +FMF2B*TEMP
      END IF
      SMAGB=SMAGB +FMGNB*GG

C--USE MAGNITUDE BASED ON ELAPSED TIME TAU.
C--USE CALC P TRAVEL TIME, NOT OBSERVED.
C--FIND STATION IN PHASE ARRAY
      DO IM=1,M
        IF ((K.EQ.IND(IM)) .AND. .NOT.MSFLAG(IM)) THEN
          TAU=MTCAL(IM)*.01 +TEMP
          GOTO 22
        END IF
      END DO
C--THIS SHOULD NEVER HAPPEN!
C--FORCE SMAGT=VNOMAG BELOW
      TAU=0.

C--GET MAGNITUDE FOR POSITIVE TAU
22    IF (TAU.GT..5) THEN
        TLOG=ALOG10 (TAU)
C--INCLUDE DEPTH, GAIN & STATION CORRECTION TERMS
        SMAGT=DMA0 +DMA1*TLOG +DMA2*TLOG**2 +DMLI*TAU
     2  +DMZ*Z1 +DMGN*G +.01*JFCOR(J)
      ELSE
C--WHEN TAU <= 0.5, SMAGT IS SET TO VNOMAG BELOW
C--SET SMAGT TO ZERO HERE TO PREVENT UNDEFINED SMAGT UNTIL THEN
        SMAGT=0.
      END IF

C--APPLY CORRECTIONS SPECIFIED FOR EACH COMPONENT
      DO I=1,NFCM
        IF (JCOMP3(J)(1:NCOMP) .EQ. CFCM(I)(1:NCOMP)) THEN
          SMAGT=SMAGT+AFCM(I)
          SMAGF=SMAGF+AFCM(I)
        END IF
      END DO

C--APPLY ADDITIONAL DISTANCE AND DEPTH CORRECTIONS
      IF (DIS(K).LT.DBRKM1) THEN
        SMAGT=SMAGT+DCOFM1*(DIS(K)-DBRKM1)
        SMAGF=SMAGF+DCOFM1*(DIS(K)-DBRKM1)
      END IF
      IF (DIS(K).GT.DBRKM2) THEN
        SMAGT=SMAGT+DCOFM2*(DIS(K)-DBRKM2)
        SMAGF=SMAGF+DCOFM2*(DIS(K)-DBRKM2)
      END IF
      IF (Z1.GT.ZBRKM) THEN
        SMAGT=SMAGT+ZCOFM*(Z1-ZBRKM)
        SMAGF=SMAGF+ZCOFM*(Z1-ZBRKM)
      END IF

C--DONT LET MAG BE LOWER THAN THE MINIMUM
      IF (SMAGF.LT.VNOMAG) SMAGF=VNOMAG
      IF (SMAGB.LT.VNOMAG) SMAGB=VNOMAG
      IF (SMAGT.LT.VNOMAG .OR. TAU.LE..5) SMAGT=VNOMAG

C--ASSIGN WEIGHT
      KFMWT=0
      IF (KFWT(K).LT.4) KFMWT=NINT(2.5*((4-KFWT(K))*JFWT(J)))

C--CHOOSE MAGNITUDE TYPE FOR FMAG1
C--NOTE THAT IF TWO FMAG TYPES ARE USED FOR THIS STATION, (I.E. ITS
C  COMPONET LETTER QUALIFIES IT FOR BOTH MAGNITUDES) EACH WILL FIGURE INTO
C  THE TWO EVENT MAGS BUT ONLY THE LAST WILL APPEAR IN THE STATION LIST.
      SMAG=VNOMAG
      KFMAG(K)=MINMAG
      IF (JFM1(J)) THEN
        IF (MAGSEL.EQ.1) THEN
          SMAG=SMAGF
        ELSE IF (MAGSEL.EQ.2) THEN
          SMAG=SMAGT
        ELSE IF (MAGSEL.EQ.3) THEN
          SMAG=SMAGB
        END IF
        KFMAG(K)=NINT(100.*SMAG)
        IF (SMAG.GT.VNOMAG .AND. KFMWT.GT.0 .AND. NMED.LT.MAXPHS) THEN
C         MFMAG=MFMAG+KFMWT
          NFMAG=NFMAG+1
          NMED=NMED+1
          IMAG(NMED)=KFMAG(K)
          IMWT(NMED)=KFMWT
        END IF
      END IF

C--CHOOSE MAGNITUDE TYPE FOR FMAG2
      IF (JFM2(J)) THEN
        IF (MAGSL2.EQ.1) THEN
          SMAG=SMAGF
        ELSE IF (MAGSL2.EQ.2) THEN
          SMAG=SMAGT
        ELSE IF (MAGSL2.EQ.3) THEN
          SMAG=SMAGB
        END IF
        KFMAG(K)=NINT(100.*SMAG)
        IF (SMAG.GT.VNOMAG .AND. KFMWT.GT.0 .AND. NMED2.LT.MAXPHS) THEN
C         MFMAG2=MFMAG2+KFMWT
          NFMAG2=NFMAG2+1
          NMED2=NMED2+1
          IMAG2(NMED2)=KFMAG(K)
          IMWT2(NMED2)=KFMWT
        END IF
      END IF

C--END OF STATION LOOP
50    CONTINUE

C-------- CALCULATE CODA MAGNITUDE FOR EVENT ---------------------------
C--CALCULATE WEIGHTED MEDIAN OF FMAG1 (OR MEAN)
      IF (NMED.GT.0) THEN
        CALL MEDWT (LMED,NMED,IMAG,IMWT,MEDIAN,MSD)
        FMAG=.01*MEDIAN
        IF (.NOT.LMED) FMMAD=.01*MSD
      END IF

C--CALC MEDIAN ABSOLUTE DIFFERENCE OF FMAG1
C  THE MAGS, WEIGHTS & MEDIAN ARE STILL IN THE ARRAYS
      IF (NMED.GT.1 .AND. LMED) THEN
        DO I=1,NMED
          IMAG(I)=ABS(IMAG(I)-MEDIAN)
        END DO
        CALL MEDWT (LMED,NMED,IMAG,IMWT,MEDIAN,MSD)
        FMMAD=.01*MEDIAN
      END IF
      NMED=0

C--CALCULATE WEIGHTED MEDIAN OF FMAG2
      IF (NMED2.GT.0) THEN
        CALL MEDWT (LMED,NMED2,IMAG2,IMWT2,MEDIAN,MSD)
        FMAG2=.01*MEDIAN
        IF (.NOT.LMED) FMMAD2=.01*MSD
      END IF

C--CALC MEDIAN ABSOLUTE DIFFERENCE OF FMAG2
C  THE MAGS, WEIGHTS & MEDIAN ARE STILL IN THE ARRAYS
      IF (NMED2.GT.1 .AND. LMED) THEN
        DO I=1,NMED2
          IMAG2(I)=ABS(IMAG2(I)-MEDIAN)
        END DO
        CALL MEDWT (LMED,NMED2,IMAG2,IMWT2,MEDIAN,MSD)
        FMMAD2=.01*MEDIAN
      END IF
      NMED2=0

C--CALCULATE AMPLITUDE MAGNITUDE --------------------------------------

C--COMPUTE XMAG AS A SUM OF 3 TERMS:
C--1: LOG(MAX PEAK-TO-PEAK AMP / 2*CAL FACTOR)
C--2: -LOG RESPONSE OF INST REL TO WOOD-ANDERSON.
C--3: -LOG(A(0)), THE LOCAL MAGNITUDE DISTANCE CORRECTION
C--ALFRQ = LOG(FREQUENCY OF SIGNAL AT MAX AMPLITUDE / 5HZ.)
C--SLDIS = HYPOCENTRAL DISTANCE IN KM
CC--ALDSQ = LOG((HYPOCENTRAL DISTANCE)**2)

C--LOOP OVER STATIONS
      DO 100 K=1,KSTA
C     KXMAG(K)=0
      KXMAG(K)=MINMAG
      J=KINDX(K)

C--DO NOT CALC A MAG IF THERE IS NO AMP OR IF AMP IS FROM THE WRONG INSTRUMENT
      KIMTYP(K)=0               !DEFAULT IS NO MAG FOR THIS STATION
      IF (AMPK(K).LE.0. .OR. KAMPTYP(K).GT.2) GOTO 100

C--DETERMINE MAGNITUDE TYPE
C--USE AMPLITUDE TYPE IF IT IS SPECIFIED FOR A WA OR VELOCITY
C--IF AMP TYPE IS UNSPECIFIED, USE INSTRUMENT TYPE
      IF (KAMPTYP(K).EQ.1) THEN
        KIMTYP(K)=1                     !WOOD ANDERSON ML
      ELSE IF (KAMPTYP(K).EQ.2) THEN
        KIMTYP(K)=2                     !VELOCITY MX
      ELSE IF (KAMPTYP(K).EQ.0) THEN
        IF (JTYPE(J).EQ.0 .OR. JTYPE(J).EQ.2) THEN
          KIMTYP(K)=1                   !WOOD ANDERSON ML
        ELSE
          KIMTYP(K)=2                   !VELOCITY MX
        END IF
      END IF

C--USE KCAL IN PREFERENCE TO JCAL IF KCAL IS PRESENT
      IF (KCAL(K).EQ.0) THEN
        CAL=.001*JCAL(J)
      ELSE
        CAL=.01*KCAL(K)
      END IF
      IF (CAL.EQ.0.) THEN
        KIMTYP(K)=0
        GO TO 100
      END IF

C--SET THE PERIOD AS STANDARD ONE FOR STATION IF NOT ON PHASE CARD
C  JPER(J) IS THE STATION PERIOD IN .1 SEC.
C  KPER(K) IS THE PHASE CARD PERIOD IN .01 SEC.
      IF (KPER(K).LE.0) THEN
        PER= .1*JPER(J)
      ELSE
        PER=.01*KPER(K)
      END IF
C--SKIP THE CALCULATIONS IF PERIOD IS OUT OF ALLOWED RANGE
      IF (PER.GT.6.3 .OR. PER.LT..02) THEN
        KIMTYP(K)=0
        GO TO 100
      END IF

C--USE AMPLITUDE IN PP MM FOR MAGNITUDE
C--THE COUNT-TO-MM CONVERSION IS USED FOR ALL DIGITAL SYSTEMS, CAL
C  IS USED TO ALLOW FOR DIFFERENCES FROM THE 12 BIT STANDARD SYSTEM.
C--USGS STANDARD SYSTEM IS 4096 COUNTS = 5.0 VOLT = 200 MM
      AMPMM=AMPK(K)                     !KAMPU(K)=0
      IF (KAMPU(K).EQ.1) THEN           !HALF AMPLITUDE
        AMPMM=AMPK(K)*2.
      ELSE IF (KAMPU(K).EQ.2) THEN      !DIGITAL COUNTS FROM EARTHWORM OR MENLO
        AMPMM=AMPK(K)*0.04883
      ELSE IF (KAMPU(K).EQ.3) THEN      !DIGITAL COUNTS FROM HVO
        AMPMM=AMPK(K)*0.012207
      ELSE IF (KAMPU(K).EQ.4) THEN      !DIGITAL COUNTS IN HUNDREDS FROM MENLO
        AMPMM=AMPK(K)*4.883
      END IF

C--COMPUTE MAGNITUDE
      SLDIS=SQRT(Z1**2 +DIS(K)**2) +.01
C--AMPK IS PEAK-TO-PEAK AMPLITUDE IN MM (.5 IS TO USE HALF-AMP)
      SMAG=ALOG10(AMPMM*.5/CAL)

C--CORRECT TO WA RESPONSE IF A VELOCITY SEISMOMETER
C--ALSO CORRECT DIGITAL TELEMETRY STATIONS WITH SAME CURVE
C  (THIS IS THE OLD APPROXIMATION FOR DOMINANT FREQ)
C      IF (JTYPE(J).EQ.1 .OR. JTYPE(J).EQ.3) SMAG=SMAG-1.3-.95*ALFRQ
C      IF (JTYPE(J).EQ.1 .OR. (JTYPE(J).GT.2 .AND.
C     2 JTYPE(J).LE.7)) THEN

C--APPLY CORRECTION FOR VELOCITY SEISMOMETERS
C--DO NOT APPLY CORRECTION FOR WA (TYPE 0) OR SPRENGNETHER (TYPE 2)
      IF (KIMTYP(K).EQ.2) THEN
        FQ=10.*ALOG10(1./PER)+9.
        IFQ=FQ
C--USE RESPONSE CORRECTION INTERPOLATED FROM TABLE
        SMAG=SMAG -(RSPA(IFQ) +(FQ-IFQ) *(RSPA(IFQ+1)-RSPA(IFQ)))
C--APPLY MOTOR CONSTANT OF SEISMOMETER
        SMAG=SMAG -ALOG10 (S (JTYPE(J)))
      END IF

C--CORRECT AN HVO TYPE SPRENGNETHER
      IF (JTYPE(J).EQ.2) SMAG=SMAG+.41+.56*ALOG10(.2/PER)

C--APPLY CORRECTIONS SPECIFIED FOR EACH COMPONENT.
      DO I=1,NXCM
        IF (JCOMP3(J)(1:NCOMP) .EQ. CXCM(I)(1:NCOMP)) SMAG=SMAG+AXCM(I)
      END DO

C--USE THE DEFAULT LOGA0 RELATION FOR MOST COMPONENTS
      LATYP=MLOGA0
C--TEST COMPONENTS TO USE A COMPONENT SPECIFIC LOGA0 RELATION
      DO I=1,NLA0
        IF (JCOMP3(J)(1:NCOMP) .EQ. CLA0(I)(1:NCOMP)) LATYP=MLA0(I)
      END DO

C--APPLY RICHTER'S LOGA0 DISTANCE CORRECTION TERM
      CALL LOGA0 (LATYP,A0MAG,DIS(K),SLDIS)
      SMAG=SMAG+A0MAG
C--OUTPUT LOG(A0) TERM IN REMARK COLUMN OF PRINTOUT (DEBUG ONLY)
C      WRITE (KRMK6(K),'(F4.2,''-'',I1)') A0MAG,LATYP

C--APPLY MAG CORRECTION & ACCUMULATE SUMS
      SMAG=SMAG+.01*JXCOR(J)
C      IF (SMAG.LT.0.) SMAG=0.
      KXMAG(K)=NINT(100.*SMAG)
      KXMWT=0
      IF (KXWT(K).LT.4) KXMWT=NINT (2.5*((4-KXWT(K))*JXWT(J)))

C--USE XMAG IN FIRST AMP MAG
C--MXMAG IS 100X THE TOTAL OF STATION WEIGHTS (NO LONGER USED)
C--USE THIS MAG IN XMAG1 IF
C  1) THE COMPONENT IS ON THE COMPONENT LIST (OR ALL COMPS ARE USED), AND
C  2) THE TYPE (MX OR ML) MATCHES THAT SELECTED FOR XMAG1 (OR ALL MAGS USED)
      IF(JXM1(J) .AND. (MAG1TYPX.EQ.0 .OR. KIMTYP(K).EQ.MAG1TYPX)) THEN
C       IF (SMAG.GT.0. .AND. KXMWT.GT.0 .AND. NMED.LT.MAXPHS) THEN
        IF (KXMWT.GT.0 .AND. NMED.LT.MAXPHS) THEN
C         MXMAG=MXMAG+KXMWT
          NXMAG=NXMAG+1
          NMED=NMED+1
          IMAG(NMED)=KXMAG(K)
          IMWT(NMED)=KXMWT
        END IF
      END IF

C--USE XMAG IN SECOND AMP MAG
      IF(JXM2(J) .AND. (MAG2TYPX.EQ.0 .OR. KIMTYP(K).EQ.MAG2TYPX)) THEN
C       IF (SMAG.GT.0. .AND. KXMWT.GT.0 .AND. NMED2.LT.MAXPHS) THEN
        IF (KXMWT.GT.0 .AND. NMED2.LT.MAXPHS) THEN
C         MXMAG2=MXMAG2+KXMWT
          NXMAG2=NXMAG2+1
          NMED2=NMED2+1
          IMAG2(NMED2)=KXMAG(K)
          IMWT2(NMED2)=KXMWT
        END IF
      END IF

100   CONTINUE

C----------- CALCULATE AMPLITUDE MAGNITUDE FOR EVENT -------------------
C--CALCULATE WEIGHTED MEDIAN OF FIRST AMP MAG
      IF (NMED.GT.0) THEN
        CALL MEDWT (LMED,NMED,IMAG,IMWT,MEDIAN,MSD)
        XMAG=.01*MEDIAN
        IF (.NOT.LMED) XMMAD=.01*MSD
      END IF

C--CALC MEDIAN ABSOLUTE DIFFERENCE OF FIRST AMP MAG
C  THE MAGS, WEIGHTS & MEDIAN ARE STILL IN THE ARRAYS
      IF (NMED.GT.1 .AND. LMED) THEN
        DO I=1,NMED
          IMAG(I)=ABS(IMAG(I)-MEDIAN)
        END DO
        CALL MEDWT (LMED,NMED,IMAG,IMWT,MEDIAN,MSD)
        XMMAD=.01*MEDIAN
      END IF

C--CALCULATE WEIGHTED MEDIAN OF SECOND AMP MAG
      IF (NMED2.GT.0) THEN
        CALL MEDWT (LMED,NMED2,IMAG2,IMWT2,MEDIAN,MSD)
        XMAG2=.01*MEDIAN
        IF (.NOT.LMED) XMMAD2=.01*MSD
      END IF

C--CALC MEDIAN ABSOLUTE DIFFERENCE OF SECOND AMP MAG
C  THE MAGS, WEIGHTS & MEDIAN ARE STILL IN THE ARRAYS
      IF (NMED2.GT.1 .AND. LMED) THEN
        DO I=1,NMED2
          IMAG2(I)=ABS(IMAG2(I)-MEDIAN)
        END DO
        CALL MEDWT (LMED,NMED2,IMAG2,IMWT2,MEDIAN,MSD)
        XMMAD2=.01*MEDIAN
      END IF

      RETURN
      END

C-----------------------------------------------------------
C--APPLY THE LOGA0 TERM TO THE AMPLITUDE MAGNITUDE
      SUBROUTINE LOGA0 (LATYP,A0MAG,DIS,SLDIS)

      IMPLICIT NONE

      REAL     FLOGA0
      EXTERNAL FLOGA0

      INTEGER  LATYP
      REAL     A0MAG, DIS, SLDIS

      REAL*8   DPRDIST, CISN_mlAo
      REAL     ALDSQ
      INTEGER  INDX

C-- -LOGA0 VALUES & THEIR MAX DISTANCES DIRECTLY FROM RICHTERS TABLE
      REAL A0RICT(71)
      SAVE A0RICT
      DATA A0RICT /2*1.4,1.5, 1.6, 1.7, 1.9, 2.1, 2.3, 2.4, 2.5, 2.6,
     2 2.7,3*2.8,2.85,2*2.9,3*3.0,2*3.1,2*3.2, 2*3.3,2*3.4,2*3.5,3.6,
     3 3.65, 2*3.7, 2*3.8, 2*3.9, 2*4.0, 2*4.1, 2*4.2, 3*4.3, 2*4.4,
     4 3*4.5, 4*4.6, 4*4.7, 5*4.8, 5*4.9/

      A0MAG=0.
      IF (LATYP.LE.0 .OR. LATYP.GT.6) RETURN
C--USE THE APPROXIMATION OF EATON
      IF (LATYP.EQ.1) THEN
        IF (SLDIS.LT.185.3) THEN
          A0MAG=.821*ALOG10(SLDIS) +.00405*SLDIS +.955
        ELSE
          A0MAG=2.55*ALOG10(SLDIS) -2.21
        END IF

C--APPLY A SMALL SINUSOIDAL CORRECTION AT DISTANCES CLOSER THAN 70 KM
        IF (DIS.LT.70.) A0MAG=A0MAG -.09*SIN(.07*(DIS-25.))

C--ADD A (SIGNIFICANT) CORRECTION FOR P AMPLITUDE MAGNITUDES
C--CORRECTION CONSISTS OF 4 LINEAR SEGMENTS WITH DISTACE
C          IF (LATYP.EQ.5) THEN
C            IF (DIS.LT.52.) THEN
C              XX= -.08 +.00942*DIS
C            ELSE IF (DIS.GE.52. .AND. DIS.LT.115.) THEN
C              XX= +.41
C            ELSE IF (DIS.GE.115. .AND. DIS.LT.280.) THEN
C              XX= +.812 -.0035*DIS
C            ELSE IF (DIS.GE.280.) THEN
C              XX= -.168
C            END IF
C          A0MAG=A0MAG +XX
C          END IF

      ELSE IF (LATYP.EQ.2) THEN
C--USE THE APPROXIMATION OF BAKUN & JOYNER
        A0MAG= ALOG10(SLDIS) +.003*SLDIS +.7

      ELSE IF (LATYP.EQ.3) THEN
C--USE RICHTERS ORIGINAL APPROXIMATION
        ALDSQ=ALOG10(SLDIS**2)
        A0MAG=-.15+.8*ALDSQ
C       IF (ALDSQ.GT.4.6) A0MAG=-3.23+.7*ALDSQ
        IF (ALDSQ.GT.4.6) A0MAG= A0MAG -3.23+.7*ALDSQ

      ELSE IF (LATYP.EQ.4) THEN
C--APPLY BERKELEY'S CODE FOR THE NORDQUIST NOMOGRAM
        A0MAG=FLOGA0(DIS)

C--USE RICHTERS ORIGINAL TABLE
      ELSE IF (LATYP.EQ.5) THEN
        A0MAG=5.0
C--INDX IS THE INDEX OF THE NEAREST ENTRY IN RICHTERS TABLE
        IF (DIS.LT.100.) THEN
          INDX=NINT(DIS/5.)+1
        ELSE
          INDX=NINT(DIS/10.)+11
        END IF
        IF (INDX.GT.71) INDX=71
        A0MAG=A0RICT(INDX)

C--USE LOGA0 FUNCTION FROM CISN NETWORKS (USED IN TRIMAG & JIGGLE)
      ELSE IF (LATYP.EQ.6) THEN
        IF (SLDIS .LE. 0.1) RETURN
        DPRDIST=SLDIS
        A0MAG=CISN_mlAo (DPRDIST)
        RETURN
      END IF
      RETURN
      END

C---------------------------------------------------------------
      REAL FUNCTION FLOGA0(DELTA)

C************************************************************************
C*   DEFINE LOG(AO) FUNCTION
C*   BASED ON LINEAR INTERPOLATION OF NORDQUIST NOMOGRAPH (BSSA, 1948)
C************************************************************************
C--THIS CODE FROM BERKELEY/BOB UHRHAMMER 5/95

      IMPLICIT NONE

      REAL    DELTA

      INTEGER I, J, J0

      REAL D(19),A(19)
      DATA A /1.4,1.4,1.5,1.7,2.3,2.7,2.85,
     1        3.0,3.65,4.2,4.6,5.0,5.3,5.7,
     2        6.0,6.2,6.3,6.4,6.5/
      DATA D /0.,5.,10.,20.,35.,55.,73.,96.,
     1        220.,334.,445.,600.,750.,1010.,
     2        1340.,1600.,1750.,1960.,2230./
      SAVE A, D

      FLOGA0=0.
      IF ((DELTA.LT.0.).OR.(DELTA.GT.2224.)) RETURN
      J=0

      DO I=1,19
        J=J+1
        IF (DELTA.LT.D(I)) GOTO 2
      END DO
      RETURN
2     J0=J-1
      FLOGA0=A(J0)+(A(J)-A(J0))*(DELTA-D(J0))/(D(J)-D(J0))

      RETURN
      END

C------------------------------------------------------------------
      real*8 function CISN_mlAo( rdist )
c--Originally from Bob Uhrhammer from mid 2007 via Al Walter
c--Incorporated in hypoinverse in 2/2011
c--rdist is the slant (hypocentral) distance
c
c ...... calculate CISN -logAo ML attenuation function
c
      implicit none

      real*8   logAo, T
      external logAo, T

      real*8 rdist

      real*8 TP(6), mlogAo, z, x, b
      integer j
c
      data TP( 1 ) / +0.056d0 /
      data TP( 2 ) / -0.031d0 /
      data TP( 3 ) / -0.053d0 /
      data TP( 4 ) / -0.080d0 /
      data TP( 5 ) / -0.028d0 /
      data TP( 6 ) / +0.015d0 /
      save TP
c
      if( rdist .le. 0.1d0 ) then
c
c ...... invalid for rdist less that 0.1 km
c        return with -9.d0
c
        mlogAo = -9.d0
c
      elseif( rdist .le. 8.d0 ) then
c
c ...... linear extrapolation of average slope between 8 km and 60 km
c
        b = ( 2.6182d0-1.5429d0)/(dlog10(60.d0)-dlog10(8.d0))
        mlogAo = 1.5429d0 + b * ( dlog10( rdist ) - dlog10( 8.d0 ) )
c
      elseif( rdist .le. 500.d0 ) then
c
c ...... Chebychev polynomial expansion
c
        x = z( rdist )
        mlogAo = logAo( rdist ) + 0.0054d0
        do j = 1 , 6
          mlogAo = mlogAo + TP( j ) * T( j , x )
        end do
c
      else
c
c ...... invalid for rdist greater than 500 km
c        return with -9.d0
c
        mlogAo = -9.d0
c
      endif
c
      CISN_mlAo = mlogAo
c
      return
      end

      real*8 function T( n , x )
c
c ...... Chebyshev Polynomial
c
      implicit none

      integer n
      real*8 x
      
      real*8 theta
c
      theta = dacos( x )
      T = dcos( dble( n ) * theta )
c
      return
      end

      real*8 function z( r )
c
c ...... translate scale from r to z
c
      implicit none

      real*8 r

      real*8 r0, r1, z0, z1, a , b, l_r0, l_r1
      integer ncall
c
      data ncall /0/
      data r0,r1 /8.d0,500.d0/
      data z0,z1 /-1.d0,+1.d0/
      save ncall, r0, r1, z0, z1
c
      if( ncall .eq. 0 ) then
        l_r0 = dlog10( r0 )
        l_r1 = dlog10( r1 )
        b = ( z1 - z0 ) / ( l_r1 - l_r0 )
        a = z0 - b * l_r0
        ncall = 1
      endif
c
      z = a + b * dlog10( r )
c
      return
      end

      real*8 function logAo( rdist )
c
c ...... -logAo attenuation function
c
      implicit none

      real*8 rdist
      
c
      logAo = 1.11d0 * dlog10( rdist ) +
     1    0.00189d0 * rdist + 0.591d0
c
      return
      end
