C 
C Program GRDSTAT to rescale gridded data to given statistic properties.
C
C Version: 5.40
C Date: 2000, February 21
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     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 to rescale the random values:
C     DSD=positive real... Desired Standard Deviation:
C             The output grid values are scaled to have standard
C             deviation DSD.
C             Default: DSD=1.
C     VMEAN=real... Desired mean value.
C             The output grid values are shifted to have the average
C             value of VMEAN.
C             Default: VMEAN=0.
C     DEVMAX=positive real... Maximum deviation from the mean value.
C             For finite DEVMAX, the grid values V with mean value VMEAN
C             and standard deviation DSD are rescaled using
C               Vnew=VMEAN+(V-VMEAN)
C                         /(1+ABS((V-VMEAN)/DEVMAX)**DEVEXP)**(1/DEVEXP)
C             This rescaling does not influence values close to mean
C             VMEAN, especially for larger exponents DEVEXP.
C             For DEVEXP=999999. (infinity), rescaling does not change
C             values up to the deviation of DEVMAX from VMEAN.
C             Default: DEVMAX=999999. (infinity. i.e. no rescaling)
C     DEVEXP=positive real... Exponent for the renormalization to the
C             maximum deviation from the mean value.
C             Has no effect if DEVMAX=999999. (infinity).
C             Default: DEVEXP=2.0
C Names of input and output formatted files with the values:
C     STATIN='string' ... Name of the input file containing the
C             gridded data.
C             No default, 'STATIN' must be specified and cannot be blank
C     STATOUT='string' ... Name of the output file containing the
C             gridded data rescaled to the given statistic properties.
C             Default: STATOUT='grdstat.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,RARRAY,WARRAY
C     ERROR ... File error.for.
C     RSEP1,RSEP3T,RSEP3R,RSEP3I ...
C             File sep.for.
C     RARRAY,WARRAY ... File forms.for.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      CHARACTER*80 FILSEP,FILIN,FILOUT
      INTEGER LU1
      PARAMETER (LU1=1)
      INTEGER N1,N2,N3,I1,N1N2N3
      REAL DSD,VMEAN,DEVMAX,DEVEXP
      REAL DEVINV,V,VMAX,RMEAN,RMS
C
C-----------------------------------------------------------------------
C
C     Reading in a name of the file with the input data:
      WRITE(*,'(A)') '+GRDSTAT: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
C
C     Reading all the data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU1,FILSEP)
      ELSE
C       GRDSTAT-01
        CALL ERROR('GRDSTAT-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
C
C     Reading filenames of the input and output files:
      CALL RSEP3T('STATIN',FILIN,' ')
      IF (FILIN.EQ.' ') THEN
C       GRDSTAT-02
        CALL ERROR('GRDSTAT-02: No input file given.')
      ENDIF
      CALL RSEP3T('STATOUT',FILOUT,'grdstat.out')
C
C     Reading the values describing the grid:
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      N1N2N3=N1*N2*N3
      IF (N1N2N3.GT.MRAM) THEN
C       GRDSTAT-03
        CALL ERROR('GRDSTAT-03: Small array RAM.')
      ENDIF
C     Data to rescale the random values
      CALL RSEP3R('DSD'   ,DSD   ,1.)
      CALL RSEP3R('VMEAN' ,VMEAN ,0.)
      CALL RSEP3R('DEVMAX',DEVMAX,999999.)
      CALL RSEP3R('DEVEXP',DEVEXP,2.)
C
C     Input gridded data:
      CALL RARRAY(LU1,FILIN,'FORMATTED',.TRUE.,0.,N1N2N3,RAM)
C
      WRITE(*,'(A)') '+GRDSTAT: Working ...           '
C
C     Demean grid:
      RMEAN=0.
      DO 12, I1=1,N1N2N3
        RMEAN=RMEAN+RAM(I1)
  12  CONTINUE
      RMEAN=RMEAN/FLOAT(N1N2N3)
      DO 14, I1=1,N1N2N3
        RAM(I1)=RAM(I1)-RMEAN
  14  CONTINUE
C
      IF (DSD.NE.0.) THEN
C       Computing RMS:
        RMS=0.
        DO 20, I1=1,N1N2N3
          RMS=RMS+RAM(I1)**2
  20    CONTINUE
        RMS=SQRT(RMS/FLOAT(N1N2N3))
C
C       Scaling to desired RMS, DSD:
        DO 30, I1=1,N1N2N3
          RAM(I1)=DSD*RAM(I1)/RMS
  30    CONTINUE
      ENDIF
C
C     Rearranging and rescaling the grid values:
      DEVINV=1./DEVEXP
      IF (DEVEXP.GT.999998.) THEN
        VMAX=DEVMAX
      ELSE
        VMAX=DEVMAX*16000000.**DEVINV
      END IF
      DO 50 I1=1,N1N2N3
C       Rescaling the grid values:
        V=RAM(I1)
        IF (DEVMAX.GT.999998.) THEN
          RAM(I1)=VMEAN+V
        ELSE
          IF (ABS(V).GT.VMAX) THEN
            RAM(I1)=VMEAN+SIGN(DEVMAX,V)
          ELSEIF (DEVEXP.GT.999998.) THEN
            RAM(I1)=VMEAN+V
          ELSE
            RAM(I1)=VMEAN+V/(1.+ABS(V/DEVMAX)**DEVEXP)**DEVINV
          ENDIF
        ENDIF
   50 CONTINUE
C
      IF (FILOUT.NE.' ') THEN
        CALL WARRAY(LU1,FILOUT,'FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *              N1N2N3,RAM)
      ENDIF
      WRITE(*,'(A)') '+GRDSTAT: Finished.             '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'sep.for'
C     sep.for
C
C=======================================================================
C