C 
C Subroutine file 'gse.for' to write and read seismograms in the GSE 1.0
C data exchange format.
C
C Version: 5.50
C Date: 2000, November 26
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
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.  Usually not used.
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     Example:
C       If we like to write real-valued data into the comment lines in
C       the SEP format, subroutine
C       WSEPR may be called to generate
C       string containing 2 spaces followed 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 values of the
C       real-valued parameters may be obtained by the invocations of
C       subroutine 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
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 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
      INTEGER NWIDTH,NFORM
      PARAMETER (FIXED=.FALSE.)
      PARAMETER (NWIDTH=5,NFORM=80/NWIDTH)
C
C     NWIDTH is the width of the output field reserved for one integer
C             value of the seismogram.  The number of output digits is
C             then NWIDTH-2.  NWIDTH must be taken from interval 3 to 9,
C             but NWIDTH=3 or 4 is strongly discouraged.
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 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
      DATA NLINES/0/
C     MLINES..Only the first MLINES comment lines are stored.
C
C.......................................................................
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
      TMS4=(TSTART+FLOAT(NSEIS1-1)*TSTEP)*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(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
C
C     Writing waveform identification section - Comments:
      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,80
            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
      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 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-02
        CALL ERROR('GSE-02: 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-03
        CALL ERROR('GSE-03: Data format not supported')
      END IF
      IF(I.NE.0) THEN
C       GSE-04
        CALL ERROR('GSE-04: 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
C         Reading waveform identification section - Comments
          IF(NLINES.LT.MLINES) THEN
            NLINES=NLINES+1
            LINES(NLINES)=LINE
          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