source: trunk/src/seismic_processing/hyp2000/hyp.for @ 3172

C--HYP IS THE VAX/SUN VERSION OF THE LOCATION PROGRAM HYPOINVERSE. THE PROGRAM
C  IS DESIGNED TO BE FAST, FILE ORIENTED AND COMMAND DRIVEN.

C--WRITTEN BY FRED KLEIN. SEE GREETING IN HYBDA.FOR FOR VERSION #

C++++++++++++++++ LIST OF REQUIRED SUBROUTINES +++++++++++++++++
C--SUBROUTINES MARKED WITH * REQUIRE THE COMMON BLOCK INCLUDE FILE 'common.inc'
C  'common.inc' INCLUDES THE FILE 'integer.for' WHERE ARRAYS ARE DECLARED
C    INTEGER*2 AND LOGICAL*2. *4 VARIABLES MAY BE USED IF NECESSARY.
C--SUBROUTINES MARKED WITH & HAVE DIFFERENT VERSIONS ON THE SUN, OS2 AND VAX.
C--USE *.FOR ON ALL SYSTEMS IF IT IS THE ONLY FILE THAT EXISTS.
C--USE *.F ON SUN AND *.FOR ON THE VAX IF BOTH EXIST.
C--USE *.OS2 ON OS2 IF IT EXISTS, OTHERWISE *.F THEN *.FOR IN THAT ORDER.
C--USE integer.os2 INSTEAD OF integer.for ON OS2.
C  HYP     * MAIN PROGRAM.
C  HYBDA   * BLOCK DATA INITIALIZATION OF COMMON.
C  HYBEG  &* INITIALIZATION OF OTHER VARIABLES. 
C  HYCMD   * GETS AND PROCESSES COMMANDS.
C  HYSTA   * READS IN STATIONS.
C  HYDEL   * READS IN STATION DELAYS (FOR MULTIPLE MODELS).
C  HYATE   * READS IN STATION ATTENUATION HISTORY.
C  HYCAL   * READS IN STATION CALIBRATION FACTOR HISTORY.
C  HYFMC   * READS IN STATION FMAG CORRECTIONS.
C  HYXMC   * READS IN STATION XMAG CORRECTIONS.
C  HYCRH   * READS IN HOMOGENOUS LAYER CRUSTAL MODELS.
C  HYCRT   * READS IN TRAVEL-TIME-TABLE CRUSTAL MODELS.
C  HYSTL   * OUTPUTS STATIONS, CRUST & PARAMETERS TO PRINT FILE.
C  HYOPEN  * OPENS FILES FOR LOCATION RUN.
C  HYINIT  * INITIALIZES SOME VARIABLES FOR LOCATION RUN.
C  HYPHS   * READS IN PHASE DATA FOR ONE EVENT.
C  HYCIN  &* INPUTS PHASE DATA FROM CUSP MEM FILES (ALTERNATE TO HYPHS)
C  HYCOUT &* OUTPUTS RESULTS TO THE CUSP MEM FILE
C            (HYCIN & HYCOUT ARE NOT USED IN UNIX/SUN VERSION)
C  HYTRL   * SETS TRIAL HYPOCENTER.
C  HYLOC   * LOCATES ONE EVENT.
C  HYSOL   * PERFORMS INVERSION FOR ONE ITERATION.
C  HYSVD  &  CANNED SINGLE-VALUE-DECOMPOSITION ROUTINE.
C            (SEE NOTES IN HYSVD.F AND FPE-TRAPS COMMENT FOR SPECIAL SUN MODS)
C  HYTRA   * MANAGE CRUST MODEL CHOICE & AVERAGING.
C  HYTRH   * CALC TRAVEL TIMES AND DERIVS FOR HOMO LAYER MODEL.
C  HYTRT   * CALC TRAVEL TIMES AND DERIVS FROM TRAV-TIME TABLE.
C  HYMAG   * COMPUTES ALL MAGNITUDES.
C  HYMAGP  * COMPUTES P AMPLITUDE MAGNITUDES.
C  HYPREF  * SELECTS THE PREFERRED MAGNITUDE FROM ALL AVAILABLE.
C  HYREP   * REPORTS A LOCATION ON THE TERMINAL.
C  HYSOU   * TABULATES MOST COMMON DATA SOURCE CODES
C  HYLST   * OUTPUTS DATA BY STATION TO PRINT & ARCHIVE FILES.
C  HYSUM   * OUTPUTS SUMMARY RECORD (FINAL LOCATION).
C  HYINP     FOR INTERACTIVE ENTRY OF PHASE DATA.
C  HYPRO   * INTERACTIVELY PROCESSES A SERIES OF EVENTS.
C  MEDWT     COMPUTES THE WEIGHTED MEDIAN OF A SERIES OF MAGNITUDES.
C  HYDELT  & DELETES FILES IN INTERACTIVE PROCESSING
C  HYEDTI  & RUNS AN EDTIOR WITHIN A PROCESS
C  HYTIME  & GETS CURRENT SYSTEM TIME FOR LABELING PRINT FILE
C  UTMCAL    COMPUTES STATION DISTANCES ON A UTM GRID

C--GENERAL PURPOSE SUBROUTINES
C  KLAS      ASSIGNS A NAME AND NUMBER TO AN EVENT BASED ON LOCATION.
C  KLASS (NET 1), BOX2 (NET2), BOX3 (NET3), KSIC - USED BY KLAS.
C  UPSTR     CONVERTS A STRING TO UPPER CASE.
C  JASK      INTERACTIVE PROMPT & ENTRY OF AN INTEGER.
C  ASKC      INTERACTIVE PROMPT AND ENTRY OF A STRING.
C  ASKR      INTERACTIVE PROMPT AND ENTRY OF A REAL VALUE.
C  LASK      INTERACTIVE PROMPT AND ENTRY OF A LOGICAL VALUE.
C  LENG      DETERMINES THE NON-BLANK LENGTH OF A STRING.
C  DAYJL     RETURNS PERPETUAL JULIAN DAY FOR A DATE.
C  JDATE     RETURNS DATE FOR A PERPETUAL JULIAN DAY.
C  OPENR  &  OPENS A FILE FOR READING.
C  OPENW  &  OPENS A FILE FOR WRITING.
C  ERRSET &  VAX SYSTEM SUBROUTINE ONLY: CONTROLS ERROR MESSAGES ON OVERFLOWS.
C               (A DUMMY ERRSET.F IS SUPPLIED WITH THE UNIX VERSION)
C  SPAWN  &  SPAWNS AN OPERATING SYSTEM COMMAND.
C  READQ  &  READS AN ASCII RECORD AND RETURNS ITS LENGTH.
C  GETENV &  ON UNIX, RETURNS ENVIRONMENT VAR W/NAME OF INI COMMAND FILE
C               (A DUMMY GETENV.VAX IS SUPPLIED FOR VAX & OS2 VERSIONS)

C--CUSP SUBROUTINES
C  MEM_DUMP      READS A CUSP MEM FILE AND PUTS DATA INTO STRUCTURES
C  OPHASE      PARSES REMARK, WEIGHT & FIRST MOT FROM PHASE DESCRIPTOR
C  MEM_EQ_UPDATE  PUTS HYPOCENTER & STA INFO INTO CUSP MEMORY (& MEM FILE)

C--DIFFERENCES BETWEEN THE VAX AND SUN/UNIX VERSIONS:
C--WHERE THEY DIFFER, THE SUBROUTINE SOURCE CODE FILES THAT END IN .F ARE
C  FOR SUN; THOSE ENDING IN .FOR ARE FOR VAX.  FILES FOR WHICH THERE IS
C  ONLY A .FOR VERSION WILL COMPILE ON EITHER MACHINE.
C--ROUTINES WITH DIFFERENT VERSIONS ARE HYBEG, HYDELT, HYEDIT, HYTIME,
C  SPAWN, INIT_EVENT, HYCIN, OPENR AND OPENW.
C--HYBEG INITIALIZES FILENAMES AND DEVICES THAT ARE SYSTEM SPECIFIC.

C--THE FOLLOWING "NON-ANSI" FORTRAN FEATURES ARE USED (THESE WERE FLAGGED
C  BY SUN'S FORTRAN COMPILER WHEN THE -ansi COMMAND FLAG WAS USED):
C    OPTIONAL INTEGER*2 AND LOGICAL*2 VARIABLES IN COMMON (SEE INTEGER.FOR)
C    INCLUDE STATEMENT
C    DO ... END DO STATEMENTS
C    ! TO BEGIN END OF LINE COMMENTS ('common.' FILE ONLY)
C    SUBROUTINE NAMES (HYPOINV) LONGER THAN 6 CHARACTERS
C    LIST-DIRECTED FORMATTING FROM AN INTERNAL STRING
C    CHARACTER*(*) IN CONCATENATION

C--FPE (FLOATING POINT EXCEPTION) TRAPS.
C  THE SUN VERSION, WHEN PRESENTED WITH UNDERDETERMINED EARTHQUAKES WITH FEW
C  READINGS, WOULD SOMETIMES ATTEMPT A ZERO / ZERO OPERATION IN HYSVD. THE
C  SOLUTION PROGRAMMED INTO HYSVD.F WAS TO TEST THE DIVIDEND AND DIVISOR
C  BEFORE DIVISION AND TO RETURN A ZERO QUOTIENT IF BOTH WERE 0. THE VAX DOES 
C  THIS AUTOMATICALLY. 0/0 ON THE SUN YIELDS AN IEEE NaN (NOT A NUMBER) WHICH
C  CONTAMINATES ALL SUCCEEDING VARIABLES THAT DEPEND ON THIS NUMBER. WHEN
C  EVENTUALLY USED AS A SUBSCRIPT, THE PROGRAM HANGS UNTIL STOPPED WITH ^C.
C--HERE ARE SUGGESTIONS TO TRAP FPE'S THAT WERE USED TO DETECT THIS 0/0 FPE.
C  ON THE SUN, COMPILE f77 WITH THE -g OPTION TO STORE LINE NUMBERS IN THE
C  SOURCE CODE. SUN DOES NOT PERMIT COMPILING THE MAIN PROGRAM WITH BOTH -g
C  AND A COMMON BLOCK, SO USE A DUMMY MAIN PROGRAM HYPMAIN.F:
C       CALL TRAPFPE
C       CALL HYPM
C       END
C--THE TRAPFPE SUBROUTINE ENABLES A IEEE HANDLER WHICH IS CALLED WHEN AN FPE
C  EXCEPTION OCCURS (0/0, OVERFLOW, ETC). THE HANDLER PRINTS THE HEX ADDRESS OF
C  THE CODE GENERATING THE EXCEPTION, THEN USE THE dbx DEBUGGER TO FIND THE
C  SOURCE CODE LINE NUMBER. ALTERNATIVELY, USE THE dbxtool WITH CODE COMPILED
C  WITH THE -g OPTION, AND GIVE THE dbx COMMAND
c       catch FPE
C       dbxenv case insensitive ALLOWS EXAMINING HI VARIABLES WITH dbxtool.
C--SEE SUN'S DEBUGGING TOOLS AND FORTRAN NUMERICAL COMPUTATION DOCUMENTATION
C  FOR MORE INFO.

C--REMOVAL OF UNNEEDED SUBROUTINES AND DATA TO SAVE MEMORY SPACE:
C--IF CUSP DATA WILL NOT BE USED (JCP 6 OR 7):
C    ELIMINATE THE CALLS TO INIT_EVENT AND HYCIN FROM HYP. ALSO ELIMINATE THE
C    SUBROUTINES CALLED BY HYCIN (SUCH AS GET_* & THE CUSP LIBRARY).
C    THE COMMAND FID WILL NOT THEN BE NEEDED. ALSO ELIMINATE THE VARIABLES
C    IRES, LCUSP & FORID FROM COMMON.
C--IF SHADOW PHASE FORMATS WILL NOT BE USED (JCP 4 & 5):
C    ELIMINATE THE VARIABLES KSHAD, KLSHA, LENSHA, SHADO, LSHA1 & SHAD1
C    FROM COMMON, AND REFERENCES TO THEM IN HYPHS AND HYLST.
C--THE NUMBER OF STATIONS ARE CONTROLED BY THE COMMON PARAMETERS
C    MAXSTA > MAXPHS > MMAX. THESE CAN BE MADE SMALLER FOR SMALL NETWORKS.

C--VERY PARTIAL VERSION HISTORY
C
C  1978 VERSION 0.x (Eclipse computer, simple & compact)
C  1985 VERSION 0.x (VAX and pro-350 computer, compact)
C  1989 VERSION 0.x (multiple crust models)
C  2002 VERSION 1.0 (Full documentation, completely Y2000 capable)
C  2/2007 VERSION 1.1 (Can fix origin time, many other changes)
C  5/2007 VERSION 1.11 (Can fix origin time, greeting)
C  12/2007 VERSION 1.2 (Can fix origin time; g77 compiler)

C++++++++++++++++ I/O DEVICE NUMBERS USED ++++++++++++++++++++++
C  5  TERMINAL INPUT.
C  6  TERMINAL OUTPUT.
C  7  ARCHIVE OUTPUT FILE.
C  8,9,10,11  INPUT COMMAND FILES.
C 12  SUMMARY OUTPUT FILE.
C 13  STATION DATA FILES (ATTEN, DELAY, FMAG & XMAG CORRECTIONS).
C 14  CRUST, STATION & PHASE INPUT FILES.
C 15  PRINT OUTPUT FILE.
C 16  MAGNITUDE DATA OUTPUT FILE.
C 17  EVENT LIST FILE FOR INTERACTIVE PROCESSING

      INCLUDE 'common.inc'
      CHARACTER XMON(12)*3,CC*1
      DATA XMON/'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP',
     2 'OCT','NOV','DEC'/

C--DONT GIVE ERROR MESSAGES WHEN DATA OVERFLOW OUTPUT FIELDS
      CALL ERRSET (63,.TRUE.,.FALSE.,.FALSE.,.FALSE.)
C--DONT GIVE ERROR MESSAGES WHEN INTEGERS OVERFLOW (BUSTED EVENTS THAT ARE
C  KICKED OUTSIDE THE NETWORK)
      CALL ERRSET (70,.TRUE.,.FALSE.,.FALSE.,.FALSE.)
C--INITIALIZE FLAG TO SUCCESS
      IRES=1

C--SEND A MESSAGE TO THE TERMINAL. THIS ALSO ASSIGNS UNIT 5 TO TERMINAL
      WRITE (6,1000)
1000  FORMAT (' HYPOINVERSE 2000 STARTING')
      WRITE (6,'(A)') GREETING

C--INITIALIZE VARIABLES NOT INITIALIZED IN BLOCK DATA
      CALL HYBEG

C--OPEN AND BEGIN READING THE OPTIONAL STARTUP COMMAND FILE HYPINST
      INP=5
      CALL OPENR (8,INFILE(1),'F',IOS)
      IF (IOS.NE.0) GOTO 5
      INP=8

C--GO GET A COMMAND AND EXECUTE IT. RETURN HERE IF CALLING A SUBROUTINE.
5     CALL HYCMD

C--ISTAT DIRECTS WHICH SUBROUTINE OR SECTION TO INVOKE.
C  ISTAT IS ONLY ASSIGNED A VALUE BY HYCMD.
C  1  CRH  READ LAYER CRUST MODEL
C  2  CRT  READ GRADIENT TRAVEL TIME TABLE
C  3  STA  READ STATION FILE
C  4  INP  INPUT PHASE DATA INTERACTIVELY
C  5  BUG  READ A PHASE FILE ONLY TO CHECK FOR ERRORS
C  7  LOC  LOCATE A PHASE FILE
C  8  STO  STOP HYPOINVERSE, OR RETURN TO MASTER CALLING PROGRAM

      GOTO (10,12,14,74,78,5,84,13), ISTAT
      GOTO 5

C--<CRH> READ A HOMOGENEOUS LAYER CRUSTAL MODEL
10    CALL HYCRH
      CLOSE (14)
      GOTO 5

C--<CRT> READ A LINEAR GRADIENT CRUSTAL MODEL
12    CALL HYCRT
      CLOSE (14)
      GOTO 5

C--<STO> STOP THE PROGRAM
13    IF (SUBMOD) THEN
        CONTINUE
      ELSE
        STOP
      END IF

C--<STA> READ STATION FILE
14    CALL HYSTA
      CLOSE (14)
      GOTO 5

C--<INP> ENTER PHASE DATA MANUALLY INTO A CONDENSED FORMAT FILE.
74    CALL HYINP
      GOTO 5

C--<BUG> READ PHASE FILE ONLY TO CHECK FOR ERRORS
C--OPEN PRINT OUTPUT FILE
78    LPRT=.TRUE.
      IF (PRTFIL(1:4).EQ.'    ') PRTFIL=TERMIN
      IF (LAPP(1)) THEN
        CALL OPENW (15,PRTFIL,'F',IOS,'A')
      ELSE
        CALL OPENW (15,PRTFIL,'F',IOS,'S')
      END IF

C--OPEN PHASE FILE
      CALL OPENR (15,PHSFIL,'F',IOS)
      IF (IOS.NE.0) THEN
        WRITE (6,1008)
1008    FORMAT (' *** ERROR - PHASE FILE NOT FOUND')
        CLOSE (15)
        IRES=-61
        GOTO 5
      END IF

C--LOOP TO READ EVENTS
80    CALL HYPHS
      IF (KEND.EQ.0) GOTO 80
      CLOSE (15)
      CLOSE (14)
      GOTO 5

C++++++++++++++++++++ EARTHQUAKE LOCATION SECTION ++++++++++++++++++

C--<LOC>ATE ALL EARTHQUAKES IN THE SPECIFIED FILE, USING PRESENT PARAMETERS
C--INITIALIZE SOME VARIABLES
84    CALL HYINIT
C--OPEN FILES
      CALL HYOPEN

C--STOP NOW IF THERE IS NO PHASE FILE
      IF (ISTAT2.EQ.0) GOTO 5
C--LIST AVAILABLE STATIONS & MODELS IN PRINT FILE
      CALL HYSTL

      LCUSP=JCP.EQ.6 .OR. JCP.EQ.7
C--READ IN PHASE DATA FOR ONE EVENT
50    IF (LCUSP) THEN
        CALL HYCIN
      ELSE
        CALL HYPHS
      END IF

C--KEND IS SET BY HYPHS DEPENDING ON END-OF-FILE STATUS
C  =-1  END OF FILE, STOP RIGHT AWAY
C  = 0  LOCATE THIS EVENT, THEN READ ANOTHER
C  = 1  END OF FILE, LOCATE THIS EVENT THEN STOP
C--CLOSE FILES & QUIT IF END OF FILE OCCURRED IN PHASE FILE
      IF (KEND.LT.0) GOTO 70

C--SET THE TRIAL HYPOCENTER
      CALL HYTRL

C--PRINT THE EVENT DATE AND TIME AS HEADING
      IF (LPRT .AND. KPRINT.GT.0) THEN
        IF (LEJCT) THEN
          CC='1'
        ELSE
          CC=' '
          WRITE (15,'(1X,21(''####''))')
        END IF

        WRITE (15,1005) CC,KDAY,XMON(KMONTH),
     2  KYEAR2,KHOUR,KMIN,INUM,IDNO
1005    FORMAT (A1,I3,1X,A3,I5,',',I3,':',I2.2,'  SEQUENCE NO.',
     2  I5,', ID NO.',I10)
      END IF

C--LOCATE THE EVENT
      CALL HYLOC

C--ASSIGN A 3-LETTER CODE AND NAME BASED ON LOCATION
C  I IS THE REGION NUMBER, PRESENTLY UNUSED
      IF (NET.GT.0) I=KLAS (NET,CLAT,-CLON,Z1,REMK,FULNAM)

C--CALCULATE THE EARTHQUAKE'S MAGNITUDE
      CALL HYMAG

C--CALCULATE THE EARTHQUAKE'S P AMPLITUDE MAGNITUDE
      CALL HYMAGP

C--SELECT PREFERRED MAGNITUDE
      CALL HYPREF

C--TABULATE DATA SOURCE CODES
      CALL HYSOU

C--GENERATE PRINTED AND ARCHIVE OUTPUT
      CALL HYLST

C--ABORT THE LOOP IF THERE ARE NOT ENOUGH READINGS
      IF (NWR.LT.MINSTA) THEN
        WRITE (6,1002) NWR,KYEAR2,KMONTH,KDAY,KHOUR,KMIN
        IF (LPRT) WRITE (15,1002) NWR,KYEAR,KMONTH,KDAY,KHOUR,KMIN
1002    FORMAT (' *** ABANDON EVENT WITH ONLY',I2,' READINGS:',I4,4I3)
        IRES=-51
        GOTO 50         !FOR HYP (MAIN PROGRAM) ONLY
      END IF

C--WRITE RESULTS TO CUSP MEM FILE
C--JCPO CONTROLS TO WHAT EXTENT RESULTS ARE WRITTEN OUT TO CUSP
C  =0 NOTHING WRITTEN ANYWHERE
C  =1 STRUTURES UPDATED
C  =2 ABOVE PLUS SHARED MEMORY UPDATED
C  =3 ABOVE PLUS MEM FILE RE-WRITTEN
      IF (LCUSP .AND. JCPO.GT.0) CALL HYCOUT

C--OUTPUT SUMMARY DATA USING UNIT NUMBER FOR SUMMARY FILE
      IF (LSUM) CALL HYSUM (12)

C--OUTPUT A MESSAGE ON THE CONSOLE FOR EACH EVENT
      IF (LREP) CALL HYREP

C--END OF LOCATION LOOP
      IF (KEND.EQ.0) GOTO 50

C--CLOSE FILES THEN GET ANOTHER COMMAND
70    CLOSE (12)
      CLOSE (7)
      CLOSE (15)
      CLOSE (14)
      CLOSE (16)
      GOTO 5

      END