C                                                            <HTML><PRE>
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                                                    <A NAME="data"></A>
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                                                     <A NAME="SEP"></A>
C Input data file 'SEP':
C     File 'SEP' has the form of the <A HREF="../forms/sep.htm">SEP</A>
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 <A HREF="../crt/mtt.for">mtt.for</A>
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 <A HREF="../forms/mgrd.for">mgrd.for</A>.
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 <A HREF="../crt/mtt.for">mtt.for</A>
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             <A HREF="../forms/mgrd.for">mgrd.for</A>.
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     <A HREF="ram.inc" TYPE="text/html">ram.inc</A>
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3T,RSEP3I,RSEP3R,RARRAY,WARRAY
C     ERROR... File <A HREF="error.for" TYPE="text/html">error.for/A>.
C     RSEP1,RSEP3T,RSEP3I,RSEP3R... File <A HREF="sep.for" TYPE="text/html">sep.for</A>.
C     RARRAY,WARRAY... File <A HREF="forms.for" TYPE="text/html">forms.for</A>.
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       <A NAME="error01"></A><FONT COLOR="RED">GRDAMPKM-01</FONT>
        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       <A NAME="error02"></A><FONT COLOR="RED">GRDAMPKM-02</FONT>
        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       <A HREF="ram.inc" TYPE="text/html">ram.inc</A>.
      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         <A NAME="error03"></A><FONT COLOR="RED">GRDAMPKM-03</FONT>
          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       <A NAME="error04"></A><FONT COLOR="RED">GRDAMPKM-04</FONT>
        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     <A HREF="error.for" TYPE="text/html">error.for</A>
      INCLUDE 'sep.for'
C     <A HREF="sep.for" TYPE="text/html">sep.for</A>
      INCLUDE 'forms.for'
C     <A HREF="forms.for" TYPE="text/html">forms.for</A>
      INCLUDE 'length.for'
C     <A HREF="length.for" TYPE="text/html">length.for</A>
      INCLUDE 'eigen.for'
C     <A HREF="eigen.for" TYPE="text/html">eigen.for</A>
C
C=======================================================================
C                                                                 </PRE>
