      PROGRAM PLOTGRP

C     Plot the group velocity curves

      INTEGER PRTFIL, FFTFIL
      PARAMETER ( PRTFIL = 6, FFTFIL = 5, PI = 3.141592,
     &    MAXM = 1500, MAXMED = 51, MAXNF = 1000, MAXN = 16001 )


      INTEGER N( MAXMED ), M( MAXNF )
      REAL FREQ( MAXNF ), Z( MAXN ), DEPTH( MAXMED ), RHO( MAXMED )
      COMPLEX CK( MAXM, MAXNF ), KTOP2, KBOT2
      CHARACTER TITLE*80, MATER( MAXMED )*8, BCTOP*1, BCBOT*1

C     *** Get modal values ***

      MACT  = 0
      NFREQ = 0

      DO IF = 1, MAXNF
         READ( FFTFIL, *, END = 9000 ) F
         NFREQ = IF

         CALL MODHDR( IF, ' ',
     &      TITLE, FREQ( IF ), NMEDIA, NTOT, NMAT,
     &      N, MATER, DEPTH, RHO,
     &      BCTOP, RHOT, DEPTHT,
     &      BCBOT, RHOB, DEPTHB,
     &      M( IF ), MAXM, Z, CK( 1, IF), KTOP2, KBOT2, MODFIL )

         CLOSE( 29 + IF )

         MACT = MAX( M( IF ), MACT )

      END DO   ! next frequency
 9000 CONTINUE

      LAST = INDEX( TITLE, '    ' ) - 1
      WRITE( PRTFIL, 4000 ) TITLE(1:LAST)
 4000 FORMAT( '@ title "', A, '"' )

      WRITE( PRTFIL, 100 )
  100 FORMAT( '@   legend on' )

C     *** Compute group velocity using a centered difference ***

      WRITE( *, 5000 )
 5000 FORMAT( '@WITH G1', /,
     &        '@G1 ON', /,
     &        '@TYPE xy' )
 
       DO MODE = 1, MACT
          WRITE( PRTFIL, 600 ) MODE - 1, MODE
  600     FORMAT( '@   legend string ', I2, ' "Mode ', I2, '"' ) 

          DO IF = 2, NFREQ - 1

C            --- do we have a mode at this frequency?
             IF ( M( IF - 1 ) .GE. MODE .AND.
     &            M( IF + 1 ) .GE. MODE ) THEN

                DELTAF =     FREQ( IF + 1 ) -     FREQ( IF - 1 )
                DELTAK = CK( MODE, IF + 1 ) - CK( MODE, IF - 1 )

                VEL = 0.0
                IF ( DELTAK .NE. 0.0 ) VEL = 2.0 * PI * DELTAF / DELTAK
                WRITE( *, * ) FREQ( IF ), VEL
             ENDIF

         END DO   ! next frequency

         WRITE( *, * ) '&'

      END DO   ! next mode

C     *** Compute phase velocity ***

      WRITE( *, 6000 )
 6000 FORMAT( '@WITH G0', /,
     &        '@G0 ON', /,
     &        '@TYPE xy' )
 
C      WRITE( *, * ) '   linestyle 2'
      DO MODE = 1, MACT
         DO IF = 1, NFREQ

C           --- do we have a mode at this frequency?
            IF ( M( IF ) .GE. MODE ) THEN
               VEL = 2.0 * PI * FREQ( IF ) / CK( MODE, IF )
               WRITE( *, * ) FREQ( IF ), VEL
            ENDIF

         END DO   ! next frequency

         WRITE( *, * ) '&'

      END DO   ! next mode

      STOP
      END
