C 
C Program GRDTRANS to transpose the coordinate axes of the gridded data
C
C Version: 5.50
C Date: 2001, May 28
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'... Names of the input ASCII file with the grid
C             values.
C             Default: GRD='grd.out'
C     GRDNEW='string'... Name of the output ASCII file containing the
C             input grid values ordered according to the transposed
C             coordinates.
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 fastest X1
C             axis (inner loop).
C             Default: N1=1
C     N2=positive integer... Number of gridpoints along the X2 axis
C             (intermediate loop).
C             Default: N2=1
C     N3=positive integer... Number of gridpoints along the slowest
C             spatial X3 axis (outer spatial loop).
C     N4=positive integer... Number of gridpoints along optional time X4
C             axis (outermost, temporal loop).
C             Default: N4=1
C Data specifying the output grid coordinates:
C     NEWX1=positive integer... Index of the input axis, corresponding
C             to the fastest output axis NEWX1 (inner output loop).
C             NEWX1=1: Fastest input axis (inner input loop),
C             NEWX1=2: Medium input axis (intermediate input loop),
C             NEWX1=3: Slowest input spatial axis (outer spatial input
C                      loop).
C             NEWX1=4: Input time axis (outermost, temporal input loop).
C             Default: NEWX1=1
C     NEWX2=positive integer... Index of the input axis, corresponding
C             to the output axis NEWX2.  Analogous to NEWX1.
C             Default: NEWX2=2
C     NEWX3=positive integer... Index of the input axis, corresponding
C             to the slowest output axis NEWX3.  Analogous to NEWX1.
C             Default: NEWX3=3
C     NEWX4=positive integer... Index of the input axis, corresponding
C             to the slowest output axis NEWX4.  Analogous to NEWX1.
C             Default: NEWX4=4
C Optional parameters specifying the form of the real quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C.......................................................................
C
C     Filenames and parameters:
      CHARACTER*80 FSEP,FGRD1,FGRD2
      INTEGER LU
      REAL UNDEF
      PARAMETER (LU=1,UNDEF=-999999.)
C     Input data:
      INTEGER N1,N2,N3,N4,NEWX1,NEWX2,NEWX3,NEWX4
C     Other variables:
      INTEGER NEWN1,NEWN2,NEWN3,NEWN4,ITRANS(4),I1,I2,I3,I4,I,J3,J4,J
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+GRDTRANS: Enter input filename: '
      FSEP=' '
      READ(*,*) FSEP
      WRITE(*,'(A)') '+GRDTRANS: Working ...           '
C
C     Reading all data from the SEP file into the memory:
      IF (FSEP.NE.' ') THEN
        CALL RSEP1(LU,FSEP)
      ELSE
C       GRDTRANS-01
        CALL ERROR('GRDTRANS-01: 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'   ,FGRD1,'grd.out'   )
      CALL RSEP3T('GRDNEW',FGRD2,'grdnew.out')
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3I('N4',N4,1)
      IF(2*N1*N2*N3*N4.GT.MRAM) THEN
C       GRDTRANS-02
        CALL ERROR('GRDTRANS-02: Too small array RAM(MRAM)')
C       Array RAM(MRAM) allocated in include file 'ram.inc' is too small
C       to contain two input grids (2*N1*N2*N3*N4 values).  You may wish
C       to increase the dimension MRAM in file 'ram.inc'.
C       ram.inc
      END IF
      CALL RSEP3I('NEWX1',NEWX1,1)
      CALL RSEP3I('NEWX2',NEWX2,2)
      CALL RSEP3I('NEWX3',NEWX3,3)
      CALL RSEP3I('NEWX4',NEWX4,4)
      IF(NEWX1.LT.1.OR.4.LT.NEWX1) THEN
C       GRDTRANS-03
        CALL ERROR('GRDTRANS-03: Incorrect value of parameter NEWX1')
      END IF
      IF(NEWX2.LT.1.OR.4.LT.NEWX2) THEN
C       GRDTRANS-04
        CALL ERROR('GRDTRANS-04: Incorrect value of parameter NEWX2')
      END IF
      IF(NEWX3.LT.1.OR.4.LT.NEWX3) THEN
C       GRDTRANS-05
        CALL ERROR('GRDTRANS-05: Incorrect value of parameter NEWX3')
      END IF
      IF(NEWX4.LT.1.OR.4.LT.NEWX4) THEN
C       GRDTRANS-06
        CALL ERROR('GRDTRANS-06: Incorrect value of parameter NEWX4')
      END IF
      IF(NEWX1.EQ.NEWX2.OR.NEWX1.EQ.NEWX3.OR.NEWX2.EQ.NEWX3.OR.
     *   NEWX1.EQ.NEWX4.OR.NEWX2.EQ.NEWX4.OR.NEWX3.EQ.NEWX4) THEN
C       GRDTRANS-07
        CALL ERROR('GRDTRANS-07: Coinciding output axes')
      END IF
C
      ITRANS(1)=N1
      ITRANS(2)=N2
      ITRANS(3)=N3
      ITRANS(4)=N4
      NEWN1=ITRANS(NEWX1)
      NEWN2=ITRANS(NEWX2)
      NEWN3=ITRANS(NEWX3)
      NEWN4=ITRANS(NEWX4)
      ITRANS(1)=1
      ITRANS(2)=N1
      ITRANS(3)=N1*N2
      ITRANS(4)=N1*N2*N3
      NEWX1=ITRANS(NEWX1)
      NEWX2=ITRANS(NEWX2)
      NEWX3=ITRANS(NEWX3)
      NEWX4=ITRANS(NEWX4)
C
C     Reading input grid:
      CALL RARAY(LU,FGRD1,'FORMATTED',.TRUE.,UNDEF,N1*N2*N3,N4,RAM)
C
C     Tranposing the grid:
      I=N1*N2*N3*N4
      DO 14 I4=0,NEWN4-1
        J4=1+I4*NEWX4
        DO 13 I3=0,NEWN3-1
          J3=J4+I3*NEWX3
          DO 12 I2=0,NEWN2-1
            J=J3+I2*NEWX2
            DO 11 I1=0,NEWN1-1
              I=I+1
              RAM(I)=RAM(J)
              J=J+NEWX1
   11       CONTINUE
   12     CONTINUE
   13   CONTINUE
   14 CONTINUE
C
C     Writing output grid:
      CALL WARAY(LU,FGRD2,'FORMATTED',.TRUE.,UNDEF,.FALSE.,0.,N1*N2*N3,
     *                                            N4,RAM(N1*N2*N3*N4+1))
      WRITE(*,'(A)') '+GRDTRANS: 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