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: 7.40 C Date: 2017, May 19 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 http://sw3d.cz/staff/klimes.htm 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 C responsibility 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 initial distribution before wavenumber filtering: C RANDIS='string'... Character specifying the distribution: C RANDIS='U' or RANDIS='u': Uniform distribution between C -0.5 and 0.5. C RANDIS='G' or RANDIS='g': Gaussian distribution. C Default: RANDIS='U' 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 Default: CTYPE='D' 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 Optional parameters specifying the form of the quantities C written in the output formatted files: C MINDIG,MAXDIG=positive integers ... See the description in file C forms.for. C NUMLIN=positive integer ... Number of reals or integers in one C line of the output file. See the description in file C forms.for. C Value of undefined quantities: C UNDEF=real... The value to be used for undefined real quantities. C Default: UNDEF=undefined value used in forms.for 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,RANDIS*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 RSEP3T('RANDIS',RANDIS,'U') 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, RANDIS, 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,RANDIS, 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, RANDIS, CFUNC C PI = 3.1415926 C C TOL is a numeric tolerance for checking that IFFT result is pure real TOL = .5E-5 C C GRDRAN2D-01 IF (ISEED .GE. 0) CALL BOMOU2 (LP, LTW, FLOAT(ISEED), . 'GRDRAN2D-01: ENTER ISEED NEGATIVE. ILLEGAL ISEED=') C 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 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 C first determine the powers of 2 to use in the 2D FFT C 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 C 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 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 C determine physical extent of grid, wavenumber interval XLEN = NX2 * DX YLEN = NY2 * DY DKX = 2 * PI / XLEN DKY = 2 * PI / YLEN C C fill work array with random numbers distributed uniformly between C -0.5 and 0.5 IF (RANDIS.EQ.'U'.OR.RANDIS.EQ.'u') THEN DO 100 I = 1, NX2*NY2 * WORK(I) = CMPLX(RAN2(ISEED)-.5, 0.) WORK(I) = CMPLX(RAN3(ISEED)-.5, 0.) 100 CONTINUE C C Gaussian distribution ELSE DO 111 I = 1, NX2*NY2, 2 * WORK(I) = CMPLX(RAN2(ISEED)-.5, 0.) 110 CONTINUE V1 = 2.*RAN3(ISEED)-1. V2 = 2.*RAN3(ISEED)-1. VV = V1*V1+V2*V2 IF (VV.GE.1.) GO TO 110 VV = SQRT(-2.*ALOG(VV)/VV) WORK(I) = CMPLX(V1*VV, 0.) IF (I.LT.NX2*NY2) THEN WORK(I+1) = CMPLX(V2*VV, 0.) END IF 111 CONTINUE ENDIF C NDIM = 2 ISIGN = 1 C two-dimensional FFT CALL FOURN (WORK, NN, NDIM, ISIGN) C 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 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. C 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 C DO 200 IY = 1, NY2 RKY = (IY-1) * DKY IF (IY .GT. IYNYQ) RKY = (IY - NY2 - 1) * DKY C 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 C IWORK = IWORK + 1 WORK(IWORK) = P * WORK(IWORK) C 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 C 200 CONTINUE C C inverse FFT ISIGN = -1 CALL FOURN (WORK, NN, NDIM, ISIGN) C C Check that the result is pure real (to numerical accuracy). C 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 C 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 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 C NG = NXG*NYG GMEAN = GMEAN / NG C 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 C RMS = SQRT (RMS / NG) C C scale grid to desired rms, DDSD IF (DSD .NE. 0) THEN DDSD = DSD ELSE DDSD = 1. RMS = 1. END IF C DO 400 IY = 1, NYG DO 400 IX = 1, NXG GRID(IX,IY) = DDSD * GRID(IX,IY) / RMS 400 CONTINUE C C ------ RETURN C ------ 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 C OF " *** FATAL ERROR ***" CONCATENATED WITH THE INPUT STRING, AND C THEN STOPS EXECUTION OF THE CALLING PROGRAM. IF LU1=LU2, ONLY ONE C MESSAGE 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 C OF " *** 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 of Numerical Recipes * INCLUDE 'ran2.for' C ran2.for of Numerical Recipes INCLUDE 'ran3.for' C ran3.for of Numerical Recipes C C======================================================================= C