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.60 C Date: 2001, September 5 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.EX.LT.0. * .AND.ABS(XK1).LT.0.1*D1 * .AND.ABS(XK2).LT.0.1*D2 * .AND.ABS(XK3).LT.0.1*D3) THEN C Nulling infinite value corresponding to zero wavenumber RAM(I4)=0. ELSE IF (EX.EQ.0.) THEN RAM(I4)=RKAPPA*EXP(-(ACORG*XK)**2/8.) 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