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 Optional parameters specifying the form of the real quantities C written in the output formatted files: C MINDIG,MAXDIG=positive integers ... See the description in file C forms.for. 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