C
C Program GRDAMPKM to read the grided amplitudes of the Green function C generated be program 'mtt.for' and to prepare the gridded amplitudes C for common-source Kirchhoff migration of vectorial three-component C seismograms by program 'grdmigr.for' C C Version: 8.10 C Date: 2022, December 15 C C Coded by: Ludek Klimes C Department of Geophysics, Charles University, 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 files: C GRDAKI='string'... Name of the output file AMPKI of program mtt.for C with the grid values of the amplitude modified for use in C the Kirchoff integral, converted from a multivalued grid C to a singlevalued grid by program mgrd.for. C No default. The filename must be specified. C GRDAR11='string',GRDAI11='string',GRDAR21='string', C GRDAI21='string',GRDAR31='string',GRDAI31='string', C GRDAR12='string',GRDAI12='string',GRDAR22='string', C GRDAI22='string',GRDAR32='string',GRDAI32='string', C GRDAR13='string',GRDAI13='string',GRDAR23='string', C GRDAI23='string',GRDAR33='string',GRDAI33='string'... C Names of the output files AMPR11,AMPI11,AMPR21,...,AMPI33 C of program mtt.for C with the grid values of the real and imaginary parts of C the tensorial amplitude of the Green function, converted C from a multivalued grid to a singlevalued grid by program C mgrd.for. C No default. The filenames must be specified. C Names of the output files: C GRDAK1='string',GRDAK2='string',GRDAK3='string'... Singlevalued C grid files with three components of vectorial amplitude C corresponding to the three seismogram components C for the Kirchhoff integral. C No default. The filenames must be specified. C Data specifying the dimensions of the input and output grids: C N1=positive integer... Number of gridpoints along the X1 axis C (inner loop). C Default: N1=1 C N2=positive integer... Number of gridpoints along the X2 axis C (intermediate loop). C Default: N2=1 C N3=positive integer... Number of gridpoints along the X3 axis C (outer loop). C Default: N3=1 C Numerical parameters specifying the output quantities: C ORIKI1=real, ORIKI2=real, ORIKI3=real... Output amplitude vector C (GRDAK1,GRDAK2,GRDAK3) is determined up to the sign. C The sign specifies the sign in the migrated image. C The sign is determined so that the scalar product of C vectors (GRDAK1,GRDAK2,GRDAK3) and (ORIKI1,ORIKI2,ORIKI3) C is positive. The three-component migration then should C resemble the one-component migration of the seismogram C component along vector (ORIKI1,ORIKI2,ORIKI3). C Default: ORIKI1=0., ORIKI2=0., ORIKI3=1. (sign rougly C corresponding to the scalar migration of the third C seismogram component) C C======================================================================= C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C Subroutines and external functions required: EXTERNAL ERROR,RSEP1,RSEP3T,RSEP3I,RSEP3R,RARRAY,WARRAY C ERROR... File error.for/A>. C RSEP1,RSEP3T,RSEP3I,RSEP3R... File sep.for. C RARRAY,WARRAY... File forms.for. C C....................................................................... C C Filenames and parameters: CHARACTER*80 FSEP,FGRD1,FGRD2,FGRD3 INTEGER LU1 PARAMETER (LU1=1) C C Input data: INTEGER N1,N2,N3 REAL ORIKI1,ORIKI2,ORIKI3 C C Memory allocation: INTEGER IB11,IB12,IB22,IB13,IB23,IB33,IAMP1,IAMP2,IAMP3,NRAM C C Other variables: INTEGER JB11,JB12,JB22,JB13,JB23,JB33,JAMP1,JAMP2,JAMP3 INTEGER N123,I,ICOMP REAL B(6),E(9),A1,A2,A3,A C C----------------------------------------------------------------------- C C Reading input SEP parameter file: WRITE(*,'(A)') '+GRDMIGRV: Enter input filename: ' FSEP=' ' READ(*,*) FSEP IF (FSEP.EQ.' ') THEN C GRDAMPKM-01 CALL ERROR('GRDAMPKM-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 CALL RSEP1(LU1,FSEP) WRITE(*,'(A)') '+GRDAMPKM: Working ... ' C C Reading input parameters from the SEP file: CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) N123=N1*N2*N3 CALL RSEP3R('ORIKI1',ORIKI1,0.) CALL RSEP3R('ORIKI2',ORIKI2,0.) CALL RSEP3R('ORIKI3',ORIKI3,1.) C C Memory allocation: IB11=0 IB12=IB11+N123 IB22=IB12+N123 IB13=IB22+N123 IB23=IB13+N123 IB33=IB23+N123 IAMP1=IB33+N123 IAMP2=IAMP1+N123 IAMP3=IAMP2+N123 NRAM=IAMP3+N123 IF(NRAM.GT.MRAM) THEN C GRDAMPKM-02 CALL ERROR('GRDAMPKM-02: Too small array RAM(MRAM)') C Array RAM(MRAM) allocated in include file 'ram.inc' is too small C to contain both input and output grids. C Dimension MRAM must be at least 9*N1*N2*N3. C You may wish to increase the dimension MRAM in file C ram.inc. END IF C C Initializing the 3*3 matrices: JB11=IB11 JB12=IB12 JB22=IB22 JB13=IB13 JB23=IB23 JB33=IB33 DO 11 I=1,N123 JB11=JB11+1 JB12=JB12+1 JB22=JB22+1 JB13=JB13+1 JB23=JB23+1 JB33=JB33+1 RAM(JB11)=0. RAM(JB12)=0. RAM(JB22)=0. RAM(JB13)=0. RAM(JB23)=0. RAM(JB33)=0. 11 CONTINUE C C Accumulating the 3*3 matrices of scalar products: DO 22 ICOMP=1,6 IF(ICOMP.EQ.1) THEN CALL RSEP3T('GRDAR11',FGRD1,' ') CALL RSEP3T('GRDAR12',FGRD2,' ') CALL RSEP3T('GRDAR13',FGRD3,' ') ELSE IF(ICOMP.EQ.2) THEN CALL RSEP3T('GRDAI11',FGRD1,' ') CALL RSEP3T('GRDAI12',FGRD2,' ') CALL RSEP3T('GRDAI13',FGRD3,' ') ELSE IF(ICOMP.EQ.3) THEN CALL RSEP3T('GRDAR21',FGRD1,' ') CALL RSEP3T('GRDAR22',FGRD2,' ') CALL RSEP3T('GRDAR23',FGRD3,' ') ELSE IF(ICOMP.EQ.4) THEN CALL RSEP3T('GRDAI21',FGRD1,' ') CALL RSEP3T('GRDAI22',FGRD2,' ') CALL RSEP3T('GRDAI23',FGRD3,' ') ELSE IF(ICOMP.EQ.5) THEN CALL RSEP3T('GRDAR31',FGRD1,' ') CALL RSEP3T('GRDAR32',FGRD2,' ') CALL RSEP3T('GRDAR33',FGRD3,' ') ELSE CALL RSEP3T('GRDAI31',FGRD1,' ') CALL RSEP3T('GRDAI32',FGRD2,' ') CALL RSEP3T('GRDAI33',FGRD3,' ') END IF IF(FGRD1.EQ.' '.OR.FGRD2.EQ.' '.OR.FGRD3.EQ.' ') THEN C GRDAMPKM-03 CALL ERROR('GRDAMPKM-03: Missing Green function file') C At least one of 18 input files with the grid values of the C real and imaginary parts of the tensorial amplitude of the C Green function has blank name. C The filenames have no default and must be specified. END IF CALL RARRAY(LU1,FGRD1,.TRUE.,0.,N123,RAM(IAMP1+1)) CALL RARRAY(LU1,FGRD2,.TRUE.,0.,N123,RAM(IAMP2+1)) CALL RARRAY(LU1,FGRD3,.TRUE.,0.,N123,RAM(IAMP3+1)) JB11=IB11 JB12=IB12 JB22=IB22 JB13=IB13 JB23=IB23 JB33=IB33 JAMP1=IAMP1 JAMP2=IAMP2 JAMP3=IAMP3 DO 21 I=1,N123 JB11=JB11+1 JB12=JB12+1 JB22=JB22+1 JB13=JB13+1 JB23=JB23+1 JB33=JB33+1 JAMP1=JAMP1+1 JAMP2=JAMP2+1 JAMP3=JAMP3+1 A1=RAM(JAMP1) A2=RAM(JAMP2) A3=RAM(JAMP3) RAM(JB11)=RAM(JB11)+A1*A1 RAM(JB12)=RAM(JB12)+A1*A2 RAM(JB22)=RAM(JB22)+A2*A2 RAM(JB13)=RAM(JB13)+A1*A3 RAM(JB23)=RAM(JB23)+A2*A3 RAM(JB33)=RAM(JB33)+A3*A3 21 CONTINUE 22 CONTINUE C C Determining the eigenvector of the largest eigenvalue: WRITE(*,'(A)') '+GRDAMPKM: Calculating ... ' JB11=IB11 JB12=IB12 JB22=IB22 JB13=IB13 JB23=IB23 JB33=IB33 JAMP1=IAMP1 JAMP2=IAMP2 JAMP3=IAMP3 DO 31 I=1,N123 JB11=JB11+1 JB12=JB12+1 JB22=JB22+1 JB13=JB13+1 JB23=JB23+1 JB33=JB33+1 JAMP1=JAMP1+1 JAMP2=JAMP2+1 JAMP3=JAMP3+1 A=RAM(JB11)+RAM(JB22)+RAM(JB33) IF(A.GT.0.) THEN B(1)=RAM(JB11)/A B(2)=RAM(JB12)/A B(3)=RAM(JB22)/A B(4)=RAM(JB13)/A B(5)=RAM(JB23)/A B(6)=RAM(JB33)/A CALL EIGEN(B,E,3,0) A1=E(1) A2=E(2) A3=E(3) IF(A1*ORIKI1+A2*ORIKI2+A3*ORIKI3.LT.0.) THEN A1=-A1 A2=-A2 A3=-A3 END IF A=SQRT(A1*A1+A2*A2+A3*A3) IF(A.GT.0.) THEN RAM(JAMP1)=A1/A RAM(JAMP2)=A2/A RAM(JAMP3)=A3/A ELSE RAM(JAMP1)=0. RAM(JAMP2)=0. RAM(JAMP3)=0. END IF ELSE RAM(JAMP1)=0. RAM(JAMP2)=0. RAM(JAMP3)=0. END IF 31 CONTINUE C C Multiplying the eigenvector by the amplitude: CALL RSEP3T('GRDAKI',FGRD1,' ') CALL RARRAY(LU1,FGRD1,.TRUE.,0.,N123,RAM(IB33+1)) IF(FGRD1.EQ.' ') THEN C GRDAMPKM-04 CALL ERROR('GRDAMPKM-04: File GRDAKI not given') C Filename GRDAKI with the input grid values is blank. C The filename has no default and must be specified. END IF JB33=IB33 JAMP1=IAMP1 JAMP2=IAMP2 JAMP3=IAMP3 DO 41 I=1,N123 JB33=JB33+1 JAMP1=JAMP1+1 JAMP2=JAMP2+1 JAMP3=JAMP3+1 A=RAM(JB33) RAM(JAMP1)=RAM(JAMP1)*A RAM(JAMP2)=RAM(JAMP2)*A RAM(JAMP3)=RAM(JAMP3)*A 41 CONTINUE C C Writing the output vectorial amplitude: CALL RSEP3T('GRDAK1',FGRD1,' ') CALL RSEP3T('GRDAK2',FGRD2,' ') CALL RSEP3T('GRDAK3',FGRD3,' ') CALL WARRAY(LU1,FGRD1,.FALSE.,0.,.FALSE.,0.,N123,RAM(IAMP1+1)) CALL WARRAY(LU1,FGRD2,.FALSE.,0.,.FALSE.,0.,N123,RAM(IAMP2+1)) CALL WARRAY(LU1,FGRD3,.FALSE.,0.,.FALSE.,0.,N123,RAM(IAMP3+1)) WRITE(*,'(A)') '+GRDAMPKM: Done. ' STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for INCLUDE 'eigen.for' C eigen.for C C======================================================================= C