      SUBROUTINE HYLOC
C--LOCATES ONE EARTHQUAKE FOR THE PROGRAM HYPOINVERSE
      LOGICAL CONERR,LZFIX,LPKTMP,LMESS,LFREEW,LSHAD
      CHARACTER C1*1,C3*3
      INCLUDE 'common.inc'

C--THIS IS THE HYPOCENTER ADJUSTMENT VECTOR
      DIMENSION Y(4)

C--NRES IS THE NUMBER OF TIMES THAT RESIDUAL WEIGHTING IS APPLIED EACH ITER.
      SAVE NRES
      DATA NRES /2/

C--S ARRIVAL FLAG
      LOGICAL LPS

C--BLANK OUT FIELDS COMPUTED AFTER EQ IS LOCATED IN CASE RESULTS ARE PRINTED 
C  EACH ITERATION
      SOUCOD=' '
      FMSOU=' '
      XMSOU=' '
      LABPR=' '
      PMAG=0.

C--INITIALIZE VARIABLES
      LSHAD=JCP.EQ.4 .OR. JCP.EQ.5
      N=3
C--WHEN HYPO IS FIXED, SOLVE FOR EVERYTHING TO GET ERZ, ETC.
      IF (ITRLIM.EQ.0) N=4
      ITR=0
      LFREEW=.TRUE.
      DONE=.FALSE.
      RMSMIN=10000.
      OLDRMS=10000.
      RR=10000.
      LMESS=KPRINT.GT.1 .AND. LPRT

C--DEPFIX INDICATES WHETHER DEPTH IS FIXED FOR THIS EVENT, & LZFIX INDICATES
C  WHETHER IT IS FIXED ON THIS ITERATION.
      IF (Z1.EQ.0.) Z1=ZTR
      LZFIX=.TRUE.

C--ANNOUNCE WHETHER WE ARE FIXING THE HYPOCENTER
      IF (ALLFIX .AND. LPRT) THEN
        WRITE (15,*) ' FIX HYPOCENTER AND ORIGIN TIME FOR THIS EVENT'
      ELSE IF (HYPOFIX .AND. LPRT) THEN
        WRITE (15,*) ' FIX HYPOCENTER FOR THIS EVENT'
      ELSE IF (DEPFIX .AND. LPRT) THEN
        WRITE (15,*) ' FIX DEPTH FOR THIS EVENT'
      END IF

C--ZERO THE DATA NOT CALCULATED EACH ITERATION IF WE ARE PRINTING
C  A STATION LIST EACH ITERATION.
      IF (KPRINT.GT.5) THEN
        DO I=1,MMAX
          IMPORT(I)=0
        END DO
C        MXMAG=0
        NXMAG=0
        XMAG=0.
        XMMAD=0.
C        MFMAG=0
        NFMAG=0
        FMAG=0.
        FMMAD=0.
        REMK=' '
        DO K=1,KSTA
          KFMAG(K)=0
          KXMAG(K)=0
        END DO
      END IF

C******************* BEGIN ITERATION LOOP ********************************

30    M=0
      NWR=0
      D1=10000.
      D2=10000.
      ITR=ITR+1
      COSLAT=COS(CLAT/RDEG)

C--PRINT HEADING ON FIRST ITERATION
      IF (ITR.EQ.1 .AND. LMESS) WRITE (15,1010)

C--SIGNAL THE BEGINNING OF RESIDUAL & DISTANCE WEIGHTING
      IF ((ITR.EQ.ITRDI1) .AND. LMESS) WRITE (15,1081)
1081  FORMAT (' BEGIN FIRST DISTANCE WEIGHTING')
      IF ((HYPOFIX .OR. ITR.EQ.ITRDIS) .AND. LMESS) WRITE (15,1001)
1001  FORMAT (' BEGIN MAIN DISTANCE WEIGHTING')
      IF ((HYPOFIX .OR. ITR.EQ.ITRRES) .AND. LMESS) WRITE (15,1002)
1002  FORMAT (' BEGIN RESIDUAL WEIGHTING')

C++++++++++++++++++++++ CALC STATION DISTANCES AND AZIMUTHS +++++++++++++

35    DO 40 K=1,KSTA
      J=KINDX(K)

C--NEW DISTANCE CALCULATION USING UTM GRID CALCULATION
      DLATU=JLATD(J)+JLATM(J)/6000.
      DLONU=JLOND(J)+JLONM(J)/6000.
      CALL UTMCAL(CLAT,CLON,DLATU,DLONU,XSUTM,YSUTM,DELTKM)
      DIS(K)=DELTKM
C--CALCULATE AZIMUTH TO STATION IN DEG BETWEEN +/- 180
      KTEMP=0
      IF (XSUTM*YSUTM .NE. 0.) KTEMP=RDEG*ATAN2(-XSUTM,YSUTM)
      KAZEM(K)=180*KTEMP

C--SKIP UNWEIGHTED STATIONS IN FINDING THE 2ND CLOSEST STATION
C--DECODE P & S WEIGHTS
      LSWT=KWT(K)/10
      LPWT=KWT(K)-10*LSWT
      IF ((KPRK(K).EQ.'  '.OR.LPWT.GT.3) .AND.
     2 (KSRK(K).EQ.'  '.OR.LSWT.GT.3)) GOTO 40

C--DETERMINE THE INDEX OF THE 2ND CLOSEST STATION
      IF (DELTKM.EQ.D1 .OR. DELTKM.EQ.D2) THEN
        GOTO 40
      ELSE IF (DELTKM.LT.D1) THEN
        IF (D2 .NE. D1) THEN
          D2=D1
          K2=K1
        END IF
        D1=DELTKM
        K1=K
      ELSE IF (DELTKM.LT.D2) THEN
        D2=DELTKM
        K2=K
      END IF
40    CONTINUE

C--TEMP3 IS THE DISTANCE WEIGHTING SCALE FACTOR
      TEMP3=DIS(K2)
      IF (TEMP3.LT.DISCUT) TEMP3=DISCUT
      TEMP=TEMP3*DISW1
      TEMP2=TEMP3*DISW2

      TEMP4=DIS(K2)
      IF (TEMP4.LT.DISCU1) TEMP4=DISCU1
      TEMP5=TEMP4*DISW11
      TEMP6=TEMP4*DISW21

C--ABORT EVENT IF DISTANCE TO SECOND CLOSEST STATION EXCEEDS D2FAR KM
      IF (D2.GT.D2FAR) THEN
        J=KINDX(K2)
        WRITE (6,1003) KYEAR2,KMONTH,KDAY,KHOUR,KMIN,
     2  STANAM(J),JNET(J),JCOMP3(J),D2,D2FAR

        IF (LPRT) WRITE (15,1003) KYEAR2,KMONTH,KDAY,KHOUR,KMIN,
     2  STANAM(J),JNET(J),JCOMP3(J),D2,D2FAR

1003    FORMAT (' *** STOP ITERATING EVENT ',I4,4I3,
     2  ' BECAUSE DISTANCE TO SECOND CLOSEST'/
     3  ' *** STATION ',A5,'-',A2,'-',A3,' IS',
     4  F6.0,' & EXCEEDS',F6.0,' KM.')
        IRES=-53
        DONE=.TRUE.
      END IF

C++++++ LOOP OVER REPORTING STATIONS TO CALC DIST WTS & LOG IN TIMES ++++++

      DO 50 K=1,KSTA
      J=KINDX(K)

C-- APPLY DISTANCE AND STATION WEIGHTS
C--CHECK MAIN DISTANCE WEIGHTING FIRST (FOR EX. ITERATIONS 4-30),
C  THEN FIRST DISTANCE WEIGHTING SECOND (FOR EX. ITERATIONS 1-3)
      WFAC=.1*JPSWT(J)

      IF (DIS(K).GE.TEMP .AND. (ITR.GE.ITRDIS .OR. HYPOFIX) .AND.
     2 LFREEW) THEN
        IF (DIS(K).GT.TEMP2) THEN
          WFAC=0.
        ELSE
          WFAC=(.5*COS(PI*(DIS(K)-TEMP)/(TEMP2-TEMP))+.5)*WFAC
        END IF

      ELSE IF (DIS(K).GE.TEMP5 .AND. (ITR.GE.ITRDI1 .OR. HYPOFIX)
     2 .AND. LFREEW) THEN
        IF (DIS(K).GT.(TEMP4*DISW21)) THEN
          WFAC=0.
        ELSE
          WFAC=(.5*COS(PI*(DIS(K)-TEMP5)/(TEMP6-TEMP5))+.5)*WFAC
        END IF
      END IF

C--DECODE P & S WEIGHTS
      LSWT=KWT(K)/10
      LPWT=KWT(K)-10*LSWT
C--IF ONLY A MAGNITUDE IS TO BE COMPUTED (BUT THERE IS NO P OR S), OR IF THERE
C  IS A SHADOW CARD, RESERVE A PLACE FOR STATION BY PRETENDING TO LOG IN A P.
      LPKTMP=KPRK(K).NE.'  '
      IF (LPKTMP .OR. LSHAD .OR. ((AMPK(K).NE.0. .OR.
     2 FMPK(K).GT.0.) .AND. KSRK(K).EQ.'  ')) THEN

C--LOG IN P ARRIVAL TIMES AND WEIGHTS
        M=M+1
        IF (M.GT.MMAX) GOTO 52
        W(M)=0.

C--USE A FLEXIBLE WEIGHT SCHEME INSTEAD OF THE TRADITIONAL 0.25 WEIGHT STEPS
        IF (LPWT.LT.4 .AND. LPKTMP) W(M)= WTVALS(LPWT+1) *WFAC
        IF (W(M).GT.0.) NWR=NWR+1
        IND(M)=K
C--FLAG P ARRIVALS
        MSFLAG(M)=.FALSE.
      END IF

      IF (KSRK(K).NE.'  ') THEN
C--LOG IN S ARRIVAL TIMES AND WEIGHTS
        M=M+1
        IF (M.GT.MMAX) GOTO 52
        W(M)=0.

C--USE A FLEXIBLE WEIGHT SCHEME INSTEAD OF THE TRADITIONAL 0.25 WEIGHT STEPS
        IF (LSWT.LT.4) W(M)= WTVALS(LSWT+1) *WFAC*SWT
        IF (W(M).GT.0.) NWR=NWR+1
        IND(M)=K
C--FLAG S ARRIVALS
        MSFLAG(M)=.TRUE.
      END IF

C--END OF STATION LOOP
50    CONTINUE
      GOTO 54

C--COME HERE IN CASE OF TOO MANY PHASES
52    M=MMAX
C--FIND STATION INDEX
      K=IND(M)
      J=KINDX(K)
      IF (LPRT .AND. ITR.EQ.1)
     2 WRITE (15,1004) STANAM(J),JNET(J),JCOMP3(J)
1004  FORMAT (' TOO MANY PHASES. LAST STATION INCLUDED IS ',
     2 A5,'-',A2,'-',A3)

C---> LMB: Should this be 4 or more stations instead of arrivals?
54    IF (M.LT.4) GO TO 500

C++++++++++++++ CALC ALL TRAVEL TIMES AND DERIVS ++++++++++++++++++++

C--PASS TRAVEL TIMES & DERIVATIVES IN COLS OF ARRAY A:
C--THIS USES ARRAY A TWICE TO SAVE SPACE
C  A(I,1)  AIN   ANGLE OF INCIDENCE
C  A(I,2)  T     TRAVEL TIME
C  A(I,3)  DTDR  TT DERIV WRT DISTANCE
C  A(I,4)  DTDZ  TT DERIV WRT DEPTH

C--GET ALL TRAVEL TIMES & DERIVATIVES, INCLUDING MULTIPLE MODELS
      CALL HYTRA

C--LOOP OVER ALL ARRIVAL TIMES TO CALC EMERGENCE ANGLES, RESIDUALS, ETC.
      DO 60 IM=1,M

C--FIND STATION INDEX AND S FLAG
      K=IND(IM)
      LPS=MSFLAG(IM)
      J=KINDX(K)

C--RECORD EMERGENCE ANGLE AND AZIMUTH
C--MAKE THEM THE SAME SIGN TO FIT BOTH IN A 16-BIT WORD
      KTEMP=KAZEM(K)/180
      AZM=KTEMP/RDEG
      TEMP=A(IM,1)
      IF (KTEMP.LT.0) TEMP=-TEMP
      KAZEM(K)=KTEMP*180+TEMP
C--CALCULATED TRAVEL TIME IN UNITS OF 0.01 SEC
C      MTCAL(IM)=100.*A(IM,2)

C--FIX FOR INFINITE LOOP IN HYLST CAUSED BY ARITHMETIC CONVERSION ERROR
C--HERE: TRUNCATE CALCULATED TRAVEL TIME TO FIT IN MTCAL(IM) (THERE WILL
C--STILL BE PROBLEMS WITH THE LOCATION, WHICH WILL HOPEFULLY BE OBVIOUS
C--WHEN THE ANALYST SEES THE MESSAGES PRINTED BELOW)
C  LMB 1.43
C--ORIGINAL CODE TRUNCATED TOWARDS ZERO; NOW WE ROUND TO NEAREST INTEGER
      LTMP = NINT( 100.* A(IM,2) )
      IF (LTMP .GT. TTMAX) THEN
        LTMP = TTOVER
      ELSE IF (LTMP .LT. TTMIN) THEN
        LTMP = TTUNDR
      END IF
      IF (LTMP .NE. NINT( 100.* A(IM,2) )) THEN
        IF (MSFLAG(IM)) THEN
          C3=KSRK(K)
        ELSE
          C3=KPRK(K)
        END IF
        IF (LPRT) WRITE (15,1036) STANAM(J), C3
        IF (LERR) WRITE (6,1036) STANAM(J), C3
1036    FORMAT (' *** ', A,
     1          ': CALCULATED ', A, ' TRAVEL TIME TRUNCATED TO FIT')
      END IF
      MTCAL(IM)= LTMP

C+++++++++++++++++ CALC DELAYS AND RESIDUALS ++++++++++++++++++

C--CALC DELAYS, OBS ARRIVAL TIME, AND VEL INDEPENDENTLY FOR P AND S
C--KP, JPD, KS, TEMP, TEMP2, DLY, PSB IN UNITS OF 0.01 SEC
C--PSFAC IS THE SLOPE OF THE P OVER S LINEAR VELOCITY RELATION FOR
C  USE WITH PARTIAL DERIVATIVES
      IF (.NOT.LPS) THEN
C--THIS IS A P WAVE
        PSFAC=1.
        PSB=0.
        TEMP2=KP(K)
        IF (LMULT) THEN
          DLY=0.
          DO I=1,NMOD
            IT=MODS(I)
            DLY=DLY+ WMOD(I)*JPD(IT,J)
          END DO
        ELSE
          DLY=JPD(MOD,J)
        END IF
      ELSE

C--THIS IS AN S WAVE
C--FOR AN S MODEL, POSM SHOULD BE 1, POSB SHOULD BE 0
        TEMP2=KS(K)
        IF (LMULT) THEN
          PSFAC=0.
          PSB=0.
          DLY=0.
          DO I=1,NMOD
            IT=MODS(I)
            PSFAC=PSFAC+ WMOD(I)*POSM(IT)
            PSB=PSB+ 100.*WMOD(I)*POSB(IT)
            DLY=DLY +(WMOD(I)*POSM(IT)*JPD(IT,J))
          END DO
        ELSE
          PSFAC=POSM(MOD)
          PSB=100.*POSB(MOD)
          DLY=POSM(MOD)*JPD(MOD,J)
        END IF
      END IF

C--CALC TRAVEL TIME RESIDUAL
C--MTCAL, TEMP2, PSB, DLY IN UNITS OF 0.01 SEC. 
C--PSFAC IS A CLOSE APPROX TO WEIGHTING DIFFERENT MODELS POS VALUES
      R(IM)=.01*(TEMP2 -(MTCAL(IM)*PSFAC+PSB) -DLY)-T1

C--SET PARTIAL DERIVATIVES
C--(1)=DTDT,(2)=DTDY,(3)=DTDX,(4)=DTDZ
      A(IM,1)=1.
      TEMP=PSFAC*A(IM,3)
C--SET LAT & LON PARTIAL DERIVATIVES
      A(IM,2)=-TEMP*COS(AZM)
      A(IM,3)= TEMP*SIN(AZM)
C--SET DT/DZ, WHICH MAY BE UNUSED IF DEPTH IS FIXED
      A(IM,4)=PSFAC*A(IM,4)

C--END OF PHASE LOOP
60    CONTINUE

C++++++++++++++++ RESIDUAL WEIGHTING ++++++++++++++++++++

C--LOOP TO APPLY RESIDUAL WEIGHTING NRES TIMES
C--NRES WILL APPROXIMATE THE NUMBER OF VERY LARGE RESIDUALS
C  TO BE DISCARDED BEFORE INVERTING FOR THE HYPOCENTER
C--IF NRES=1, ONLY 1 VERY LARGE RESIDUAL WILL BE DISCARDED,
C  LEAVING ERRORS IN THE REMAINING DATA TO CONTAMINATE THE SOLUTION

      DO 80 I=1,NRES
      RMSWT=0.
      WNORM2=0.
      NWR=0

C--ACCUMULATE RMSWT AND WEIGHTS
      DO IM=1,M
        WNORM2=WNORM2+W(IM)**2
        RMSWT=RMSWT+(W(IM)*R(IM))**2
      END DO
      IF (WNORM2.EQ.0.) GO TO 500
      RMSWT=SQRT(RMSWT/WNORM2)
C--DON'T SHRINK RESIDUAL WEIGHTING LIMITS IF RMSWT BECOMES TOO SMALL
      TEMP2=RMSWT
      IF (TEMP2.LT.RMSCUT) TEMP2=RMSCUT
      TEMP=(RMSW2-RMSW1)*TEMP2
      WNORM=0.

C--LOOP TO CALC AND APPLY RESIDUAL WEIGHTING FOR ALL ARRIVALS
      DO 70 IM=1,M
      IF ((ITR.GE.ITRRES .OR. HYPOFIX) .AND. LFREEW) THEN 
        IF (W(IM).EQ.0.) GO TO 70
        RES=(ABS(R(IM))-RMSW1*TEMP2)/TEMP
        IF (RES.GT.1.) THEN
          W(IM)=0.
          GOTO 70
        END IF
        IF (RES.GT.0.) W(IM)=W(IM)*(.5+.5*COS(PI*RES))
      END IF
      WNORM=WNORM+W(IM)
      NWR=NWR+1
70    CONTINUE
80    CONTINUE

      WFAC=NWR/WNORM
      WNORM2=0.
      RMS=0.
      NWS=0
      NWR=0

C--NOW NORMALIZE WEIGHTS AND APPLY THEM TO PARTIAL DERIVATIVES
      DO 95 IM=1,M

C--NORMALIZE WEIGHTS
        W(IM)=W(IM)*WFAC
        WNORM2=WNORM2+W(IM)**2
        RMS=RMS+(R(IM)*W(IM))**2
C--COUNT WEIGHTED ARRIVALS
        IF (W(IM).GT..1) THEN
          NWR=NWR+1
C--FIND S FLAG
          LPS=MSFLAG(IM)
          IF (LPS) NWS=NWS+1
        END IF

C--APPLY WEIGHTS TO PARTIAL DERIVS
        DO 90 I=1,N
90      A(IM,I)=A(IM,I)*W(IM)

C--LOAD THE DATA VECTOR R INTO THE N+1ST COL OF A FOR INVERSION IN SOLVE
95    A(IM,N+1)=R(IM)*W(IM)

      RMS=SQRT(RMS/WNORM2)

C--TEST TO SEE IF WE STILL HAVE ENOUGH DATA AFTER ALL WEIGHTING
      IF (NWR.LT.4 .AND. LFREEW) THEN
        LFREEW=.NOT.LJUNK
C--OPTIONALLY FORCE A SOLUTION WITHOUT DIST & RESIDUAL WEIGHTING
        IF (.NOT.LFREEW) THEN
          IF (LMESS) WRITE (15,1005)
1005      FORMAT (' DISTANCE & RESIDUAL WEIGHTING CANCELLED')
          IF (LREP) WRITE (6,1008) KYEAR2,KMONTH,KDAY,KHOUR,KMIN
1008      FORMAT (' *** DISTANCE & RESIDUAL WEIGHTING CANCELLED',
     2    ' FOR EVENT: ',I4,4I3)
        ELSE
          GOTO 500
        END IF
      END IF

C--SET DONE FLAG IF HYPOCENTER IS FIXED 
      IF (HYPOFIX) DONE=.TRUE. 

C--MAKE THIS THE LAST ITER IF 1) CHANGE IN RMS OR ADJ VECTOR IS SMALL
C  (ONLY IF THE DEPTH HAS BEEN FREED) OR 2) THE MAX NO. OF ITERATIONS 
C  IS REACHED (AND THE RMS DID NOT INCREASE LAST TIME)
C--BUT CONTINUE ITERATING IF DIS & RES WEIGHTS HAVENT BEEN APPLIED.
      IF (ITR.GT.ITRDIS .AND. ITR.GT.ITRRES .AND. .NOT.DONE)
     2 DONE=(((ABS(OLDRMS-RMS).LT.DRQT) .OR.
     3 RR.LT.DQUIT) .AND. (.NOT.LZFIX .OR. DEPFIX))
     4 .OR. (ITR.GT.ITRLIM .AND. RMS.LT.RMSMIN+RBACK)
     
C--GO BACK TOWARD OLD SOLUTION IF RMS INCREASED LAST TIME
      IF (RMS.GT.RMSMIN+RBACK .AND. .NOT.DONE) THEN
        JBAC=JBAC+1
        OLDRMS=RMS
        IF (RMS.LT.RMSMIN) RMSMIN=RMS
        IF (JBAC.EQ.1) THEN
C--BACKUP TOWARD OLD SOLUTION FOR FIRST TIME
          DO 105 I=1,4
105       Y(I)=-BACFAC*Y(I)
          RR=RR*BACFAC
        ELSE
C--BACKUP TOWARD OLD SOLUTION FOR SECOND OR LATER TIME
          DO 100 I=1,4
100       Y(I)=(1.-BACFAC)*Y(I)
          RR=RR*(1.-BACFAC)
        END IF
        IF (LMESS) WRITE (15,1006) BACFAC
1006    FORMAT (' RMS INCREASE - MOVE HYPO ',F4.2,' BACK')
        CONERR=.TRUE.
        GOTO 130
      END IF

C+++++++++++++++ INVERT FOR HYPO ADJUSTMENT ++++++++++++++++++++++++++

C--DO INVERSION FOR 1 ITER, FINDING OTHER USEFUL MATRICES IF REQUESTED
C--IF DEPTH HELD FIXED, N=3. N=4 OTHERWISE.
      CALL HYSOL (N,Y,NFREE,LZFIX)

C--SET CONVERGENCE ERROR FLAG (USED TO POST A '#' REMARK)
      CONERR=.FALSE.

C--ALLOW DEPTHS USING HYPOELLIPSE GRADIENT MODELS TO BECOME NEGATIVE, BUT
C  ONLY UP TO THE REFERENCE (MAX) ELEVATION
      IF (MODTYP(MOD).EQ.2 .OR. MODTYP(MOD).EQ.3
     2 .OR. MODTYP(MOD).EQ.4) THEN

C--PREVENT HYPO FROM ITERATING ABOVE REFERENCE ELEVATION
C        IF (Z1+Y(4) .LT. -ELEVMX(MOD)) THEN
        IF (Z1+Y(4) .LT. -ELEVMAX) THEN		!FWK 1.37

          IF (LMESS) WRITE (15,*)
     2    ' KEEP QUAKE BELOW REFERENCE ELEVATION'
          CONERR=.TRUE.
C          TEMP=-DZAIR*(ELEVMX(MOD)+Z1)/Y(4)
C          Y(4)=-DZAIR*(ELEVMX(MOD)+Z1)

          TEMP=-DZAIR*(ELEVMAX+Z1)/Y(4)	!FWK 1.37
          Y(4)=-DZAIR*(ELEVMAX+Z1)	!FWK 1.37
          DO I=1,3
            Y(I)=Y(I)*TEMP
          END DO
        END IF

      ELSE
C--ABORT PROGRAM IF A NEGATIVE DEPTH IS FOUND 
C  (BY TRIAL DEPTH OR MIGRATION FROM A CRE MODEL)

        IF (Z1.LT.0.) THEN
          IF (MODTYP(MOD).EQ.0) C1='T'
          IF (MODTYP(MOD).EQ.1) C1='H'
          WRITE (6,1060) Z1,C1,MOD,KYEAR2,KMONTH,KDAY,KHOUR,KMIN,T1
1060      FORMAT('*** ILLEGAL NEGATIVE DEPTH',F8.2,
     2    ' USED FOR CR',A1,' MODEL # ',I3,/
     3    ' FOR EVENT',I5,4I3,' SEC=',F7.2,'. STOP PROGRAM.'/
     4    ' NOTE THAT TRIAL DEPTH IS MODEL DEPTH, NOT GEOID DEPTH.')
          IF (LPRT) WRITE (15,1060) 
     2    Z1,C1,MOD,KYEAR2,KMONTH,KDAY,KHOUR,KMIN,T1
          STOP
        END IF

C--HANDLE LAYER & GRADIENT MODELS AS BEFORE PREVENTING AIRQUAKES
C--PREVENT HYPO FROM ITERATING INTO THE AIR
C        IF (Z1+Y(4).LT.0.) THEN
        IF (Z1+Y(4) .LT. -ELEVMAX) THEN	!FWK 1.37

          IF (LMESS) WRITE (15,'(" AIRQUAKE PREVENTED")')
          CONERR=.TRUE.
C          TEMP=-DZAIR*Z1/Y(4)
C          Y(4)=-DZAIR*Z1

          TEMP=-DZAIR*(ELEVMAX+Z1)/Y(4)	!FWK 1.37
          Y(4)=-DZAIR*(ELEVMAX+Z1)	!FWK 1.37
          DO I=1,3
            Y(I)=Y(I)*TEMP
          END DO
        END IF
      END IF

C--DAMP DEPTH ADJUSTMENT IF TOO LARGE
      IF (ABS(Y(4)).GT.DZMAX) THEN
        IF (LMESS) WRITE (15,1220)
1220    FORMAT (' DEPTH ADJUSTMENT DAMPED')
        CONERR=.TRUE.
        Y(4)=Y(4)*DZMAX/(ABS(Y(4))+DZMAX)
      END IF

C--SAVE RMS AND COMPUTE ADJUSTMENT VECTOR LENGTH
      OLDRMS=RMS
      IF (RMS.LT.RMSMIN) RMSMIN=RMS
      JBAC=0
      RR=0.
      DO 120 I=2,4
120   RR=RR+Y(I)**2
      RR=SQRT(RR)

C--DAMP DISTANCE ADJUSTMENT
      IF (RR.GT.DXMAX) THEN
        CONERR=.TRUE.
        TEMP=DXMAX/RR
        IF (LMESS) WRITE (15,1233) TEMP
1233    FORMAT (' CURB DISTANCE STEP. ALL ADJUSTMENTS DAMPED BY',F6.3)
        DO I=1,4
          Y(I)=Y(I)*TEMP
        END DO

C--RE-COMPUTE ADJUSTMENT VECTOR LENGTH
        RR=0.
        DO I=2,4
          RR=RR+Y(I)**2
        END DO
        RR=SQRT(RR)
      END IF

C+++++++++++++++ PRINT INVERSION INFO +++++++++++++++++++++

C--GET COORDS IN DEG AND MIN, AND WRITE LOC FOR THIS ITERATION
C--LOCATION IS THAT PRIOR TO ADJUSTMENTS GIVEN
130   IS=' '
      IE=' '
C--LATITUDE
      IF (CLAT.LT.0.) IS='S'
      TEMP=ABS(CLAT)
      LAT=TEMP
      XLTM=60.*(TEMP-LAT)

C--LONGITUDE
C--BRING EQS WEST OF WEST HEMISPHERE INTO EAST HEMISPHERE V.1.38
      TEMP7=CLON
      IF (TEMP7.GT.180.) TEMP7=TEMP7-360.
      IE='W'
      IF (TEMP7.LT.0.) IE='E'
      TEMP=ABS(TEMP7)      
      LON=TEMP
      XLNM=60.*(TEMP-LON)

C--OUTPUT ADJUSTMENTS TO PRINT FILE
      IF (LMESS) THEN
C--PRINT HEADING
        IF (KPRINT.GT.4 .AND. ITR.NE.1) WRITE (15,1010)
1010    FORMAT (56X,'ADJUSTMENTS (KM)'/'  I ORIGIN   LAT N',5X,'LON W',
     2  5X,'Z   NWR  RMS    DT    DLAT   DLON    DZ    RR NF MOD')
C--PRINT ADJUSTMENT
        WRITE (15,1020) ITR,T1, LAT,IS,XLTM, LON,IE,XLNM,
     2  Z1,NWR,RMS, Y,RR, NFREE,CRODE(MOD)
1020    FORMAT (1X,I2,F6.2, I4,A1,F5.2, I5,A1,F5.2,
     2  F6.2,I4,F5.2, 4F7.3,F6.3, I2,1X,A3)
      END IF

C--IF HYPOCENTER IS FIXED 
      IF (ITRLIM.EQ.0 .OR. HYPOFIX) THEN 
        WNORM=0.
        TEMP=0.
C--FIND AVERAGE RESIDUAL
        DO IM=1,M
          WNORM=WNORM+W(IM)
          TEMP=TEMP+R(IM)*W(IM)
        END DO

C--ADJUST ORIGIN TIME IF OT IS NOT FIXED
        IF (.NOT.ALLFIX) THEN
          IF (LMESS) WRITE (15,1221)
1221      FORMAT (' AVERAGE RESIDUAL REMOVED; O.T. ADJUSTED')

          DT=TEMP/WNORM
          T1=T1+DT
C--ADJUST RESIDUALS BY SAME TIME SHIFT AND RECOMPUTE RMS.
C  RESIDUALS WILL BE RECALCULATED IN HYLST.
          RMS=0.
          WNORM2=0.
          DO IM=1,M
            R(IM)=R(IM)-DT
            WNORM2=WNORM2+W(IM)**2
            RMS=RMS+(R(IM)*W(IM))**2
          END DO
C--CORRECT VALUE OF RMS
          RMS=SQRT(RMS/WNORM2)
        END IF
      END IF

      IF (DONE) GOTO 510

C--OPTIONALLY LIST ALL STATIONS NOW
      IF (KPRINT.GT.5 .AND. LPRT) CALL HYLST

C+++++++++++++++++++++ ADJUST HYPOCENTER +++++++++++++++++++++

      T1=T1+Y(1)
      CLAT=CLAT+Y(2)/111.19
      CLON=CLON+Y(3)/(111.19*COS(CLAT/RDEG))
      Z1=Z1+Y(4)

C--FREE THE FOCAL DEPTH ON NEXT ITER IF ADJUSTMENT IS SMALL ENOUGH
C  AND THE DEPTH IS NOT PERMANENTLY FIXED
C      NLAST=N
      N=3
      IF (RR.LT.DXFIX .AND. .NOT.DEPFIX) THEN
        IF (LMESS .AND. LZFIX) WRITE (15,1230)
1230    FORMAT (' FREE DEPTH')
        LZFIX=.FALSE.
      END IF
C--SOLVE FOR DEPTH IF IT IS NOT FIXED
      IF (.NOT.LZFIX) N=N+1

      GOTO 30

C*************** END ITERATION LOOP ***************************

C--ABANDON SOLUTION WITH INSUFFICIENT DATA
C--OUTPUT STATIONS HERE, SINCE IT WILL NOT BE DONE FROM MAIN PROG
500   IF (LPRT) WRITE (15,1052) NWR,INUM,IDNO
      IF (LERR) WRITE (6,1052) NWR,INUM,IDNO
1052  FORMAT (' *** ',I1,' PHASES CANT SOLVE SEQUENCE NO.',
     2 I5,', ID NO.',I10)
      IRES=-51

C--FLAG THIS EVENT IF IT HAD CONVERGENCE PROBLEMS
510   IF (ITR.GT.ITRLIM .OR. CONERR .OR. .NOT.LFREEW) RMK2='#'
C--FLAG EVENT IF DEPTH HELD FIXED
      IF (NFREE.LT.4) RMK2='-'
C--FLAG EVENT IF HYPOCENTER HELD FIXED
      IF (HYPOFIX) RMK2='X' 
      IF (ALLFIX) RMK2='O' 
      RETURN
      END
