C
C Subroutine file 'gse.for' to write and read seismograms in the GSE 1.0 C data exchange format. C C Version: 6.60 C Date: 2012, February 11 C C Coded by: Ludek Klimes C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/klimes.htm C C Reference to GSE 1.0: C GSETT-2 (1990): Instructions for the conduct of the preparatory C test of phase 3 of GSETT-2. C Ad hoc Group of Scientific Experts (GSE) to consider C international cooperative measures to detect and identify C seismic events - Conference Room Paper 190/Rev.4, C Geneva, Switzerland, September 1990. C Reference to GSE 2.0: C GSETT-3 (1995): Formats. C Group of Scientific Experts - Conference Room Paper 243, C Conference on Disarment, United Nations, Geneva, C Switzerland, July 1995. C Reference to GSE 2.1: C GSETT-3 (1997): Provisional GSE 2.1 Message Formats and Protocols. C Operation Annex 3. May 1997. C C....................................................................... C C This FORTRAN77 file consists of the following external procedures: C WGSE1...Subroutine designed to write the header section. It C should be called once after opening the output GSE file. C WGSE1 C WGSE2...Subroutine designed to write one seismogram. C WGSE2 C WGSE2C..Entry of subroutine WGSE2 designed to accumulate C individual lines to be written to the Comment sections of C the next seismograms written by WGSE2. Entry WGSE2C need C not be considered unless using the Comment sections of the C Waveform identification sections for special purposes. C WGSE2C C WGSE2D..Entry of subroutine WGSE2 designed to delete the lines C accumulated by previous invocations of WGSE2C. Usually C not needed. C WGSE2D C WGSE3...Subroutine designed to write the end of data file C indicator. It should be called once before closing the C output GSE file. C WGSE3 C RGSE1...Subroutine designed to read the text from the header C section. It may be called once after opening the input C GSE file. C RGSE1 C RGSE2...Subroutine designed to read one seismogram. The next C seismogram will be read during the next invocation of this C subroutine. C RGSE2 C RGSE2C..Entry of subroutine RGSE2 designed to pass individual C lines of the Comment section of the last seismogram read C by RGSE2. Entry RGSE2C need not be considered unless C using the Comment sections of the Waveform identification C sections for special purposes. C RGSE2C C C....................................................................... C C C Structure of the GSE 1.0 waveform data exchange file: C (1) Header section: C Position: Format: Item: C 01-04 A4 'XW01'... Obligatory string, an abbreviation for C the waveform data exchange file. C 05-45 41X Reserved for the information on the data C transmission, see the GSE 1.0 standard (blank C for synthetic data). Left blank during writing, C skipped during reading. C 46-80 1X,A34 Text describing the file or the data. Left C justified. The GSE 1.0 standard requires the C first 6 characters of the text (positions 47-52) C to identify the data day in the form: YYMMDD, C see parameter text of subroutine WGSE1. C (2) For each seismogram (2.1), (2.2), (2.3), (2.4), (2.5), and (2.6): C (2.1) Waveform identification section, Line 1: C Position: Format: Item: C 01-04 A4 'WID1'... Obligatory string, an abbreviation for C the waveform identification section. C 05-10 1X,I5 Year (often zero for synthetic data). Left C blank during writing, skipped during reading. C 11-13 I3 Day of the year (often zero for synthetic data). C If start time (argument TSTART of subroutine C WGSE2) is positive, WGSE2 writes day 001. If C start time TSTART is negative, WGSE2 writes day C 000 and time TSTART increased by 86400 seconds C before conversion to hours, minutes and seconds. C Subroutine RGSE2 deems a positive day or a blank C day to be 001, and a non-positive day to be 000. C 14-16 1X,I2 Start hours (often zero for synthetic data). C 17-19 1X,I2 Start minutes (often zero for synthetic data). C 20-22 1X,I2 Start seconds. C 23-26 1X,I3 Start miliseconds. C 27 I1 Start tenths of miliseconds (in the GSE 1.0 C standard, start time is rounded to miliseconds C and this position should be left free). C Written by subroutine WGSE2 if non-zero. C This position makes likely the greatest C difference between the GSE 1.0 standard and this C implementation, but need not necessarily C influence the portability of the data files. C 28-35 I8 Number of samples. C 36-42 1X,A6 Station name (sometimes used for synthetic C data, sometimes blank or '-'). Left justified. C 43-51 1X,A8 Channel name (occasionally used for synthetic C data, often blank or '-' even for real data). C Left justified. C 52-54 1X,A2 Component, recommendation for synthetic data: C '- ' OR ' 0'... Undefined, C '-E' OR ' 1'... X1 component, C '-N' OR ' 2'... X2 component, C '-Z' OR ' 3'... X3 component, etc. C The component (channel identifier) is written to C the GSE file by subroutine WGSE2 in the form of C integer (format I2), which is probably unusual C for the field data. However, GSE 1.0 standard C admits an arbitrary string composed of two C uppercase characters. C 55-66 1X,F11.7 Samples per second. C 67-73 1X,A6 Instrument type ('-' for synthetic data). C '- ' written by subroutine WGSE2, skipped C during reading. C 74-78 1X,A4 Data format, recommendation for synthetic data: C 'INT5'... Format(16(1X,I4), C 'INTV'... Integers, list directed input (free C format), linewidth 80 characters. C 79-80 1X,I1 Differencing flag, recommendation for synthetic C data: C '0' or ' '... Actual seismogram stored, no C differences. Differences are not C assumed by the subroutines of this C file. C (2.2) Waveform identification section, Line 2: C Position: Format: Item: C 01-09 F9.6 Amplitude corresponding to 1 in the seismogram. C 10 A1 Units of the amplitude: C '0' or ' '... Displacement in nanometers, C '1'... Particle velocity in nanometers/second, C '2'... Acceleration in nanometers/second**2. C Left blank by subroutine WGSE2, not read by C subroutine RGSE2. C 11-17 F7.4 Calibration period (-1 for synthetic data). C -1. written by subroutine WGSE2, not read by C subroutine RGSE2. C 18-27 1X,F9.4 X2 coordinate (only geographical coordinates C allowed by the GSE 1.0 standard for field data). C 28-37 1X,F9.4 X1 coordinate (only geographical coordinates C allowed by the GSE 1.0 standard for field data). C 38-47 1X,F9.4 X3 coordinate (only elevation above sea level in C metres allowed by the GSE standard for field C data). C 48-57 1X,F9.4 Depth of sensor (often -999 (N/A) for synthetic C data). C -999. written by subroutine WGSE2, not read by C subroutine RGSE2. C 58-65 1X,F7.2 Beam azimuth (-1 (N/A) for synthetic data). C -1. written by subroutine WGSE2, not read by C subroutine RGSE2. C 66-73 1X,F7.2 Beam slowness (-1 (N/A) for synthetic data). C -1. written by subroutine WGSE2, not read by C subroutine RGSE2. C 74-80 1X,F6.1 Orientation of horizontal sensors, measured C clockwise from the X2-axis (mostly -1 (N/A) for C synthetic data). C -1. written by subroutine WGSE2, not read by C subroutine RGSE2. C (2.3) Waveform identification section - Comments: C None to several lines with character strings, FORMAT(A80), C containing application specific information. C First four characters of a line must not match any of the C following strings: 'X***', 'WID*', 'DAT*', 'STOP', where * stands C for a wild character. To use columns 3 to 80 only seems to be the C best way of contingently introducing this specific information. C The comment lines are accumulated by means of invocation of entry C WGSE2C of subroutine WGSE2. C Note: Subroutine WGSE2 writes the start time to the Comments C as a SEP parameter OT. Subroutine RGSE2 reads the start time C from the Comments, and substitutes it for the GSE start time C from the Line 1 if the GSE start time is lower than 10 seconds C in order to reduce the GSE start time rounding errors. C Example: C If we like to write real-valued data into the comment lines in C the SEP format, C subroutine WSEPR C may be called to generate string containing 2 spaces followed C by the NAME=RVAL couple. C Here NAME is the name of the parameter and RVAL is the value. C For instance, string ' X1SRC=2.5' indicates that parameter C X1SRC, used by some programs for the X1 source coordinate, has C the value of 2.5. C The string may then be stored in the memory by the invocation of C entry WGSE2C. All strings stored by invocations of WGSE2C are C written to individual comment lines of the corresponding C seismogram by each invocation of subroutine WGSE2. C After reading a seismogram by the invocation of subroutine C RGSE2, corresponding individual comment lines may be retrieved C by the respective invocations of entry RGSE2C. The comment C lines may be collected in the memory by subroutine C RSEP2 and the C values of the real-valued parameters may be obtained by the C invocations of subroutine C RSEP3R. C (2.4) Waveform data section - Beginning: C Position: Format: Item: C 01-04 A4 'DAT1'... Obligatory string, an abbreviation for C the waveform data section. C (2.5) Waveform data section: C Seismogram formatted according to the specified data format, C assumed to be readable by list directed input into an integer or C real array. C (2.6) Waveform data section - End: C Position: Format: Item: C 01-04 A4 'CHK1'... Obligatory string, an abbreviation for C the checksum. C 05-09 5X C 10-24 I or F Sum of all data values (before compression or C differencing). C (3) End of data file indicator: C Position: Format: Item: C 01-04 A4 'STOP'... Obligatory string at the end. C Note: GSE file may contain many other sections or lines than described C above, but they are not written and read by the subroutines of this C file. C C======================================================================= C C C SUBROUTINE WGSE1(LU,TEXT) CHARACTER*(*) TEXT INTEGER LU C C Subroutine designed to write the header section. It should be called C once after opening the output GSE file. C C Input: C LU... Logical unit number of the external output device already C open for formatted sequential output. C TEXT... String containing the text to be stored in the header C section. Maximum of 34 characters will be stored. C Note that the GSE standard requires the first 6 characters C of the text to identify the data day in the form: YYMMDD. C Such a restriction is not required by the subroutines of C this file. C The input parameters are not altered. C C No output. C C----------------------------------------------------------------------- C C No temporary storage locations. C C....................................................................... C C Writing header section: WRITE(LU,'(A4,42X,A)') 'XW01',TEXT(1:MIN0(LEN(TEXT),34)) RETURN END C C======================================================================= C C C SUBROUTINE WGSE2(LU,STAT,CHAN,KOMP, * X1,X2,X3,TSTART,TSTEP,NSEIS,SEIS) INTEGER LU,KOMP,NSEIS CHARACTER*(*) STAT,CHAN REAL X1,X2,X3,TSTART,TSTEP,SEIS(NSEIS) C C Subroutine designed to write one seismogram. C C Input: C LU... Logical unit number of the external output device already C open for formatted sequential output. C STAT... String with the station name (not required, may be ' ' or C '-'). C CHAN... String with the channel name (often ' ' or '-'). C KOMP... Component (nonnegative integer not exceeding 99): C 0... Undefined component (e.g., scalar time function), C 1... X1 component, C 2... X2 component, C 3... X3 component, etc. C The component (channel identifier) is written to the GSE C file in the form of integer (format I2), which is probably C unusual for the field data. However, GSE standard admits C an arbitrary string composed of 2 uppercase characters. C For example, the components of the 3*3 Green function may C be numbered 1,2,3,4,5,6,7,8,9. In the case of complex C functions, the imaginary parts may be distinguished from C the corresponding real parts by KOMP increased by 10. C X1,X2,X3... Coordinates of the receiver. C -999.0000 stands for the undefined value, that may be C applicable especially for, e.g., X2 in 2-D synthetic C calculations. C TSTART..Start time, i.e. the time corresponding to the first C sample in seconds. In the GSE file, the time is converted C to hours, minutes, seconds, and days from the beginning of C an unspecified year. Times greater than 86400 seconds in C absolute value are not assumed. Formal day corresponding C to negative times in GSE 1.0 format is January 0. C TSTEP...Time step between samples. C NSEIS...Number of samples. C (SEIS(I),I=1,NSEIS)... Seismogram. The leading and trailing zeros C will be removed before writing. C The input parameters are not altered. C C No output. C C C Input SEP parameter: C GSEWIDTH=positive integer... Width of the output field reserved C for one integer value of the seismogram. C The number of output digits is then GSEWIDTH-2. C GSEWIDTH must be taken from interval 5 to 8. C Note that GSEWIDTH=9 might generate incorrect output C caused by rounding errors in single precision, and that C GSEWIDTH=3 or 4 would output too inaccurate seismograms. C Default: GSEWIDTH=5 C C....................................................................... C C C C ENTRY WGSE2C(COMLIN) CHARACTER*(*) COMLIN C C Entry designed to accumulate individual lines to be written to the C Comment sections of the next seismograms written by WGSE2. Entry C WGSE2C need not be considered unless using the Comment sections of the C Waveform identification sections for special purposes. C C Input: C COMLIN..Line to be written to the comment lines of the waveform C identification section. The lines passed by the first C MLINES invocations of WGSE2C are accumulated and written C by each next invocation of WGSE2, until the lines are C deleted by invocation of WGSE2D. Parameter MLINES may be C adjusted below. C Only the first 80 characters of COMLIN are considered. C The input parameter is not altered. C C No output. C C....................................................................... C C C C ENTRY WGSE2D C C Entry designed to delete the lines accumulated by previous invocations C of WGSE2C. The accumulation by invocations of WGSE2C then begins from C the scratch. WGSE2D is not required too often. C C No Input, no output. C C----------------------------------------------------------------------- C C Parameters: LOGICAL FIXED PARAMETER (FIXED=.FALSE.) C FIXED...FIXED=.TRUE.: Fixed-length format, C FIXED=.FALSE.: Variable-length format. C C Temporary storage locations: CHARACTER*6 FORMAT CHARACTER*4 FORM CHARACTER*80 LINE CHARACTER*6 STAT6 CHARACTER*8 CHAN8 INTEGER IDAY,IH,IM,IS,IMS,IMS4,I,L,NSEIS1,NSEIS2,ISEIS(0:25) REAL OT,TMS4,FSTEP,AMPL,AUX C C Storage location for Comments of Waveform identification section: INTEGER MLINES,NLINES PARAMETER (MLINES=20) CHARACTER*80 LINES(MLINES) SAVE NLINES,LINES C MLINES..Only the first MLINES comment lines are stored. C C Width of the output field reserved for one integer value: INTEGER NWIDTH,NFORM SAVE NWIDTH,NFORM C DATA NLINES/0/ DATA NWIDTH/0/ C C....................................................................... C C Reading input SEP parameter GSEWIDTH: IF(NWIDTH.EQ.0) THEN CALL RSEP3I('GSEWIDTH',NWIDTH,5) IF(NWIDTH.LT.5.OR.NWIDTH.GT.8) THEN C GSE-08 CALL ERROR('GSE-08: Incorrect SEP parameter GSEWIDTH') C SEP parameter GSEWIDTH must be taken from interval 5 to 8. END IF NFORM=80/NWIDTH END IF C C Output format: FORMAT='(16I5)' IF(NWIDTH.NE.5) THEN WRITE(FORMAT(5:5),'(I1)') NWIDTH END IF FORM='INTV' IF(FIXED) THEN FORM(4:4)=FORMAT(5:5) END IF C C Nonzero part of the seismogram and the maximum amplitude: DO 11 NSEIS1=1,NSEIS IF(SEIS(NSEIS1).NE.0.) THEN GO TO 12 END IF 11 CONTINUE 12 CONTINUE DO 13 NSEIS2=NSEIS,NSEIS1,-1 IF(SEIS(NSEIS2).NE.0.) THEN GO TO 14 END IF 13 CONTINUE 14 CONTINUE AMPL=0. DO 15 I=NSEIS1,NSEIS2 AMPL=AMAX1(ABS(SEIS(I)),AMPL) 15 CONTINUE AMPL=AMPL/FLOAT(10**(NWIDTH-2)-1) AUX=AMPL/2. DO 16 NSEIS1=NSEIS1,NSEIS2 IF(ABS(SEIS(NSEIS1)).GT.AUX) THEN GO TO 17 END IF 16 CONTINUE 17 CONTINUE DO 18 NSEIS2=NSEIS2,NSEIS1,-1 IF(ABS(SEIS(NSEIS2)).GT.AUX) THEN GO TO 19 END IF 18 CONTINUE 19 CONTINUE C C Writing waveform identification section: IF(ABS(TSTART).GT.86400.) THEN C GSE-01 CALL ERROR('GSE-01: Starting time greater than 86400') END IF OT=TSTART+FLOAT(NSEIS1-1)*TSTEP TMS4=OT*10000. IF(TMS4.LT.-.5) THEN IDAY=000 IMS4=INT(TMS4-.5)+864000000 ELSE IDAY=001 IMS4=INT(TMS4+.5) END IF IMS=IMS4/10 IMS4=IMS4-IMS*10 IS=IMS/1000 IMS=IMS-IS*1000 IM=IS/60 IS=IS-IM*60 IH=IM/60 IM=IM-IH*60 I=MAX0(NSEIS2-NSEIS1+1,0) STAT6=STAT CHAN8=CHAN FSTEP=1./TSTEP IF(IMS4.EQ.0) THEN IF(FSTEP.LT.999.99999995) THEN WRITE(LU,'(A4,6X,I3.3,3(1X,I2.2),I4.3,1X, * I8,1X,A6,1X,A8,1X,I2,1X,F11.7,1X,A1,6X,A4)') * 'WID1',IDAY, IH,IM,IS,IMS, * I,STAT6,CHAN8, KOMP, FSTEP, '-', FORM ELSE IF(FSTEP.LT.9999.9999995) THEN WRITE(LU,'(A4,6X,I3.3,3(1X,I2.2),I4.3,1X, * I8,1X,A6,1X,A8,1X,I2,1X,F11.6,1X,A1,6X,A4)') * 'WID1',IDAY, IH,IM,IS,IMS, * I,STAT6,CHAN8, KOMP, FSTEP, '-', FORM ELSE IF(FSTEP.LT.9999999.9995) THEN WRITE(LU,'(A4,6X,I3.3,3(1X,I2.2),I4.3,1X, * I8,1X,A6,1X,A8,1X,I2,1X,F11.3,1X,A1,6X,A4)') * 'WID1',IDAY, IH,IM,IS,IMS, * I,STAT6,CHAN8, KOMP, FSTEP, '-', FORM ELSE WRITE(LU,'(A4,6X,I3.3,3(1X,I2.2),I4.3,1X, * I8,1X,A6,1X,A8,1X,I2,1X,F11.0,1X,A1,6X,A4)') * 'WID1',IDAY, IH,IM,IS,IMS, * I,STAT6,CHAN8, KOMP, FSTEP, '-', FORM END IF ELSE IF(FSTEP.LT.999.99999995) THEN WRITE(LU,'(A4,6X,I3.3,3(1X,I2.2),I4.3,I1, * I8,1X,A6,1X,A8,1X,I2,1X,F11.7,1X,A1,6X,A4)') * 'WID1',IDAY, IH,IM,IS,IMS,IMS4, * I,STAT6,CHAN8, KOMP, FSTEP, '-', FORM ELSE IF(FSTEP.LT.9999.9999995) THEN WRITE(LU,'(A4,6X,I3.3,3(1X,I2.2),I4.3,I1, * I8,1X,A6,1X,A8,1X,I2,1X,F11.6,1X,A1,6X,A4)') * 'WID1',IDAY, IH,IM,IS,IMS,IMS4, * I,STAT6,CHAN8, KOMP, FSTEP, '-', FORM ELSE IF(FSTEP.LT.9999999.9995) THEN WRITE(LU,'(A4,6X,I3.3,3(1X,I2.2),I4.3,I1, * I8,1X,A6,1X,A8,1X,I2,1X,F11.3,1X,A1,6X,A4)') * 'WID1',IDAY, IH,IM,IS,IMS,IMS4, * I,STAT6,CHAN8, KOMP, FSTEP, '-', FORM ELSE WRITE(LU,'(A4,6X,I3.3,3(1X,I2.2),I4.3,I1, * I8,1X,A6,1X,A8,1X,I2,1X,F11.0,1X,A1,6X,A4)') * 'WID1',IDAY, IH,IM,IS,IMS,IMS4, * I,STAT6,CHAN8, KOMP, FSTEP, '-', FORM END IF END IF IF(-999.9999.LE.X1.AND.X1.LE.9999.9999.AND. * -999.9999.LE.X2.AND.X2.LE.9999.9999.AND. * -999.9999.LE.X3.AND.X3.LE.9999.9999) THEN IF(0.0099999.LT.AMPL.AND.AMPL.LE.99.999999) THEN WRITE(LU, * '( F9.6, 1X, F7.0,3(1X, F9.4 ), F10.0,F8.0,F8.0,F7.0)') * AMPL, -1.0, X2,X1,X3, -999.,-1.0,-1.0,-1.0 ELSE WRITE(LU, * '(1PE9.3E2,1X,0PF7.0,3(1X, F9.4 ), F10.0,F8.0,F8.0,F7.0)') * AMPL, -1.0, X2,X1,X3, -999.,-1.0,-1.0,-1.0 END IF ELSE IF(0.0099999.LT.AMPL.AND.AMPL.LE.99.999999) THEN WRITE(LU, * '( F9.6, 1X, F7.0,3(1X,1PE9.3E1),0PF10.0,F8.0,F8.0,F7.0)') * AMPL, -1.0, X2,X1,X3, -999.,-1.0,-1.0,-1.0 ELSE WRITE(LU, * '(1PE9.3E2,1X,0PF7.0,3(1X,1PE9.3E1),0PF10.0,F8.0,F8.0,F7.0)') * AMPL, -1.0, X2,X1,X3, -999.,-1.0,-1.0,-1.0 END IF END IF C C Writing waveform identification section - Comments: CALL WSEPR(LINE,'OT',OT) WRITE(LU,'(A)') LINE(3:LENGTH(LINE)) DO 23 I=1,NLINES DO 21 L=LEN(LINES(I)),2,-1 IF(LINES(I)(L:L).NE.' ') THEN GO TO 22 END IF 21 CONTINUE L=1 22 CONTINUE WRITE(LU,'(A)') LINES(I)(1:L) 23 CONTINUE C C Writing waveform data section: WRITE(LU,'(A4)') 'DAT1' IH=0 DO 59 IS=NSEIS1,NSEIS2,NFORM IM=MIN0(NSEIS2-IS,NFORM-1) DO 51 I=0,IM AUX=SEIS(IS+I)/AMPL AUX=SIGN(ABS(AUX)+.5,AUX) ISEIS(I)=INT(AUX) IH=IH+ISEIS(I) 51 CONTINUE IF(FIXED) THEN WRITE(LU,FORMAT) (ISEIS(I),I=0,IM) ELSE WRITE(LINE,FORMAT) (ISEIS(I),I=0,IM) IM=0 DO 52 I=2,NWIDTH*NFORM IF(LINE(I-1:I).NE.' ') THEN IM=IM+1 LINE(IM:IM)=LINE(I:I) END IF 52 CONTINUE WRITE(LU,'(A)') LINE(1:IM) END IF 59 CONTINUE WRITE(LU,'(A4,5X,I15)') 'CHK1',IH C RETURN C C----------------------------------------------------------------------- C ENTRY WGSE2C(COMLIN) C C....................................................................... C IF(NLINES.LT.MLINES) THEN NLINES=NLINES+1 LINES(NLINES)=COMLIN ELSE C GSE-02 CALL ERROR('GSE-02: Small parameter MLINES in subroutine WGSE2') END IF RETURN C C----------------------------------------------------------------------- C ENTRY WGSE2D C C....................................................................... C NLINES=0 RETURN END C C======================================================================= C C C SUBROUTINE WGSE3(LU) INTEGER LU C C Subroutine designed to write the end of data file indicator. It C should be called once before closing the output GSE file. C C Input: C LU... Logical unit number of the external output device already C open for formatted sequential output. The logical unit C is not closed in this subroutine. C The input parameter is not altered. C C No output. C C----------------------------------------------------------------------- C C No temporary storage locations. C C....................................................................... C C Writing end of data file indicator: WRITE(LU,'(A4)') 'STOP' RETURN END C C======================================================================= C C C SUBROUTINE RGSE1(LU,TEXT) CHARACTER*(*) TEXT INTEGER LU C C Subroutine designed to read the text from the header section. It may C be called once after opening the input GSE file. C C Input: C LU... Logical unit number of the external input device C containing the input data, already open for formatted C sequential input. C The input parameter is not altered. C C Output: C TEXT... String containing the text to read from the header C section. Maximum of 34 characters are stored in the file. C C----------------------------------------------------------------------- C C Temporary storage locations: CHARACTER*4 XW01 C C....................................................................... C C Searching for header section: 10 CONTINUE READ(LU,'(A4,42X,A)',END=80) XW01,TEXT IF(XW01.EQ.'XW01') THEN RETURN END IF GO TO 10 C C End of file: 80 CONTINUE C GSE-51 CALL WARN('GSE-51: No header section in input GSE file') RETURN END C C======================================================================= C C C SUBROUTINE RGSE2(LU,STAT,CHAN,KOMP, * X1,X2,X3,TSTART,TSTEP,NSEIS,MSEIS,SEIS) INTEGER LU,KOMP,NSEIS,MSEIS CHARACTER*(*) STAT,CHAN REAL X1,X2,X3,TSTART,TSTEP,SEIS(MSEIS) C C Subroutine designed to read one seismogram. The next seismogram will C be read during the next invocation of this subroutine. C C Input: C LU... Logical unit number of the external input device C containing the input data, already open for formatted C sequential input. C MSEIS...Dimension of array SEIS, i.e. the maximum number of C samples. C The input parameters are not altered. C C Output: C STAT... String with the station name. 6 characters at the most, C but may be declared even as CHARACTER*1 if not required. C See also subroutine WGSE2. C CHAN... String with the channel name. 8 characters at the most, C but may be declared even as CHARACTER*1 if not required. C See also subroutine WGSE2. C KOMP... Component (nonnegative integer): C 0... Undefined or unrecognized component, C 1... X1 component, C 2... X2 component, C 3... X3 component, C or another positive integer written by subroutine WGSE2. C In addition to the nonnegative integers written by WGSE2, C channel identifiers in the form of strings containing C characters 'E','N','Z' are recognized as components 1,2,3, C respectively. C X1,X2,X3... Coordinates of the receiver. C -999.0000 means an undefined value, that may apply C especially to one of the coordinates in 2-D synthetic C calculations. C TSTART..Start time, i.e. the time corresponding to the first C sample. Days are almost ignored: the time is measured C since the last midnight if the day of the year in the GSE C file is positive or blank, and from the next midnight C (i.e. TSTART is negative) if the day of the year is not C positive (e.g. January 0). C TSTEP...Time step between samples. C NSEIS...Number of output samples, or -1 if there is no more C seismogram to be read. C (SEIS(I),I=1,NSEIS)... Seismogram without the leading and trailing C zeros. C C....................................................................... C C C C ENTRY RGSE2C(COMLIN,*) CHARACTER*(*) COMLIN C C Entry designed to pass individual lines of the Comment section of the C last seismogram read by RGSE2. Entry RGSE2C need not be considered C unless using the Comment sections of the Waveform identification C sections for special purposes. C C Input: C *... Label of the statement in the calling subroutine for the C alternate return when no comment line remains for the C output. The label should be preceded by *. For example, C if the comment lines are to be processed by subroutine C RSEP2, the calling code may read C 10 CONTINUE C CALL RGSE2C(COMLIN,*20) C CALL RSEP2(COMLIN) C GO TO 10 C 20 CONTINUE C C Output: C COMLIN..One line of the Comment section of the last seismogram C read by RGSE2. The first invocation of RGSE2 returns the C first comment line, the second invocation of RGSE2 returns C the second comment line and so on. At most the first C MLINES comment lines are recorded and may be output. C Parameter MLINES may be adjusted below. C C----------------------------------------------------------------------- C C Storage location for Comments of Waveform identification section: INTEGER MLINES,NLINES PARAMETER (MLINES=20) CHARACTER*80 LINES(MLINES) SAVE NLINES,LINES C MLINES..Only the first MLINES comment lines are stored. C C Temporary storage locations: CHARACTER*80 LINE CHARACTER*1 COMP1,COMP2 CHARACTER*4 FORM INTEGER IDAY,IH,IM,IS,IMS,IMS4,I REAL OT,FSTEP,AMPL C C....................................................................... C C Searching for waveform identification section: 10 CONTINUE READ(LU,'(A)',END=80) LINE IF(LINE(1:4).EQ.'WID1') THEN GO TO 20 ELSE IF(LINE(1:4).EQ.'STOP') THEN GO TO 90 END IF GO TO 10 C C Reading waveform identification section: 20 CONTINUE READ(LINE,'(13X, * 3(1X,I2),1X,I3,I1,I8,1X,A6,1X,A8,1X,2A1,1X,F11.7,8X,A4,1X,I1)') * IH,IM,IS,IMS,IMS4,NSEIS,STAT, CHAN,COMP1,COMP2,FSTEP,FORM,I IF(LINE(11:13).NE.' ') THEN READ(LINE(11:13),'(I3)') IDAY IF(IDAY.LE.0) THEN IH=IH-24 END IF END IF TSTART=FLOAT((IH*60+IM)*60+IS)+FLOAT(10*IMS+IMS4)/10000. IF(NSEIS.GT.MSEIS) THEN C GSE-03 CALL ERROR('GSE-03: Seismogram longer than array') END IF IF(LLE('0',COMP1).AND.LLE(COMP1,'9')) THEN KOMP=ICHAR(COMP1)-ICHAR('0') ELSE KOMP=0 END IF IF(LLE('0',COMP2).AND.LLE(COMP2,'9')) THEN KOMP=ICHAR(COMP2)-ICHAR('0')+10*KOMP END IF IF(KOMP.EQ.0) THEN IF(COMP1.EQ.'Z'.OR.COMP2.EQ.'Z') THEN KOMP=3 ELSE IF(COMP1.EQ.'N'.OR.COMP2.EQ.'N') THEN KOMP=2 ELSE IF(COMP1.EQ.'E'.OR.COMP2.EQ.'E') THEN KOMP=1 END IF END IF TSTEP=1./FSTEP IF(FORM(1:3).NE.'INT') THEN C GSE-04 CALL ERROR('GSE-04: Data format not supported') END IF IF(I.NE.0) THEN C GSE-05 CALL ERROR('GSE-05: Data differences not supported') END IF READ(LU,'(F9.6,8X,3(1X,F9.4))') AMPL,X2,X1,X3 C C Searching for waveform data section: NLINES=0 30 CONTINUE READ(LU,'(A)') LINE IF(LINE(1:4).EQ.'DAT1') THEN GO TO 50 ELSE IF (LINE(1:3).EQ.'OT=') THEN IF(ABS(TSTART).LT.10.) THEN READ(LINE,'(3X,F77.0)') OT IF(ABS(TSTART-OT).LT.0.00006) THEN TSTART=OT ELSE C GSE-06 CALL ERROR('GSE-06: Ambiguous start time in the GSE file') END IF END IF ELSE C Reading waveform identification section - Comments IF(NLINES.LT.MLINES) THEN NLINES=NLINES+1 LINES(NLINES)=LINE ELSE C GSE-07 CALL ERROR('GSE-07: Small parameter MLINES in subr. RGSE2') END IF END IF GO TO 30 C C Reading waveform data section: 50 CONTINUE IF(NSEIS.GT.0) THEN READ(LU,*) (SEIS(I),I=1,NSEIS) END IF C Checksum is not read here. C C Removing leading and trailing zeros, and renormalizing: DO 51 NSEIS=NSEIS,1,-1 IF(SEIS(NSEIS).NE.0.) THEN GO TO 52 END IF 51 CONTINUE 52 CONTINUE DO 53 IS=1,NSEIS IF(SEIS(IS).NE.0.) THEN GO TO 54 END IF 53 CONTINUE 54 CONTINUE IF(IS.GT.1) THEN IS=IS-1 TSTART=TSTART+FLOAT(IS)*TSTEP NSEIS=NSEIS-IS DO 55 I=1,NSEIS-IS SEIS(I)=SEIS(I+IS) 55 CONTINUE END IF DO 59 I=1,NSEIS SEIS(I)=SEIS(I)*AMPL 59 CONTINUE C RETURN C C End of file: 80 CONTINUE C GSE-52 CALL WARN('GSE-52: End of input GSE file encountered') C GSE file is not terminated by string STOP in the first 4 columns, C or NSEIS=-1 has not been checked during previous invocation and C subroutine RGSE2 was called again. 90 CONTINUE NSEIS=-1 RETURN C C----------------------------------------------------------------------- C ENTRY RGSE2C(COMLIN,*) C C....................................................................... C IF(NLINES.LE.0) THEN C No comment line left - return to the label given by argument *: RETURN 1 END IF C COMLIN=LINES(1) NLINES=NLINES-1 DO 99 I=1,NLINES LINES(I)=LINES(I+1) 99 CONTINUE RETURN END C C======================================================================= C