C
C Program GSECAL calculates a linear combination of the given C seismograms stored in the GSE data exchange format. C C Version: 6.60 C Date: 2012, February 13 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....................................................................... 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 GSE files and corresponding coefficients: C Here, # represents the value of integer constant MFILSS specified C in the this program. C E.g., if MFILSS=12, SS# stands for SS12, COEF# stands for C COEF12, etc. C SS1='string', SS2='string', ..., SS#='string'... Strings with the C names of the input data files in the GSE data exchange C format, containing the given seismograms. C For the GSE data exchange format, refer to the description C in file 'gse.for'. C Just nonblank filenames are considered. C Defaults: SS1=' ', SS2=' ', SS3=' ', ..., SS#=' ' C COEF1=real, COEF2=real, COEF3=real, ..., COEF#=real: C Coefficients of the linear combination to multiply C seismograms from the corresponding input GSE files. C Defaults: COEF1=1., COEF2=1., COEF3=1., ..., COEF#=1. C Name of the output GSE file: C SS='string'... String with the name of the output file in the GSE C data exchange format, containing the specified linear C combination of the input seismograms. C For the GSE data exchange format, refer to the description C in file 'gse.for'. C Default: SS='ss.gse' C Optional parameters: C TMIN=real... Times of the samples in the output seismograms C will differ from TMIN by integer multiples of DT. C Default: TMIN=0. (mostly sufficient) C DT=real... Time interval to digitize the output seismograms. C If DT=0. (default), time TMIN and time interval DT are C taken, for each component at each receiver, from the C first zero or nonzero input seismogram. C Default: DT=0. C GSEWIDTH=positive integer... Width of the output field reserved C for one integer value of the seismogram. Refer to the C description in file 'gse.for'. C C======================================================================= C C External functions and subroutines: EXTERNAL ERROR C ERROR EXTERNAL RSEP1,RSEP2,RSEP3T,RSEP3R,WSEPT,WSEPR,SSEP C RSEP1,RSEP2,RSEP3T,RSEP3R,WSEPT,WSEPR,SSEP... Subroutines of the C Fortran 77 file sep.for (package FORMS). EXTERNAL RGSE1,RGSE2,RGSE2C,WGSE1,WGSE2D,WGSE2C,WGSE2,WGSE3 C RGSE1,RGSE2,RGSE2C,WGSE1,WGSE2D,WGSE2C,WGSE2,WGSE3... Subroutines C of the Fortran 77 file gse.for (package FORMS) C to store seismograms in the GSE data exchange format. EXTERNAL LENGTH INTEGER LENGTH C LENGTH EXTERNAL STRIND C STRIND C The length of character function STRIND is declared later on. EXTERNAL UARRAY REAL UARRAY C UARRAY C C....................................................................... C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C....................................................................... C C Parameters: INTEGER LU,MFILSS,MDIGSS,MREC REAL UNDEF PARAMETER (LU=1,MFILSS=29,MDIGSS=2,MREC=2048) C LU... Logical unit number C MFILSS..Maximum number of input files with synthetic seismograms. C MDIGSS..Number of digits of MFILSS. C MREC... Maximum number of receivers. C UNDEF...Undefined value. Determined later on. C C Length of the external character function: CHARACTER*(4+MDIGSS) STRIND C C Filenames: CHARACTER*80 FSEP,FILESS C C Coefficient of the linear combination corresponding to the file: REAL COEF C C Storing seismograms in memory CHARACTER*6 REC(MREC),SRC(MREC) COMMON/RECC/REC INTEGER KSS(0:3*MREC+1) REAL TMIN(0:3*MREC),DT(0:3*MREC) REAL X1(MREC),X2(MREC),X3(MREC) REAL X1SRC(MREC),X2SRC(MREC),X3SRC(MREC) EQUIVALENCE (KSS(0),RAM(1)) EQUIVALENCE (TMIN(0),RAM(3*MREC+3)) EQUIVALENCE (DT(0),RAM(6*MREC+4)) EQUIVALENCE (X1(1),RAM(9*MREC+5)) EQUIVALENCE (X2(1),RAM(10*MREC+5)) EQUIVALENCE (X3(1),RAM(11*MREC+5)) EQUIVALENCE (X1SRC(1),RAM(12*MREC+5)) EQUIVALENCE (X2SRC(1),RAM(13*MREC+5)) EQUIVALENCE (X3SRC(1),RAM(14*MREC+5)) C REC(I)... Name of the I-th receiver. C KSS(3*I-3+K)... Index of the last sample of K-th seismogram C component at the I-th receiver. C KSS(3*MREC+1)... Index of the last sample of the current input C seismogram C TMIN(3*I-3+K)... Time of the first seismogram sample at the I-th C receiver. C DT(3*I-3+K)... Digitization interval of the seismogram. C X1(I),X2(I),X3(I)... Coordinates of the I-th receiver. C X1SRC(I),X2SRC(I),X3SRC(I)... Coordinates of the point source C corresponding to the I-th receiver. C C Parameters of the seismogram: CHARACTER*1 CHAN,TEXT CHARACTER*6 NAMREC,NAMSRC INTEGER KOMP,NSS REAL X1R,X2R,X3R,X1S,X2S,X3S,T0,TD 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 Other variables: INTEGER IFILSS,NREC,IREC,J,ILEFT,IRIGHT,ISHIFT C IFILSS..Index of the input file. C NREC... Number of already identified receivers. C IREC... Index of the receiver. C J... Index of the seismogram. C ILEFT,IRIGHT,ISHIFT... Extending memory for the new seismogram. C C Temporary storage locations: INTEGER IT,ITRAM,I REAL T,RT,S C C Occupied array RAM: KSS(0)=15*MREC+4 C C Undefined value: UNDEF=UARRAY() C C----------------------------------------------------------------------- C C Reading name of SEP file with input data: WRITE(*,'(A)') '+GSECAL: Enter input filename: ' FSEP=' ' READ(*,*) FSEP WRITE(*,'(A)') '+GSECAL: Working... ' C C Reading all data from the SEP file into the memory: IF (FSEP.NE.' ') THEN CALL RSEP1(LU,FSEP) ELSE C GSECAL-01 CALL ERROR('GSECAL-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. END IF C C Defining index ISEP of the working set of SEP parameters: ISEP=0 CALL SSEP(ISEP,IOLD) CALL SSEP(IOLD,ISEP) C C Reading optional input data: CALL RSEP3R('TMIN',TMIN(0),0.) CALL RSEP3R('DT',DT(0),0.) C C Initial values: C No receiver NREC=0 C Zero seismograms DO 11 I=1,3*MREC+1 KSS(I)=KSS(0) 11 CONTINUE C No seismogram defined DO 12 I=1,3*MREC DT(I)=0. 12 CONTINUE C No parameters from the seismogram comment lines DO 13 I=1,MREC SRC(I)=' ' X1SRC(I)=UNDEF X2SRC(I)=UNDEF X3SRC(I)=UNDEF 13 CONTINUE C C....................................................................... C C Loop for input GSE files: DO 88 IFILSS=1,MFILSS CALL RSEP3T(STRIND('SS',IFILSS),FILESS,' ') CALL RSEP3R(STRIND('COEF',IFILSS),COEF,1.) IF(FILESS.NE.' ') THEN C C Opening the input GSE file WRITE(*,'(2A)') '+GSECAL: Reading ',FILESS(1:63) OPEN(LU,FILE=FILESS,STATUS='OLD') CALL RGSE1(LU,TEXT) C C Loop for input seismograms 20 CONTINUE I=KSS(3*MREC) CALL RGSE2(LU,NAMREC,CHAN,KOMP,X1R,X2R,X3R,T0,TD, * NSS,MRAM-I,RAM(I+1)) IF(NSS.LE.-1) THEN C End of the GSE file GO TO 87 END IF KSS(3*MREC+1)=KSS(3*MREC)+NSS IF(KOMP.LT.1.OR.KOMP.GT.3) THEN C GSECAL-02 CALL ERROR('GSECAL-02: Component other than 1, 2 or 3') C Input GSE file contains seismogram component other than C 1, 2 or 3, which is not allowed by this version of this C program. END IF C C Determining the index of the receiver DO 31 IREC=1,NREC IF(REC(IREC).EQ.NAMREC) THEN C Already identified receiver IF(ABS(X1(IREC)+999.).LT.0.001) THEN X1(IREC)=X1R END IF IF(ABS(X2(IREC)+999.).LT.0.001) THEN X2(IREC)=X2R END IF IF(ABS(X3(IREC)+999.).LT.0.001) THEN X3(IREC)=X3R END IF IF(ABS(X1R+999.).LT.0.001) THEN X1R=X1(IREC) END IF IF(ABS(X2R+999.).LT.0.001) THEN X2R=X2(IREC) END IF IF(ABS(X3R+999.).LT.0.001) THEN X3R=X3(IREC) END IF IF(X1(IREC).NE.X1R.OR. * X2(IREC).NE.X2R.OR. * X3(IREC).NE.X3R) THEN C GSECAL-03 CALL ERROR('GSECAL-03: Different coordinates') C The same receiver has different coordinates in C different input GSE files. END IF GO TO 32 END IF 31 CONTINUE IF(NREC.GE.MREC) THEN C GSECAL-04 CALL ERROR('GSECAL-04: Too many receivers') C Number of receivers exceeds integer constant MREC C specified in the this program. C You may wish to increase parameter MREC. END IF NREC=NREC+1 REC(NREC)=NAMREC IREC=NREC X1(IREC)=X1R X2(IREC)=X2R X3(IREC)=X3R 32 CONTINUE C Index of the seismogram J=3*IREC-3+KOMP C C Reading the comment lines corresponding to the seismogram ISEP=-ISEP CALL SSEP(ISEP,IOLD) 21 CONTINUE CALL RGSE2C(LINE,*22) CALL RSEP2(LINE) GO TO 21 22 CONTINUE CALL RSEP3T('NAMESRC',NAMSRC,' ') CALL RSEP3R('X1SRC',X1S,UNDEF) CALL RSEP3R('X2SRC',X2S,UNDEF) CALL RSEP3R('X3SRC',X3S,UNDEF) CALL SSEP(IOLD,ISEP) IF(NAMSRC.NE.' ') THEN IF(SRC(IREC).EQ.' ') THEN SRC(IREC)=NAMSRC ELSE IF(SRC(IREC).NE.NAMSRC) THEN C GSECAL-05 CALL ERROR('GSECAL-05: Different source name') C The source corresponding to the same receiver has C different name in different input GSE files. END IF END IF IF(X1S.NE.UNDEF) THEN IF(X1SRC(IREC).EQ.UNDEF) THEN X1SRC(IREC)=X1S ELSE IF(X1SRC(IREC).NE.X1S) THEN C GSECAL-06 CALL ERROR('GSECAL-06: Different source coordinates') C The source corresponding to the same receiver has C different X1 coordinate in different input GSE files. END IF END IF IF(X2S.NE.UNDEF) THEN IF(X2SRC(IREC).EQ.UNDEF) THEN X2SRC(IREC)=X2S ELSE IF(X2SRC(IREC).NE.X2S) THEN C GSECAL-07 CALL ERROR('GSECAL-07: Different source coordinates') C The source corresponding to the same receiver has C different X2 coordinate in different input GSE files. END IF END IF IF(X3S.NE.UNDEF) THEN IF(X3SRC(IREC).EQ.UNDEF) THEN X3SRC(IREC)=X3S ELSE IF(X3SRC(IREC).NE.X3S) THEN C GSECAL-08 CALL ERROR('GSECAL-08: Different source coordinates') C The source corresponding to the same receiver has C different X3 coordinate in different input GSE files. END IF END IF C C Defining a new seismogram IF(DT(J).EQ.0.) THEN IF(DT(0).EQ.0.) THEN TMIN(J)=T0 DT(J)=TD ELSE TMIN(J)=TMIN(0) DT(J)=DT(0) END IF END IF C IF(NSS.GE.1) THEN C C Determining the memory for the new seismogram ILEFT=NINT((TMIN(J)-T0+TD)/DT(J)-0.501) IF(KSS(J).EQ.KSS(J-1)) THEN C Zero stored seismogram TMIN(J)=TMIN(J)-FLOAT(ILEFT)*DT(J) ILEFT=0 IRIGHT=NINT((T0+FLOAT(NSS)*TD-TMIN(J))/DT(J)+0.499) ELSE C Non-zero stored seismogram IRIGHT=NINT((T0+FLOAT(NSS)*TD-TMIN(J))/DT(J)+0.499) * -(KSS(J)-KSS(J-1)) ILEFT=MAX0(0,ILEFT) IRIGHT=MAX0(0,IRIGHT) END IF C ILEFT is the number of new samples before the seismogram C IRIGHT is the number of new samples after the seismogram C C Preparing the memory for the new seismogram TMIN(J)=TMIN(J)-FLOAT(ILEFT)*DT(J) ISHIFT=ILEFT+IRIGHT IF(KSS(3*MREC+1)+ISHIFT.GT.MRAM) THEN C GSECAL-09 CALL ERROR('GSECAL-09: Too small array RAM(MRAM)') C Array RAM(MRAM) allocated in include file 'ram.inc' C is too small to contain the samples of the calculated C linear combination of the given seismograms and the C samples of the currently read seismogram. C You may wish to increase the dimension MRAM in file C ram.inc. END IF DO 41 I=J,3*MREC+1 KSS(I)=KSS(I)+ISHIFT 41 CONTINUE DO 42 I=KSS(3*MREC+1),KSS(J)+1,-1 RAM(I)=RAM(I-ISHIFT) 42 CONTINUE DO 43 I=KSS(J),KSS(J)-IRIGHT+1,-1 RAM(I)=0. 43 CONTINUE DO 44 I=KSS(J)-IRIGHT,KSS(J-1)+ILEFT+1,-1 RAM(I)=RAM(I-ILEFT) 44 CONTINUE DO 45 I=KSS(J-1)+ILEFT,KSS(J-1)+1,-1 RAM(I)=0. 45 CONTINUE C C Loop over the samples of the stored seismogram DO 60 I=KSS(J-1)+1,KSS(J) C Time of the stored sample T=TMIN(J)+FLOAT(I-KSS(J-1)-1)*DT(J) C Number of intervals from the first input sample RT=(T-T0)/TD C Integer part of the number of intervals IT=INT(RT+2.)-2 C Fractional part of the number of intervals RT=RT-FLOAT(IT) C Index of the left-hand sample of the interval IT=IT+1 C Index of the left-hand sample of the interval in RAM ITRAM=KSS(3*MREC)+IT IF(0.LE.IT.AND.IT.LE.NSS) THEN C Linear interpolation of the input seismogram IF(IT.EQ.0) THEN S=RT*RAM(ITRAM+1) ELSE IF(IT.EQ.NSS) THEN S=(1.-RT)*RAM(ITRAM) ELSE S=(1.-RT)*RAM(ITRAM)+RT*RAM(ITRAM+1) END IF C Adding the seismogram to the memory RAM(I)=RAM(I)+COEF*S END IF 60 CONTINUE C END IF GO TO 20 C C End of loop for input seismograms 87 CONTINUE CLOSE(LU) C END IF 88 CONTINUE C End of loop for input GSE files C C....................................................................... C C Writing the output GSE file: CALL RSEP3T('SS',FILESS,'ss.gse') OPEN(LU,FILE=FILESS) CALL WGSE1(LU,' ') DO 92 IREC=1,NREC DO 91 KOMP=1,3 J=3*IREC-3+KOMP IF(DT(J).NE.0.) THEN CALL WGSE2D IF(SRC(IREC).NE.' ') THEN CALL WSEPT(LINE,'NAMESRC',SRC(IREC)) CALL WGSE2C(LINE) END IF IF(X1SRC(IREC).NE.UNDEF) THEN CALL WSEPR(LINE,'X1SRC',X1SRC(IREC)) CALL WGSE2C(LINE) END IF IF(X2SRC(IREC).NE.UNDEF) THEN CALL WSEPR(LINE,'X2SRC',X2SRC(IREC)) CALL WGSE2C(LINE) END IF IF(X3SRC(IREC).NE.UNDEF) THEN CALL WSEPR(LINE,'X3SRC',X3SRC(IREC)) CALL WGSE2C(LINE) END IF CALL WGSE2(LU,REC(IREC),' ',KOMP,X1(IREC),X2(IREC),X3(IREC), * TMIN(J),DT(J),KSS(J)-KSS(J-1),RAM(KSS(J-1)+1)) END IF 91 CONTINUE 92 CONTINUE CALL WGSE3(LU) CLOSE(LU) C WRITE(*,'(A)') '+GSECAL: 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 INCLUDE 'forms.for' C forms.for C C======================================================================= C