C
C Program GRDCROS to calculate autocorrelations and crosscorrelations
C of given grid values
C
C Version: 6.00
C Date: 2006, March 3
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                                                   
C Names of input and output files:
C     GRD1='string',GRD2='string'... Names of two input ASCII files with
C             the grid values.
C             Default: GRD1='grd1.out', GRD2=GRD1
C     GRD3='string'... Name of the output ASCII file containing the
C             crosscorrelations of grid values read from files GRD1 and
C             GRD2.
C             Default: GRD3='grd3.out'
C     GRD4='string'... Name of the output ASCII file containing the
C             numbers of points used to determinine individual
C             crosscorrelations written to file GRD3, multiplied by the
C             given cosine window.
C             Not written if GRD4=' '.
C             Default: GRD4=' '
C     For general description of the files with gridded data refer
C     to file forms.htm.
C                                                   
C Data specifying dimensions of the input grid:
C     D1=real... Grid interval in the direction of the first coordinate
C             axis in files GRD1, GRD2 and GRD3.
C             Default: D1=1.
C     D2=real... Grid interval in the direction of the second coordinate
C             axis in files GRD1, GRD2 and GRD3.
C             Default: D2=1.
C     D3=real... Grid interval in the direction of the third coordinate
C             axis in files GRD1, GRD2 and GRD3.
C             Default: D3=1.
C     N11=positive integer... Number of gridpoints along the X1 axis in
C             input file GRD1.
C             Default: N11=1
C     N21=positive integer... Number of gridpoints along the X2 axis in
C             input file GRD1.
C             Default: N21=1
C     N31=positive integer... Number of gridpoints along the X3 axis in
C             input file GRD1.
C             Default: N31=1
C     O11=real... First coordinate of the grid origin (first point of
C             the grid) in input file GRD1.
C             Default: O11=0.
C     O21=real... Second coordinate of the grid origin in input file
C             GRD1.
C             Default: O21=0.
C     O31=real... Third coordinate of the grid origin in input file
C             GRD1.
C             Default: O31=0.
C     N12=positive integer... Number of gridpoints along the X1 axis in
C             input file GRD2.
C             Default: N12=N11
C     N22=positive integer... Number of gridpoints along the X2 axis in
C             input file GRD2.
C             Default: N22=N21
C     N32=positive integer... Number of gridpoints along the X3 axis in
C             input file GRD2.
C             Default: N32=N31
C     O12=real... First coordinate of the grid origin (first point of
C             the grid) in input file GRD2.
C             Default: O12=O11
C     O22=real... Second coordinate of the grid origin in input file
C             GRD2.
C             Default: O22=O21
C     O32=real... Third coordinate of the grid origin in input file
C             GRD2.
C             Default: O32=O31
C                                                   
C Data specifying dimensions of the output grid:
C     N1=positive integer... Number of crosscorrelations (discretized
C             with step D1) in the direction of the X1 axis.
C             N1 cannot be greater than N11+N12-1.
C             Default: N1=N11+N12-1
C     N2=positive integer... Number of crosscorrelations (discretized
C             with step D2) in the direction of the X2 axis.
C             N2 cannot be greater than N21+N22-1.
C             Default: N2=N21+N22-1
C     N3=positive integer... Number of crosscorrelations (discretized
C             with step D3) in the direction of the X3 axis.
C             N3 cannot be greater than N31+N32-1.
C             Default: N3=N31+N32-1
C     O1=real... Shift (in the direction of the X1 axis) of GRD2 with
C             respect to GRD1 corresponding to the first output
C             crosscorrelation.
C             (O1+O12-O11) should equal an integer multiple of D1.
C             O1 cannot be smaller than O11-O12-(N12-1)*D1 and cannot
C             be greater than O11+(N11-1)*D1-O12-(N1-1)*D1.
C             Default: O1=O11-O12-(N12-1)*D1
C     O2=real... Shift (in the direction of the X2 axis) of GRD2 with
C             respect to GRD1 corresponding to the first output
C             crosscorrelation.
C             (O2+O22-O21) should equal an integer multiple of D2.
C             O2 cannot be smaller than O21-O22-(N22-1)*D2 and cannot
C             be greater than O21+(N21-1)*D2-O22-(N2-1)*D2.
C             Default: O2=O21-O22-(N22-1)*D2
C     O3=real... Shift (in the direction of the X3 axis) of GRD2 with
C             respect to GRD1 corresponding to the first output
C             crosscorrelation.
C             (O3+O32-O31) should equal an integer multiple of D3.
C             O3 cannot be smaller than O31-O32-(N32-1)*D3 and cannot
C             be greater than O31+(N31-1)*D3-O32-(N3-1)*D3.
C             Default: O3=O31-O32-(N32-1)*D3
C                                                   
C Data specifying the cosine window:
C     X1MIN=real, X1LOW=real, X1HIGH=real, X1MAX=real... Parameters of
C             the cosine window in the X1 direction to be applied to
C             the values in output file GRD4.
C             If X1MIN.LT.X1MAX, the window is zero for X1 smaller than
C             X1MIN or greater than X1MAX, and the window equals 1
C             between X1LOW and X1HIGH.
C             If X1MAX.LT.X1MIN, the window is zero for X1 greater than
C             X1MAX and smaller than X1MIN.  Then the window equals 1
C             for X1 smaller than X1HIGH or greater than X1LOW.
C             Between X1MIN and X1LOW, cosine tapering
C               ( 0.5-0.5*cos(pi*(X1-X1MIN)/(X1LOW-X1MIN)) )**KEXP
C             is used.
C             Between X1MIN and X1LOW, cosine tapering
C               ( 0.5-0.5*cos(pi*(X1-X1MAX)/(X1HIGH-X1MAX)) )**KEXP
C             is used.
C             Default: X1MIN=0., X1LOW=0., X1HIGH=0., X1MAX=0.
C     KEXP=real... Exponent controlling the cosine window.
C             Usually need not be specified because the default is the
C             most common option.
C             Default: KEXP=1 (mostly sufficient)
C                                                   
C Optional parameters specifying the form of the real-valued quantities
C written to the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C.......................................................................
C
C     External functions and subroutines:
      EXTERNAL ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,UARRAY,WARRAY,RARRAY
      REAL UARRAY
C
C     Filenames and parameters:
      CHARACTER*80 FILE1,FILE2,FILE3,FILE4
      INTEGER LU
      REAL UNDEF
      PARAMETER (LU=1)
C     Input data:
      INTEGER N1,N2,N3,N11,N21,N31,N12,N22,N32
      REAL D1,D2,D3,O1,O2,O3,O11,O21,O31,O12,O22,O32
      REAL X1MIN,X1LOW,X1HIGH,X1MAX
      REAL X2MIN,X2LOW,X2HIGH,X2MAX
      REAL X3MIN,X3LOW,X3HIGH,X3MAX
C     Other variables:
      INTEGER I1,I2,I3,I01,I11,I21,I31,I02,I12,I22,I32
      INTEGER I11MIN,I21MIN,I31MIN,I11MAX,I21MAX,I31MAX
      INTEGER NO1,NO2,NO3,IGRD1,IGRD2,IGRD3,IGRD4,I,IAUX
      REAL X1,X2,X3,X1WIN,X2WIN,X3WIN,AUX
C
      UNDEF=UARRAY()
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+GRDCROS: Enter input filename: '
      FILE1=' '
      READ(*,*) FILE1
      WRITE(*,'(A)') '+GRDCROS: Working ...           '
C
C     Reading all data from the SEP file into the memory:
      IF (FILE1.EQ.' ') THEN
C       GRDCROS-01
        CALL ERROR('GRDCROS-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(LU,FILE1)
C
C     Reading input parameters from the SEP file:
      CALL RSEP3T('GRD1',FILE1,'grd1.out')
      CALL RSEP3T('GRD2',FILE2,FILE2)
      CALL RSEP3T('GRD3',FILE3,'grd3.out')
      CALL RSEP3T('GRD4',FILE4,' ')
C
C     Reading grid dimensions:
C     Input grids:
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      CALL RSEP3I('N11',N11,1)
      CALL RSEP3I('N21',N21,1)
      CALL RSEP3I('N31',N31,1)
      CALL RSEP3R('O11',O11,0.)
      CALL RSEP3R('O21',O21,0.)
      CALL RSEP3R('O31',O31,0.)
      CALL RSEP3I('N12',N12,N11)
      CALL RSEP3I('N22',N22,N21)
      CALL RSEP3I('N32',N32,N31)
      CALL RSEP3R('O12',O12,O11)
      CALL RSEP3R('O22',O22,O21)
      CALL RSEP3R('O32',O32,O31)
C     Output grid:
      CALL RSEP3I('N1',N1,N11+N12-1)
      CALL RSEP3I('N2',N2,N21+N22-1)
      CALL RSEP3I('N3',N3,N31+N32-1)
      CALL RSEP3R('O1',O1,O11-O12-(N12-1)*D1)
      CALL RSEP3R('O2',O2,O21-O22-(N22-1)*D2)
      CALL RSEP3R('O3',O3,O31-O32-(N32-1)*D3)
      NO1=NINT((O1+O12-O11)/D1)
      NO2=NINT((O2+O22-O21)/D2)
      NO3=NINT((O3+O32-O31)/D3)
      WRITE(*,'(3(A,I4),3(A,F8.3))')
     *                           '+GRDCROS: N1=',N1,' N2=',N2,' N3=',N3,
     *                                    ' O1=',O1,' O2=',O2,' O3=',O3
      WRITE(*,'(A)') ' GRDCROS: Working ...           '
      IF(NO1.LT.1-N12.OR.NO1.GT.N11-N1.OR.
     *   NO2.LT.1-N22.OR.NO2.GT.N21-N2.OR.
     *   NO3.LT.1-N32.OR.NO3.GT.N31-N3) THEN
C       GRDCROS-02
        CALL ERROR('GRDCROS-02: Incorrect N1, N2, N3, O1, O2 or O3')
C       N1, N2, N3, O1, O2 or O3 do not satisfy the inequalities set
C       in the description of the input data.
      END IF
C
C     Reading the parameters of the cosine window:
      CALL RSEP3R('X1MIN' ,X1MIN ,0.)
      CALL RSEP3R('X1LOW' ,X1LOW ,0.)
      CALL RSEP3R('X1HIGH',X1HIGH,0.)
      CALL RSEP3R('X1MAX' ,X1MAX ,0.)
      CALL RSEP3R('X2MIN' ,X2MIN ,0.)
      CALL RSEP3R('X2LOW' ,X2LOW ,0.)
      CALL RSEP3R('X2HIGH',X2HIGH,0.)
      CALL RSEP3R('X2MAX' ,X2MAX ,0.)
      CALL RSEP3R('X3MIN' ,X3MIN ,0.)
      CALL RSEP3R('X3LOW' ,X3LOW ,0.)
      CALL RSEP3R('X3HIGH',X3HIGH,0.)
      CALL RSEP3R('X3MAX' ,X3MAX ,0.)
C
      IGRD1=N1*N2*N3
      IGRD2=IGRD1+N11*N21*N31
      IGRD3=IGRD2+N12*N22*N32
      IGRD4=IGRD3
      IF(FILE4.NE.' ') THEN
        IGRD4=IGRD3+N1*N2*N3
      END IF
      IF(IGRD4.GT.MRAM) THEN
C       GRDCROS-03
        CALL ERROR('GRDCROS-03: Too small array RAM(MRAM)')
C       Too small array RAM(MRAM) to allocate both input and output
C       grid values.  If possible, increase dimension MRAM in include
C       file ram.inc.
      END IF
C
C     Reading input grid values:
      CALL RARRAY(LU,FILE1,'FORMATTED',.TRUE.,UNDEF,
     *                                         N11*N21*N31,RAM(IGRD1+1))
      CALL RARRAY(LU,FILE2,'FORMATTED',.TRUE.,UNDEF,
     *                                         N12*N22*N32,RAM(IGRD2+1))
C
C     Calculating crossvariances:
      I=0
      DO 23 I3=NO3,NO3+N3-1
        I31MIN=MAX0(0,I3)
        I31MAX=MIN0(N31-1,N32-1+I3)
        X3=O3+D3*FLOAT(I3-NO3)
        X3WIN=1.
        IF(X3MAX.LE.X3MIN) THEN
          X3WIN=1.-X3WIN
        END IF
        IF(X3.LT.X3MIN) THEN
          X3WIN=1.-X3WIN
        END IF
        IF(X3MAX.LE.X3) THEN
          X3WIN=1.-X3WIN
        END IF
        DO 22 I2=NO2,NO2+N2-1
          I21MIN=MAX0(0,I2)
          I21MAX=MIN0(N21-1,N22-1+I2)
          X2=O2+D2*FLOAT(I2-NO2)
          X2WIN=1.
          IF(X2MAX.LE.X2MIN) THEN
            X2WIN=1.-X2WIN
          END IF
          IF(X2.LT.X2MIN) THEN
            X2WIN=1.-X2WIN
          END IF
          IF(X2MAX.LE.X2) THEN
            X2WIN=1.-X2WIN
          END IF
          DO 21 I1=NO1,NO1+N1-1
            I11MIN=MAX0(0,I1)+1
            I11MAX=MIN0(N11-1,N12-1+I1)+1
            X1=O1+D1*FLOAT(I1-NO1)
            X1WIN=1.
            IF(X1MAX.LE.X1MIN) THEN
              X1WIN=1.-X1WIN
            END IF
            IF(X1.LT.X1MIN) THEN
              X1WIN=1.-X1WIN
            END IF
            IF(X1MAX.LE.X1) THEN
              X1WIN=1.-X1WIN
            END IF
            AUX=0.
            DO 13 I31=I31MIN,I31MAX
              I32=I31-I3
              DO 12 I21=I21MIN,I21MAX
                I22=I21-I2
                I01=IGRD1+N11*(N21*I31+I21)
                I02=IGRD2+N12*(N22*I32+I22)
                DO 11 I11=I11MIN,I11MAX
                  I12=I11-I1
                  AUX=AUX+RAM(I01+I11)*RAM(I02+I12)
   11           CONTINUE
   12         CONTINUE
   13       CONTINUE
            I=I+1
            IAUX=(I11MAX-I11MIN+1)*(I21MAX-I21MIN+1)*(I31MAX-I31MIN+1)
            RAM(I)=AUX/FLOAT(IAUX)
            IF(FILE4.NE.' ') THEN
              RAM(IGRD3+I)=FLOAT(IAUX)*X1WIN*X2WIN*X3WIN
            END IF
   21     CONTINUE
   22   CONTINUE
   23 CONTINUE
C
C     Writing output grid values:
      CALL WARRAY(LU,FILE3,'FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                                                     N1*N2*N3,RAM)
      IF(FILE4.NE.' ') THEN
        CALL WARRAY(LU,FILE4,'FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                                            N1*N2*N3,RAM(IGRD3+1))
      END IF
      WRITE(*,'(A)') '+GRDCROS: 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
C
C=======================================================================
C