C 
C Program GRDNORM to calculate the spatial density of the Lebesgue norm
C Ln of gridded values
C
C Version: 5.40
C Date: 2000, May 19
C
C Coded 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.......................................................................
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
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C.......................................................................
C
C     Filenames and parameters:
      CHARACTER*80 FILE1,FILE2,FILE3
      INTEGER LU
      REAL UNDEF
      PARAMETER (LU=1,UNDEF=-999999999.)
C     Input data:
      INTEGER N1,N2,N3,N4,N1NEW,N2NEW,N3NEW,N4NEW
      REAL O1,O2,O3,O4,D1,D2,D3,D4
      REAL O1NEW,O2NEW,O3NEW,O4NEW,D1NEW,D2NEW,D3NEW,D4NEW,GNORM
C     Other variables:
      INTEGER N123,N1234N,NSUB,NSUBA,NALL,I,I1,I2,I3,I4
      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
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
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:
      OPEN(LU,FILE=FILE2,FORM='FORMATTED',STATUS='OLD')
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,' ','FORMATTED',.TRUE.,UNDEF,
     *                                       N123,RAM(N1234N+N123*I4+1))
   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,'FORMATTED',.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