      SUBROUTINE GBSOL(A,M,NA,B,N,NRHS,ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,
     .                 EPS,INFO,JERR)

c*********************************************************************72
c
cc GBSOL solves a large banded system by Gauss elimination using disk storage.
c
C     NAME: GBSOL            FORTRAN 77                 JANUARY 15, 1988
C
C     PURPOSE:
C     THIS ALGORITHM IS A 1977 AMERICAN NATIONAL STANDARD FORTRAN
C     IMPLEMENTATION OF A GAUSSIAN ELIMINATION PROCEDURE WITH PARTIAL
C     PIVOTING TO SOLVE SYSTEMS WITH LARGE BANDED MATRICES AND SEVERAL
C     RIGHT-HAND SIDES.
C
C     THE ALGORITHM KEEPS ONLY THE DATA NECESSARY FOR THE ACTUAL COMPUTATION
C     IN THE MEMORY AND SWAPS THE OTHER DATA ON A DISK. IN ORDER
C     TO MINIMIZE I/O OPERATIONS, THE ELEMENTS OF THE MATRIX ARE CALCULATED
C     WHEN THEY ARE NEEDED FOR THE COMPUTATION OF THE UPPER TRIANGULAR
C     MATRIX. HENCE, THE USER HAS TO PROVIDE A SUBROUTINE WHICH
C     CALCULATES THE ROWS OF THE MATRIX.
C
C     THE WHOLE PROGRAM CONFORMS TO STANDARD, BUT IT IS RECOMMENDED THAT
C     THE SUBROUTINES FOR HANDLING RANDOM-ACCESS STORAGE (GOPEN, GCLOSE,
C     AND GRWRAN) BE REPLACED BY INSTALLATION-DEPENDENT VERSIONS.
C
C     THE ELEMENTS OF THE BANDMATRIX ARE STORED IN ROWS AS ILLUSTRATED
C     IN THE FOLLOWING EXAMPLE:
C
C
C
C                *  *  31 41 51  0  0  0  0  0  0
C                   *  22 32 42 52  0  0  0  0  0
C                      13 23 33 43 53  0  0  0  0
C                       0 14 24 34 44 54  0  0  0
C                       0  0 15 25 35 45 55  0  0
C                       0  0  0 16 26 36 46 56  0
C                       0  0  0  0 17 27 37 47 57
C                       0  0  0  0  0 18 28 38 48  *
C                       0  0  0  0  0  0 19 29 39  *  *
C
C
C     THE FIRST INDEX SHOWS THE POSITION OF THE ELEMENT IN THE BAND,
C     AND THE SECOND ONE INDICATES THE ROW. ONLY THE BAND IS STORED,
C     AND THE USER HAS TO SUPPLY ZEROS FOR THE POSITIONS INDICATED WITH
C     AN ASTERISK IN THE ABOVE EXAMPLE. THE MATRIX A DECLARED IN GBSOL
C     CONTAINS THE FOLLOWING ELEMENTS:
C
C                             0  0  31 41 51
C                             0  22 32 42 52
C                            13  23 33 43 53
C                            14  24 34 44 54
C                            15  25 35 45 55
C                            16  26 36 46 56
C                            17  27 37 47 57
C                            18  28 38 48  0
C                            19  29 39  0  0
C
C     THE USER HAS TO PROVIDE A SUBROUTINE "ROW" THAT CALCULATES THE
C     ROWS OF A. THIS SUBROUTINE IS CALLED BY GBSOL IN THE FORM
C
C                       CALL ROW( A(1,L),L )
C
C     WHERE A(1,L) POINTS TO THE ADDRESS OF THE FIRST ELEMENT IN THE
C     L-TH ROW.
C
C     THE ELEMENTS OF THE ROWS OF A ARE STORED IN CONSECUTIVE ADDRESSES
C     SO THAT THE PROGRAM CAN EASILY BE MODIFIED FOR MACHINES WITH PIPE-
C     LINE ARCHITECTURE, SUCH AS THE CYBER 205 OR THE FPS164. THE MODI-
C     FICATION CONSISTS OF REPLACING THE LOOPS MARKED IN THE CODE BY
C     VECTOR OPERATIONS.
C
C     CALLING SEQUENCE:
C
C     CALL GBSOL(A,M,NA,B,N,NRHS,ALOGDT,SIGNDT,PIVMAX,PIVMIN,IPIV,
C                                                   EPS,INFO,JERR)
C     PARAMETERS:
C     A:          DOUBLE PRECISION. A(M,NA), IS THE ARRAY DECLARED IN
C                 THE SUBROUTINE TO PROVIDE MEMORY FOR THE COMPUTATION.
C     M:          INTEGER, INPUT, NUMBER OF ELEMENTS IN EACH ROW OF THE
C                 BAND.
C     NA:         INTEGER, INPUT, SECOND DIMENSION OF A.
C                 BY SPECIFYING NA, THE USER CHOOSES THE SIZE OF THE
C                 PART OF THE MATRIX TO BE KEPT IN THE MEMORY:
C                 NA>=N: THE WHOLE MATRIX IS KEPT IN THE MEMORY.
C                        THERE ARE NO I/O OPERATIONS.
C                 NA <N: ONLY PART OF THE MATRIX ( A "BLOCK") IS KEPT
C                        IN THE MEMORY. THE ONLY REQUIREMENT IS
C                                  NA >= (M+1)/2 + 2.
C     B:          DOUBLE PRECISION. B(N,NRHS), CONTAINS THE RIGHT-HAND
C                 SIDES AS INPUT AND THE SOLUTION VECTORS AS OUTPUT.
C     N:          INTEGER, INPUT, NUMBER OF EQUATIONS.
C     NRHS:       INTEGER, INPUT, NUMBER OF RIGHT-HAND SIDES.
C     ALOGDT:     DOUBLE PRECISION, OUTPUT, LOGARITHM OF THE ABSOLUTE
C                 VALUE OF THE DETERMINANT.
C     SIGNDT:     DOUBLE PRECISION, OUTPUT, SIGN OF THE DETERMINANT.
C                 THE DETERMINANT ITSELF CAN BE CALCULATED BY
C                     DET( A )   =   SIGNDT * 10.0D0**ALOGDT.
C     PIVMAX:     DOUBLE PRECISION, OUTPUT, MAXIMUM OF THE ABSOLUTE
C                 VALUES OF THE PIVOT ELEMENTS.
C     PIVMIN:     DOUBLE PRECISION, OUTPUT, MINIMUM OF THE ABSOLUTE
C                 VALUES OF THE PIVOT ELEMENTS.
C     IPIV:       INTEGER, OUTPUT, POSITION OF THE PIVOT WITH MINIMAL
C                 ABSOLUTE VALUE.
C     EPS:        DOUBLE PRECISION, INPUT. ALL PIVOT ELEMENTS WHOSE
C                 ABSOLUTE VALUES ARE LESS THAN EPS ARE REPLACED BY
C                                SIGN( PIVOT ) * EPS.
C     INFO:       INTEGER,OUTPUT,
C                 = 0, NO PIVOT ELEMENT IS REPLACED.
C                 = L, POSITION OF THE LAST ELEMENT WHICH IS REPLACED.
C
C                 FOR NORMAL USE WE RECOMMEND SETTING
C                    EPS = DSQRT( SMALLEST DOUBLE PRECISION CONSTANT )
C                 IN ORDER TO AVOID OVERFLOW IN CASE OF A SINGULAR
C                 MATRIX. WITH THIS SETTING AND A PROPERLY SCALED MATRIX
C                 (I.E. PIVMAX = O(1)), A RETURNED VALUE OF INFO DIFFER-
C                 ENT FROM ZERO INDICATES THAT THE SOLUTIONS HAVE NO
C                 SIGNIFICANT FIGURES.
C
C     JERR:       INTEGER, OUTPUT,ERROR CODE:
C                 JERR=0: NA >= (M+1)/2+2.
C                 JERR=1: NA <  (M+1)/2+2. CHOOSE A LARGER VALUE FOR NA.
C
C     USER SUPPLIED SUBROUTINES CALLED:
C
C     ROW         CALCULATES A ROW OF THE BAND OF MATRIX AS DESCRIBED
C                 ABOVE.
C
C     GBSOL CONTAINS THE FOLLOWING CALLS IN ORDER TO PERFORM THE
C     INSTALLATION-DEPENDENT DIRECT-ACCESS IN- AND OUTPUT OPERATIONS:
C
C
C     CALL GOPEN(M)    TO OPEN A TEMPORARY FILE FOR DIRECT-ACCESS I/O.
C
C     CALL GCLOSE      TO CLOSE AND FREE THE FILE OPENED WITH GOPEN.
C
C     CALL GRWRAN(IOP,A,NREAL8,NSTART,M)
C
C                      IOP=1: TO READ NREAL8 DOUBLE PRECISION NUMBERS
C                             FROM A DISK TO THE ARRAY A, STARTING WITH
C                             THE NSTART-TH NUMBER ON THE DISK, OR
C
C                      IOP=2: TO WRITE THE FIRST NREAL8 DOUBLE PRECISION
C                             NUMBERS OF THE ARRAY A ON A DISK, STARTING
C                             WITH THE POSITION NSTART.
C
C     THE USER SHOULD APPLY UTILITY PROGRAMS WRITTEN IN ASSEMBLER IN
C     ORDER TO GUARANTEE RAPID I/O OPERATIONS. STANDARD FORTRAN VERSIONS
C     OF THE ABOVE SUBROUTINES ARE SUPPLIED, BUT FORTRAN SUPPLIED I/O IS
C     TOO SLOW.
C
C     AUTHOR:     GEZA SCHRAUF, CAL TECH
C
      DOUBLE PRECISION A(M,NA),B(N,NRHS)
      DOUBLE PRECISION ABSPIV,AH,ALOGDT,EPS,PIV,PIVMAX,PIVMIN,RA1L,RA1N,
     .                 RPIV,S,SIGNDT
      INTEGER I,IH,INFO,IPIV,ISTART,J,JA,JB,JEND,JERR,JH,JPREL,JSTART,
     .        K,L,LABS, LABS1,LAMLR,LB,LBLK,LBLKM1,LDONE,LEND,LENGTH,
     .        LREL,LRELP1,LREL1,LSTABS,M,MD,MD1,M1,N,NA,NBLKM1,
     .        NBLKS,NEWABS,NEWREL,NREAL8,NRHS,NRLB,NRPB,NRPBPJ,NRPS,
     .        NRPSM1,NSTART,NTEMP,N1

      ALOGDT = 0.D0
      SIGNDT = 1.D0
      PIVMAX = 0.D0
      INFO = 0
      IPIV = 1
      JERR = 0

      MD = (M+1)/2
      MD1 = MD - 1
      M1 = M - 1
      N1 = N - 1

      IF (NA.LT.N) GO TO 10

      NRPB = N
      NRLB = N
      NBLKS = 1
      NBLKM1 = NBLKS - 1
      GO TO 30

   10 NRPB = NA - MD
      NBLKS = N/NRPB
      NTEMP = NBLKS*NRPB
      IF(N.NE.NTEMP) NBLKS = NBLKS + 1

      NBLKM1 = NBLKS - 1
      NRLB = N - NBLKM1*NRPB

      IF (NRPB.GE.2) GO TO 20
      JERR = 1
      RETURN
   20 CONTINUE

      NREAL8 = M * NRPB
C
C  OPEN TEMPORARY FILE FOR DIRECT ACCESS.
C
      CALL GOPEN(M)
   30 CONTINUE
C
C  CALCULATE THE FIRST MD ROWS OF THE MATRIX A.
C
      DO 40 L = 1,MD
         CALL ROW(A(1,L),L)
   40 CONTINUE
C
C  SHIFT THE MATRIX ELEMENTS OF THE FIRST MD1 ROWS TO THE LEFT.
C
      IH = MD
      ISTART = IH

      DO 70 J = 1,MD1
         DO 50 I = 1,IH
            A(I,J) = A(ISTART+I-1,J)
   50    CONTINUE
         IH = IH + 1
         DO 60 I = IH,M
            A(I,J) = 0.D0
   60    CONTINUE
         ISTART = ISTART - 1
   70 CONTINUE
C
C  INITIALIZE PIVMIN.
C
      PIVMIN = 0.D0

      DO 75 J=1,MD
         AH = DABS(A(1,J))
         IF (AH.GT.PIVMIN) PIVMIN = AH
   75 CONTINUE
C
C  CALCULATION OF THE UPPER TRIANGULAR MATRIX.
C
      DO 200 LBLK = 1,NBLKS
         LEND = NRPB
         IF (LBLK.EQ.NBLKS) LEND = NRLB - 1

         LDONE = (LBLK-1)*NRPB

         DO 170 L = 1,LEND
            LREL = L
            LABS = LDONE + LREL
            LREL1 = LREL - 1
            LABS1 = LABS - 1
            LSTABS = MIN0(LABS+MD1,N)
            NRPS = LSTABS - LABS1
C
C  PIVOT SEARCH.
C
            ABSPIV = 0.D0
            JPREL = 1
            DO 80 J = 1,NRPS
               AH = DABS(A(1,J+LREL1))
               IF (AH.LE.ABSPIV) GO TO 80
               JPREL = J
               ABSPIV = AH
   80       CONTINUE
C
C  INTERCHANGE THE ROWS OF A AND THE CORRESPONDING RHS'S.
C
            IF (JPREL.EQ.1) GO TO 110

            JA = JPREL + LREL1
            JB = JPREL + LABS1

            DO 90 I = 1,M
               S = A(I,JA)
               A(I,JA) = A(I,LREL)
               A(I,LREL) = S
   90       CONTINUE
            DO 100 I = 1,NRHS
               S = B(JB,I)
               B(JB,I) = B(LABS,I)
               B(LABS,I) = S
  100       CONTINUE

  110       PIV = A(1,LREL)
            IF (DABS(PIV).GE.EPS) GO TO 120
            IF (PIV.NE.0.D0) PIV = DSIGN(EPS,PIV)
            IF (PIV.EQ.0.D0) PIV = EPS
            A(1,LREL) = PIV
            INFO = LABS
  120       CONTINUE

            IF (JPREL.NE.1) SIGNDT = -SIGNDT
            IF (PIV.LT.0.D0) SIGNDT = -SIGNDT
            ABSPIV = DABS(PIV)
            ALOGDT = ALOGDT + DLOG10(ABSPIV)
            IF (ABSPIV.GT.PIVMAX) PIVMAX = ABSPIV
            IF (ABSPIV.GT.PIVMIN) GO TO 130
            PIVMIN = ABSPIV
            IPIV = LABS
  130       CONTINUE
C
C  ADD A SUITABLE MULTIPLE OF THE LREL-TH ROW OF A TO THE FOLLOWING
C  ROWS OF A, AND SHIFT THE MATRIX ROWS TO THE LEFT, SO THAT THE
C  NEXT PIVOT IS TO BE SEARCHED AGAIN IN THE FIRST COLUMN OF A.
C  ADD THE SAME MULTIPLE OF THE LABS-TH ROW OF B TO THE FOLLOWING
C  ROWS OF B.
C
            NRPSM1 = NRPS - 1
            RPIV = 1.D0/PIV

            JSTART = LABS + 1
            JEND = LABS + NRPSM1
            LAMLR = LABS - LREL

            DO 150 I = 1,NRHS
               S = -B(LABS,I)*RPIV
c
C  VECTORIZE THE FOLLOWING LOOP, IF POSSIBL.
c
               DO 161 J = JSTART,JEND
                  B(J,I) = A(1,J-LAMLR)*S + B(J,I)
  161          CONTINUE
  150       CONTINUE

            DO 160 JH = 1,NRPSM1
               JA = LREL + JH
               JB = LABS + JH
               S = -A(1,JA)*RPIV
c
C  VECTORIZE THE FOLLOWING LOOP, IF POSSIBLE.
c
               DO 140 I = 1,M1
                  A(I,JA) = A(I+1,LREL)*S + A(I+1,JA)
  140          CONTINUE
               A(M,JA) = 0.D0
  160       CONTINUE

            NEWABS = LABS + MD
            NEWREL = NEWABS - LDONE

            IF (NEWABS.LE.N) CALL ROW(A(1,NEWREL),NEWABS)
  170    CONTINUE

         IF (LBLK.EQ.NBLKS) GO TO 200
C
C  WRITE THE NREAL8 = M*NRPB NUMBERS JUST CALCULATED ON THE TEMPORARY FILE.
C
         NSTART = LDONE*M + 1
         CALL GRWRAN(2,A,NREAL8,NSTART,M)
C
C  SHIFT THE REMAINING MD ROWS UP.
C
         DO 190 J = 1,MD
            NRPBPJ = NRPB + J
            DO 180 I = 1,M
               A(I,J) = A(I,NRPBPJ)
  180       CONTINUE
  190    CONTINUE

  200 CONTINUE
c
C  THE LAST ELEMENT.
c
      LRELP1 = LREL + 1
      PIV = A(1,LRELP1)
      IF (DABS(PIV).GE.EPS) GO TO 210
      IF (PIV.NE.0.D0) PIV = DSIGN(EPS,PIV)
      IF (PIV.EQ.0.D0) PIV = EPS
      A(1,LRELP1) = PIV
      INFO = N
  210 CONTINUE

      IF (PIV.LT.0.D0) SIGNDT = -SIGNDT
      ABSPIV = DABS(PIV)
      ALOGDT = ALOGDT + DLOG10(ABSPIV)
      IF (ABSPIV.GT.PIVMAX) PIVMAX = ABSPIV
      IF (ABSPIV.GT.PIVMIN) GO TO 220
      PIVMIN = ABSPIV
      IPIV = N
  220 CONTINUE
C
C  SOLVE THE SYSTEMS WITH THE UPPER TRIANGULAR MATRIX.
C
      RA1N = 1.D0/A(1,LRELP1)
      DO 230 I = 1,NRHS
         B(N,I) = B(N,I)*RA1N
  230 CONTINUE

      DO 290 LB = 1,NBLKS
         LBLK = NBLKS + 1 - LB
         LBLKM1 = LBLK - 1
         LEND = NRPB

         IF (LBLK.NE.NBLKS) GO TO 240
         LEND = NRLB - 1
         GO TO 250

  240    CONTINUE
c
C  READ NREAL8 = M*NRPB NUMBERS FROM THE TEMPORARY FILE.
c
         NSTART = LBLKM1*NREAL8 + 1
         CALL GRWRAN(1,A,NREAL8,NSTART,M)

  250    CONTINUE
c
C  BACKSOLVE.
c
         DO 280 L = 1,LEND
            LREL = LEND + 1 - L
            LABS = LBLKM1*NRPB + LREL
            LENGTH = M1
            IF (LBLK.EQ.NBLKS) LENGTH = MIN0(L,M1)

            RA1L = 1.D0/A(1,LREL)
            DO 270 K = 1,NRHS
               S = 0.D0
c
C  VECTORIZE THE FOLLOWING LOOP, IF POSSIBLE.
c
               DO 260 I = 1,LENGTH
                  S = A(I+1,LREL)*B(I+LABS,K) + S
  260          CONTINUE
               B(LABS,K) = (B(LABS,K)-S)*RA1L
  270       CONTINUE
  280    CONTINUE

  290 CONTINUE
c
C  CLOSE AND FREE TEMPORARY FILE.
c
      IF ( NA .EQ. N ) GO TO 300
      CALL GCLOSE
  300 RETURN
      END
      SUBROUTINE GOPEN(M)

c*********************************************************************72
c
cc GOPEN opens a temporary file for direct access.
C
C     CALLING SEQUENCE: CALL GOPEN(M)
C
C     PARAMETERS:
C     M:          INTEGER, INPUT, NUMBER OF DOUBLE-PRECISION NUMBERS IN
C                 EACH DIRECT-ACCESS RECORD. USED ONLY FOR THE STANDARD
C                 FORTRAN 77 VERSION.
C
C     AUTHOR:     GEZA SCHRAUF, CAL TECH
C

      INTEGER IUNIT,IRECL,M
      IUNIT = 9
C
C     STANDARD FORTRAN 77 VERSION. THIS VERSION SHOULD BE REPLACED BY
C     A FASTER MACHINE-DEPENDENT ROUTINE.
C
C  CHOOSE THE RECORD LENGTH "IRECL" SO THAT ONE
C  RECORD CONTAINS  M  DOUBLE PRECISION NUMBERS.
C
      IRECL = M * 8
C
      OPEN(IUNIT,STATUS='SCRATCH',ACCESS='DIRECT',RECL=IRECL)
C
C     INSTALLATION FOR THE IBM 4331 OF THE DEPARTMENT OF APPLIED MATHE-
C     MATICS AT THE CALIFORNIA INSTITUTE OF TECHNOLOGY.
C
C     IBM FORTRAN UTILITIES FOR VM/370 *** (PROGRAM NUMBER: 5798-DFH)
C
C     UTILITY USED: UOPEN IN FORTRAN 77 LEVEL LANGUAGE (MANUAL P. 34).
C
C-    CALL UOPEN(IUNIT,'TEMP DATA A',5,3,0,0,' ',0,IERR)
C
C-    IF (IERR.EQ.0) GO TO 10
C
C-    WRITE (6,9001) IERR
C9001 FORMAT (///' AN ERROR OCCURRED USING THE IBM FORTRAN UTILITY',
C-   .  ' "UOPEN" WHICH WAS '/' CALLED BY "GOPEN":'//' IBM ERROR CODE:',
C-   .  I4//' 0: OPENED AS REQUESTED.'/
C-   .  ' 1: FILE ALREADY OPENED WITH ANOTHER UNIT NUMBER.'/
C-   .  ' 2: SHOULD NOT OCCUR.'/
C-   .  ' 3: UNIT NUMBER ALREADY ASSOCIATED WITH A FILE.'/
C-   .  ' 8: DATA INPUT ERROR. FILE NOT OPENED.'///)
C
C- 10 CONTINUE
C
      RETURN
      END
      SUBROUTINE GCLOSE

c*********************************************************************72
c
cc GCLOSE closes and frees a temporary file opened by GOPEN.
C
C     AUTHOR:     GEZA SCHRAUF, CAL TECH
C
      INTEGER IUNIT
      IUNIT = 9
C
C     STANDARD FORTRAN 77 VERSION. THIS VERSION SHOULD BE REPLACED BY
C     A MACHINE-DEPENDENT ROUTINE.
C
      CLOSE(IUNIT)
C
C     INSTALLATION FOR THE IBM 4331 OF THE DEPARTMENT OF APPLIED MATHE-
C     MATICS AT THE CALIFORNIA INSTITUTE OF TECHNOLOGY.
C
C     IBM FORTRAN UTILITIES FOR VM/370 *** (PROGRAM NUMBER: 5798-DFH)
C
C     UTILITY USED: UCLOSE IN FORTRAN 77 LEVEL LANGUAGE (MANUAL P. 32).
C
C-    CALL UCLOSE(IUNIT,2,IERR)
C
C-    IF (IERR.EQ.0) GO TO 10
C
C-    WRITE (6,9001) IERR
C9001 FORMAT (///' AN ERROR OCCURRED USING THE IBM FORTRAN UTILITY',
C-   .  ' "UCLOSE" WHICH WAS '/' CALLED BY "GCLOSE":'//
C-   .  ' IBM ERROR CODE:',I4//' 0: CLOSED AS REQUESTED.'/
C-   .  ' 4: UNIT IS NOT OPENED.'/' 8: DATA INPUT ERROR.'///)
C
C- 10 CONTINUE
C
      RETURN
      END
      SUBROUTINE GRWRAN(IOP,ARRAY,NREAL8,NSTART,M)

c*********************************************************************72
c
cc GRWRAN reads from or writes to a direct access file.
c
C     CALLING SEQUENCE: CALL GRWRAN(IOP,ARRAY,NREAL8,NSTART,M)
C
C     PARAMETERS:
C     IOP:        INTEGER, OPERATION MODE:
C                 IOP = 1: READ FROM THE FILE TO THE ARRAY "ARRAY."
C                 IOP = 2: WRITE FROM THE ARRAY "ARRAY" TO THE FILE.
C     ARRAY:      DOUBLE PRECISION, ARRAY TO OR FROM WHICH DATA IS
C                 TRANSFERRED.
C     NREAL8:     INTEGER, NUMBER OF DOUBLE-PRECISION NUMBERS TO BE
C                 TRANSFERRED.
C     NSTART:     INTEGER, THE FIRST DOUBLE-PRECISION NUMBER ON THE DISK
C                 WHICH IS TO BE TRANSFERRED.
C     M:          INTEGER, INPUT, NUMBER OF DOUBLE PRECISION NUMBERS IN
C                 EACH DIRECT-ACCESS RECORD. ONLY USED FOR THE STANDARD
C                 FORTRAN 77 VERSION.
C
C     AUTHOR:     GEZA SCHRAUF, CAL TECH
C

      DOUBLE PRECISION ARRAY(NREAL8)
      INTEGER IOP,IUNIT,K,KFIRST,KLAST,L,LREC,LSTART,
     .        M,NRPB,NREAL8,NSTART
      IUNIT = 9
C
C     STANDARD FORTRAN 77 VERSION. THIS VERSION SHOULD BE REPLACED BY
C     A FASTER MACHINE-DEPENDENT ROUTINE.
C
      LSTART = (NSTART - 1)/M
      NRPB   = NREAL8/M
C
      IF   (IOP.EQ.1)   THEN
         DO 10 L=1,NRPB
            LREC   = LSTART + L
            KFIRST = (L-1)*M + 1
            KLAST  = KFIRST + M - 1
            READ(IUNIT,REC=LREC)    (ARRAY(K),K=KFIRST,KLAST)
   10    CONTINUE
      ELSE IF (IOP.EQ.2)   THEN
         DO 20 L=1,NRPB
            LREC   = LSTART + L
            KFIRST = (L-1)*M + 1
            KLAST  = KFIRST + M - 1
            WRITE(IUNIT,REC=LREC)   (ARRAY(K),K=KFIRST,KLAST)
   20    CONTINUE
      ENDIF
C
C     ACCELERATED STANDARD FORTRAN 77 VERSION WITH LARGER RECORD LENGTH.
C     NREAL8 HAS TO BE  S M A L L E R  THAN THE MAXIMAL POSSIBLE RECORD
C     LENGTH. TO USE THIS VERSION REPLACE THE STATEMENT "CALL GOPEN(M)"
C     THE BY STATEMENT "CALL GOPEN(NREAL8)" IN SUBROUTINE GBSOL.
C
C-    LREC = (NSTART - 1)/NREAL8 + 1
C
C-    IF   (IOP.EQ.1)   THEN
C-       READ(IUNIT,REC=LREC)  ARRAY
C-    ELSE IF (IOP.EQ.2)   THEN
C-       WRITE(IUNIT,REC=LREC) ARRAY
C-    ENDIF
C
C     INSTALLATION FOR THE IBM 4331 OF THE DEPARTMENT OF APPLIED MATHE-
C     MATICS AT THE CALIFORNIA INSTITUTE OF TECHNOLOGY.
C
C     IBM FORTRAN UTILITIES FOR VM/370 *** (PROGRAM NUMBER: 5798-DFH)
C
C     UTILITY USED: RWFIL IN FORTRAN 77 LEVEL LANGUAGE (MANUAL P. 25).
C
C-    IFIRST = 2*NSTART - 1
C-    IWORDS = 2*NREAL8
C
C-    CALL RWFIL(IUNIT,IOP,IFIRST,ARRAY(1),IWORDS,IERR)
C
C-    IF (IERR.EQ.0) GO TO 10
C
C-    WRITE (6,9001) IERR,IOP,NREAL8,NSTART
C9001 FORMAT (///' AN ERROR OCCURRED USING THE IBM FORTRAN UTILITY',
C-   .  ' "RWFIL" WHICH WAS '/' CALLED BY "GRWRAN":'//
C-   .  ' IBM ERROR CODE:',I4//' 0: NO ERROR.'/
C-   .  ' 5: END OF FILE ENCOUNTERED.'/
C-   .  ' 6: UNIT HAS NOT BEEN OPENED AS A DIRECT-ACCESS FILE.'/
C-   .  ' 8: DATA INPUT ERROR.'//' PARAMETER VALUES OF "GRWRAN":'//
C-   .  '    OPERATION CODE             IOP    =',
C-   .  I4/'    NUMBER OF REAL*8 NUMBERS   NREAL8 =',
C-   .  I10/'    FIRST REAL*8 NUMBER        NSTART =',I10///)
C
C- 10 CONTINUE
C
      RETURN
      END
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Discussion:
c
c    This FORTRAN77 version is made available for cases where the
c    FORTRAN90 version cannot be used.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    None
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, month(m), y, h, ':', n, ':', s, '.', mm, ampm

      return
      end
