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 Program GREENSS can also read frequency-dependent coupling-ray-theory C Green functions corresponding to S waves in weakly anisotropic models. C C Version GREENRM of program GREENSS may include the response of the C transmission through a stack of fine horizontal layers at each C receiver location. C The frequency-dependent transmission coefficiens are calculated by C subroutines of package RMATRIX written by C. J. Thomson (Copyright (c) C 1996 C. J. Thomson). Please, refer to package C RMATRIX for copyright conditions C and details on compilation. C Version GREENRM is derived from the basic version of program GREENSS C by changing the asterisks '*' in the first column of the INCLUDE C statements at the end of this file. C C Version: 5.40 C Date: 2000, January 30 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 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 input files: C GREEN='string'... Name of the input formatted file with the Green C tensor. The file is generated by program 'green.for'. C Description of file GREEN C Default: GREEN='green.out' (mostly sufficient) C SOURCE='string'... Name of the input formatted file containing the C complex-valued seismic force or moment. C Description of file SOURCE C Default: SOURCE='source.dat' C RMATRIX='string'... Name of optional formatted input file with the C parameters of fine layered structures at the receiver C sites. If non-blank, the frequency-dependent transmission C coefficiens, calculated by subroutines of package RMATRIX, C are considered at each receiver. C Available just for version GREENRM of program GREENSS, see C the comment at the beginning of this file. C Default: RMATRIX=' ' C SIGGSE='string'... Name of the optional input file in the GSE C format containing the source time function and its Hilbert C transform (for seismic force), or their derivatives (for C seismic moment). C If blank (default), the frequency-domain response C functions are generated. C If non-blank, the time-domain synthetic seismograms are C generated instead of frequency-domain response functions. C It is recommended not to specify SIGGSE (use default). C Description of GSE C Default: SIGGSE=' ' (recommended) C Names of output files: C RF='string'... Name of the output file containing the C frequency-domain Response Function. C Generated just if SIGGSE=' ' (default). C Description of file RF C Default: RF='rf.out' (mostly sufficient) C SS='string'... Name of the optional output file in the GSE format C containing the ray-theory seismograms calculated in the C time domain. Generated just if SIGGSE.NE.' '. C Description of file SS C Default: SS='ss.gse' C Data describing the frequency domain common with program 'ss.for': C These data are used to generate the optimum default values of C input parameters DF, OF and NF compatible with data for C 'ss.for'. C DT=real... Time step to digitize the source time function and C seismograms. C Default: DT=1. C NFFT=integer... Number of the time samples for the fast Fourier C transform. Will be used to convert the time step to the C frequency step. C NFFT must be an integer power of 2. C Default: NFFT=1 C FMIN=real... The lowest frequency of the cosine band-pass filter. C Default: FMIN=0 (mostly sufficient) C FMAX=real... The highest frequency of the cosine band-pass filter. C Default: FMAX=0.5/DT (Nyquist frequency, mostly C sufficient) C Data describing the frequency domain specific to this program: C DF=real ... Frequency step. C Default: DF=1/(NPTS*DT) (recommended) C OF=real ... Response functions are calculated for NF frequencies C OF,OF+DF,OF+2*DF,...,OF+(NF-1)*DF. C Default: OF=DF*NINT(FMIN/DF) (i.e. FMIN rounded to the C nearest multiple of the frequency step, recommended) C NF=integer ... Number of frequencies. C Default: NF=NINT((FMAX-OF)/DF)+1 (recommended) C Hint: Do not specify DF, OF and NF. Use the data for program C 'ss.for' instead and leave the default values of DF, OF and NF C if you are going to generate synthetic seismograms by 'ss.for'. C Specify explicitly DF, OF and NF only if you generate the response C function for another purpose than the generation of synthetic C seismograms by program 'ss.for'. C Data specific to this program: C TINT=non-negative real... Maximum time interval. The contribution C of the ray to the seismogram is taken into account only if C the travel time does not exceed TMIN+TINT, where TMIN is C the minimum travel time over the preceding rays arriving C at the receiver. Useful to remove alias in the time C domain. C TINT=0: Infinite time interval. C Default: TINT=0. (mostly sufficient for reasonable NPTS) C TIMUL=non-negative real... Multiplication factor for the imaginary C part TI of the travel time. It may be used to globally C decrease or increase attenuation in the whole model. C Default: TIMUL=1. (mostly sufficient) C FREF=non-negative real... Reference frequency for the Futterman's C (1962) quasi-causal attenuation. The travel times TR and C TI 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 Default: FREF=(FMIN+FMAX)/2. 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 'SIGGSE'. 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 'SIGGSE'. 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 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 GREEN(33:*)... Analogy of GREEN(15:32) for other frequencies. C Present just for frequency-dependent elementary Green C functions. C /... A slash after at the end of line. It may also stand for C zeros at the end of line. C (3) / (a slash). C C C Optional input formatted file RMATRIX: C Parameters of fine layered structures at the receiver sites. C Once or for each receiver the following data: C (1) 'REC',/ C 'REC'...Blank string if a single fine layered structure applies to C all receivers. Otherwise the name of the receiver. C The order and names of the receivers must be the same as C in file 'GREEN'. The order is given by the receiver file. C If the receiver file is specified for the GREEN program, C it also determines the receiver names. C Default: 'REC'=' '. C (2) NLAY,IFREE C NLAY... Number of layers, including the bottom and top halfspaces. C IFREE...IFREE=1: The thickness THICKN(1) of the first (top) layer C is used for the calculation of the travel time C correction at the geophone situated in the top halfspace C at the elevation THICKN(1) above the uppermost C interface. C IFREE=0: Has the same effect as THICKN(1)=0. C (3) ISOFLG,THICKN,RHO C ISOFLG..ISOFLG=0: Anisotropic layer. C ISOFLG=1: Isotropic elastic layer. C ISOFLG=2: Fluid layer. Does not work in the current C version. ISOFLG=1 with a very small S wave velocity is C suggested instead. C THICKN..Thickness of the layer. THICKN(1) and THICKN(NLAY) are C used to position the layer with respect to the coordinates C of the receiver. The thicknesses may also be negative. C The lowermost interface (top of the bottom halfspace) is C situated at the elevation THICKN(NLAY) above the receiver C coordinates. The vertical component of the slowness C vector is used to linearly extrapolate travel time from C the receiver coordinates to the bottom of the layer stack. C The geophone is considered to be situated in the elevation C THICKN(1)+THICKN(2)+...+THICKN(NLAY) above the receiver. C Example: C For a receiver at the Earth surface, layer 1 has C THICK(1)=0 (geophone exactly at the surface), RHO(1) C equal to the density of the air (about 0.0012kg/l) or C smaller, VP(1) equal to the sound speed in the air C (about 0.34km/s), and VS(1) very small. Layers 2 to C NLAY-1 correspond to the properties of the fine layered C structure. Layer NLAY usually takes negative thickness C THICKN(NLAY)=-THICKN(1)-THICKN(2)-...-THICKN(NLAY) C in order to situate geophone exactly at the receiver C coordinates. C RHO... Density of the layer. C (4) For ISOFLG=0 (anisotropic layer) (4.1), otherwise (4.2): C (4.1) (A(I,J,K,L),I=1,3,J=1,3,K=1,3,L=1,3) C A(I,J,K,L)... Matrix of 81 elastic parameters C(I,J,K,L) divided C by density RHO. C (4.2) VP,VS C VP... P wave velocity. C VS... S wave velocity. May be very small for fluids but cannot C be zero in this version. C C C Output formatted file RF: C 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 Optional output formatted file SS: C Ray-theory time-domain synthetic seismograms (without C attenuation) in the GSE format. They 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 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 MGREEN,MSEIS PARAMETER (MGREEN=MRAM/4*3) PARAMETER (MSEIS=(MRAM-MGREEN)/6) REAL SEIS1(MSEIS),SEIS2(MSEIS),SEIS3(MSEIS) REAL SEIS4(MSEIS),SEIS5(MSEIS),SEIS6(MSEIS) REAL GREEN(MGREEN) 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)) EQUIVALENCE (GREEN,RAM(6*MSEIS+1)) C Frequency-dependent Green function requires 14+18*NF storage C locations. C C----------------------------------------------------------------------- C C Input and output files: INTEGER LU1,LU2,LU3 PARAMETER (LU1=1,LU2=2,LU3=3) CHARACTER*80 FILE1,FILE2,FILSEP,FILSIG,FILOUT,FILRM,LINE CHARACTER*80 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 Strings to identify elementary Green functions: CHARACTER*12 TXTOLD,TXTREC,TXTSRC,TXTRM C Green function: INTEGER NGREEN C Input and output arguments for RMATRIX: INTEGER NLAY,IFREE 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 REAL AR1,AI1,AR2,AI2,AR3,AI3,BR1,BI1,BR2,BI2,BR3,BI3 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,J 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 Reading name of SEP file with input data: WRITE(*,'(A)') '+GREENSS: Enter input filename: ' FILSEP=' ' READ(*,*) FILSEP C C Reading all data from the SEP file into the memory: IF (FILSEP.NE.' ') THEN CALL RSEP1(LU1,FILSEP) ELSE C GREENSS-03 CALL ERROR('GREENSS-03: 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 WRITE(*,'(A)') '+GREENSS: Working... ' C C Reading seismic force or seismic moment: CALL RSEP3T('SOURCE',FILE2 ,'source.dat') OPEN(LU1,FILE=FILE2,STATUS='OLD') SFR1=0. SFI1=0. SFR2=0. SFI2=0. SFR3=0. SFI3=0. READ(LU1,*) 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(LU1,*) SMR11,SMI11,SMR12,SMI12,SMR13,SMI13, * SMR21,SMI21,SMR22,SMI22,SMR23,SMI23, * SMR31,SMI31,SMR32,SMI32,SMR33,SMI33 CLOSE(LU1) C CALL RSEP3T('SIGGSE',FILSIG,' ') IF(FILSIG.EQ.' ') THEN C C Reading the data for the frequency domain: CALL RSEP3I('NFFT ',NPTS ,1) CALL RSEP3R('DT ',TD ,1.) CALL RSEP3R('FMIN ',FMIN ,0.) CALL RSEP3R('FMAX ',FMAX ,0.5/TD) CALL RSEP3R('TINT ',TINT ,0.) CALL RSEP3R('TIMUL',TIMUL ,1.) CALL RSEP3R('FREF ',FREF ,(FMIN+FMAX)/2.) CALL RSEP3R('DF ',FSTEP ,1./(FLOAT(NPTS)*TD)) CALL RSEP3R('OF ',FMIN ,FSTEP*NINT(FMIN/FSTEP)) CALL RSEP3I('NF ',NF ,NINT((FMAX-FMIN)/FSTEP)+1) IF(14+18*NF.GT.MGREEN) THEN C GREENSS-01 CALL ERROR('GREENSS-01: Small array for Green functions') C Frequency-dependent Green function requires 14+18*NF storage C locations. Array GREEN located within array RAM declared in C include file ram.inc has C MGREEN=MRAM/4*3 locations. END IF CALL RSEP3T('RMATRIX',FILRM,' ') C Parameters for the response functions: FMIN,FSTEP,NF,TIMUL,FREF. NCOM=MCOM CALL RSEP3T('RF',FILOUT,'rf.out') C ELSE C C Reading the source time function and its Hilbert transform: OPEN(LU1,FILE=FILSIG,STATUS='OLD') CALL RGSE2(LU1,STAT,CHAN, * I,XS1,XS2,XS3,TSTAR0,TSTEP0,NSEIS4,MSEIS,SEIS4) CALL RGSE2(LU1,STAT,CHAN, * I,XS1,XS2,XS3,TSTARH,TSTEP ,NSEIS5,MSEIS,SEIS5) CLOSE(LU1) IF(TSTEP0.NE.TSTEP) THEN C GREENSS-02 CALL ERROR('GREENSS-02: Different time steps.') END IF FILRM=' ' IF(FILRM.NE.' ') THEN C GREENSS-04 CALL ERROR('GREENSS-04: Layer stack in time domain') C The response of the transmission through a stack of fine C horizontal layers at the receiver locations cannot be C calculated in the time domain. END IF NCOM=3 CALL RSEP3T('SS',FILOUT,'ss.gse') C END IF C C....................................................................... C C Opening input file with the Green function: CALL RSEP3T('GREEN' ,FILE1 ,'green.out' ) OPEN(LU1,FILE=FILE1,STATUS='OLD') READ(LU1,*) (TEXT,I=1,20) C C Opening the output file: OPEN(LU3,FILE=FILOUT) IF(NCOM.NE.MCOM) THEN C Time domain CALL WGSE1(LU3,' ') ELSE C Frequency domain IF(LGSE) THEN CALL WGSE1(LU3,' ') ELSE LWRITE=.TRUE. END IF END IF C C Removed from the input data in version 5.20: C CALL RSEP3I('KOMP',KOMP,0) KOMP=0 C KOMP=integer... The default value is recommended. C 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: KOMP=0 (recommended) C C Loop over rays: TXTREC='$' TXTRM='$' 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 C Reading one elementary Green function: IF(NCOM.NE.MCOM) THEN C Time domain NGREEN=33 ELSE C Frequency domain NGREEN=14+18*NF END IF DO 32 I=1,NGREEN GREEN(I)=0. 32 CONTINUE GREEN(33)=UNDEF READ(LU1,*) TXTREC,TXTSRC,(GREEN(I),I=1,NGREEN) IF(GREEN(33).EQ.UNDEF) THEN C Frequency-independent Green function: NGREEN=32 ELSE C Frequency-dependent Green function: IF(NCOM.NE.MCOM) THEN C Time domain calculation C GREENSS-05 CALL ERROR('GREENSS-05: Frequency dependent Green function') C Frequency-dependent Green function has been identified in C the input file during the calculation of the seismograms in C the time domain. END IF NGREEN=14+18*NF END IF IF(TXTOLD.NE.'$'.AND.TXTREC.NE.TXTOLD) THEN C Writing the previous seismogram: C -------------------------------- IF(NCOM.NE.MCOM) THEN C Time domain IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN CALL WGSE2(LU3,TXTOLD,' ',1,XR1,XR2,XR3, * TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS1) END IF IF(KOMP.LE.0.OR.KOMP.EQ.2) THEN CALL WGSE2(LU3,TXTOLD,' ',2,XR1,XR2,XR3, * TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS2) END IF IF(KOMP.LE.0.OR.KOMP.EQ.3) THEN CALL WGSE2(LU3,TXTOLD,' ',3,XR1,XR2,XR3, * TSTART,TSTEP,MIN0(MSEIS,NSEIS),SEIS3) END IF IF(NSEIS.GT.MSEIS) THEN C GREENSS-51 WRITE(TEXT,'(A,I12,2A)') 'GREENSS-51:',NSEIS, * ' non-zero samples at receiver ',TXTOLD CALL WARN(TEXT) END IF ELSE C Frequency domain IF(LGSE) THEN IF(KOMP.LE.0.OR.KOMP.EQ.1) THEN CALL WGSE2(LU3,TXTOLD,' ', 1,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS1) CALL WGSE2(LU3,TXTOLD,' ',11,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS4) END IF IF(KOMP.LE.0.OR.KOMP.EQ.2) THEN CALL WGSE2(LU3,TXTOLD,' ', 2,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS2) CALL WGSE2(LU3,TXTOLD,' ',12,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS5) END IF IF(KOMP.LE.0.OR.KOMP.EQ.3) THEN CALL WGSE2(LU3,TXTOLD,' ', 3,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS3) CALL WGSE2(LU3,TXTOLD,' ',13,XR1,XR2,XR3, * FMIN,FSTEP,NF,SEIS6) END IF ELSE IF(LWRITE)THEN WRITE(LU3,'(A)') '/' WRITE(LU3,'(3(F7.3,1X),A)') VCOM(1),VCOM(2),VCOM(3),'/' WRITE(LU3,'(2(E11.5,1X),I8,1X,A)') FMIN,FSTEP,NF,'/' LWRITE=.FALSE. END IF AMAX=0. DO 33 I=1,NF AMAX=AMAX1(ABS(SEIS1(I)),ABS(SEIS2(I)),ABS(SEIS3(I)), * ABS(SEIS4(I)),ABS(SEIS5(I)),ABS(SEIS6(I)), * AMAX) 33 CONTINUE IF(AMAX.LE.0.) THEN IF(VCOM(4).GT.0.) THEN VCOM(5)=0. ELSE VCOM(5)=VCOM(4)-1. END IF END IF IF(XR1.LT.-99.999.OR.999.999.LT.XR1.OR. * XR2.LT.-99.999.OR.999.999.LT.XR2.OR. * XR3.LT.-99.999.OR.999.999.LT.XR3) THEN WRITE(LU3,'(3A,3(F7.0,1X),3(E11.5,1X),A)') * '''',TXTOLD,''' ',XR1,XR2,XR3,VCOM(4),VCOM(5),AMAX,'/' ELSE WRITE(LU3,'(3A,3(F7.3,1X),3(E11.5,1X),A)') * '''',TXTOLD,''' ',XR1,XR2,XR3,VCOM(4),VCOM(5),AMAX,'/' END IF IF(VCOM(4).LE.VCOM(5)) THEN DO 34 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(LU3,'(12I6)') I01,I02,I03,I04,I05,I06, * I07,I08,I09,I10,I11,I12 ELSE WRITE(LU3,'(12I6)') I01,I02,I03,I04,I05,I06 END IF 34 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 ------------- IF(NCOM.NE.MCOM) THEN C Time domain ISHIFT=INT(GREEN(1)/TSTEP) TSTART=TSTAR0+FLOAT(ISHIFT)*TSTEP NSEIS=0 ELSE C Frequency domain VCOM(4)= 999999. VCOM(5)=-999999. DO 41 I=1,NF SEIS1(I)=0. SEIS2(I)=0. SEIS3(I)=0. SEIS4(I)=0. SEIS5(I)=0. SEIS6(I)=0. 41 CONTINUE END IF C C Reading fine layered structure at the receiver location IF(FILRM.NE.' ') THEN IF(TXTRM.EQ.'$') THEN C Opening file RMATRIX with fine layered structures at the C receiver sites: OPEN(11 ,STATUS='SCRATCH') OPEN(LU2,FILE=FILRM,STATUS='OLD') END IF IF(TXTRM.NE.' ') THEN C Searching for the receiver in file RMATRIX 43 CONTINUE TXTRM=' ' READ(LU2,*,ERR=44,END=44) TXTRM CALL INPUT(LU2,NLAY,IFREE) IF(TXTRM.NE.' '.AND.TXTRM.NE.TXTREC) GO TO 43 GO TO 45 C .......................................................... C The receiver has not been found 44 CONTINUE C GREENSS-06 CALL ERROR('GREENSS-06: Receiver is not in file RMATRIX') C The receiver has not been found in file C RMATRIX C .......................................................... 45 CONTINUE END IF END IF END IF C DO 49 I=15,NGREEN GREEN(I)=GREEN(I)/1000000. 49 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 Time domain 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=MIN0(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=MIN0(MSEIS,-ISHIFT),1,-1 SEIS1(I)=0. SEIS2(I)=0. SEIS3(I)=0. 52 CONTINUE ISHIFT=0 END IF DO 53 I=NSEIS+1,MIN0(NSEIS4+ISHIFT+1,MSEIS) SEIS1(I)=0. SEIS2(I)=0. SEIS3(I)=0. 53 CONTINUE DO 54 I=1,MIN0(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) 54 CONTINUE ISHIFT=ISHIFT+1 NSEIS=MAX0(NSEIS,NSEIS4+ISHIFT) DO 55 I=1,MIN0(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) 55 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 DO 63 I=NSEIS+1,MIN0(NSEIS5+ISHIFT+1,MSEIS) SEIS1(I)=0. SEIS2(I)=0. SEIS3(I)=0. 63 CONTINUE NSEIS=MAX0(NSEIS,NSEIS5+ISHIFT) DO 64 I=1,MIN0(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) 64 CONTINUE ISHIFT=ISHIFT+1 NSEIS=MAX0(NSEIS,NSEIS5+ISHIFT) DO 65 I=1,MIN0(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) 65 CONTINUE ELSE C Frequency domain 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.LE.0..OR.TI.LE.0.) THEN TR=GREEN(1) ELSE IF(F.GT.0.) THEN TR=GREEN(1)-TI*ALOG(F/FREF)/1.5707963 ELSE TR=GREEN(1)-TI*ALOG(FSTEP/FREF)/1.5707963 END IF END IF C C Frequency-dependent green function IF(NGREEN.GT.32) THEN J=18*(I-1) AR1=GREEN(J+15)*FR1+GREEN(J+21)*FR2+GREEN(J+27)*FR3 * -GREEN(J+16)*FI1-GREEN(J+22)*FI2-GREEN(J+28)*FI3 AR2=GREEN(J+17)*FR1+GREEN(J+23)*FR2+GREEN(J+29)*FR3 * -GREEN(J+18)*FI1-GREEN(J+24)*FI2-GREEN(J+30)*FI3 AR3=GREEN(J+19)*FR1+GREEN(J+25)*FR2+GREEN(J+31)*FR3 * -GREEN(J+20)*FI1-GREEN(J+26)*FI2-GREEN(J+32)*FI3 AI1=GREEN(J+15)*FI1+GREEN(J+21)*FI2+GREEN(J+27)*FI3 * +GREEN(J+16)*FR1+GREEN(J+22)*FR2+GREEN(J+28)*FR3 AI2=GREEN(J+17)*FI1+GREEN(J+23)*FI2+GREEN(J+29)*FI3 * +GREEN(J+18)*FR1+GREEN(J+24)*FR2+GREEN(J+30)*FR3 AI3=GREEN(J+19)*FI1+GREEN(J+25)*FI2+GREEN(J+31)*FI3 * +GREEN(J+20)*FR1+GREEN(J+26)*FR2+GREEN(J+32)*FR3 END IF C C Adding the response of the fine layered structure at the C receiver location to the complex-valued amplitude: IF(FILRM.EQ.' ') THEN BR1=AR1 BI1=AI1 BR2=AR2 BI2=AI2 BR3=AR3 BI3=AI3 ELSE CALL RM(NLAY,IFREE,I, * OMEGA,TR,GREEN(9),GREEN(10),GREEN(11), * AR1,AI1,AR2,AI2,AR3,AI3,BR1,BI1,BR2,BI2,BR3,BI3) END IF C C Complex-valued phase factor (AR,AI) TI= EXP(-OMEGA*TI) AR=TI*COS(OMEGA*TR) AI=TI*SIN(OMEGA*TR) C C Contribution of the ray to the response function: SEIS1(I)=SEIS1(I)+BR1*AR-BI1*AI SEIS2(I)=SEIS2(I)+BR2*AR-BI2*AI SEIS3(I)=SEIS3(I)+BR3*AR-BI3*AI SEIS4(I)=SEIS4(I)+BI1*AR+BR1*AI SEIS5(I)=SEIS5(I)+BI2*AR+BR2*AI SEIS6(I)=SEIS6(I)+BI3*AR+BR3*AI 69 CONTINUE END IF END IF GO TO 30 C 80 CONTINUE IF(NCOM.NE.MCOM) THEN C Time domain CALL WGSE3(LU3) ELSE C Frequency domain IF(LGSE)THEN CALL WGSE3(LU3) ELSE WRITE(LU3,'(A)') '/' END IF END IF CLOSE(LU3) IF(FILRM.NE.' ') THEN CLOSE(LU2) CLOSE(11) END IF CLOSE(LU1) WRITE(*,'(A)') '+GREENSS: 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 Response of fine layers: C Include just one set of the following files: C (a) Version GREENSS: No response of a stack of fine horizontal layers: INCLUDE 'rmnul.for' C rmnul.for C (b) Version GREENRM: Stacks of fine horizontal layers at receivers: * INCLUDE 'rm.for' C rm.for * INCLUDE 'rmatrix1.f' C rmatrix1.f * INCLUDE 'rmatrix2.f' C rmatrix2.f * INCLUDE 'isotroc.f' C isotroc.f * INCLUDE 'aniso6cg.f' C aniso6cg.f * INCLUDE 'aniso6c1.f' C aniso6c1.f * INCLUDE 'aniso3c.f' C aniso3c.f * INCLUDE 'deigv.f' C deigv.f * INCLUDE 'cg.f' C cg.f C C======================================================================= C