C
C Program GREENSS to read a formatted file containing the ray-theory C elastodynamic Green function and to generate ray-theory time-domain C synthetic seismograms (without attenuation) or frequency-domain C response functions (including causal Futterman's attenuation). C C Version: 5.10 C Date: 1997, September 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....................................................................... C C C Description of data files: C C Main input data read from external interactive device (*): C The data consist of character strings read by list directed (free C format) input. The strings have thus to be enclosed in C apostrophes. The interactive * external unit may be redirected to C the file containing the data. C (1) 'GREEN','SOURCE','FRDAT','SIGNAL','SS',KOMP C 'GREEN'... Name of the input formatted file with the Green tensor. C Description of file GREEN C 'SOURCE'... Name of the input formatted file containing the C complex-valued seismic force or moment. C Description of file SOURCE C 'FRDAT'... Name of the input formatted file containing the C specification of the frequency domain for the calculation C of response functions. C If blank, the time-domain synthetic seismograms are C generated instead of frequency-domain response functions, C see 'SIGNAL'. C Description of file FRDAT C 'SIGNAL'... Name of the input file in the GSE format containing C the source time function and its Hilbert transform (for C seismic force), or their derivatives (for seismic moment). C To be submitted if 'SIGNAL'=' ', otherwise has no meaning. C (One of filenames 'FRDAT' and 'SIGNAL' must be nonblank.) C Description of GSE C 'RF'... Name of the output file in the GSE format containing the C ray-theory seismograms or the frequency-domain Response C Function. C Description of file RF C KOMP... KOMP=0: All 3 components of the synthetic seismograms are C stored in the output GSE file. C KOMP=1: The 1st component of the synthetic seismograms is C stored in the output GSE file. C KOMP=2: The 2nd component of the synthetic seismograms is C stored in the output GSE file. C KOMP=3: The 3rd component of the synthetic seismograms is C stored in the output GSE file. C Default: 'GREEN'='green.out', 'SOURCE'='source.dat', C 'FRDAT'='fr.dat', 'SIGNAL'='source.gse', 'RF'='rf.out'. C Note: File 'SIGNAL'='source.gse' need not exist (i.e. is C not read) if 'SOURCE' is not blank. C C C Input formatted file GREEN: C (1) / (a slash). C (2) For each two-point ray (2.1): C (2.1) 'R','S',(GREEN(I),I=1,32),/ C 'R'... String in apostrophes identifying the receiver. Only C the first 6 characters are written to the output GSE C file. The strings corresponding to different receivers C thus should, if possible, differ in the first 6 C characters. If this is not the case, at least in the C first 12 characters. All lines corresponding to the same C receiver must be consecutive. C 'S'... String in apostrophes describing the source. Not taken C into account. C GREEN(1)... Travel time between receiver and source. C GREEN(2)... Imaginary part of the complex-valued travel time C between receiver and source due to attenuation. C GREEN(3:8)... Coordinates of the receiver and coordinates of the C source. C GREEN(9:14)... Derivatives of the travel time with respect to the C coordinates of the receiver and coordinates of the source. C GREEN(15:32)... 1000000 times enlarged amplitude of the Green C function: contravariant components of the complex-valued C 3*3 matrix Gij in model coordinates, where the first C subscript corresponds to the receiver and the second C subscript corresponds to the source. The components are C ordered as C Re(G11),Im(G11),Re(G21),Im(G21),Re(G31),Im(G31), C Re(G12),Im(G12),Re(G22),Im(G22),Re(G32),Im(G32), C Re(G13),Im(G13),Re(G23),Im(G23),Re(G33),Im(G33). C /... An obligatory slash after at the end of line, in place C where the slowness vector components could be written. C (3) / (a slash). C C C Input formatted file SOURCE: C (1) SFR1,SFI1,SFR2,SFI2,SFR3,SFI3,/ C Components of the complex-valued vectorial seismic force. The C seismic force is assumed to be the product of this vector and the C source time function submitted in file 'SIGNAL'. C Note: The 'unit' radiation pattern corresponds to C SF = 4 pi rho v**2 C (2) SMR11,SMI11,SMR12,SMI12,SMR13,SMI13, C SMR21,SMI21,SMR22,SMI22,SMR23,SMI23, C SMR31,SMI31,SMR32,SMI32,SMR33,SMI33,/ C Components of the transposed complex-valued 3*3 seismic-moment C tensor. The tensor is transposed in order not to look transposed C in the input data. The time derivative of the seismic moment is C assumed to be the product of this vector and the time function C submitted in file 'SIGNAL'. C Note: The 'unit' radiation pattern corresponds to C SM = 4 pi rho v**3 C Example of data file SOURCE C C C Input formatted file FRDAT: C (1) NPTS,FMIN,FMAX,TD,TINT,TIMUL,FREF,/ C NPTS... Number of time steps for the fast Fourier transform. C Will be used to convert the time step to the frequency C step. C NPTS=0: Frequency step is specified instead of time step. C FMIN,FMAX... Response functions are calculated for frequencies C from FMIN to FMAX. FMIN and FMAX are rounded to the C nearest multiples of the time step. C TD... NPTS=0: Frequency step. C Otherwise: Time step. Frequency step=1/(NPTS*TD). C TINT... Maximum time interval. The contribution of the ray to C the seismogram is taken into account only if the travel C time does not exceed TMIN+TINT, where TMIN is the minimum C travel time over the preceding rays arriving at the ' C receiver. Useful to remove alias in the time domain. C TINT=0: Infinite time interval. C TIMUL...Multiplication factor for the imaginary part TI of the C travel time. It may be used to globally decrease or C increase attenuation in the whole model. C FREF... Reference frequency for the Futterman's (1962) C quasi-causal attenuation. The travel times TR and TI C describing the Green function are assumed to correspond C to this frequency. Frequency-dependent travel times are C then given by C Re TT(F)=TR-TI*ln(F/FREF)*2/pi, C Im TT(F)=TI. C FREF=0: Noncausal attenuation, just for test purposes, C Re TT(F)=TR, C Im TT(F)=TI. C Defaults: NPTS=0, FMIN=0, FMAX=1/(NPTS*TD) if NPTS.NE.0, TINT=0, C TIMUL=1, FREF=(FMIN+FMAX)/2. C C C Output formatted file RF: C If 'FRDAT'=' ': Ray-theory time-domain synthetic seismograms (without C attenuation) in the GSE format. The may be, e.g., plotted by C program 'sp.for'. C Program 'greenss.for' stores in the comment lines of the waveform C identification section the hypocentral coordinates identified by C strings 'XS1 ', 'XS2 ' and 'XS3 '. C Description of the GSE format C Otherwise: Ray-theory frequency-domain response functions (including C causal Futterman's attenuation), saved in the format 'RF' C described in 'ss.for'. C Description of file RF C Program 'greenss.for' is prepared to generate the response C functions also in the GSE format in future versions. In such a C case program 'greenss.for' would store in the comment lines of the C waveform identification section the hypocentral coordinates C identified by strings 'XS1 ', 'XS2 ' and 'XS3 ', and times of the C first and last considered arrivals to each receiver, identified by C 'TMIN' and 'TMAX'. C Description of the GSE format C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL WGSE1,WGSE2,WGSE3,RGSE2 C WGSE1,WGSE2,WGSE3,RGSE2... File 'gse.for'. C C----------------------------------------------------------------------- C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C Allocation of working arrays: INTEGER MSEIS PARAMETER (MSEIS=MRAM/6) REAL SEIS1(MSEIS),SEIS2(MSEIS),SEIS3(MSEIS) REAL SEIS4(MSEIS),SEIS5(MSEIS),SEIS6(MSEIS) EQUIVALENCE (SEIS1,RAM ) EQUIVALENCE (SEIS2,RAM( MSEIS+1)) EQUIVALENCE (SEIS3,RAM(2*MSEIS+1)) EQUIVALENCE (SEIS4,RAM(3*MSEIS+1)) EQUIVALENCE (SEIS5,RAM(4*MSEIS+1)) EQUIVALENCE (SEIS6,RAM(5*MSEIS+1)) C C----------------------------------------------------------------------- C C Input and output files: INTEGER LU1,LU2 PARAMETER (LU1=1,LU2=2) CHARACTER*80 FILE1,FILE2,FILE3,FILE4,FILE5,LINE CHARACTER*1 TEXT C Undefined value: REAL UNDEF PARAMETER (UNDEF=-999999.) C Seismic force: REAL SFR1,SFI1,SFR2,SFI2,SFR3,SFI3 C Seismic moment: REAL SMR11,SMI11,SMR12,SMI12,SMR13,SMI13 REAL SMR21,SMI21,SMR22,SMI22,SMR23,SMI23 REAL SMR31,SMI31,SMR32,SMI32,SMR33,SMI33 C Green function: CHARACTER*12 TXTOLD,TXTREC,TXTSRC REAL GREEN(32) C Data for the frequency band: INTEGER NPTS,NF REAL FMIN,FMAX,TD,TINT,TIMUL,FREF,FSTEP C Seismograms: LOGICAL LGSE,LWRITE INTEGER NSEIS,NSEIS4,NSEIS5 REAL TSTAR0,TSTEP0,TSTARH,TSTART,TSTEP C Force at the source and displacement at the receiver: REAL FR1,FI1,FR2,FI2,FR3,FI3,AR1,AI1,AR2,AI2,AR3,AI3 C Temporary storage locations: CHARACTER*1 STAT,CHAN INTEGER I01,I02,I03,I04,I05,I06,I07,I08,I09,I10,I11,I12 INTEGER KOMP,ISHIFT,I REAL AMAX,AR,AI,TR,TI,F,OMEGA REAL XR1,XR2,XR3,XS1,XS2,XS3,TSHIFT,W0,W1 C Source coordinates transferred through the GSE file: INTEGER MCOM,NCOM PARAMETER (MCOM=5) CHARACTER*4 TCOM(MCOM) REAL VCOM(MCOM) DATA TCOM/'XS1 ','XS2 ','XS3 ','TMIN','TMAX'/ C C....................................................................... C C Format of the output response function (.TRUE.'GSE', .FALSE.'RF'): LGSE=.FALSE. C C Names of input and output files: FILE1='green.out' FILE2='source.dat' FILE3='fr.dat' FILE4='source.gse' FILE5='rf.gse' KOMP=0 WRITE(*,'(A)') * ' Enter green.out,source.dat,fr.dat,source.gse,rf.gse,KOMP: ' READ(*,*) FILE1,FILE2,FILE3,FILE4,FILE5,KOMP C C Reading seismic force or seismic moment: OPEN(LU2,FILE=FILE2,STATUS='OLD') SFR1=0. SFI1=0. SFR2=0. SFI2=0. SFR3=0. SFI3=0. READ(LU2,*) SFR1,SFI1,SFR2,SFI2,SFR3,SFI3 SMR11=0. SMI11=0. SMR12=0. SMI12=0. SMR22=0. SMI22=0. SMR13=0. SMI13=0. SMR23=0. SMI23=0. SMR33=0. SMI33=0. READ(LU2,*) SMR11,SMI11,SMR12,SMI12,SMR13,SMI13, * SMR21,SMI21,SMR22,SMI22,SMR23,SMI23, * SMR31,SMI31,SMR32,SMI32,SMR33,SMI33 CLOSE(LU2) C C Reading the data for the frequency domain: IF(FILE3.NE.' ') THEN OPEN(LU2,FILE=FILE3,STATUS='OLD') NPTS=0 FMIN=0. FMAX=UNDEF TINT=0. TIMUL=1. FREF=UNDEF READ(LU2,*) NPTS,FMIN,FMAX,TD,TINT,TIMUL,FREF IF(FMAX.EQ.UNDEF) THEN FMAX=1/TD END IF IF(FREF.EQ.UNDEF) THEN FREF=(FMIN+FMAX)/2. END IF CLOSE(LU2) IF(NPTS.EQ.0) THEN FSTEP=TD ELSE FSTEP=1./(FLOAT(NPTS)*TD) END IF FMIN=FSTEP*INT(FMIN/FSTEP+.5) NF= INT((FMAX-FMIN)/FSTEP+.5)+1 C Parameters for the response functions: FMIN,FSTEP,NF,TIMUL,FREF. NCOM=MCOM ELSE NCOM=3 C C Reading the source time function and its Hilbert transform: IF(FILE4.NE.' ') THEN OPEN(LU2,FILE=FILE4,STATUS='OLD') CALL RGSE2(LU2,STAT,CHAN, * I,XS1,XS2,XS3,TSTAR0,TSTEP0,NSEIS4,MSEIS,SEIS4) CALL RGSE2(LU2,STAT,CHAN, * I,XS1,XS2,XS3,TSTARH,TSTEP ,NSEIS5,MSEIS,SEIS5) CLOSE(LU2) IF(TSTEP0.NE.TSTEP) THEN C GREENSS-01 PAUSE 'Error GREENSS-01: Different time steps.' STOP END IF ELSE C GREENSS-02 PAUSE * 'Error GREENSS-02: One of files FRDAT or SIGNAL must be given' STOP END IF C END IF C C....................................................................... C OPEN(LU1,FILE=FILE1,STATUS='OLD') READ(LU1,*) (TEXT,I=1,20) OPEN(LU2,FILE=FILE5) IF(NCOM.NE.MCOM) THEN CALL WGSE1(LU2,' ') ELSE IF(LGSE) THEN CALL WGSE1(LU2,' ') ELSE LWRITE=.TRUE. END IF END IF C C Loop over rays: TXTREC='$' 30 CONTINUE TXTOLD=TXTREC TXTREC='$' XR1=GREEN(3) XR2=GREEN(4) XR3=GREEN(5) VCOM(1)=GREEN(6) VCOM(2)=GREEN(7) VCOM(3)=GREEN(8) C Preparing source coordinates for output in the GSE file: CALL WGSE2D() DO 31 I=1,NCOM CALL WSEPR(LINE,TCOM(I),VCOM(I)) CALL WGSE2C(LINE) 31 CONTINUE C READ(LU1,*) TXTREC,TXTSRC,GREEN IF(TXTOLD.NE.'$'.AND.TXTREC.NE.TXTOLD) THEN C Writing the previous seismogram: C -------------------------------- IF(NCOM.NE.MCOM) THEN IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN CALL WGSE2(LU2,TXTOLD,' ',1,XR1,XR2,XR3, * TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS1) END IF IF(KOMP.LE.0.OR.KOMP.EQ.2) THEN CALL WGSE2(LU2,TXTOLD,' ',2,XR1,XR2,XR3, * TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS2) END IF IF(KOMP.LE.0.OR.KOMP.EQ.3) THEN CALL WGSE2(LU2,TXTOLD,' ',3,XR1,XR2,XR3, * TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS3) END IF IF(NSEIS.GT.MSEIS) THEN C GREENSS-51 WRITE(*,'(A,I12,2A)') * ' Warning GREENSS-51:',NSEIS, * ' non-zero samples at receiver ',TXTOLD PAUSE END IF ELSE IF(LGSE) THEN IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN CALL WGSE2(LU2,TXTOLD,' ', 1,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS1) CALL WGSE2(LU2,TXTOLD,' ',11,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS4) END IF IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN CALL WGSE2(LU2,TXTOLD,' ', 2,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS2) CALL WGSE2(LU2,TXTOLD,' ',12,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS5) END IF IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN CALL WGSE2(LU2,TXTOLD,' ', 3,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS3) CALL WGSE2(LU2,TXTOLD,' ',13,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS6) END IF ELSE IF(LWRITE)THEN WRITE(LU2,'(A)') '/' WRITE(LU2,'(3(F7.3,1X),A)') VCOM(1),VCOM(2),VCOM(3),'/' WRITE(LU2,'(2(E11.5,1X),I8,1X,A)') FMIN,FSTEP,NF,'/' LWRITE=.FALSE. END IF AMAX=0. DO 32 I=1,NF AMAX=AMAX1(ABS(SEIS1(I)),ABS(SEIS2(I)),ABS(SEIS3(I)), * ABS(SEIS4(I)),ABS(SEIS5(I)),ABS(SEIS6(I)), * AMAX) 32 CONTINUE WRITE(LU2,'(3(F7.3,1X),3(E11.5,1X),A)') * XR1,XR2,XR3,VCOM(4),VCOM(5),AMAX,'/' IF(VCOM(4).LE.VCOM(5)) THEN DO 33 I=1,NF,2 I01=IFIX(99999.1*SEIS1(I)/AMAX) I02=IFIX(99999.1*SEIS4(I)/AMAX) I03=IFIX(99999.1*SEIS2(I)/AMAX) I04=IFIX(99999.1*SEIS5(I)/AMAX) I05=IFIX(99999.1*SEIS3(I)/AMAX) I06=IFIX(99999.1*SEIS6(I)/AMAX) IF(I.LT.NF) THEN I07=IFIX(99999.1*SEIS1(I+1)/AMAX) I08=IFIX(99999.1*SEIS4(I+1)/AMAX) I09=IFIX(99999.1*SEIS2(I+1)/AMAX) I10=IFIX(99999.1*SEIS5(I+1)/AMAX) I11=IFIX(99999.1*SEIS3(I+1)/AMAX) I12=IFIX(99999.1*SEIS6(I+1)/AMAX) WRITE(LU2,'(12I6)') I01,I02,I03,I04,I05,I06, * I07,I08,I09,I10,I11,I12 ELSE WRITE(LU2,'(12I6)') I01,I02,I03,I04,I05,I06 END IF 33 CONTINUE END IF END IF END IF END IF IF(TXTREC.EQ.'$') THEN C No more two-point rays: C ----------------------- GO TO 80 END IF IF(TXTREC.NE.TXTOLD) THEN C New receiver: C ------------- DO 41 I=1,MSEIS SEIS1(I)=0. SEIS2(I)=0. SEIS3(I)=0. 41 CONTINUE IF(NCOM.NE.MCOM) THEN ISHIFT=INT(GREEN(1)/TSTEP) TSTART=TSTAR0+FLOAT(ISHIFT)*TSTEP NSEIS=0 ELSE VCOM(4)= 999999. VCOM(5)=-999999. DO 42 I=1,MSEIS SEIS4(I)=0. SEIS5(I)=0. SEIS6(I)=0. 42 CONTINUE END IF END IF DO 43 I=15,32 GREEN(I)=GREEN(I)/1000000. 43 CONTINUE C C Adding the contribution from a new two-point ray: C ------------------------------------------------- C C Complex-valued amplitude corresponding to the given source: FR1=SFR1-SMR11*GREEN(12)-SMR12*GREEN(13)-SMR13*GREEN(14) FI1=SFI1-SMI11*GREEN(12)-SMI12*GREEN(13)-SMI13*GREEN(14) FR2=SFR2-SMR21*GREEN(12)-SMR22*GREEN(13)-SMR23*GREEN(14) FI2=SFI2-SMI21*GREEN(12)-SMI22*GREEN(13)-SMI23*GREEN(14) FR3=SFR3-SMR31*GREEN(12)-SMR32*GREEN(13)-SMR33*GREEN(14) FI3=SFI3-SMI31*GREEN(12)-SMI32*GREEN(13)-SMI33*GREEN(14) AR1=GREEN(15)*FR1+GREEN(21)*FR2+GREEN(27)*FR3 * -GREEN(16)*FI1-GREEN(22)*FI2-GREEN(28)*FI3 AR2=GREEN(17)*FR1+GREEN(23)*FR2+GREEN(29)*FR3 * -GREEN(18)*FI1-GREEN(24)*FI2-GREEN(30)*FI3 AR3=GREEN(19)*FR1+GREEN(25)*FR2+GREEN(31)*FR3 * -GREEN(20)*FI1-GREEN(26)*FI2-GREEN(32)*FI3 AI1=GREEN(15)*FI1+GREEN(21)*FI2+GREEN(27)*FI3 * +GREEN(16)*FR1+GREEN(22)*FR2+GREEN(28)*FR3 AI2=GREEN(17)*FI1+GREEN(23)*FI2+GREEN(29)*FI3 * +GREEN(18)*FR1+GREEN(24)*FR2+GREEN(30)*FR3 AI3=GREEN(19)*FI1+GREEN(25)*FI2+GREEN(31)*FI3 * +GREEN(20)*FR1+GREEN(26)*FR2+GREEN(32)*FR3 C C Time domain or frequency domain: IF(NCOM.NE.MCOM) THEN C C Adding the multiple of the source time function: TSHIFT=(TSTAR0+GREEN(1)-TSTART)/TSTEP ISHIFT=INT(TSHIFT) W1=TSHIFT-FLOAT(ISHIFT) W0=1.-W1 C W0,W1 are the weights of shifts ISHIFT,ISHIFT+1 IF(ISHIFT.LT.0) THEN C Shifting the start time to the new position: TSTART=TSTART+FLOAT(ISHIFT)*TSTEP NSEIS=NSEIS-ISHIFT DO 51 I=MAX0(MSEIS,NSEIS),-ISHIFT+1,-1 SEIS1(I)=SEIS1(I+ISHIFT) SEIS2(I)=SEIS2(I+ISHIFT) SEIS3(I)=SEIS3(I+ISHIFT) 51 CONTINUE DO 52 I=MAX0(MSEIS,-ISHIFT),1,-1 SEIS1(I)=0. SEIS2(I)=0. SEIS3(I)=0. 52 CONTINUE ISHIFT=0 END IF NSEIS=MAX0(NSEIS,NSEIS4+ISHIFT) DO 53 I=1,MIN(NSEIS4,MSEIS-ISHIFT) SEIS1(I+ISHIFT)=SEIS1(I+ISHIFT)+W0*AR1*SEIS4(I) SEIS2(I+ISHIFT)=SEIS2(I+ISHIFT)+W0*AR2*SEIS4(I) SEIS3(I+ISHIFT)=SEIS3(I+ISHIFT)+W0*AR3*SEIS4(I) 53 CONTINUE ISHIFT=ISHIFT+1 NSEIS=MAX0(NSEIS,NSEIS4+ISHIFT) DO 54 I=1,MIN(NSEIS4,MSEIS-ISHIFT) SEIS1(I+ISHIFT)=SEIS1(I+ISHIFT)+W1*AR1*SEIS4(I) SEIS2(I+ISHIFT)=SEIS2(I+ISHIFT)+W1*AR2*SEIS4(I) SEIS3(I+ISHIFT)=SEIS3(I+ISHIFT)+W1*AR3*SEIS4(I) 54 CONTINUE C C Adding the multiple of the Hilbert transform: TSHIFT=(TSTARH+GREEN(1)-TSTART)/TSTEP ISHIFT=INT(TSHIFT) W1=TSHIFT-FLOAT(ISHIFT) W0=1.-W1 C W0,W1 are the weights of shifts ISHIFT,ISHIFT+1 IF(ISHIFT.LT.0) THEN C Shifting the start time to the new position: TSTART=TSTART+FLOAT(ISHIFT)*TSTEP NSEIS=NSEIS-ISHIFT DO 61 I=MIN0(MSEIS,NSEIS),-ISHIFT+1,-1 SEIS1(I)=SEIS1(I+ISHIFT) SEIS2(I)=SEIS2(I+ISHIFT) SEIS3(I)=SEIS3(I+ISHIFT) 61 CONTINUE DO 62 I=MIN0(MSEIS,-ISHIFT),1,-1 SEIS1(I)=0. SEIS2(I)=0. SEIS3(I)=0. 62 CONTINUE ISHIFT=0 END IF NSEIS=MAX0(NSEIS,NSEIS5+ISHIFT) DO 63 I=1,MIN(NSEIS5,MSEIS-ISHIFT) SEIS1(I+ISHIFT)=SEIS1(I+ISHIFT)-W0*AI1*SEIS5(I) SEIS2(I+ISHIFT)=SEIS2(I+ISHIFT)-W0*AI2*SEIS5(I) SEIS3(I+ISHIFT)=SEIS3(I+ISHIFT)-W0*AI3*SEIS5(I) 63 CONTINUE ISHIFT=ISHIFT+1 NSEIS=MAX0(NSEIS,NSEIS5+ISHIFT) DO 64 I=1,MIN(NSEIS5,MSEIS-ISHIFT) SEIS1(I+ISHIFT)=SEIS1(I+ISHIFT)-W1*AI1*SEIS5(I) SEIS2(I+ISHIFT)=SEIS2(I+ISHIFT)-W1*AI2*SEIS5(I) SEIS3(I+ISHIFT)=SEIS3(I+ISHIFT)-W1*AI3*SEIS5(I) 64 CONTINUE ELSE C C Response functions: VCOM(4)=AMIN1(GREEN(1),VCOM(4)) IF(TINT.EQ.0..OR.GREEN(1).LE.VCOM(4)+TINT) THEN VCOM(5)=AMAX1(GREEN(1),VCOM(5)) DO 69 I=1,NF F=FMIN+FSTEP*FLOAT(I-1) OMEGA=6.2831853*F TI=GREEN(2)*TIMUL IF(FREF.EQ.0.) THEN TR=GREEN(1) ELSE TR=GREEN(1)-TI*ALOG(F/FREF)/1.5707963 END IF TI= EXP(-OMEGA*TI) AR=TI*COS(OMEGA*TR) AI=TI*SIN(OMEGA*TR) SEIS1(I)=SEIS1(I)+AR1*AR-AI1*AI SEIS2(I)=SEIS2(I)+AR2*AR-AI2*AI SEIS3(I)=SEIS3(I)+AR3*AR-AI3*AI SEIS4(I)=SEIS4(I)+AI1*AR+AR1*AI SEIS5(I)=SEIS5(I)+AI2*AR+AR2*AI SEIS6(I)=SEIS6(I)+AI3*AR+AR3*AI 69 CONTINUE END IF END IF GO TO 30 C 80 CONTINUE IF(NCOM.NE.MCOM) THEN CALL WGSE3(LU2) ELSE IF(LGSE)THEN CALL WGSE3(LU2) ELSE WRITE(LU2,'(A)') '/' END IF END IF CLOSE(LU2) CLOSE(LU1) STOP END C C======================================================================= C INCLUDE 'sep.for' C sep.for INCLUDE 'gse.for' C gse.for INCLUDE 'length.for' C length.for C C======================================================================= C