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: 6.60 C Date: 2012, June 7 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 http://sw3d.cz/staff/bucha.htm 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 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 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 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 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 Indices of the sets of SEP parameters INTEGER ISEP,IOLD 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 Defining index ISEP of the working set of SEP parameters: ISEP=0 CALL SSEP(ISEP,IOLD) CALL SSEP(IOLD,ISEP) 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: ISEP=-ISEP CALL SSEP(ISEP,IOLD) 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.) CALL SSEP(IOLD,ISEP) 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