      SUBROUTINE HYTRA
C--MANAGE MULTIPLE CRUSTAL MODELS & CALC ALL TRAVEL TIMES & DERIVATIVES
      INCLUDE 'common.inc'
C--AA ACCUMULATES INDIVIDUAL TRAVEL TIMES & DERIVS FOR THE MODELS USED
      DIMENSION AA (MMAX,4)
C--NODE HOLDS THE CURRENT NODE NUMBER & IS TESTED FIRST FOR PROGRAM SPEED
      SAVE NODE
      DATA NODE /0/
      LOGICAL LMESS
      LMESS=KPRINT.GT.1 .AND. LPRT

C---------------- SEE IF ONLY 1 MODEL IS NEEDED -----------------------
C--USE MODEL 1 EXCLUSIVELY IF MULTIPLE MODELS ARE NOT IN USE
      IF (.NOT.LMULT) THEN
        MOD=1
        GOTO 90
      END IF
      COSLAT=COS(CLAT/RDEG)

C--IF THIS IS NOT THE FIRST ITERATION, TRY THE SAME NODE AS LAST TIME.
C  MOD SHOULD BE SET TO THE SAME MODEL AS LAST TIME. LEAVE NODE UNCHANGED.
      IF (ITR.GT.1 .AND. NODE.GT.0) THEN
        DIST=111.195* SQRT((CLAT-HLAT(NODE))**2
     2  +(COSLAT*(CLON-HLON(NODE)))**2)
        IF (DIST.LE.RAD1(NODE)) GOTO 90
      END IF

C--FOR EACH NODE IN TURN, CALC EPICENTRAL DIST & TEST IF EPI IS IN INNER CIRCLE
      DO NODE=1,NNODE
        HDIS(NODE)=111.195* SQRT((CLAT-HLAT(NODE))**2
     2  +(COSLAT* (CLON-HLON(NODE)))**2)
        IF (HDIS(NODE) .LE. RAD1(NODE)) THEN
          MOD=MODH(NODE)
          GOTO 90
        END IF
      END DO

C----------------------- TEST FOR A MIXTURE OF MODELS --------------------
C--WE NOW HAVE A MIXTURE OF MODELS OR THE DEFAULT MODEL. RECORD THIS FOR
C  THE NEXT ENTRY INTO HYTRA
      NODE=0
C--ZERO THE NODE WEIGHTS & SET THE NODE COUNTER
      MNODE=0
      DO I=1,3
        WMOD(I)=0.
      END DO

C--TEST ALL NODES TO SEE IF EPICENTER IS IN OUTER CIRCLE.
C  IF SO, RECORD IT BUT STOP AT 3 NODES
      WNORM=0.
      DO NO=1,NNODE
        IF (HDIS(NO) .LT. RAD2(NO)) THEN
          MNODE=MNODE+1
C--RECORD THE CRUST MODEL NUMBER FOR THIS NODE
          MODS(MNODE)=MODH(NO)
C--SET THE WEIGHT OF THIS MODEL FOR AVERAGING USING A COSINE TAPER
          DIST=(HDIS(NO)-RAD1(NO))/DRAD(NO)
          WMOD(MNODE)=.5+.5*COS(PI*DIST)
          WNORM=WNORM+WMOD(MNODE)

C--SKIP REMAINING NODES IF THIS MAKES 3
          IF (MNODE.EQ.3) GOTO 10
        END IF
      END DO

C--WE NOW HAVE MNODE MODELS. USE A WEIGHTED AVERAGE. THE LOGIC DEPENDS ON MNODE:
C  MNODE=0  USE 100% OF THE DEFAULT MODEL
C  MNODE=1  USE "WMOD" WEIGHT OF MODEL "MODS", 1-WMOD OF THE DEFAULT
C  MNODE=2  & WNORM>=1  USE ONLY THE 2 MODELS BUT RENORMALIZE
C    "      & WNORM <1  USE THE 2 MODELS PLUS THE DEFAULT MODEL
C  MNODE=3  USE THE 3 MODELS BUT RENORMALIZE

C--IF WE WERE OUTSIDE ALL OUTER CIRCLES, USE ONLY THE DEFAULT MODEL
10    IF (MNODE.EQ.0) THEN
        MOD=MODDEF
        GOTO 90

C--IF IN ONLY 1 OUTER CIRCLE, MIX IN SOME OF THE DEFAULT
      ELSE IF (MNODE.EQ.1) THEN
        NMOD=2
        MODS(2)=MODDEF
        WMOD(2)=1.-WMOD(1)

C--IF IN 2 CIRCLES, MIX IN SOME DEFAULT IF WE ARE NOT FULLY WEIGHTED
      ELSE IF (MNODE.EQ.2) THEN
        IF (WNORM.LT.1.) THEN
C--MIX IN SOME DEFAULT AS THE 3RD MODEL
          NMOD=3
          MODS(3)=MODDEF
          WMOD(3)=1.-WNORM
        ELSE
C--USE ONLY THE 2 MODELS, BUT RENORMALIZE THE WEIGHTS
          NMOD=2
          WMOD(1)=WMOD(1)/WNORM
          WMOD(2)=WMOD(2)/WNORM
        END IF

C--IF IN 3 CIRCLES, KEEP THEM ALL BUT RENORMALIZE WEIGHTS
      ELSE IF (MNODE.EQ.3) THEN
        NMOD=3
        DO I=1,3
          WMOD(I)=WMOD(I)/WNORM
        END DO

C------------------- CHECK THE MODELS FOR REDUNDANCY --------------------
C--CHECK THE 2ND & 3RD MODELS FOR REDUNDANCY NOW WHILE WE KNOW THERE ARE 3
        IF (MODS(2).EQ.MODS(3)) THEN
          NMOD=2
          WMOD(2)=WMOD(2)+WMOD(3)
        END IF
      END IF

C--CHECK THE 1ST & 2ND MODELS FOR REDUNDANCY. THERE MAY BE 2 OR 3 DEFINED
      IF (MODS(1).EQ.MODS(2)) THEN
        NMOD=NMOD-1
        WMOD(1)=WMOD(1)+WMOD(2)
        MODS(2)=MODS(3)
        WMOD(2)=WMOD(3)
      END IF

C--IF WE ARE LEFT WITH ONE MODEL, GO USE SINGLE MODEL LOGIC
      IF (NMOD.EQ.1) THEN
        MOD=MODS(1)
        GOTO 90
      END IF

C------------------- SORT THE MODELS BY WEIGHT -------------------------
C--MAKE SURE THE 3RD MODEL HAS THE SMALLEST WEIGHT
      IF (NMOD.EQ.3) THEN
        I=2
        IF (WMOD(1).LT.WMOD(2)) I=1      !I HAS THE LEAST WEIGHT
C--COMPARE WEIGHTS I & 3
        IF (WMOD(3).GT.WMOD(I)) THEN
          TEMP=WMOD(3)
          WMOD(3)=WMOD(I)
          WMOD(I)=TEMP
          IT=MODS(3)
          MODS(3)=MODS(I)
          MODS(I)=IT
        END IF
      END IF

C--COMPARE WEIGHTS 1 & 2 & INTERCHANGE IF NECESSARY
      IF (WMOD(2).GT.WMOD(1)) THEN
        TEMP=WMOD(2)
        WMOD(2)=WMOD(1)
        WMOD(1)=TEMP
        IT=MODS(2)
        MODS(2)=MODS(1)
        MODS(1)=IT
      END IF

C---------------- GET WEIGHTED AVERAGE OF 2 OR 3 MODELS -------------------
C--USE THE ARRAY AA TO ACCUMULATE THE TRAVEL TIME DATA
      DO 40 I=1,M
      DO 40 IT=1,4
40    AA(I,IT)=0.

C--SET THE LIMIT ON NEGATIVE EQ DEPTH TO THE MIN OF MIN ELEVS OF ALL 3 MODELS
C  FWK 1.37
      ELEVMAX=ELEVMX(MODS(1))
      DO NO=1,NMOD
        IF (ELEVMX(MODS(NO)) .LT. ELEVMAX) ELEVMAX=ELEVMX(MODS(NO))
      END DO

C--LIMIT DEPTH TO BE BELOW THE MAXIMUM ELEVATION TO AVOID PROGRAM CRASH
C  FWK 1.37
C      IF (Z1.LT.ELEVMAX) THEN
C        Z1=ELEVMAX
      IF (Z1.LT.-ELEVMAX) THEN	!FWK 1/2/2018 JP bug correction
        Z1=-ELEVMAX
        IF (LMESS) WRITE (15,*)
     2  ' LIMIT QUAKE BELOW MINIMUM REFERENCE ELEVATION OF MODELS USED'
      END IF

C--LOOP OVER THE 2 OR 3 MODELS
      DO 50 NO=1,NMOD
        WNORM=WMOD(NO)
        MOD=MODS(NO)

C--CALC TRAVEL TIMES & DERIVS FOR ALL STATIONS. HYTRT, HYTRH, HYLINV USE ONLY
C  MOD TO DECIDE WHAT MODEL TO USE, AND KNOW NOTHING OF MULTIPLE MODELS.
        IF (MODTYP(MOD).EQ.0) THEN
          IF (Z1.LT.0.) THEN
            WRITE (6,1000) Z1,MOD,KYEAR2,KMONTH,KDAY,KHOUR,KMIN,T1
1000        FORMAT('*** ILLEGAL NEGATIVE DEPTH',F8.2,
     2      ' USED FOR CRT MODEL # ',I3,/
     3      ' FOR EVENT',I5,4I3,' SEC=',F7.2,'. STOP PROGRAM.')
            IF (LPRT) WRITE (15,1000) 
     2      Z1,MOD,KYEAR2,KMONTH,KDAY,KHOUR,KMIN,T1
            STOP
          END IF
          CALL HYTRT (.FALSE.)	!CALL ONCE FOR P MODEL, AGAIN FOR S IF NEEDED
          IF (MODSAL(MOD).GT.0) CALL HYTRT (.TRUE.)

        ELSE IF (MODTYP(MOD).EQ.1) THEN
          IF (Z1.LT.0.) THEN
            WRITE (6,1001) Z1,MOD,KYEAR2,KMONTH,KDAY,KHOUR,KMIN,T1
1001        FORMAT('*** ILLEGAL NEGATIVE DEPTH',F8.2,
     2      ' USED FOR CRH MODEL # ',I3,/
     3      ' FOR EVENT',I5,4I3,' SEC=',F7.2,'. STOP PROGRAM.')
            IF (LPRT) WRITE (15,1001) 
     2      Z1,MOD,KYEAR2,KMONTH,KDAY,KHOUR,KMIN,T1
            STOP
          END IF
          CALL HYTRH

        ELSE IF (MODTYP(MOD).EQ.2) THEN
          CALL HYLINV
        ELSE IF (MODTYP(MOD).EQ.3) THEN
          CALL HYLIN
        ELSE IF (MODTYP(MOD).EQ.4) THEN
          CALL HYTRE
        END IF

C--WEIGHT THE RESULTS & ACCUMULATE IN AA
        DO 45 I=1,M
        DO 45 IT=1,4
45      AA(I,IT)=AA(I,IT)+A(I,IT)*WNORM
50    CONTINUE

C--TRANSFER THE TRAVEL TIME DATA BACK TO A, WHICH HYLOC USES
      DO 55 I=1,M
      DO 55 IT=1,4
55    A(I,IT)=AA(I,IT)
C--BE SURE TO LEAVE MOD SET TO THE DOMINANT MODEL
      MOD=MODS(1)
      RETURN

C------------------- USE ONLY ONE MODEL -----------------------
90    NMOD=1
      MODS(1)=MOD
      WMOD(1)=1.
      ELEVMAX=ELEVMX(MOD)	!FWK 1.37

C--LIMIT DEPTH TO BE BELOW THE MAXIMUM ELEVATION TO AVOID PROGRAM CRASH
C  LMB 1.42
      IF (Z1.LT.-ELEVMAX) THEN
        Z1=-ELEVMAX
        IF (LMESS) WRITE (15,*)
     2  ' LIMIT QUAKE BELOW MINIMUM REFERENCE ELEVATION OF MODELS USED'
      END IF

      IF (MODTYP(MOD).EQ.0) THEN
          CALL HYTRT (.FALSE.)	!CALL ONCE FOR P MODEL, AGAIN FOR S IF NEEDED
          IF (MODSAL(MOD).GT.0) CALL HYTRT (.TRUE.)
      ELSE IF (MODTYP(MOD).EQ.1) THEN
        CALL HYTRH
      ELSE IF (MODTYP(MOD).EQ.2) THEN
        CALL HYLINV
      ELSE IF (MODTYP(MOD).EQ.3) THEN
        CALL HYLIN
      ELSE IF (MODTYP(MOD).EQ.4) THEN
        CALL HYTRE
      END IF
      RETURN
      END
