      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/

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--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 AND HYTRH USE ONLY MOD
C  TO DECIDE WHAT TO DO, AND KNOW NOTHING OF MULTIPLE MODELS.
      IF (MODTYP(MOD).EQ.0) THEN
        CALL HYTRT
      ELSE
        CALL HYTRH
      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.
      IF (MODTYP(MOD).EQ.0) THEN
        CALL HYTRT
      ELSE
        CALL HYTRH
      END IF
      RETURN
      END
