C 
C Program GRDCOR to compute the values of the spectral filter
C according to the formula (3.1) of the paper Klimes, L.:
C Correlation functions of random media. In: Seismic waves in 3-D
C structures, Report 6, Department of Geophysics, Charles University,
C Prague (1997).
C
C Version: 5.40
C Date: 2000, May 22
C
C Coded by Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bulant@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 Data specifying the parameters of the grid:
C     O1=real, O2=real, O3=real ... Coordinates of the origin of the
C             grid.
C             Default: O1=0.  O2=0.  O3=0.
C     D1=real... Grid interval along the X1 axis.
C             Default: D1=0.
C     D2=real... Grid interval along the X2 axis.
C             Default: D2=0.
C     D3=real... Grid interval along the X3 axis.
C             Default: D3=0.
C     N1=positive integer... Number of gridpoints along the X1 axis.
C             Default: N1=1
C     N2=positive integer... Number of gridpoints along the X2 axis.
C             Default: N2=1
C     N3=positive integer... Number of gridpoints along the X3 axis.
C             Default: N3=1
C Data for the spectral filter:
C     NDIM=positive integer ... Spatial dimension.
C             Default: NDIM equal to the dimension of the above grid.
C     KAPPA=real ... Multiplicative factor.
C             Default: KAPPA=1.
C     POWERN=real ... Exponent or index related to fractal dimension:
C             Medium is self-affine at distances L:
C                                            ACORG .LT. L .LT. ACOR
C             Reasonable values for geology: -0.5 .LT. POWERN .LT. 0.0
C             Default: POWERN=0.0
C     ACORG=nonnegative real ... Gaussian (small-scale) correlation
C             length:
C             Removes small details (smaller than ACORG).
C             Default: ACORG=0.0
C     ACOR=positive real .. Von Karman (large-scale) correlation length:
C             Suppresses large heterogeneities (larger than ACOR).
C             Default: ACOR=999999. (infinity)
C Data for the cosine filter:
C     AMIN=nonnegative real... Minimum wavelength.  All wavelengths
C             shorter than AMIN are removed.
C             Default: AMIN=0.
C     AMAX=nonnegative real... Maximum wavelength.
C             Default: AMAX=999999. (infinity)
C     ASMALL=nonnegative real... Refer to ABIG.
C             Default: ASMALL=0.
C     ABIG=nonnegative real
C             Default: ABIG=999999. (infinity)
C             The cosine filter has value 0. at wavenumbers smaller than
C             2*PI/AMAX (i.e. at wavelengths longer than AMAX),
C             rises between 2*PI/AMAX and 2*PI/ABIG,
C             equals 1 between 2*PI/ABIG and 2*PI/ASMALL,
C             tapers off between 2*PI/ASMALL and 2*PI/AMIN,
C             and is zero at wavenumbers greater than 2*PI/AMIN
C             (i.e. at wavelengths shorter than AMIN).
C Name of output formatted file with the computed values:
C     COROUT='string'
C             Default: COROUT='grdcor.out'
C     For general description of the files with gridded data refer to
C     file forms.htm.
C
C-----------------------------------------------------------------------
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,WARRAY
C     ERROR ... File error.for.
C     RSEP1,RSEP3T,RSEP3R,RSEP3I ...
C             File sep.for.
C     WARRAY ... File forms.for.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER NDIM
      REAL DIM,RKAPPA,POWERN,ACORG,ACOR,ACOR2,AMIN,ASMALL,ABIG,AMAX
      CHARACTER*80 FILSEP,FILOUT
      REAL PI
      PARAMETER (PI=3.141592653589793)
      INTEGER LU1
      PARAMETER (LU1=1)
      INTEGER N1,N2,N3,I1,I2,I3,I4
      REAL O1,O2,O3,D1,D2,D3
      REAL XK1,XK2,XK3,XK,EX,CF
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+GRDCOR: 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       GRDCOR-05
        CALL ERROR('GRDCOR-05: 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)') '+GRDCOR: Working ...           '
C
C     Reading filename of the output file:
      CALL RSEP3T('COROUT',FILOUT,'grdcor.out')
C     Reading the values describing the grid:
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3R('D1',D1,0.)
      CALL RSEP3R('D2',D2,0.)
      CALL RSEP3R('D3',D3,0.)
      IF (((D1.EQ.0.).AND.(N1.NE.1)).OR.
     *    ((D2.EQ.0.).AND.(N2.NE.1)).OR.
     *    ((D3.EQ.0.).AND.(N3.NE.1))) THEN
C       GRDCOR-01
        CALL ERROR('GRDCOR-01: Wrong input grid.')
      ENDIF
      IF ((D1.EQ.0.).AND.(N1.EQ.1)) D1=1.
      IF ((D2.EQ.0.).AND.(N2.EQ.1)) D2=1.
      IF ((D3.EQ.0.).AND.(N3.EQ.1)) D3=1.
      IF (N1*N2*N3.GT.MRAM) THEN
C       GRDCOR-02
        CALL ERROR('GRDCOR-02: Small array RAM.')
      ENDIF
C     Reading numerical constants:
      NDIM=3
      IF (N1.EQ.1) NDIM=NDIM-1
      IF (N2.EQ.1) NDIM=NDIM-1
      IF (N3.EQ.1) NDIM=NDIM-1
      IF (NDIM.EQ.0) NDIM=NDIM+1
      I1=NDIM
      CALL RSEP3I('NDIM',NDIM,I1)
      DIM=FLOAT(NDIM)
      CALL RSEP3R('KAPPA',RKAPPA,1.)
      CALL RSEP3R('POWERN',POWERN,0.)
      CALL RSEP3R('ACORG',ACORG,0.)
      IF (ACORG.LT.0.) THEN
C       GRDCOR-03
        CALL ERROR('GRDCOR-03: ACORG less than zero.')
      ENDIF
      CALL RSEP3R('ACOR',ACOR,999999.)
      IF (ACOR.LE.0.) THEN
C       GRDCOR-04
        CALL ERROR('GRDCOR-04: ACOR less than, or equal to zero.')
      ENDIF
      IF (ACOR.GT.999998.) THEN
        ACOR2=0.
      ELSE
        ACOR2=1./ACOR**2
      ENDIF
      CALL RSEP3R('AMIN',AMIN,0.)
      CALL RSEP3R('ASMALL',ASMALL,0.)
      CALL RSEP3R('ABIG',ABIG,999999.)
      CALL RSEP3R('AMAX',AMAX,999999.)
      IF (AMIN.EQ.0.) THEN
        AMIN=999999.
      ELSE
        AMIN=2.*PI/AMIN
      ENDIF
      IF (ASMALL.EQ.0.) THEN
        ASMALL=999999.
      ELSE
        ASMALL=2.*PI/ASMALL
      ENDIF
      IF (ABIG.EQ.999999.) THEN
        ABIG=0.
      ELSE
        ABIG=2.*PI/ABIG
      ENDIF
      IF (AMAX.EQ.999999.) THEN
        AMAX=0.
      ELSE
        AMAX=2.*PI/AMAX
      ENDIF
C
      EX=(-(DIM/2.+POWERN)/2.)
      I4=0
      DO 30, I3=1,N3
        DO 20, I2=1,N2
          DO 10, I1=1,N1
            I4=I4+1
            XK1=O1+(I1-1)*D1
            XK2=O2+(I2-1)*D2
            XK3=O3+(I3-1)*D3
            XK=SQRT(XK1**2+XK2**2+XK3**2)
            IF (ACOR.GT.999998..AND.ABS(XK1).LT.0.1*D1
     *                         .AND.ABS(XK2).LT.0.1*D2
     *                         .AND.ABS(XK3).LT.0.1*D3) THEN
              RAM(I4)=0.
            ELSE
              RAM(I4)=RKAPPA*EXP(-(ACORG*XK)**2/8.)*(ACOR2+XK**2)**EX
            ENDIF
C           Cosine filter:
            IF ((XK.LE.AMAX).OR.(XK.GE.AMIN)) THEN
              CF=0.
            ELSEIF ((XK.GE.ABIG).AND.(XK.LE.ASMALL)) THEN
              CF=1.
            ELSEIF ((XK.GT.AMAX).AND.(XK.LT.ABIG)) THEN
              CF=.5 - .5*COS((XK-AMAX)/(ABIG-AMAX)*PI)
            ELSEIF ((XK.GT.ASMALL).AND.(XK.LT.AMIN)) THEN
              CF=.5 + .5*COS((XK-ASMALL)/(AMIN-ASMALL)*PI)
            ENDIF
            RAM(I4)=RAM(I4)*CF
  10      CONTINUE
  20    CONTINUE
  30  CONTINUE
      IF (FILOUT.NE.' ') THEN
        CALL WARRAY(LU1,FILOUT,'FORMATTED',.FALSE.,0.,.FALSE.,0.,I4,RAM)
      ENDIF
      WRITE(*,'(A)') '+GRDCOR: 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