C 
C Subroutine file 'forms.for' to facilitate writing and reading data.
C
C Version: 5.50
C Date: 2000, November 25
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 This file consists of the following external procedures:
C     OMAT... Subroutine designed to set the form of the file for reading
C             or writing a matrix, and to open the file.
C             OMAT
C     WMAT... Subroutine designed to write a given matrix into the given
C             file.
C             WMAT
C     RMAT... Subroutine designed to read a matrix from the given file.
C             RMAT
C     UARRAY..Function returning the undefined value used in the
C             unformatted files with real arrays.
C             UARRAY
C     WARRAY..Subroutine designed to write a given real array into the
C             given formatted or unformatted file.
C             WARRAY
C     WARRAI..Subroutine designed to write a given integer array into
C             the given formatted or unformatted file.
C             WARRAI
C     RARRAY..Subroutine designed to read the real array from the given
C             formatted or unformatted file.
C             RARRAY
C     RARRAI..Subroutine designed to read the integer array from the
C             given formatted or unformatted file.
C             RARRAI
C     WARAY...Subroutine calling WARRAY for N4 individual time levels.
C             WARAY
C     WARAI...Subroutine calling WARRAI for N4 individual time levels.
C             WARAI
C     RARAY...Subroutine calling RARRAY for N4 individual time levels.
C             RARAY
C     RARAI...Subroutine calling RARRAI for N4 individual time levels.
C             RARAI
C     FORM1...Subroutine designed to determine the best output format
C             for reals.
C             FORM1
C     FORM2...Subroutine designed to determine the best output format
C             for multiples of real numbers.
C             FORM2
C
C=======================================================================
C
C     
C
      SUBROUTINE OMAT(LU,FILE,IRW,FORMM)
      CHARACTER*(*) FILE,FORMM
      INTEGER LU,IRW
C
C Subroutine designed to set the form FORMM of the file with matrix
C to be read or written, and to open the file, if FILE is specified.
C
C Input:
C     LU...   Logical unit number to be used for the output.
C     FILE... Destination filename.  If not blank, the file will be
C             opened.
C     IRW ... Identifies, whether the file will be read or written:
C             IRW=1 ... reading
C             IRW=2 ... writing
C Output:
C     FORMM...Form of the file to be read or written.
C
C Date: 2000, October 20
C Coded by Petr Bulant
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
C
      CHARACTER*13 FORMMR,FORMMW
      SAVE FORMMR,FORMMW
      DATA  FORMMR/'undefined'/
C
C     FORMMR ..Form of the files with matrices to be read.
C     FORMMW ..Form of the files with matrices to be written.
C
C.......................................................................
C
      IF (FORMMR.EQ.'undefined') THEN
C       Reading the forms of the files with matrices:
        CALL RSEP3T('FORMM',FORMM,'formatted')
        CALL RSEP3T('FORMMR',FORMMR,FORMM)
        CALL RSEP3T('FORMMW',FORMMW,FORMM)
        CALL LOWER(FORMM)
        CALL LOWER(FORMMR)
        CALL LOWER(FORMMW)
        IF ((FORMM.NE.'formatted').AND.(FORMM.NE.'unformatted').OR.
     *      (FORMMR.NE.'formatted').AND.(FORMMR.NE.'unformatted').OR.
     *      (FORMMW.NE.'formatted').AND.(FORMMW.NE.'unformatted')) THEN
C         FORMS-01
          CALL ERROR('FORMS-01: Wrong value of FORMM, FORMMR or FORMMW')
C         Input parameters FORMM, FORMMR and FORMMW, if specified,
C         must equal either 'formatted' or 'unformatted'.
        ENDIF
      ENDIF
C
C     Setting the form FORMM of the file to be opened:
      IF (IRW.EQ.1) THEN
        FORMM=FORMMR
      ELSEIF (IRW.EQ.2) THEN
        FORMM=FORMMW
      ELSE
C       FORMS-02
        CALL ERROR('FORMS-02: Wrong value of IRW')
C       Dumy argument IRW must equal either 1 or 2.
      ENDIF
C
      IF (FILE.NE.' ') THEN
C       Opening the file for reading or writing:
        IF (IRW.EQ.1) THEN
          WRITE(*,'(2A)') '+Reading: ',FILE(1:MIN0(LEN(FILE),70))
          OPEN(LU,FILE=FILE,FORM=FORMM,STATUS='OLD')
        ELSEIF (IRW.EQ.2) THEN
          WRITE(*,'(2A)') '+Writing: ',FILE(1:MIN0(LEN(FILE),70))
          OPEN(LU,FILE=FILE,FORM=FORMM)
        ENDIF
      ENDIF
C
      RETURN
      END
C=======================================================================
C
C     
C
      SUBROUTINE WMAT(LU,FILE,M1,M2,OUT)
      CHARACTER*(*) FILE
      INTEGER LU,M1,M2
      REAL OUT(*)
C
C Subroutine designed to write a given matrix into the file.
C
C Input:
C     LU...   Logical unit number to be used for the output.
C     FILE... Destination filename.  If not blank, the file will be
C             opened and closed.  If blank, the file is assumed to be
C             already open, and will not be closed in this subroutine.
C     M1...   Number of rows of the given matrix.
C     M2...   M2=0 for a symmetric matrix,
C             M2=1 for a diagonal matrix,
C             M2=number of columns for a general matrix.
C     OUT...  Components of the given matrix stored columnwise.
C             For a symmetric matrix, just components from the first row
C             to the diagonal are stored for each column, i.e., array
C             OUT has M1*(M1+1)/2 matrix components.
C             For a diagonal matrix, just M1 diagonal components are
C             stored.
C
C No output.
C
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
C
      CHARACTER*13 FORMAT,FORMM
      INTEGER I1,I2
C
C     FORMAT..String containing the output format.
C     FORMM ..Form of the files with matrices.
C     I1,I2.. Loop variables.
C
C.......................................................................
C
C     Setting output format:
      FORMAT='(5(G13.7,1X))'
C
C     Form of the file with the matrix, opening the file:
      CALL OMAT(LU,FILE,2,FORMM)
C
C     Writing the matrix:
      IF(M2.LE.0) THEN
C       Symmetric matrix
          IF (FORMM.EQ.'formatted') THEN
            DO 11 I2=1,M1
              WRITE(LU,FORMAT) (OUT(I1),I1=I2*(I2-1)/2+1,I2*(I2+1)/2)
  11        CONTINUE
          ELSE
            WRITE(LU) (OUT(I1),I1=1,M1*(M1+1)/2)
          ENDIF
      ELSE
C       Diagonal or general matrix
        IF (FORMM.EQ.'formatted') THEN
          DO 12 I2=M1,M1*M2,M1
            WRITE(LU,FORMAT) (OUT(I1),I1=I2-M1+1,I2)
  12      CONTINUE
        ELSE
          WRITE(LU) (OUT(I1),I1=1,M1*M2)
        ENDIF
      END IF
C
C     Closing output file:
      IF(FILE.NE.' ') THEN
        CLOSE(LU)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RMAT(LU,FILE,M1,M2,ARRAY)
      CHARACTER*(*) FILE
      INTEGER LU,M1,M2
      REAL ARRAY(*)
C
C Subroutine designed to read a matrix from the file.
C
C Input:
C     LU...   Logical unit number to be used for the input.
C     FILE... Destination filename.  If not blank, the file will be
C             opened and closed.  If blank, the file is assumed to be
C             already open, and will not be closed in this subroutine.
C     M1...   Number of rows of the matrix.
C     M2...   M2=0 for a symmetric matrix,
C             M2=1 for a diagonal matrix,
C             M2=number of columns for a general matrix.
C
C Output:
C     ARRAY...Components of the given matrix stored columnwise.
C             For a symmetric matrix, just components from the first row
C             to the diagonal are stored for each column, i.e., array
C             out has M1*(M1+1)/2 matrix components.
C             For a diagonal matrix, just M1 diagonal components are
C             stored.
C
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Local storage location:
      CHARACTER*13 FORMM
      INTEGER I
C     I...    Loop variable.
C
C.......................................................................
C
C     Form of the file with the matrix, opening the file:
      CALL OMAT(LU,FILE,1,FORMM)
C
C     Reading the matrix:
      IF(M2.LE.0) THEN
C       Symmetric matrix
        IF (FORMM.EQ.'formatted') THEN
          READ(LU,*) (ARRAY(I),I=1,M1*(M1+1)/2)
        ELSE
          READ(LU)   (ARRAY(I),I=1,M1*(M1+1)/2)
        ENDIF
      ELSE
C       Diagonal or general matrix
        IF (FORMM.EQ.'formatted') THEN
          READ(LU,*) (ARRAY(I),I=1,M1*M2)
        ELSE
          READ(LU)   (ARRAY(I),I=1,M1*M2)
        ENDIF
      END IF
C
C     Closing input file:
      IF(FILE.NE.' ') THEN
        CLOSE(LU)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      REAL FUNCTION UARRAY()
C
C Function returning the undefined value used in the unformatted files
C with real-valued arrays.
C
C No input.
C
C Output:
C     UARRAY..The value used as "undefined value" in the unformatted
C             files with real-valued arrays by subroutines WARRAY and
C             RARRAY.
C
C Date: 2000, November 25
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Parameters:
      REAL UNDEF
      PARAMETER (UNDEF=-999999999.)
C
      UARRAY=UNDEF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WARRAY(LU,FILE,FORM,LMIN,VMIN,LMAX,VMAX,NOUT,OUT)
      CHARACTER*(*) FILE,FORM
      LOGICAL LMIN,LMAX
      INTEGER LU,NOUT
      REAL VMIN,VMAX,OUT(NOUT)
C
C Subroutine designed to write a given real array into the file.
C
C Input:
C     LU...   Logical unit number to be used for the output.
C     FILE... Destination filename.  If not blank, the file will be
C             opened and closed.  If blank, the file is assumed to be
C             already open, and will not be closed in this subroutine.
C     FORM... Form of the output file: either 'FORMATTED' or
C             'UNFORMATTED'.
C     LMIN... TRUE if the null values are to be written in place of
C             array elements less than or equal to VMIN, otherwise
C             FALSE.
C             Formatted output:
C               The null values are treated as default values when read
C               by list-directed input (free format).
C               Example: 124 null values are written as '    124*'.
C             Unformatted output:
C               The values of -999999999 are written in place of the
C               null values.
C     VMIN... Trade-off limit.
C     LMAX... TRUE if the null values are to be written in place of
C             array elements greater than or equal to VMAX, otherwise
C             FALSE.
C     VMAX... Trade-off limit.
C     NOUT... Dimension of the array OUT.
C     OUT...  Array to be written.
C
C No output.
C
C Date: 2000, November 25
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Parameters:
      EXTERNAL UARRAY
      REAL     UARRAY
      REAL     UNDEF
C
C     Local storage locations:
      CHARACTER*11 FORML
      CHARACTER*14 FORMAT
      INTEGER IMIN,IADR
      REAL OUTMIN,OUTMAX,VMINA,VMAXA
C     FORMAT..String containing the output format, e.g. like (10F8.3).
C     IMIN... Loop lower bound, locally also loop variable.
C     IADR... Loop variable.
C     OUTMIN,OUTMAX... Minimum and maximum defined element to determine
C             the best format for printing.
C     VMINA,VMAXA... Local storage locations for VMIN, VMAX.
C
C.......................................................................
C
      UNDEF=UARRAY()
C
      IF(FILE.NE.' ') THEN
        WRITE(*,'(''+'',79('' ''))')
        WRITE(*,'(2A)') '+Writing: ',FILE(1:MIN0(LEN(FILE),70))
        OPEN(LU,FILE=FILE,FORM=FORM)
      END IF
C
C     Formatted or unformatted output:
      FORML=FORM
      CALL LOWER(FORML)
      IF(FORML.EQ.'formatted') THEN
C
C       Minimum and maximum elements:
        OUTMIN=0.
        IF(LMIN) THEN
          VMINA=VMIN
          DO 11 IADR=1,NOUT
            IF(OUTMIN.GT.OUT(IADR)) THEN
              IF(OUT(IADR).GT.VMINA) THEN
                OUTMIN=OUT(IADR)
              END IF
            END IF
   11     CONTINUE
        ELSE
          DO 12 IADR=1,NOUT
            IF(OUTMIN.GT.OUT(IADR)) THEN
              OUTMIN=OUT(IADR)
            END IF
   12     CONTINUE
        END IF
        OUTMAX=0.
        IF(LMAX) THEN
          VMAXA=VMAX
          DO 13 IADR=1,NOUT
            IF(OUTMAX.LT.OUT(IADR)) THEN
              IF(OUT(IADR).LT.VMAXA) THEN
                OUTMAX=OUT(IADR)
              END IF
            END IF
   13     CONTINUE
        ELSE
          DO 14 IADR=1,NOUT
            IF(OUTMAX.LT.OUT(IADR)) THEN
              OUTMAX=OUT(IADR)
            END IF
   14     CONTINUE
        END IF
C
C       Setting output format for the array:
        FORMAT='(10(F00.0,1X))'
        CALL FORM1(OUTMIN,OUTMAX,FORMAT(5:12))
        FORMAT(11:14)=   '1X))'
C       Output format is set.
C
C       Printing loop:
C       Initial value of the first element to print
        IADR=1
C       Beginning of the loop
   20   CONTINUE
C
C         Trade off (searching for undefined elements):
          IMIN=IADR
          IF(LMIN) THEN
            IF(LMAX) THEN
              DO 21 IADR=IMIN,NOUT
                IF(OUT(IADR).LE.VMINA.OR.OUT(IADR).GE.VMAXA) THEN
                  GO TO 29
                END IF
   21         CONTINUE
            ELSE
              DO 22 IADR=IMIN,NOUT
                IF(OUT(IADR).LE.VMINA) THEN
                  GO TO 29
                END IF
   22         CONTINUE
            END IF
          ELSE
            IF(LMAX) THEN
              DO 23 IADR=IMIN,NOUT
                IF(OUT(IADR).GE.VMAXA) THEN
                  GO TO 29
                END IF
   23         CONTINUE
            ELSE
              IADR=NOUT+1
            END IF
          END IF
   29     CONTINUE
C         IADR is the first undefined element.
C
C         Writing the array (defined elements):
          IF(IMIN.EQ.1.AND.IADR.GT.NOUT) THEN
            WRITE(LU,FORMAT) OUT
            GO TO 90
          ELSE
            WRITE(LU,FORMAT) (OUT(IMIN),IMIN=IMIN,IADR-1)
            IF(IADR.GT.NOUT) THEN
              GO TO 90
            END IF
          END IF
C
C         Searching for the next defined elements:
          IMIN=IADR
          IF(LMIN) THEN
            IF(LMAX) THEN
              DO 31 IADR=IADR,NOUT
                IF(OUT(IADR).GT.VMINA.AND.OUT(IADR).LT.VMAXA) THEN
                  GO TO 39
                END IF
   31         CONTINUE
            ELSE
              DO 32 IADR=IADR,NOUT
                IF(OUT(IADR).GT.VMINA) THEN
                  GO TO 39
                END IF
   32         CONTINUE
            END IF
          ELSE
            IF(LMAX) THEN
              DO 33 IADR=IADR,NOUT
                IF(OUT(IADR).LT.VMAXA) THEN
                  GO TO 39
                END IF
   33         CONTINUE
            ELSE
              IADR=NOUT+1
            END IF
          END IF
   39     CONTINUE
C         IADR is the first defined element.
C
C         Writing the array (undefined elements):
          WRITE(LU,'(I7,A)') IADR-IMIN,'*'
          IF(NOUT.LT.IADR) THEN
            GO TO 90
          END IF
C
        GO TO 20
      ELSE
C
C       Null values:
        IF(LMIN) THEN
          VMINA=VMIN
          IF(LMAX) THEN
            VMAXA=VMAX
            DO 51 IADR=1,NOUT
              IF(OUT(IADR).LE.VMINA.OR.VMAXA.LE.OUT(IADR)) THEN
                OUT(IADR)=UNDEF
              END IF
   51       CONTINUE
          ELSE
            DO 52 IADR=1,NOUT
              IF(OUT(IADR).LE.VMINA) THEN
                OUT(IADR)=UNDEF
              END IF
   52       CONTINUE
          END IF
        ELSE
          IF(LMAX) THEN
            VMAXA=VMAX
            DO 53 IADR=1,NOUT
              IF(VMAXA.LE.OUT(IADR)) THEN
                OUT(IADR)=UNDEF
              END IF
   53       CONTINUE
          END IF
        END IF
C
C       Writing the array:
        WRITE(LU) OUT
C
      END IF
   90 CONTINUE
      IF(FILE.NE.' ') THEN
        CLOSE(LU)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WARRAI(LU,FILE,FORM,LMIN,IVMIN,LMAX,IVMAX,NOUT,IOUT)
      CHARACTER*(*) FILE,FORM
      LOGICAL LMIN,LMAX
      INTEGER LU,NOUT,IVMIN,IVMAX,IOUT(NOUT)
C
C Subroutine designed to write a given integer array into the file.
C
C Input:
C     LU...   Logical unit number to be used for the output.
C     FILE... Destination filename.  If not blank, the file will be
C             opened and closed.  If blank, the file is assumed to be
C             already open, and will not be closed in this subroutine.
C     FORM... Form of the output file: either 'FORMATTED' or
C             'UNFORMATTED'.
C     LMIN... TRUE if the null values are to be written in place of
C             array elements less than or equal to IVMIN, otherwise
C             FALSE.
C             Formatted output:
C               The null values are treated as default values when read
C               by list-directed input (free format).
C               Example: 124 null values are written as '    124*'.
C             Unformatted output:
C               The values of -999999999 are written in place of the
C               null values.
C     IVMIN . Trade-off limit.
C     LMAX... TRUE if the null values are to be written in place of
C             array elements greater than or equal to IVMAX, otherwise
C             FALSE.
C     IVMAX...Trade-off limit.
C     NOUT... Dimension of the array IOUT.
C     IOUT... Array to be written.
C
C No output.
C
C Date: 2000, November 25
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Parameters:
      INTEGER IUNDEF
      PARAMETER (IUNDEF=-999999999)
C
C     Local storage locations:
      CHARACTER*11 FORML
      CHARACTER*12 FORMAT
      INTEGER IMIN,IADR,MINOUT,MAXOUT,IVMINA,IVMAXA
C     FORMAT..String containing the output format, e.g. like (10I08).
C     IMIN... Loop lower bound, locally also loop variable.
C     IADR... Loop variable.
C     MINOUT,MAXOUT... Minimum and maximum defined element to determine
C             the best format for printing.
C     IVMINA,IVMAXA... Local storage locations for IVMIN, IVMAX.
C
C.......................................................................
C
      IF(FILE.NE.' ') THEN
        WRITE(*,'(''+'',79('' ''))')
        WRITE(*,'(2A)') '+Writing: ',FILE(1:MIN0(LEN(FILE),70))
        OPEN(LU,FILE=FILE,FORM=FORM)
      END IF
C
C     Formatted or unformatted output:
      FORML=FORM
      CALL LOWER(FORML)
      IF(FORML.EQ.'formatted') THEN
C
C       Minimum and maximum elements:
        MINOUT=0
        IF(LMIN) THEN
          IVMINA=IVMIN
          DO 11 IADR=1,NOUT
            IF(MINOUT.GT.IOUT(IADR)) THEN
              IF(IOUT(IADR).GT.IVMINA) THEN
                MINOUT=IOUT(IADR)
              END IF
            END IF
   11     CONTINUE
        ELSE
          DO 12 IADR=1,NOUT
            IF(MINOUT.GT.IOUT(IADR)) THEN
              MINOUT=IOUT(IADR)
            END IF
   12     CONTINUE
        END IF
        MAXOUT=0
        IF(LMAX) THEN
          IVMAXA=IVMAX
          DO 13 IADR=1,NOUT
            IF(MAXOUT.LT.IOUT(IADR)) THEN
              IF(IOUT(IADR).LT.IVMAXA) THEN
                MAXOUT=IOUT(IADR)
              END IF
            END IF
   13     CONTINUE
        ELSE
          DO 14 IADR=1,NOUT
            IF(MAXOUT.LT.IOUT(IADR)) THEN
              MAXOUT=IOUT(IADR)
            END IF
   14     CONTINUE
        END IF
C
C       Setting output format for the array:
        FORMAT='(10(I00,1X))'
        IMIN=MAXOUT
        IF(MINOUT.LT.0.) THEN
          IMIN=MAX0(IMIN,-10*MINOUT)
        END IF
        DO 15 IADR=1,99
          IMIN=IMIN/10
          IF(IMIN.LT.1) THEN
            FORMAT(6:6)=CHAR(ICHAR('0')+IADR/10)
            FORMAT(7:7)=CHAR(ICHAR('0')+MOD(IADR,10))
            GO TO 16
          END IF
   15   CONTINUE
   16   CONTINUE
C       Output format is set.
C
C       Printing loop:
C       Initial value of the first element to print
        IADR=1
C       Beginning of the loop
   20   CONTINUE
C
C         Trade off (searching for undefined elements):
          IMIN=IADR
          IF(LMIN) THEN
            IF(LMAX) THEN
              DO 21 IADR=IMIN,NOUT
                IF(IOUT(IADR).LE.IVMINA.OR.IOUT(IADR).GE.IVMAXA) THEN
                  GO TO 29
                END IF
   21         CONTINUE
            ELSE
              DO 22 IADR=IMIN,NOUT
                IF(IOUT(IADR).LE.IVMINA) THEN
                  GO TO 29
                END IF
   22         CONTINUE
            END IF
          ELSE
            IF(LMAX) THEN
              DO 23 IADR=IMIN,NOUT
                IF(IOUT(IADR).GE.IVMAXA) THEN
                  GO TO 29
                END IF
   23         CONTINUE
            ELSE
              IADR=NOUT+1
            END IF
          END IF
   29     CONTINUE
C         IADR is the first undefined element.
C
C         Writing the array (defined elements):
          IF(IMIN.EQ.1.AND.IADR.GT.NOUT) THEN
            WRITE(LU,FORMAT) IOUT
            GO TO 90
          ELSE
            WRITE(LU,FORMAT) (IOUT(IMIN),IMIN=IMIN,IADR-1)
            IF(IADR.GT.NOUT) THEN
              GO TO 90
            END IF
          END IF
C
C         Searching for the next defined elements:
          IMIN=IADR
          IF(LMIN) THEN
            IF(LMAX) THEN
              DO 31 IADR=IADR,NOUT
                IF(IOUT(IADR).GT.IVMINA.AND.IOUT(IADR).LT.IVMAXA) THEN
                  GO TO 39
                END IF
   31         CONTINUE
            ELSE
              DO 32 IADR=IADR,NOUT
                IF(IOUT(IADR).GT.IVMINA) THEN
                  GO TO 39
                END IF
   32         CONTINUE
            END IF
          ELSE
            IF(LMAX) THEN
              DO 33 IADR=IADR,NOUT
                IF(IOUT(IADR).LT.IVMAXA) THEN
                  GO TO 39
                END IF
   33         CONTINUE
            ELSE
              IADR=NOUT+1
            END IF
          END IF
   39     CONTINUE
C         IADR is the first defined element.
C
C         Writing the array (undefined elements):
          WRITE(LU,'(I7,A)') IADR-IMIN,'*'
          IF(NOUT.LT.IADR) THEN
            GO TO 90
          END IF
C
        GO TO 20
      ELSE
C
C       Null values:
        IF(LMIN) THEN
          IVMINA=IVMIN
          IF(LMAX) THEN
            IVMAXA=IVMAX
            DO 51 IADR=1,NOUT
              IF(IOUT(IADR).LE.IVMINA.OR.IVMAXA.LE.IOUT(IADR)) THEN
                IOUT(IADR)=IUNDEF
              END IF
   51       CONTINUE
          ELSE
            DO 52 IADR=1,NOUT
              IF(IOUT(IADR).LE.IVMINA) THEN
                IOUT(IADR)=IUNDEF
              END IF
   52       CONTINUE
          END IF
        ELSE
          IF(LMAX) THEN
            IVMAXA=IVMAX
            DO 53 IADR=1,NOUT
              IF(IVMAXA.LE.IOUT(IADR)) THEN
                IOUT(IADR)=IUNDEF
              END IF
   53       CONTINUE
          END IF
        END IF
C
C       Writing the array:
        WRITE(LU) IOUT
C
      END IF
   90 CONTINUE
      IF(FILE.NE.' ') THEN
        CLOSE(LU)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RARRAY(LU,FILE,FORM,LDEF,DEF,N,ARRAY)
      CHARACTER*(*) FILE,FORM
      LOGICAL LDEF
      INTEGER LU,N
      REAL DEF,ARRAY(N)
C
C Subroutine designed to read the real array from the disk.
C
C Input:
C     LU...   Logical unit number to be used.
C     FILE... Source filename.  If not blank, the file will be
C             opened and closed.  If blank, the file is assumed to be
C             already open, and will not be closed in this subroutine.
C     FORM... Form of the output file: either 'FORMATTED' or
C             'UNFORMATTED'.
C     LDEF... True if the null values are to be replaced by the given
C             default value DEF.
C             If FORM='FORMATTED' and LDEF=.FALSE., the array elements
C             corresponding to null values remain unchanged.
C     DEF...  Default value.
C     N...    Array dimension (number of elements to read).
C
C Output:
C     ARRAY.. Array having been read.
C
C Date: 2000, November 25
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Parameters:
      EXTERNAL UARRAY
      REAL     UARRAY
      REAL UNDEF
C
      CHARACTER*11 FORML
      INTEGER I
      REAL AUX
C
      UNDEF=UARRAY()
C
      IF(FILE.NE.' ') THEN
        WRITE(*,'(''+'',79('' ''))')
        WRITE(*,'(2A)') '+Reading: ',FILE(1:MIN0(LEN(FILE),70))
        OPEN(LU,FILE=FILE,FORM=FORM,STATUS='OLD')
      END IF
C
      FORML=FORM
      CALL LOWER(FORML)
      IF(FORML.EQ.'formatted') THEN
        IF(LDEF) THEN
          AUX=DEF
          DO 10 I=1,N
            ARRAY(I)=AUX
   10     CONTINUE
        END IF
        READ(LU,*) ARRAY
      ELSE
        READ(LU) ARRAY
        IF(LDEF) THEN
          IF(DEF.NE.UNDEF) THEN
            AUX=DEF
            DO 20 I=1,N
              IF(ARRAY(I).EQ.UNDEF) THEN
                ARRAY(I)=AUX
              END IF
   20       CONTINUE
          END IF
        END IF
      END IF
C
      IF(FILE.NE.' ') THEN
        CLOSE(LU)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RARRAI(LU,FILE,FORM,LDEF,IDEF,N,IARRAY)
      CHARACTER*(*) FILE,FORM
      LOGICAL LDEF
      INTEGER LU,IDEF,N,IARRAY(N)
C
C Subroutine designed to read the integer array from the disk.
C
C Input:
C     LU...   Logical unit number to be used.
C     FILE... Source filename.  If not blank, the file will be
C             opened and closed.  If blank, the file is assumed to be
C             already open, and will not be closed in this subroutine.
C     FORM... Form of the output file: either 'FORMATTED' or
C             'UNFORMATTED'.
C     LDEF... True if the null values are to be replaced by the given
C             default value IDEF.
C             If FORM='FORMATTED' and LDEF=.FALSE., the array elements
C             corresponding to null values remain unchanged.
C     IDEF... Default value.
C     N...    Array dimension (number of elements to read).
C
C Output:
C     IARRAY..Array having been read.
C
C Date: 2000, November 25
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Parameters:
      INTEGER IUNDEF
      PARAMETER (IUNDEF=-999999999)
C
      CHARACTER*11 FORML
      INTEGER I,IAUX
C
      IF(FILE.NE.' ') THEN
        WRITE(*,'(''+'',79('' ''))')
        WRITE(*,'(2A)') '+Reading: ',FILE(1:MIN0(LEN(FILE),70))
        OPEN(LU,FILE=FILE,FORM=FORM,STATUS='OLD')
      END IF
C
      FORML=FORM
      CALL LOWER(FORML)
      IF(FORML.EQ.'formatted') THEN
        IF(LDEF) THEN
          IAUX=IDEF
          DO 10 I=1,N
            IARRAY(I)=IAUX
   10     CONTINUE
        END IF
        READ(LU,*) IARRAY
      ELSE
        READ(LU) IARRAY
        IF(LDEF) THEN
C         IF(IDEF.NE.IUNDEF) THEN
            IAUX=IDEF
            DO 20 I=1,N
              IF(IARRAY(I).EQ.IUNDEF) THEN
                IARRAY(I)=IAUX
              END IF
   20       CONTINUE
C         END IF
        END IF
      END IF
C
      IF(FILE.NE.' ') THEN
        CLOSE(LU)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WARAY(LU,FILE,FORM,LMIN,VMIN,LMAX,VMAX,NOUT,N4,OUT)
      CHARACTER*(*) FILE,FORM
      LOGICAL LMIN,LMAX
      INTEGER LU,NOUT,N4
      REAL VMIN,VMAX,OUT(NOUT,N4)
C
C Subroutine designed to N4 times call subroutine WARRAY, for individual
C time levels.
C
C Input:
C     LU,FILE,FORM,LMIN,VMIN,LMAX,VMAX,NOUT... Refer to subroutine
C             WARRAY
C     N4...   Number of time levels.  NOUT values corresponding to each
C             level are written through an individual invocation of
C             subroutine WARRAY.
C     OUT...  Array of dimension NOUT*N4 to be written.
C
C No output.
C
C Date: 1998, March 21
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      INTEGER I4
C
C.......................................................................
C
      IF(FILE.NE.' ') THEN
        WRITE(*,'(''+'',79('' ''))')
        WRITE(*,'(2A)') '+Writing: ',FILE(1:MIN0(LEN(FILE),70))
        OPEN(LU,FILE=FILE,FORM=FORM)
      END IF
C
      DO 10 I4=1,N4
        CALL WARRAY(LU,' ',FORM,LMIN,VMIN,LMAX,VMAX,NOUT,OUT(1,I4))
   10 CONTINUE
C
      IF(FILE.NE.' ') THEN
        CLOSE(LU)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WARAI(LU,FILE,FORM,LMIN,IVMIN,LMAX,IVMAX,NOUT,N4,IOUT)
      CHARACTER*(*) FILE,FORM
      LOGICAL LMIN,LMAX
      INTEGER LU,IVMIN,IVMAX,NOUT,N4,IOUT(NOUT,N4)
C
C Subroutine designed to N4 times call subroutine WARRAI, for individual
C time levels.
C
C Input:
C     LU,FILE,FORM,LMIN,IVMIN,LMAX,IVMAX,NOUT... Refer to subroutine
C             WARRAI
C     N4...   Number of time levels.  NOUT values corresponding to each
C             level are written through an individual invocation of
C             subroutine WARRAI.
C     IOUT... Array of dimension NOUT*N4 to be written.
C
C No output.
C
C Date: 1998, May 28
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      INTEGER I4
C
C.......................................................................
C
      IF(FILE.NE.' ') THEN
        WRITE(*,'(''+'',79('' ''))')
        WRITE(*,'(2A)') '+Writing: ',FILE(1:MIN0(LEN(FILE),70))
        OPEN(LU,FILE=FILE,FORM=FORM)
      END IF
C
      DO 10 I4=1,N4
        CALL WARRAI(LU,' ',FORM,LMIN,IVMIN,LMAX,IVMAX,NOUT,IOUT(1,I4))
   10 CONTINUE
C
      IF(FILE.NE.' ') THEN
        CLOSE(LU)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RARAY(LU,FILE,FORM,LDEF,DEF,N,N4,ARRAY)
      CHARACTER*(*) FILE,FORM
      LOGICAL LDEF
      INTEGER LU,N,N4
      REAL DEF,ARRAY(N,N4)
C
C Subroutine designed to N4 times call subroutine RARRAY, for individual
C time levels.
C
C Input:
C     LU,FILE,FORM,LDEF,DEF,N... Refer to subroutine
C             RARRAY
C     N4...   Number of time levels.  N values corresponding to each
C             level are read by an individual invocation of subroutine
C             RARRAY.
C
C Output:
C     ARRAY...Array of dimension N*N4 having been read.
C
C Date: 2000, July 31
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      INTEGER I4
C
C.......................................................................
C
      IF(FILE.NE.' ') THEN
        WRITE(*,'(''+'',79('' ''))')
        WRITE(*,'(2A)') '+Reading: ',FILE(1:MIN0(LEN(FILE),70))
        OPEN(LU,FILE=FILE,FORM=FORM,STATUS='OLD')
      END IF
C
      DO 10 I4=1,N4
        CALL RARRAY(LU,' ',FORM,LDEF,DEF,N,ARRAY(1,I4))
   10 CONTINUE
C
      IF(FILE.NE.' ') THEN
        CLOSE(LU)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RARAI(LU,FILE,FORM,LDEF,IDEF,N,N4,IARRAY)
      CHARACTER*(*) FILE,FORM
      LOGICAL LDEF
      INTEGER LU,N,IDEF,N4,IARRAY(N,N4)
C
C Subroutine designed to N4 times call subroutine RARRAI, for individual
C time levels.
C
C Input:
C     LU,FILE,FORM,LDEF,IDEF,N,N4... Refer to subroutine
C             RARRAI
C
C Output:
C     IARRAY..Array of dimension N*N4 having been read.
C
C Date: 2000, July 31
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      INTEGER I4
C
C.......................................................................
C
      IF(FILE.NE.' ') THEN
        WRITE(*,'(''+'',79('' ''))')
        WRITE(*,'(2A)') '+Reading: ',FILE(1:MIN0(LEN(FILE),70))
        OPEN(LU,FILE=FILE,FORM=FORM,STATUS='OLD')
      END IF
C
      DO 10 I4=1,N4
        CALL RARRAI(LU,' ',FORM,LDEF,IDEF,N,IARRAY(1,I4))
   10 CONTINUE
C
      IF(FILE.NE.' ') THEN
        CLOSE(LU)
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE FORM1(OUTMIN,OUTMAX,FORMAT)
      REAL OUTMIN,OUTMAX
      CHARACTER*8 FORMAT
C
C Subroutine designed to determine the best output format for reals.
C
C Input:
C     OUTMIN,OUTMAX... Minimum and maximum real number to be written.
C
C Output:
C     FORMAT..String containing the output format e.g. like 'F07.3,A,'.
C             The width of the defined string is 8 characters.
C             It has the form of 'F00.0,A,', where zeros are replaced
C             by reasonable values.  The subroutine attempts to output
C             at least MAXDIG digits (including all zeros after the
C             decimal point) of the largest positive number OUTMAX and
C             MAXDIG-1 digits of the most negative number if OUTMIN is
C             negative, and to adjust the width of the output field to
C             MAXDIG+1 columns, if possible.  The (MAXDIG+2)th column is
C             reserved for a space or another separator.
C             If OUTMIN=0 and OUTMAX=0, the width of the output field is
C             adjusted to 2 columns.
C             If the number of digits (without leading zeros) is smaller
C             than MINDIG, format Fnn.d with nn=MAXDIG+1 is changed to
C             Gmm.d with mm=MAXDIG+5.  The (MAXDIG+6)th column is
C             reserved for a space or another separator.
C             ---------------------
              INTEGER MAXDIG,MINDIG
              PARAMETER (MINDIG=4)
              PARAMETER (MAXDIG=6)
C             ---------------------
C             MAXDIG must be less than 10,
C             MINDIG should be less than MAXDIG.
C
C Date: 2000, January 8
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      INTEGER IFORM1,IFORM2
      REAL SMALL
C     IFORM1,IFORM2... Define format to write the travel times.
C             Limits: 0.LE.IFORM2.LE.9, IFORM2+1.LE.IFORM1.LE.99.
C
C.......................................................................
C
C     Setting output format:
      IFORM1=MAX0(INT(ALOG10(AMAX1(OUTMAX,0.001))+0.3*0.1**MAXDIG+1.),0)
      IF(OUTMIN.LT.0.) THEN
        IFORM1=MAX0(INT(ALOG10(AMAX1(-OUTMIN,0.001))+3.0*0.1**MAXDIG+2.)
     *                                                        ,1,IFORM1)
      END IF
C     Here, IFORM1 is the number of digits left to the decimal point.
      IFORM2=MAX0(MAXDIG-IFORM1,0)
C     IFORM2 is the number of decimal places.
      IFORM1=IFORM1+IFORM2+1
C     IFORM1 is the width of the output field for one element.
      FORMAT='F02.0,A,'
      IF(OUTMIN.NE.0..OR.OUTMAX.NE.0.) THEN
        SMALL=10.**(MINDIG-IFORM2)-0.5*10.**(-IFORM2)
        IF(-SMALL.LE.OUTMIN.AND.OUTMAX.LE.SMALL) THEN
          FORMAT='G00.0,A,'
          IFORM1=MAXDIG+5
          IFORM2=MAXDIG
          IF(OUTMIN.LT.0.) THEN
            IFORM2=MAXDIG-1
          ELSE
            IFORM2=MAXDIG
          END IF
        ELSE IF(IFORM1.GT.MAXDIG+5) THEN
          FORMAT='G00.0,A,'
          IFORM1=MAXDIG+5
          IF(OUTMIN.LT.0.) THEN
            IFORM2=MAXDIG-1
          ELSE
            IFORM2=MAXDIG
          END IF
        END IF
        FORMAT(2:2)=CHAR(ICHAR('0')+IFORM1/10)
        FORMAT(3:3)=CHAR(ICHAR('0')+MOD(IFORM1,10))
        FORMAT(5:5)=CHAR(ICHAR('0')+IFORM2)
      END IF
C     Output format is set.
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE FORM2(NQ,OUTMIN,OUTMAX,FORMAT)
      INTEGER NQ
      REAL OUTMIN(NQ),OUTMAX(NQ)
      CHARACTER*(*) FORMAT
C
C Subroutine designed to determine the best output format for multiples
C of real numbers.
C
C Input:
C     NQ...   Number of reals in each output line.
C     OUTMIN,OUTMAX... Minimum and maximum real numbers to be written.
C     FORMAT..String of at least 8*NQ characters.
C
C Output:
C     FORMAT..String containing the output format, e.g. like
C             'F07.3,A,F07.3,A,F07.3,A,F07.6,A,F07.4,A)'.  The width of
C             the defined string is 8*NQ characters.  It has the above
C             form, where digits are replaced by reasonable values.
C             Note ')' at the end instead of ','.   The subroutine
C             attempts to output at least MAXDIG digits (including all
C             zeros after the decimal point) of the largest positive
C             number and MAXDIG-1 digits of the most negative number,
C             and to adjust the width of the output field to MAXDIG+2
C             columns including the space after the number, if possible.
C             If the number of digits (without leading zeros) is smaller
C             than MINDIG, format Fnn.d with nn=MAXDIG+1 is changed to
C             Gmm.d with mm=MAXDIG+5.  The (MAXDIG+6)th column is
C             reserved for a space or another separator.
C
C Date: 1999, August 16
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      INTEGER I
C
C.......................................................................
C
      DO 10 I=1,NQ
        CALL FORM1(OUTMIN(I),OUTMAX(I),FORMAT(8*I-7:8*I))
   10 CONTINUE
      FORMAT(8*NQ:8*NQ)=')'
C
      RETURN
      END
C
C=======================================================================
C