C
C Program GRDNORM to calculate the spatial density of the Lebesgue norm C Ln of gridded values C C Version: 8.10 C Date: 2022, December 7 C C Coded 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....................................................................... 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 input and output files: C GRD='string'... Name of the input ASCII file with the grid values. C Default: GRD='grd.out' C GRDNEW='string'... Name of the output ASCII file containing the C grid values of the calculated norm. C Default: GRDNEW='grdnew.out' C For general description of the files with gridded data refer C to file forms.htm. C Data specifying dimensions of the input 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 N4=positive integer... Number of spatial grids (number of time C levels). C Default: N4=1 C O1=real... First coordinate of the grid origin (first point of the C grid). C Default: O1=0. C O2=real... Second coordinate of the grid origin. C Default: O2=0. C O3=real... Third coordinate of the grid origin. C Default: O3=0. C O4=real... Time corresponding to the first spatial grid. C Default: O4=0. C D1=real... Grid interval in the direction of the first coordinate C axis. C Default: D1=1. C D2=real... Grid interval in the direction of the second coordinate C axis. C Default: D2=1. C D3=real... Grid interval in the direction of the second coordinate C axis. C Default: D3=1. C D4=real... Time interval. C Default: D4=1. C Data specifying dimensions of the output grid: C N1NEW=positive integer C N2NEW=positive integer C N3NEW=positive integer C N4NEW=positive integer C O1NEW=real C O2NEW=real C O3NEW=real C O4NEW=real C D1NEW=real C D2NEW=real C D3NEW=real C D4NEW=real... Analogous to N1, N2, N3, N4, O1, O2, O3, O4, D1, C D2, D3 and D4, respectively, but for the output grid. C The output grid should be a coarser grid than the input C grid. The input grid is then divided into the subgrids C corresponding to individual gridpoints of the coarser C output grid. Each subgrid is composed of gridpoints of C the input grid, which are closer to the corresponding C gridpoint than to other gridpoints of the output grid. C Empty subgrids are not allowed. The Lebesgue norm Ln, C where n is given by parameter GNORM, is calculated over C each subgrid. Output file GRDNEW thus contains the grid C of N1NEW*N2NEW*N3NEW*N4NEW norms of individual subgrids. C Examples: C (a) If N1NEW=1, N2NEW=1, N3NEW=1 and N4NEW=1, the norm of C the whole grid is calculated. C (b) If N1NEW=1 N2NEW=1 N3NEW=1, whereas N4NEW is not C specified and defaults to N4, the norm of the spatial C grid is calculated separately at each time level. C Output than consists of N4 norms of spatial grids. C (c) If N2NEW=1, whereas N1NEW, N3NEW and N4NEW are not C specified and default to the dimensions of the input C grid, input grid contains a probability density C function, and GNORM=1, the output grid contains the C gridded marginal probability density in the X1X3 plane C at each time level. Output than consists of N1*N3*N4 C values. C Defaults: N1NEW=N1, N2NEW=N2, N3NEW=N3, N4NEW=N4, C O1NEW=O1, O2NEW=O2, O3NEW=O3, O4NEW=O4, C D1NEW=D1, D2NEW=D2, D3NEW=D3, D4NEW=D4. C Type of the norm: C GNORM=real... The norm is [sum( ABS(G)**GNORM )/NSUB]**(1/GNORM). C The summation is performed over each subgrid. NSUB is the C number of points with defined values within the subgrid. C If GNORM.EQ.0., the harmonic average is calculated. C If GNORM.EQ.1., the arithmetic average is calculated. C If GNORM.GT.998., the maximum is calculated. C If GNORM.LT.-998., the minimum is calculated. C Default: GNORM=1. C Range of considered values: C VMIN=real, VMAX=real... Values from input file GRD less than or C equal to VMIN, or greater than or equal to VMAX will be C deemed undefined. C Defaults: VMIN=-999999, VMAX=999999000000. 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 Form of the input/output files: C FORM='string'... Form of the input/output files: either C 'FORMATTED' or 'UNFORMATTED' (case insensitive). C Default: FORM='FORMATTED' C FORMW='string'... Form of the output file. If both FORM and FORMW C are specified, FORMW is used for output files. C Default: FORMW=FORM C FORMR='string'... Form of the input file. If both FORM and FORMR C are specified, FORMR is used for output files. C Default: FORMR=FORM 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 C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C....................................................................... C EXTERNAL UARRAY REAL UARRAY C C Filenames and parameters: CHARACTER*80 FILE1,FILE2,FILE3 INTEGER LU REAL UNDEF PARAMETER (LU=1) C Input data: INTEGER N1,N2,N3,N4,N1NEW,N2NEW,N3NEW,N4NEW REAL O1,O2,O3,O4,D1,D2,D3,D4,VMIN,VMAX REAL O1NEW,O2NEW,O3NEW,O4NEW,D1NEW,D2NEW,D3NEW,D4NEW,GNORM C Other variables: INTEGER N123,N1234N,NSUB,NSUBA,NALL,I,I1,I2,I3,I4,I5 INTEGER INEW,I1NEW,I2NEW,I3NEW,I4NEW INTEGER I1MIN,I2MIN,I3MIN,I4MIN,I1MAX,I2MAX,I3MAX,I4MAX REAL X1,X2,X3,X4,RAMNEW C C----------------------------------------------------------------------- C C Reading name of SEP file with input data: WRITE(*,'(A)') '+GRDNORM: Enter input filename: ' FILE1=' ' READ(*,*) FILE1 WRITE(*,'(A)') '+GRDNORM: Working ... ' C C Reading all data from the SEP file into the memory: IF (FILE1.NE.' ') THEN CALL RSEP1(LU,FILE1) ELSE C GRDNORM-06 CALL ERROR('GRDNORM-06: 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 UNDEF=UARRAY() C C Reading input parameters from the SEP file: CALL RSEP3T('GRD',FILE2,'grd.out') CALL RSEP3T('GRDNEW',FILE3,'grdnew.out') C C Reading grid dimensions: C Original grid: CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3I('N4',N4,1) CALL RSEP3R('O1',O1,0.) CALL RSEP3R('O2',O2,0.) CALL RSEP3R('O3',O3,0.) CALL RSEP3R('O4',O4,0.) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) CALL RSEP3R('D4',D4,1.) C New grid: CALL RSEP3I('N1NEW',N1NEW,N1) CALL RSEP3I('N2NEW',N2NEW,N2) CALL RSEP3I('N3NEW',N3NEW,N3) CALL RSEP3I('N4NEW',N4NEW,N4) CALL RSEP3R('O1NEW',O1NEW,O1) CALL RSEP3R('O2NEW',O2NEW,O2) CALL RSEP3R('O3NEW',O3NEW,O3) CALL RSEP3R('O4NEW',O4NEW,O4) CALL RSEP3R('D1NEW',D1NEW,D1) CALL RSEP3R('D2NEW',D2NEW,D2) CALL RSEP3R('D3NEW',D3NEW,D3) CALL RSEP3R('D4NEW',D4NEW,D4) C Type of the norm: CALL RSEP3R('GNORM',GNORM,1.) C Range of considered values: CALL RSEP3R('VMIN',VMIN,-999999.) CALL RSEP3R('VMAX',VMAX, 999999000000.) C C Dimensions of the new grid related to the old one: O1=(O1NEW+0.5*D1NEW-O1)/D1 O2=(O2NEW+0.5*D2NEW-O2)/D2 O3=(O3NEW+0.5*D3NEW-O3)/D3 O4=(O4NEW+0.5*D4NEW-O4)/D4 D1=D1NEW/D1 D2=D2NEW/D2 D3=D3NEW/D3 D4=D4NEW/D4 N123 =N1*N2*N3 N1234N=N1NEW*N2NEW*N3NEW*N4NEW C IF((N1NEW.GT.0.AND.D1.LE.0.).OR. * (N2NEW.GT.0.AND.D2.LE.0.).OR. * (N3NEW.GT.0.AND.D3.LE.0.).OR. * (N4NEW.GT.0.AND.D4.LE.0.)) THEN C GRDNORM-04 CALL ERROR('GRDNORM-04: Oposite signs of grid steps)') C D1NEW must have the same sign as D1, C D2NEW must have the same sign as D2, C D3NEW must have the same sign as D3 and C D4NEW must have the same sign as D4 in this version. END IF IF(N1234N+N123.GT.MRAM) THEN C GRDNORM-01 CALL ERROR('GRDNORM-01: Too small array RAM(MRAM)') C Too small array RAM(MRAM) to allocate both input and output C grid values. If possible, increase dimension MRAM in include C file ram.inc. END IF C C Reading input grid values: CALL OARRAY(LU,FILE2,'R') C C Calculating the spatial densities of the Lebesgue norm: NALL=0 I4MIN=0 DO 24 I4NEW=0,N4NEW-1 IF(I4NEW.LT.N4NEW-1) THEN X4=O4+FLOAT(I4NEW)*D4 I4MAX=MIN0(INT(X4),N4-1) ELSE I4MAX=N4-1 END IF IF(N1234N+N123*(I4MAX-I4MIN).GT.MRAM) THEN C GRDNORM-05 CALL ERROR('GRDNORM-05: Too small array RAM(MRAM)') C Too small array RAM(MRAM) to allocate both input and output C grid values. If possible, increase dimension MRAM in include C file ram.inc. END IF DO 10 I4=0,I4MAX-I4MIN CALL RARRAY(LU,' ',.TRUE.,UNDEF, * N123,RAM(N1234N+N123*I4+1)) DO 9, I5=N1234N+N123*I4+1,N1234N+N123*I4+N123 IF(RAM(I5).LE.VMIN.OR.VMAX.LE.RAM(I5)) THEN RAM(I5)=UNDEF ENDIF 9 CONTINUE 10 CONTINUE I3MIN=0 DO 23 I3NEW=0,N3NEW-1 IF(I3NEW.LT.N3NEW-1) THEN X3=O3+FLOAT(I3NEW)*D3 I3MAX=MIN0(INT(X3),N3-1) ELSE I3MAX=N3-1 END IF I2MIN=0 DO 22 I2NEW=0,N2NEW-1 IF(I2NEW.LT.N2NEW-1) THEN X2=O2+FLOAT(I2NEW)*D2 I2MAX=MIN0(INT(X2),N2-1) ELSE I2MAX=N2-1 END IF INEW=N1NEW*(I2NEW+N2NEW*(I3NEW+N3NEW*I4NEW)) I1MIN=0 DO 21 I1NEW=0,N1NEW-1 IF(I1NEW.LT.N1NEW-1) THEN X1=O1+FLOAT(I1NEW)*D1 I1MAX=MIN0(INT(X1),N1-1) ELSE I1MAX=N1-1 END IF RAMNEW=0. NSUB=0 NSUBA=0 DO 14 I4=0,I4MAX-I4MIN DO 13 I3=I3MIN,I3MAX DO 12 I2=I2MIN,I2MAX I=N1234N+I1MIN+N1*(I2+N2*(I3+N3*I4)) DO 11 I1=I1MIN,I1MAX I=I+1 NSUBA=NSUBA+1 IF(RAM(I).NE.UNDEF) THEN NSUB=NSUB+1 IF(GNORM.EQ.0.) THEN RAMNEW=RAMNEW+ALOG(ABS(RAM(I))) ELSE IF(GNORM.EQ.1.) THEN RAMNEW=RAMNEW+RAM(I) ELSE IF(GNORM.GT.998.) THEN IF(NSUB.LE.1) THEN RAMNEW= RAM(I) ELSE RAMNEW=AMAX1(RAM(I),RAMNEW) END IF ELSE IF(GNORM.LT.-998.) THEN IF(NSUB.LE.1) THEN RAMNEW= RAM(I) ELSE RAMNEW=AMIN1(RAM(I),RAMNEW) END IF ELSE RAMNEW=RAMNEW+ABS(RAM(I))**GNORM END IF END IF 11 CONTINUE 12 CONTINUE 13 CONTINUE 14 CONTINUE NALL=NALL+NSUBA IF(NSUBA.EQ.0) THEN C C GRDNORM-02 CALL ERROR('GRDNORM-02: Empty subgrid') C The Lebesgue norm cannot be calculated and averaged C over an empty subgrid consisting of no gridpoint of C the given grid. Check the data for the 'new' grid. END IF IF(NSUB.EQ.0) THEN RAMNEW=UNDEF ELSE IF(GNORM.GE.-998..AND.GNORM.LE.998.) THEN RAMNEW=RAMNEW/FLOAT(NSUB) IF(GNORM.EQ.0.) THEN RAMNEW=EXP(RAMNEW) ELSE IF(GNORM.NE.1.) THEN RAMNEW=RAMNEW**(1./GNORM) END IF END IF INEW=INEW+1 RAM(INEW)=RAMNEW I1MIN=I1MAX+1 21 CONTINUE I2MIN=I2MAX+1 22 CONTINUE I3MIN=I3MAX+1 23 CONTINUE I4MIN=I4MAX+1 24 CONTINUE IF(NALL.NE.N1*N2*N3*N4) THEN WRITE(*,*) ' Subgrids cover',NALL,' gridpoints of',N1*N2*N3*N4 C GRDNORM-03 CALL ERROR('GRDNORM-03: Gaps between subgrids') C This error should not apperar. Contact the author. END IF C C Writing output grid values: CALL WARAY(LU,FILE3,.TRUE.,UNDEF,.FALSE.,0., * N1NEW*N2NEW*N3NEW,N4NEW,RAM(1)) WRITE(*,'(A)') * '+GRDNORM: Done. ' STOP 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 C C======================================================================= C