C 
C Program GSE2SEGY converts seismograms stored in the GSE data exchange
C format to SEGY format. Trace Data sample values are written in
C 32-bit IEEE floating-point format. This format was added by
C SEGY revision 1 (Digital tape standards (2004). Society of
C Exploration Geophysicists).
C This version of program GSE2SEGY writes only limited number
C of SEGY reel and trace header parameters.
C
C Part of code is taken from sp.for.
C
C Version: 5.90
C Date: 2005, June 15
C
C Assembled and coded by: Vaclav Bucha
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bucha@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C Attention:  Functionality of program GSE2SEGY is strongly affected by
C the Fortran compiler and by the options of the compiler.
C Program GSE2SEGY can work only if the compiler supports unformatted
C direct-access files "without headers".
C Program GSE2SEGY uses INTEGER*2 statement, which violates Fotran 77
C standard.
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Names of the input files:
C     SS='string'... String with the name of the input data file in the
C             GSE data exchange format, containing the seismograms to be
C             converted.
C             The source and receiver names cannot be longer than
C               6 characters.
C             Description of GSE file SS
C             Default: SS='ss.gse'
C     SRC='string'... String with the name of the input data file
C             with the name(s) and coordinates of the source(s).
C             If SRC is not blank and the source names are specified
C               comment lines of the GSE file, only seismograms
C               corresponding to the sources listed in file SRC will be
C               plotted, otherwise no selection according to the source
C               names is performed.
C             The source names cannot be longer than 6 characters.
C             File SRC usually need not
C             be specified if the seismograms are generated by programs
C             'greenss.for' or 'ss.for', because those programs enter
C             the source coordinates directly to the comments of
C             the data section in GSE data exchange file.
C             Description of file SRC
C             Default: SRC=' '
C     REC='string'... String with the name of the input data file
C             containing the list of receivers.
C             Only the first 6 characters of receiver names are
C               significant.
C             The receiver names cannot be longer than 6 characters.
C             In most cases, file REC will be the same as for the
C             calculation of synthetic seismograms.
C             Description of file REC
C             Default: REC=' '
C Names of output files:
C     FISEGY='string'... String with the name of the output SEGY file.
C             Default: FISEGY='ss.sgy'
C Component selection:
C     KOMP=integer... Selected component of the GSE seismograms
C             that will be converted to output SEGY file.
C             Default: KOMP=1
C SEGY parameters:
C     NTRACE=integer... Number of data traces per record (inludes dummy
C             and zero traces inserted to fill out the record
C             or common depth point).
C             Default: NTRACE=1
C     NSAMPL=integer... Number of samples per data trace (for this reel
C             of data).
C             Default: NSAMPL=1
C     SINTER=real... Sample interval in microseconds (for this reel
C             of data).
C             Default: SINTER=1.
C     ISFORM=integer... Data sample format code.
C             Default: ISFORM=1 (floating point (4 bytes))
C     ICDPF=integer... CDP fold (expected number of data traces per CDP
C             ensemble).
C             Default: ICDPF=1
C.......................................................................
C
C SEGY format specification according to
C Barry, K.M, Cavers, D.A. and Kneale, C.W. (1975):
C   Recommended standards for digital tape formats.
C   Geophysics, 40, no. 02, 344-352.
C (1) Reel identification header (3600 bytes)
C (1.1) Part 1, the EBCDIC card image block (3200 bytes - 40 cards)
C       The EBCDIC part of the reel header describes the data
C       from a line of shotpoints in a fixed specified format
C       consisting of 40 card images with each image containing
C       80 bytes. All unused card image characters are EBCDIC
C       Blank. Card image numbers 23 through 39 are unassigned
C       for optional use. Each card image should contain the
C       character C in the first card column. Each 80 bytes would
C       yield one line of format print.
C (1.2) Part 2, the binary coded block (400 bytes)
C       The binary coded section of the reel header consists
C       of 400 bytes of information common to the seismic data
C       on the related reel. There are 60 bytes assigned;
C       340 are unassigned for optional use.
C     BINARY CODE - Right Justified
C     Byte Numbers: Description:
C     3201-3204     Job identification number.
C     3205-3208  *  Line number (only one line per reel).
C     3209-3212  *  Reel number.
C     3213-3214  *  Number of data traces per record (inludes dummy
C                   and zero traces inserted to fill out the record
C                   or common depth point).
C     3215-3216  *  Number of auxiliary traces per record (includes
C                   sweep, timing, gain, sync and all other non-data
C                   traces).
C     3217-3218  *  Sample interval in microseconds (for this reel
C                   of data).
C     3219-3220     Sample interval in microseconds (for original
C                   field recording).
C     3221-3222  *  Number of samples per data trace (for this reel
C                   of data).
C     3223-3224     Number of samples per data trace (for original
C                   field recording).
C     3225-3226  *  Data sample format code:
C                   1=floating point (4 bytes)  3=fixed point (2 bytes)
C                   2=fixed point (4 bytes)     4=fixed point w/gain
C                                                 code (4 bytes)
C                   Auxiliary traces use the same number of bytes
C                   per sample.
C     3227-3228  *  CDP fold (expected number of data traces per CDP
C                   ensemble).
C
C     3229-3230     Trace sorting code:
C                   1=as recorded (no sorting)  3=single fold continuous
C                                                 profile
C                   2=CDP ensemble              4=horizontally stacked
C     3231-3232     Vertical sum code: 1=no sum, 2=two sum,...,
C                   N=N sum (N=32767)
C     3233-3234     Sweep frequency at start.
C     3235-3236     Sweep frequency at end.
C     3237-3238     Sweep length (ms).
C     3239-3240     Sweep type code: 1=linear    3=exponential
C                                    2=parabolic 4other
C     3241-3242     Trace number of sweep channel.
C     3243-3244     Sweep trace taper length in ms at start if tapered
C                   (the taper starts at zero time and is effective
C                   for this length).
C     3245-3246     Sweep trace taper length in ms at end (the ending
C                   taper starts at sweep length minus the taper
C                   length at end).
C     3247-3248     Taper type: 1=linear    3=other
C                               2=cos(exp2)
C     3249-3250     Correlated data traces: 1=no    2=yes
C     3251-3252     Binary gain recovered:  1=yes   2=no
C     3253-3254     Amplitude recovery method:
C                   1=none                   3=AGC
C                   2=spherical divergence   4=other
C     3255-3256  *  Measurement system: 1=meters  2=feet
C     3257-3258     Impulse signal polarity:
C                   1=Increase in pressure or upward geophone case
C                     movement gives negative number on tape.
C                   2=Increase in pressure or upward geophone case
C                     movement gives positive number on tape.
C     3259-3260     Vibratory polarity code.
C     3261-3600     Unassigned - for optional information.
C
C * Strongly recommended that this information always be recorded.
C
C (2) Trace data block
C (2.1) Trace identification header (240 bytes)
C       Trace header is written in binary code
C     Byte numbers: Description:
C        1-4     *  Trace sequence number within line - numbers
C                   continue to increase if additional reels are
C                   required on same line.
C        5-8        Trace sequence number within reel - each reel
C                   starts with trace number one.
C        9-12    *  Original field record number.
C       13-16    *  Trace number within the original field record.
C       17-20       Energy source point number - used when more than
C                   one record occurs at the same effective surface
C                   location.
C       21-24       CDP ensemble number.
C       25-28       Trace number within the CDP ensemble - each
C                   ensemble starts with trace number one.
C       29-30    *  Trace identification code:
C                   1=seismic data 4=time break 7=timing
C                   2=dead         5=uphole     8=water break
C                   3=dummy        6=sweep      9..N=optional use
C       31-32       Number of vertically summed traces yielding this
C                   trace. (1 is one trace, 2 is two stacked
C                   traces, etc.)
C       33-34       Number of horizontally summed traces yielding this
C                   trace. (1 is one trace, 2 is two stacked
C                   traces, etc.)
C       35-36       Data use: 1=production. 2=test.
C       37-40       Distance from source point to receiver group
C                   (negative if opposite to direction in which line
C                   is shot).
C       41-44       Receiver group elevation; all elevations above sea
C                   level are positive and below sea level are negative.
C       45-48       Surface elevation at source.
C       49-52       Source depth below surface (a positive number).
C       53-56       Datum elevation at receiver group.
C       57-60       Datum elevation at source.
C       61-64       Water depth at source.
C       65-68       Water depth at group.
C       69-70       Scaler to be applied to all elevations and depths
C                   specified in bytes 41-68 to give the real value.
C                   Scaler=1,+/-10,+/-100,+/-1000, or +/-10000.
C                   If positive, scaler is used as a multiplier;
C                   if negative, scaler is used as a divisor.
C       71-72       Scaler to be applied to all coordinates
C                   specified in bytes 73-88 to give the real value.
C                   Scaler=1,+/-10,+/-100,+/-1000, or +/-10000.
C                   If positive, scaler is used as a multiplier;
C                   if negative, scaler is used as a divisor.
C       73-76       Source coordinate -X. | If the coordinate units are
C                                         | in seconds of arc, the X
C                                         | values represent longitude
C       77-80       Source coordinate -Y. | and the Y values latitude.
C                                         | A positive value designates
C                                         | the number of seconds east
C       81-84       Group coordinate -X.  | of Greenwich Meridian or
C                                         | north of the equator and
C                                         | a negative value designates
C       85-88       Group coordinate -Y.  | the number of seconds south
C                                         | or west.
C       89-90       Coordinate units:1=length (meters or feet).
C                   2=seconds of arc.
C       91-92       Weathering velocity.
C       93-94       Subweathering velocity.
C       95-96       Uphole time at source.
C       97-98       Uphole time at group.
C       99-100      Source static correction.
C      101-102      Group static correction.
C      103-104      Total static applied. (Zero if no static has been
C                   applied.)
C      105-106      Lag time A. Time in ms between end of 240-byte trace
C                   identification header and time break. Positive if
C                   time break occurs before end of header. Time break
C                   is defined as the initiation pulse which may be
C                   recorded on an auxiliary trace or as otherwise
C                   specified by the recording system.
C      107-108      Lag time B. Time in ms between time break and the
C                   initiation time of the energy source. May be
C                   positive or negative.
C      109-110      Delay recording time. Time in ms between initiation
C                   time of energy source and time when recording of
C                   data samples begins. (For deep water work if data
C                   recording does not start at zero time.)
C      111-112      Mute time--start.
C      113-114      Mute time--end.
C      115-116   *  Number of samples in this trace.
C      117-118   *  Sample interval in microseconds for this trace.
C      119-120      Gain type of field instruments: 1=fixed, 2=binary,
C                   3=floating point, 4...N=optional use.
C      121-122      Instrument gain constant.
C      123-124      Instrument early or initial gain (db).
C      125-126      Correlated: 1=no, 2=yes.
C      127-128      Sweep frequency at start.
C      129-130      Sweep frequency at end.
C      131-132      Sweep length in ms.
C      133-134      Sweep type: 1=linear, 2=parabolic, 3=exponential,
C                   4=other.
C      135-136      Sweep trace taper length at start in ms.
C      137-138      Sweep trace taper length at end in ms.
C      139-140      Taper type: 1=linear, 2=cos(exp2), 3=other.
C      141-142      Alias filter frequency, if used.
C      143-144      Alias filter slope.
C      145-146      Notch filter frequency, if used.
C      147-148      Notch filter slope.
C      149-150      Low cut frequency, if used.
C      151-152      High cut frequency, if used.
C      153-154      Low cut slope.
C      155-156      High cut slope.
C      157-158      Year data recorded.
C      159-160      Day of year.
C      161-162      Hour of day (24 hour clock).
C      163-164      Minute of hour.
C      165-166      Second of minute.
C      167-168      Time basis code: 1=local, 2=GMT, 3=other.
C      169-170      Trace weighting factor--defined as 2(exp-N) volts
C                   for the least significant bit. (N=0, 1, ... 32767.)
C      171-172      Geophone group number of roll switch position one.
C      173-174      Geophone group number of trace number one within
C                   original field record.
C      175-176      Geophone group number of last trace within original
C                   field record.
C      177-178      Gap size (total number of groups dropped).
C      179-180      Overtravel associated with taper at beginning or
C                   end of line: 1=down (or behind), 2=up (or ahead).
C      181-240      Unassigned--for optional information.
C
C * Strongly recommended that this information always be recorded.
C
C (2.2) Trace data samples
C       Trace data samples can be written in one of four data sample
C       formats:
C       32 bit floating point format - in which each data value of
C        a seismic channel is recorded in four successive bytes,
C        in IBM compatible floating point notation as defined in IBM
C        Form GA 22-6821.
C       32 bit fixed point - each data value of a seismic channel is
C        recorded in four successive bytes.
C       16 bit fixed point - each data value of a seismic channel is
C        recorded in two successive bytes.
C       32 bit fixed point format with gain values.
C       In all four data formats, the channel or trace data should
C       represent the absolute input voltage at the recording
C       instrument.
C
C.......................................................................
C                                                
C This Fortran 77 file consists of the following external procedures:
C     GSE2SEGY...   Main program to read and plot the seismograms.
C             GSE2SEGY
C
C Other external procedures required:
C     RGSE1,RGSE2... Subroutines of the Fortran 77 file 'gse.for'
C             (package MODEL), designed to read seismograms in the GSE
C             data exchange format.
C
C=======================================================================
C
C     
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (RAM,IRAM)
C
C     Allocation of working arrays:
      INTEGER MPTS,MSS
      PARAMETER (MPTS=3000,MSS=MRAM-3*MPTS)
      CHARACTER*6 PTS(MPTS)
      COMMON/PTSC/PTS
C
C-----------------------------------------------------------------------
C
C     External functions and subroutines:
      EXTERNAL LENGTH
      INTEGER  LENGTH
C
C     Input and output data filenames and logical unit numbers:
      INTEGER LU
      PARAMETER (LU=1)
      CHARACTER*80 FILSEP,FILSRC,FILREC,FISEGY
      CHARACTER*80 FILESS(0:9)
C
C     Parameters and small working arrays:
      INTEGER ISS,ISP
C
      CHARACTER*1  CHAN,TEXT
      CHARACTER*6  NAMSRC,NAMREC
      CHARACTER*80 NAMPTS
C
C     Line of optional SEP parameter file or comments in the GSE file
      CHARACTER*80 LINE
C
C     Lists of point coordinates, sources and receivers:
C     INTEGER NSRC,NREC
C
      DATA FILESS/10*' '/
C
C     SEGY
      INTEGER NUMS,LUSEGY
      PARAMETER (NUMS=9000,LUSEGY=3)
      CHARACTER*80 RIHE(40)
      INTEGER*2 RIHB(200),TIH2(120)
      INTEGER TIH(60)
      INTEGER NTRACE,NSAMPL,ICDPF,ISFORM
      REAL TDS(NUMS)
      REAL SINTER
      EQUIVALENCE (TIH,TIH2)
C
C.......................................................................
C
C     Reading name of SEP file with input data:
      FILSEP=' '
      WRITE(*,'(A)') '+GSE2SEGY: Enter input filename: '
      READ(*,*) FILSEP
      IF (FILSEP.EQ.' ') THEN
C       GSE2SEGY-01
        CALL ERROR('GSE2SEGY-01: SEP file not given')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
C
C     Reading all data from the SEP file into the memory:
      CALL RSEP1(LU,FILSEP)
      WRITE(*,'(A)') '+GSE2SEGY: Working...            '
C
C     Output SEGY filename, seismogram component:
      CALL RSEP3T('FISEGY',FISEGY,'ss.sgy')
      CALL RSEP3I('KOMP',KOMP,1)
C
C     SEGY reel identification header
C     Part 1 - the EBCDIC card image block (3200 bytes - 40 cards)
      RIHE(1)(1:42)='SEGY file written by program gse2segy.for'
C     Part 2 - the binary coded block (400 bytes)
      CALL RSEP3I('NTRACE',NTRACE,1)
      CALL RSEP3I('NSAMPL',NSAMPL,1)
      CALL RSEP3R('SINTER',SINTER,1.)
      CALL RSEP3I('ISFORM',ISFORM,1)
      CALL RSEP3I('ICDPF',ICDPF,1)
      IF(NSAMPL.GT.NUMS) THEN
      WRITE(*,*)'NSAMPL=',NSAMPL,'  NUMS=',NUMS
C       GSE2SEGY-02
        CALL ERROR
     *  ('GSE2SEGY-02: NSAMPL is greater then NUMS. Increase the value
     *of NUMS in file gse2segy.for.')
      END IF
      RIHB(7)=NTRACE
      RIHB(9)=SINTER*1000.
      RIHB(10)=SINTER*1000.
      RIHB(11)=NSAMPL
      RIHB(12)=NSAMPL
      RIHB(13)=ISFORM
      RIHB(14)=ICDPF
      OPEN(LUSEGY,FILE=FISEGY,FORM='UNFORMATTED',ACCESS='DIRECT',
     &     RECL=1)
C     Write EBCDIC card image block (3200 bytes - 40 cards)
      DO 1 I=1,40
        DO 2 J=1,80
          WRITE(LUSEGY,REC=(I-1)*80+J) RIHE(I)(J:J)
    2   CONTINUE
    1 CONTINUE
      CLOSE(LUSEGY)
C     Write binary coded block (400 bytes)
      OPEN(LUSEGY,FILE=FISEGY,FORM='UNFORMATTED',ACCESS='DIRECT',
     &     RECL=2,STATUS='OLD')
      DO 3 J=1,200
        WRITE(LUSEGY,REC=1600+J) RIHB(J)
    3 CONTINUE
      CLOSE(LUSEGY)
C
C     Input and output filenames:
      CALL RSEP3T('SS' ,FILESS(0),'ss.gse')
      ISS0=0
C^^^
C     Reading lists of sources and receivers:
C     Sources
      NSRC=0
      CALL RSEP3T('SRC',FILSRC,' ')
      IF(FILSRC.NE.' ') THEN
        OPEN(LU,FILE=FILSRC,STATUS='OLD')
        READ(LU,*) (TEXT,I=1,20)
        DO 13 I=1,MPTS
          I0=MSS+3*I
          NAMPTS='$$$$$$$'
          RAM(I0-2)=0.
          RAM(I0-1)=0.
          RAM(I0  )=0.
          READ(LU,*,END=14) NAMPTS,RAM(I0-2),RAM(I0-1),RAM(I0)
          IF(NAMPTS(1:7).EQ.'$$$$$$$') THEN
            GO TO 14
          END IF
          IF(NAMPTS(7:).NE.' ') THEN
C           GSE2SEGY-05
            CALL ERROR('GSE2SEGY-05: Source name exceeds 6 characters')
C           Names of sources in file given by input parameter SRC cannot
C           be longer than 6 characters.  This is the limitation imposed
C           by the GSE standard on the receiver names and applied here
C           also to the source names.
          END IF
          PTS(I)=NAMPTS(1:6)
   13   CONTINUE
C         GSE2SEGY-03
          CALL ERROR
     *    ('GSE2SEGY-03: Array dimension MPTS small for sources')
   14   CONTINUE
        NSRC=I-1
        CLOSE(LU)
      END IF
C     Receivers
      NREC=0
      RECNUM=999.
      CALL RSEP3T('REC',FILREC,' ')
      IF(FILREC.NE.' ') THEN
        OPEN(LU,FILE=FILREC,STATUS='OLD')
        READ(LU,*) (TEXT,I=1,20)
        DO 15 I=NSRC+1,MPTS
          I0=MSS+3*I
          NAMPTS='$$$$$$$'
          RAM(I0-2)=0.
          RAM(I0-1)=0.
          RAM(I0  )=0.
          READ(LU,*,END=16) NAMPTS,RAM(I0-2),RAM(I0-1),RAM(I0)
          IF(NAMPTS(1:7).EQ.'$$$$$$$') THEN
            GO TO 16
          END IF
          IF(NAMPTS(7:).NE.' ') THEN
C           GSE2SEGY-06
            CALL ERROR
     *      ('GSE2SEGY-06: Receiver name exceeds 6 characters')
C           Names of receivers in file given by input parameter REC
C           cannot be longer than 6 characters.  This limitation is
C           imposed by the GSE standard.
          END IF
          PTS(I)=NAMPTS(1:6)
   15   CONTINUE
C         GSE2SEGY-04
          CALL ERROR
     *    ('GSE2SEGY-04: Array dimension MPTS small for points')
   16   CONTINUE
        NREC=I-NSRC-1
        RECNUM=FLOAT(NREC)
        CLOSE(LU)
      END IF
C
C     Loop for GSE files
      ISS=0
      IF(FILESS(ISS).NE.' ') THEN
C     Opening input GSE file with the seismograms:
        OPEN(LU,FILE=FILESS(ISS),STATUS='OLD')
        CALL RGSE1(LU,TEXT)
C^^^
C       Loop for seismograms
        ITR=1
        IREC=0
   40   CONTINUE
   41   CONTINUE
C###
C       Selecting the component:
   42   CONTINUE
        ISS1=ISS0
        CALL RGSE2
     *       (LU,NAMREC,CHAN,I,X1R,X2R,X3R,T0,TD,NSS,MSS,RAM)
        IF(NSS.LE.-1) THEN
C       End of the GSE file
          GO TO 80
        END IF
        IF(I.NE.KOMP) GO TO 42
C^^^
C       Selecting the receiver:
        IF(FILREC.EQ.' ') THEN
          IREC=IREC+1
        ELSE
C         Loop for receivers
          DO 51 I=NSRC+1,NSRC+NREC
            IF(PTS(I).EQ.NAMREC) THEN
              I0=MSS+3*I
              X1R=RAM(I0-2)
              X2R=RAM(I0-1)
              X3R=RAM(I0)
              IREC=I-NSRC
              GO TO 52
            END IF
   51     CONTINUE
          GO TO 41
        END IF
   52   CONTINUE
C###
C       Reading the source information:
   53   CONTINUE
          CALL RGSE2C(LINE,*54)
          CALL RSEP2(LINE)
        GO TO 53
   54   CONTINUE
        CALL RSEP3T('NAMESRC',NAMSRC,' ')
        CALL RSEP3R('X1SRC',X1S,0.)
        CALL RSEP3R('X2SRC',X2S,0.)
        CALL RSEP3R('X3SRC',X3S,0.)
C^^^
C       Selecting the source:
        IF(FILSRC.NE.' ') THEN
          IF(NAMSRC.EQ.' ') THEN
            I0=MSS
            X1S=RAM(I0+1)
            X2S=RAM(I0+2)
            X3S=RAM(I0+3)
          ELSE
C           Loop for sources
            DO 55 I=1,NSRC
              IF(PTS(I).EQ.NAMSRC) THEN
                I0=MSS+3*I
                X1S=RAM(I0-2)
                X2S=RAM(I0-1)
                X3S=RAM(I0)
                GO TO 56
              END IF
   55       CONTINUE
            GO TO 41
          END IF
        END IF
   56   CONTINUE
C
C       Determine first sample of seismogram
        IF(T0.GT.0.) THEN
          IZERO=T0/TD
C         Add leading zeros
          DO 70 I=1,IZERO
            TDS(I)=0.
   70     CONTINUE
        ELSE
          IZERO=T0/TD
        END IF
C       Add seismogram
        IF(IZERO+NSS.GT.NSAMPL) THEN
          IAUX=IZERO+NSS
          WRITE(*,*)'GSE samples=',IAUX,'  NSAMPL=',NSAMPL
C         GSE2SEGY-07
          CALL ERROR
     *    ('GSE2SEGY-07: Number of GSE samples is greater then NSAMPL.
     *Increase the value of NSAMPL in history file.')
        END IF
        IF(T0.GT.0.) THEN
          DO 72 I=1,NSS
            TDS(IZERO+I)=RAM(I)
   72     CONTINUE
C         Add trailing zeros
          DO 74 I=IZERO+NSS+1,NSAMPL
            TDS(I)=0.
   74     CONTINUE
        ELSE
C         Add seismogram
          DO 76 I=1,NSS
            TDS(I)=RAM(I-IZERO)
   76     CONTINUE
C         Add trailing zeros
          DO 78 I=NSS+1,NSAMPL
            TDS(I)=0.
   78     CONTINUE
        END IF
C
C       SEGY trace data block
        OPEN(LUSEGY,FILE=FISEGY,FORM='UNFORMATTED',
     &       ACCESS='DIRECT',RECL=4,STATUS='OLD')
C       Trace identification header (240 bytes)
        TIH(3)=1
        TIH(4)=ITR
        TIH2(58)=NSAMPL
        TIH2(59)=SINTER*1000.
        DO 82 K=1,60
          WRITE(LUSEGY,REC=900+(ITR-1)*(NSAMPL+60)+K) TIH(K)
   82   CONTINUE
C       Seismic trace data
        DO 84 J=1,NSAMPL
          WRITE(LUSEGY,REC=960+(ITR-1)*(NSAMPL+60)+J) TDS(J)
   84   CONTINUE
        CLOSE(LUSEGY)
C
        ITR=ITR+1
        GO TO 40
   80   CONTINUE
C       End of loop for receivers
C###
C       Closing GSE file
        CLOSE(LU)
      END IF
C     End of loop for GSE files
C
      WRITE(*,'(A)') '+GSE2SEGY: Done.                                 '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'gse.for'
C     gse.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C