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.10 C Date: 1997, October 11 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 C Main input data read from the * input device: C The data are read in by the list directed input (free format), by C a single READ statement. The following items are read: C 'SEP'...String in apostrophes containing the name of the input C file with the data specifying grid dimensions. C 'DAT'...String in apostrophes containing the name of the input C file with the data specifying the statistics. C Note that all input data may be collected in a single file C ('SEP' or 'DAT'), leaving the filename of the other one C blank. C 'FOUT'..String in apostrophes containing the name of the output C ASCII file with the generated grid values. C 'FLOG'..String in apostrophes containing the name of the output C log file which may contain error messages. C ISEED...Nonzero integer seed for the generation of pseudorandom C numbers. Its sign is ignored. C /... Input data line must be terminated by a slash. C Defaults: 'SEP'='grd.h', 'DAT'='grdran2d.dat', 'FOUT'='grd.out', C FLOG='grdran2d.prt', ISEED=-1. C C Both data files 'SEP' and 'DAT' have the same form of the SEP C (Stanford Exploration Project) parameter files: C All the data are specified in the form of PARAMETER=VALUE, e.g. C N1=50, with PARAMETER directly preceding = without intervening C spaces and with VALUE directly following = without intervening C spaces. The PARAMETER=VALUE couple must be delimited by a space C or comma from both sides. C The PARAMETER string is not case-sensitive. C PARAMETER= followed by a space resets the default parameter value. C All other text in the input files is ignored. The file thus may C contain unused data or comments without leading comment character. C Everything between comment character # and the end of the C respective line is ignored, too. C The PARAMETER=VALUE couples may be specified in any order. C The last appearance takes precedence. C Note that the division of the input data into two files 'SEP' and C 'DAT' is, in principle, arbitrary. The data may be arbitrarily mixed C or concentrated to a single file, specifying the name of the other C input file blank. Data of file 'DAT' are read after 'SEP' and take C thus precedence. C C Data 'SEP' specifying grid dimensions (SEP parameter file form): C These data are common to nearly all programs dealing with SEP-like C gridded data formats. 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 C Data 'DAT' specifying the statistics (SEP parameter file form): C These data are specific to this program. 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: ACORG << L << ACOR C Reasonable values for geology: -0.5 < POWERN < 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,FDAT*80,FOUT*80,FLOG*80,CTYPE*1 C C....................................................................... C C Main input data: C Setting the defaults FGRD='grd.h' FDAT='grdran2d.dat' FOUT='grd.out' FLOG='grdran2d.prt' ISEED=-1 C Reading filenames and seed WRITE (*,'(2A)') . ' Enter 4 filenames (Grid data + Statistics data + Output ', . ' grid + Log file), Seed and /: ' READ (*,*) FGRD,FDAT,FOUT,FLOG,ISEED C C Reading all the data from files FGRD and FDAT to the memory C (SEP parameter file form): CALL RSEP1(LU,FGRD) CALL RSEP1(LU,FDAT) 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: 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=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 Closing output log file: CLOSE(LP) 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) 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 IF (ISEED .GE. 0) CALL BOMOU2 (LP, LTW, FLOAT(ISEED), . 'GEN2D1: 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' IF (CFUNC .EQ. ' ') CALL BOMOUT (LP, LTW, . 'GEN2D1: ILLEGAL COR FUNCTION REQUESTED, CTYPE='//CTYPE) IF (RLMAX .LT. RLMIN) CALL BOMOU2(LP, LTW, 0., . 'GEN2D1: RLMIN >= 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 IF (NX2 .LE. 2) CALL BOMOU2 (LP, LTW, FLOAT(NX2), . 'GEN2D1: X FFT HAS BIZARRELY FEW POINTS, NX2=') IF (NY2 .LE. 2) CALL BOMOU2 (LP, LTW, FLOAT(NY2), . 'GEN2D1: 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 CALL BOMOUT (LP, LTW, ' 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 IF (RATIM .GT. TOL) CALL BOMOU2 (LP, LTW, RATIM, . 'GEN2D1: 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 IF (LU1 .NE. 0) WRITE (LU1, 1000) STRING(1:NBLEN(STRING)) 1000 FORMAT (' *******************' / * ' * FATAL ERROR ***** - ', A, / * ' *******************' ) IF (LU2 .NE. LU1 .AND. LU2 .NE. 0) WRITE (LU2 ,1000) * STRING(1:NBLEN(STRING)) STOP 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 IF (LU1 .NE. 0) WRITE (LU1, 1000) STRING(1:NBLEN(STRING)), R 1000 FORMAT (' ******************' / * ' * FATAL ERROR **** - ', * A, G10.3 / * ' ******************' ) IF (LU2 .NE. 0 .AND. LU2 .NE. LU1) WRITE (LU2, 1000) * STRING(1:NBLEN(STRING)), R STOP 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 '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