C 
C Program GRDRAN2D to generate a 2-D rectangular grid of real numbers
C having a desired spatial correlation function, specified variance and
C mean and, if required, falling into the given interval of the
C functional values.
C
C Version: 5.50
C Date: 2001, June 7
C
C Coded by Paul Spudich
C             U.S. Geological Survey
C             345 Middlefield Road,  Menlo Park,  CA 94025,  U.S.A.
C             E-mail: spudich@samoa.wr.usgs.gov
C Input/output modified 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 This code has no known defects, but it is the sole responsibility of
C the user to test the code to ensure that it gives accurate answers in
C the user's application.  The authors of the code take no responsibility
C for damage resulting from errors in this code.
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 Names of the output files:
C     GRD='string' ... Name of the output ASCII file with the generated
C             grid values.
C             For general description of files with gridded data refer
C             to file forms.htm.
C             Defaults: GRD='grd.out'
C     LOG='string' ... Name of the output log file which may contain
C             error messages.
C             Defaults: LOG='grdran2d.out'
C Data specifying grid dimensions:
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     D1=positive real... Grid spacing along the X1 axis.
C             Default: D1=1.
C     D2=positive real... Grid spacing along the X2 axis.
C             Default: D2=1.
C Data specifying the statistics:
C     These data are specific to this program.
C Selection of a particular pseudorandom representation:
C     ISEED=integer ... Nonzero integer seed for the generation of
C             pseudorandom numbers.  Its sign is ignored.
C     Default: ISEED=-1.
C Data for the correlation function:
C     CTYPE='string' ... String specifying the correlation function:
C          ='D'... Default.  The most general option (generalization of
C             all the subsequently listed options).
C             Spectrum of white noise is filtered by function
C               EXP(-(ACORG*k)**2/8.)
C               *(ACOR**(-2)+k**2)**(-(dim/2+POWERN)/2.)
C             Here dim=2 is the spatial dimension.
C             Parameters of the correlation function
C               ACORG,ACOR,POWERN
C          ='L'... Laguerr.
C             Spectrum of white noise is filtered by function
C               EXP(-(ACORG*k)**2/8.)*k**(-(dim/2+POWERN))
C             Parameters of the correlation function
C               ACORG,POWERN
C             Equivalent to: 'D' ACOR=999999.
C             Here 999999. represents infinity.
C          ='K'... Von Karman.
C             Spectrum of white noise is filtered by function
C               (ACOR**(-2)+k**2)**(-(dim/2+POWERN)/2.)
C             Parameters of the correlation function
C               ACOR,POWERN
C             Equivalent to: 'D' ACORG=0.
C          ='S'... Self-affine.
C             Spectrum of white noise is filtered by function
C               k**(-(dim/2+POWERN))
C             Parameter of the correlation function
C               POWERN
C             Equivalent to: 'D' ACORG=0. ACOR=999999.
C                            'L' ACORG=0.
C                            'K' ACOR=999999.
C          ='G'... Gaussian.
C             Spectrum of white noise is filtered by function
C               EXP(-(ACORG*k)**2/8.)
C             Parameter of the correlation function
C               ACORG
C             Equivalent to: 'D' POWERN=-dim/2 ACOR=999999.
C                            'L' POWERN=-dim/2
C             Here dim=2 is the spatial dimension.
C          ='E'... Exponential.
C             Spectrum of white noise is filtered by function
C               (ACOR**(-2)+k**2)**(-(dim/2+0.5)/2.)
C             Parameter of the correlation function
C               ACOR
C             Equivalent to 'D' POWERN=0.5 ACORG=0.
C                           'K' POWERN=0.5
C          ='W'... White noise.
C             No parameters of the correlation function.
C             Equivalent to 'D' POWERN=-dim/2 ACORG=0.
C                           'L' POWERN=-dim/2 ACORG=0.
C                           'K' POWERN=-dim/2
C                           'S' POWERN=-dim/2
C                           'G' ACORG=0.
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 low-pass filter:
C             Cosine filter may be applied to wavenumber power spectrum.
C     RLMIN=nonnegative real... Minimum wavelength.  All wavelengths
C             shorter than RLMIN are removed.
C             Default: RLMIN=0.
C     RLMAX=nonnegative real...
C             Low-pass filter has value 1 at wavenumbers smaller than
C             1/RLMAX (i.e. at wavelengths longer than RLMAX), tapers
C             down between 1/RLMAX and 1/RLMIN, and is zero at
C             wavenumbers greater than 1/RLMIN (i.e. at wavelengths
C             shorter than RLMIN).
C             Default: RLMAX=0. (no cosine filter)
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
C=======================================================================
C
      PARAMETER (IXGDIM=1024, IYGDIM=1024, IWDIM=IXGDIM*IYGDIM)
      PARAMETER (LU=8, LP=7, LTW=0)
C
C     LU...   Logical unit number used to read all input files and write
C             the output file with grid values.
C
      REAL GRID(IXGDIM,IYGDIM)
      COMPLEX WORK(IWDIM)
      INTEGER ISEED
      CHARACTER FGRD*80,FOUT*80,FLOG*80,CTYPE*1
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+GRDRAN2D: Enter input filename: '
      FGRD=' '
      READ(*,*) FGRD
      WRITE(*,'(A)') '+GRDRAN2D: Working ...           '
C
C     Reading all data from the SEP file into the memory:
      IF (FGRD.NE.' ') THEN
        CALL RSEP1(LU,FGRD)
      ELSE
C       GRDRAN2D-08
        CALL ERROR('GRDRAN2D-08: 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 input parameters from the SEP file:
      CALL RSEP3T('GRD',FOUT,'grd.out')
      CALL RSEP3T('LOG',FLOG,'grdran2d.out')
C
C     Recalling the data specifying grid dimensions
C     (arguments: Name of value in input data, Variable, Default):
      CALL RSEP3I('N1',NXG,1)
      CALL RSEP3I('N2',NYG,1)
      CALL RSEP3R('D1',DX ,1.)
      CALL RSEP3R('D2',DY ,1.)
C
C     Recalling the data specifying the statistics:
      CALL RSEP3I('ISEED',ISEED,-1)
C     Data for the correlation function
      CALL RSEP3T('CTYPE' ,CTYPE ,'D')
      CALL RSEP3R('POWERN',POWERN,0.)
      CALL RSEP3R('ACORG' ,ACORG ,0.)
      CALL RSEP3R('ACOR'  ,ACOR  ,999999.)
C     Data for the cosine low-pass filter
      CALL RSEP3R('RLMIN' ,RLMIN ,0.)
      CALL RSEP3R('RLMAX' ,RLMAX ,0.)
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     Opening output log file:
      OPEN (UNIT=LUWARN(LP), FILE=FLOG, FORM='FORMATTED')
C
C     Generating the random grid values with zero mean:
C     The seed must be negative
      JSEED = -IABS(ISEED)
      CALL GEN2D1 (CTYPE, DX, DY, POWERN, ACORG, ACOR, DSD,
     .             RLMIN, RLMAX, GRID, IXGDIM, IYGDIM, NXG, NYG,
     .             WORK, IWDIM, JSEED , LP, LTW)
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
      I=0
      DO 52 IYG=1,NYG
        DO 51 IXG=1,NXG
C         Rearranging 2-D array into 1-D array (new indices IX,IY):
          IY=I/IXGDIM
          I=I+1
          IX=I-IY*IXGDIM
          IY=IY+1
C         Rescaling the grid values:
          V=GRID(IXG,IYG)
          IF(DEVMAX.GT.999998.) THEN
            GRID(IX,IY)=VMEAN+V
          ELSE
            IF(ABS(V).GT.VMAX) THEN
              GRID(IX,IY)=VMEAN+SIGN(DEVMAX,V)
            ELSE IF(DEVEXP.GT.999998.) THEN
              GRID(IX,IY)=VMEAN+V
            ELSE
              GRID(IX,IY)=VMEAN+V/(1+ABS(V/DEVMAX)**DEVEXP)**DEVINV
            END IF
          END IF
   51   CONTINUE
   52 CONTINUE
C
C     Writing output grid values:
      CALL WARRAY(LU,FOUT,'FORMATTED',.FALSE.,0.,
     *                                .FALSE.,0.,NXG*NYG,GRID)
C
C     Closing output log file:
      CLOSE(LP)
C
      WRITE(*,'(A)') '+GRDRAN2D: Done.                 '
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE GEN2D1 (CTYPE, DX, DY, POWERN, ACORG, ACOR, DSD,
     .                   RLMIN, RLMAX, GRID, IXGDIM, IYGDIM, NXG, NYG,
     .                   WORK, IWDIM, ISEED, LP, LTW)
C-----------------------------------------------------------------------
C
C Subroutine to generate a rectangular grid of real numbers having a
C desired spatial correlation function and [# an approximately Gaussian
C distribution function having] a specified variance and zero mean.  The
C application this routine is specifically written for is to generate
C a 2D array of seismic velocity perturbations having either Gaussian,
C exponential, or self-similar autocorrelation functions, following
C Frankel and Clayton, J. Geophys. Res, 1986, pp. 6465-6489.  Note that
C in this implementation of their functions, I have anti-aliased the
C wavenumber power spectrum using a cosine taper.
C
C # Note: Present version seems to generate rectangular distribution
C         subsequently modified by the filtration (L.K.).
C
C Inputs:
C   CTYPE - CHARACTER*1 variable specifying which correlation function
C               to use:
C               = 'D' The most general option (generalization of all the
C                     subsequently listed options):
C                     Spectrum of white noise is filtered by function
C                       EXP( -(ACORG * k)**2 / 8. )
C                       * ( ACOR**(-2) + k**2 )**( -(1.+POWERN)/2. )
C               = 'L' Laguerr (common generalization of self-affine and
C                     Gaussian),  equivalent to 'D' with ACOR=infinity.
C               = 'K' Von Karman (Frankel and Clayton eqn B12 with
C                     a=ACOR and exponent -1 generalized to -1-2*POWERN,
C                     see eqn 4 with m=POWERN),  equivalent to 'D' with
C                     ACORG=0.
C               = 'S' self-affine (Von Karman with ACOR=infinity),
C                     equivalent to 'D' with ACORG=0 and ACOR=infinity.
C               = 'G' Gaussian (Frankel and Clayton eqn B3 with
C                     a=ACORG),  equivalent to 'D' with POWERN=-1 and
C                     ACOR=infinity.
C               = 'E' exponential (Frankel and Clayton eqn B6 with
C                     a=ACOR),  equivalent to 'D' with POWERN=0.5 and
C                     ACORG=0.
C               = 'W' white noise,  equivalent to 'D' with POWERN=-1,
C                     ACORG=0 and ACOR=infinity.
C   DX - spatial grid spacing corresponding to left index of grid array
C               (typically km)
C   DY - spatial grid spacing corresponding to right index of grid array
C               (typically km)
C   ACORG - correlation distance of Gaussian media, in same units as DX
C               and DY.  ACORG corresponds to the variable 'a' in
C               Frankel and Clayton's equation B3.
C   ACOR - correlation distance of Von Karman media, in same units as DX
C               and DY.  ACOR corresponds to the variable 'a' in Frankel
C               and Clayton's equations B6 and B12.
C   DSD - desired standard deviation of the perturbations calculated.
C               = 0 means do not alter standard deviation as calculated
C               The output grid is scaled to have standard
C               deviation = DSD if DSD is input nonzero.  Note that
C               DSD does not correspond to sigma-c in F&C eqn B3.
C   RLMIN, RLMAX -cosine filter is applied to wavenumber power spectrum.
C               Filter has value 1 at wavenumbers corresponding to
C               a wavelength of RLMAX, filter has value 0 at k
C               corresponding to wavelength of RLMIN.  I.e. if you
C               input RLMIN = 2, all wavelengths shorter than 2
C               will be filtered out of the result.
C               Units = units of DX, DY, ACORG, ACOR.
C   IXGDIM, IYGDIM - grid is dimensioned (IXGDIM, IYGDIM) outside this
C               routine
C   NXG, NYG - grid (1 to NXG, 1 to NYG) used
C   WORK - 1D complex array for working storage; contents on input
C               ignored.
C   IWDIM - WORK is dimensioned iwdim outside this routine.  Note -
C               IWDIM must equal or exceed (next power of 2 larger than
C               NXG) times (next power of 2 larger than NYG), i.e. if
C               NXG = 10 and NYG = 20, IWDIM must be .GE. 16*32
C   ISEED - integer seed for generating random numbers.
C               ISEED must be negative.
C   LP, LTW - logical unit numbers to which to write error messages
C               ( LP = LTW is okay)
C
C Outputs:
C   GRID - real array of velocities having desired statistical
C               properties.
C               GRID(1 to NXG, 1 to NYG) are set by this program.
C------------------------------------------------------------------------
C
      DIMENSION GRID(IXGDIM, IYGDIM), NN(2)
      COMPLEX WORK (IWDIM)
      CHARACTER*1 CTYPE, CFUNC

      PI = 3.1415926

C TOL is a numeric tolerance for checking that IFFT result is pure real
      TOL = .5E-5

C     GRDRAN2D-01
      IF (ISEED .GE. 0) CALL BOMOU2 (LP, LTW, FLOAT(ISEED),
     .     'GRDRAN2D-01: ENTER ISEED NEGATIVE.  ILLEGAL ISEED=')

      CFUNC = ' '
      IF (CTYPE .EQ. 'D' .OR. CTYPE .EQ. 'd') CFUNC = 'D'
      IF (CTYPE .EQ. 'L' .OR. CTYPE .EQ. 'l') CFUNC = 'L'
      IF (CTYPE .EQ. 'K' .OR. CTYPE .EQ. 'k') CFUNC = 'K'
      IF (CTYPE .EQ. 'S' .OR. CTYPE .EQ. 's') CFUNC = 'S'
      IF (CTYPE .EQ. 'G' .OR. CTYPE .EQ. 'g') CFUNC = 'G'
      IF (CTYPE .EQ. 'E' .OR. CTYPE .EQ. 'e') CFUNC = 'E'
      IF (CTYPE .EQ. 'W' .OR. CTYPE .EQ. 'w') CFUNC = 'W'
C     GRDRAN2D-02
      IF (CFUNC .EQ. ' ') CALL BOMOUT (LP, LTW,
     .    'GRDRAN2D-02: ILLEGAL COR FUNCTION REQUESTED, CTYPE='//CTYPE)

C     GRDRAN2D-03
      IF (RLMAX .LT. RLMIN) CALL BOMOU2(LP, LTW, 0.,
     .    'GRDRAN2D-03: RLMIN .GE. RLMAX =')
      IF(RLMIN.GT.0.) THEN
        RKMIN = 2*PI/RLMAX
      ELSE
        RKMIN = 999999.
      END IF
      IF(RLMAX.GT.0.) THEN
        RKMAX = 2*PI/RLMIN
      ELSE
        RKMAX = 999999.
      END IF

C first determine the powers of 2 to use in the 2D FFT

      EXP2 = LOG10(FLOAT(NXG))/LOG10(2.)
      DIF = EXP2 - IFIX(EXP2)
      IF (DIF .EQ. 0) THEN
         NX2 = 2**IFIX(EXP2)
      ELSE
         NX2 = 2** (IFIX(EXP2)+1)
      END IF

      EXP2 = LOG10(FLOAT(NYG))/LOG10(2.)
      DIF = EXP2 - IFIX(EXP2)
      IF (DIF .EQ. 0) THEN
         NY2 = 2**IFIX(EXP2)
      ELSE
         NY2 = 2** (IFIX(EXP2)+1)
      END IF
      NN(1) = NX2
      NN(2) = NY2

C     GRDRAN2D-04
      IF (NX2 .LE. 2) CALL BOMOU2 (LP, LTW, FLOAT(NX2),
     .      'GRDRAN2D-04: X FFT HAS BIZARRELY FEW POINTS, NX2=')
C     GRDRAN2D-05
      IF (NY2 .LE. 2) CALL BOMOU2 (LP, LTW, FLOAT(NY2),
     .      'GRDRAN2D-05: Y FFT HAS BIZARRELY FEW POINTS, NY2=')

C determine physical extent of grid, wavenumber interval
      XLEN = NX2 * DX
      YLEN = NY2 * DY
      DKX = 2 * PI / XLEN
      DKY = 2 * PI / YLEN

C fill work array with random numbers distributed uniformly between
C -0.5 and 0.5
      DO 100 I = 1, NX2*NY2
*       WORK(I) = CMPLX(RAN2(ISEED)-.5, 0.)
        WORK(I) = CMPLX(RAN3(ISEED)-.5, 0.)
  100   CONTINUE

      NDIM = 2
      ISIGN = 1
C two-dimensional FFT
      CALL FOURN (WORK, NN, NDIM, ISIGN)

C set kx=ky=0 component of WORK to zero so that its inverse
C transform has zero mean
      WORK(1) = CMPLX(0.,0.)

C now multiply all elements of the WORK array by the desired wavenumber
C filter.  We must be careful to multiply each quadrant properly to
C preserve the Hermitian property of the WORK array.  Recall that the
C discrete 2D fft puts the RKX=0, RKY=0 spectral component at (1,1), and
C the results for negative RKX and negative RKY appear at high values
C of the array indices.

      IXNYQ = NX2/2 + 1
      IYNYQ = NY2/2 + 1
      IWORK = 0
      POWER = -(1.+POWERN) / 2.
      AGSQ8 = ACORG**2 / 8.
      IF(ACOR.LE.999998.) THEN
        ASQINV = 1./ACOR**2
      ELSE
        ASQINV = 0.
      END IF

      DO 200 IY = 1, NY2
        RKY = (IY-1) * DKY
        IF (IY .GT. IYNYQ) RKY = (IY - NY2 - 1) * DKY

        DO 200 IX = 1, NX2
          RKX = (IX-1) * DKX
          IF (IX .GT. IXNYQ) RKX = (IX - NX2 -1) * DKX
          RKSQ = RKX**2 + RKY**2
          RK = SQRT (RKSQ)
C
          IF(IX.EQ.1.AND.IY.EQ.1) THEN
C           Zero average value of each realization:
            P=0.
          ELSE
C wavenumber amplitude spectra of the different correlation functions.
C Note that I will
C rescale the result to the desired variance at the end.
            IF (CFUNC .EQ. 'D') THEN
C             General:
              P = EXP(-AGSQ8*RKSQ) * (ASQINV+RKSQ)**POWER
            ELSE IF (CFUNC .EQ. 'L') THEN
C             Laguerr:
              P = EXP(-AGSQ8*RKSQ) * RKSQ**POWER
            ELSE IF (CFUNC .EQ. 'K') THEN
C             Von Karman:
              P = (ASQINV+RKSQ)**POWER
            ELSE IF (CFUNC .EQ. 'S') THEN
C             Self-affine:
              P = RKSQ**POWER
            ELSE IF (CFUNC .EQ. 'G') THEN
C             Gaussian:
              P = EXP(-AGSQ8*RKSQ)
            ELSE IF (CFUNC .EQ. 'E') THEN
C             Exponential:
              P = (ASQINV+RKSQ)**(-0.75)
            ELSE IF (CFUNC .EQ. 'W') THEN
C             White noise:
              P = 1.
            ELSE
C             GRDRAN2D-06
              CALL BOMOUT(LP,LTW,'GRDRAN2D-06: ILLEGAL CFUNC='//CFUNC)
            END IF
          END IF

          IWORK = IWORK + 1
          WORK(IWORK) = P * WORK(IWORK)

C apply cosine taper
          IF (RK .GE. RKMAX) THEN
            COSTAP = 0.
          ELSE IF (RK .GT. RKMIN .AND. RK .LT. RKMAX) THEN
            COSTAP = .5 + .5 * COS ( (RK-RKMIN)/(RKMAX-RKMIN)*PI )
          ELSE
            COSTAP = 1.
          END IF
C now apply filter to work array
          WORK(IWORK) = WORK(IWORK) * COSTAP

  200   CONTINUE

C inverse FFT
      ISIGN = -1
      CALL FOURN (WORK, NN, NDIM, ISIGN)

C Check that the result is pure real (to numerical accuracy).

      ABMXR = 0.
      ABMXI = 0.
      NSQ = NX2*NY2
      DO 300 I = 1, NSQ
C normalize for the 2D FFT
        WORK(I) = WORK(I) / NSQ
        ABMXR = MAX(ABMXR, ABS(REAL(WORK(I))))
        ABMXI = MAX(ABMXI, ABS(AIMAG(WORK(I))))
  300   CONTINUE

      RATIM = ABMXI/ABMXR
C     GRDRAN2D-07
      IF (RATIM .GT. TOL) CALL BOMOU2 (LP, LTW, RATIM,
     .     'GRDRAN2D-07: IFFT RESULT NOT PURE REAL, MAX(IM)/MAX(REAL)=')

C normalize to desired variance and load in output grid array
      GMEAN = 0.
      DO 401 IY = 1, NYG
        DO 401 IX = 1, NXG
          IWORK = NX2 * (IY-1) + IX
          GRID(IX,IY) = REAL(WORK(IWORK))
          GMEAN = GMEAN + GRID(IX,IY)
  401     CONTINUE

      NG = NXG*NYG
      GMEAN = GMEAN / NG

C demean grid and calculate rms value
      RMS = 0.
      DO 402 IY = 1, NYG
        DO 402 IX = 1, NXG
          GRID(IX,IY) = GRID(IX,IY) - GMEAN
          RMS = RMS + GRID(IX,IY)**2
  402     CONTINUE

      RMS = SQRT (RMS / NG)

C scale grid to desired rms, DDSD
      IF (DSD .NE. 0) THEN
        DDSD = DSD
      ELSE
        DDSD = 1.
        RMS = 1.
      END IF

      DO 400 IY = 1, NYG
        DO 400 IX = 1, NXG
          GRID(IX,IY) = DDSD * GRID(IX,IY) / RMS
  400     CONTINUE

C     ------
      RETURN
C     ------

      END
C
C=======================================================================
C
      SUBROUTINE BOMOUT (LU1, LU2, STRING)
C----------------------------------------------------------------------
C
C  THIS ROUTINE WRITES TO UNITS LU1 AND LU2 AN ERROR MESSAGE CONSISTING OF
C  " *** FATAL ERROR ***" CONCATENATED WITH THE INPUT STRING, AND THEN
C  STOPS EXECUTION OF THE CALLING PROGRAM.  IF LU1=LU2, ONLY ONE MESSAGE
C  IS WRITTEN TO LU1.
C
      CHARACTER*(*) STRING
c     IF (LU1 .NE. 0) WRITE (LU1, 1000) STRING(1:NBLEN(STRING))
c1000 FORMAT (' *******************' /
c    *        ' * FATAL ERROR ***** - ', A, /
c    *        ' *******************' )
c     IF (LU2 .NE. LU1 .AND. LU2 .NE. 0) WRITE (LU2 ,1000)
c    *                 STRING(1:NBLEN(STRING))
c     STOP
      CALL ERROR(STRING)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE BOMOU2 (LU1, LU2, R, STRING)
C----------------------------------------------------------------------
C
C  THIS ROUTINE WRITES TO UNITS LU1 AND LU2 AN ERROR MESSAGE CONSISTING OF
C  " *** FATAL ERROR ***" CONCATENATED WITH THE INPUT STRING,
C  CONCATENATED WITH THE REAL NUMBER R.  BOMOU2 THEN
C  STOPS EXECUTION OF THE CALLING PROGRAM.  IF LU1=LU2, ONLY ONE MESSAGE
C  IS WRITTEN TO LU1.
C
      CHARACTER*(*) STRING
c     IF (LU1 .NE. 0) WRITE (LU1, 1000) STRING(1:NBLEN(STRING)), R
c1000 FORMAT (' ******************' /
c    *        ' * FATAL ERROR **** - ',
c    *            A, G10.3 /
c    *        ' ******************' )
c     IF (LU2 .NE. 0 .AND. LU2 .NE. LU1) WRITE (LU2, 1000)
c    *              STRING(1:NBLEN(STRING)), R
c     STOP
      CHARACTER*200 TEXT
      WRITE(TEXT,'(A,G10.3)') STRING(1:NBLEN(STRING)),R
      CALL ERROR(TEXT(1:NBLEN(TEXT)))
      RETURN
      END
C
C=======================================================================
C
      FUNCTION NBLEN (STRING)
C--------------------------------------------------------------------
C
C  GIVEN A CHARACTER STRING, NBLEN RETURNS THE LENGTH OF THE STRING
C TO THE LAST NON-BLANK CHARACTER, PRESUMING THE STRING IS LEFT-
C JUSTIFIED, I.E. IF STRING = '   XS  J   ', NBLEN = 8.
C
C CALLED NON-LIBRARY ROUTINES: NONE
C LANGUAGE: STANDARD FORTRAN 77
C
      CHARACTER*(*) STRING, BLANK*1, NULL*1
      DATA BLANK /' '/
C
      NULL = CHAR(0)
      NBLEN = 0
      LS = LEN(STRING)
      IF (LS .EQ. 0) RETURN
      DO 1 I = LS, 1, -1
         IF (STRING(I:I) .NE. BLANK .AND. STRING(I:I) .NE. NULL) GO TO 2
    1    CONTINUE
      RETURN
    2 NBLEN = I
      RETURN
      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
      INCLUDE 'fourn.for'
C     fourn.for
*     INCLUDE 'ran2.for'
C     ran2.for
      INCLUDE 'ran3.for'
C     ran3.for
C
C=======================================================================
C