      PROGRAM PLOTMODE

C     Plots a specified set of modes

      INTEGER PLPFIL
      PARAMETER ( PLPFIL = 5, MAXN = 16000, MAXMED = 51,
     &   MAXM = 3000 )

      INTEGER   N( MAXMED ), IRD( MAXN )
      REAL      Z( MAXN ), PHIR( MAXN ), PHII( MAXN ), RD( MAXN ),
     &          W( MAXN ), DEPTH( MAXMED ), RHO( MAXMED )
      COMPLEX   PHI( MAXN ), CK( MAXM ), KTOP2, KBOT2
      CHARACTER PLTITL*80, COCY*8, MATER( MAXMED )*8, BCTOP*1, BCBOT*1

      OPEN( FILE = 'PLPFIL', UNIT = PLPFIL, STATUS = 'OLD' )
      READ( PLPFIL, * ) COCY
      READ( PLPFIL, * ) ZMIN, ZMAX

C     *** Read the header data from the mode file ***

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

C     *** Set-up vector of depth points ***

      NRD     = 1001
      RD( 1 ) = ZMIN
      RD( 2 ) = ZMAX
      RD( 3 ) = -999.9
      CALL SUBTAB( RD, NRD )
      CALL WEIGHT( Z, NTOT, RD, NRD, W, IRD )

C     *** Read in the mode ***

      DO I = 1, 50
         READ( PLPFIL, *, END = 4000 ) MODE

C        --- Is the mode in the file?
         IF ( MODE .LE. 0 .OR. MODE .GT. M ) THEN
            WRITE( *, * ) 'ERROR for mode = ', MODE
            WRITE( *, * ) 'Highest order mode in file is = ', M
            STOP
         ELSE
C           --- Cycle through COCY string
            L = MOD( I, LEN(COCY) ) + 1

            CALL GETONE( MODE, MODFIL, NTOT, NMAT, W, IRD,
     &         N, MATER, NMEDIA, COCY(L:L),
     &         KTOP2, DEPTHT, BCTOP,
     &         KBOT2, DEPTHB, BCBOT,
     &         RD, NRD, CK, PHI )

C           --- Scale the mode

            T = 0.0
            DO J = 1, NRD
               T = AMAX1( T, ABS( PHI( J ) ) )
            END DO

            WRITE( *, 1000 ) I-1, I-1
 1000       FORMAT( '@WITH G', I1.1, /,
     &              '@G', I1.1, ' ON', /,
     &              '@TYPE xy' )

            WRITE( *, 2000 ) MODE
 2000       FORMAT( '@ subtitle "Mode ', I3, '"' )

            WRITE( *, 3000 ) INT( ZMIN  ), INT( ZMAX  )
 3000       FORMAT( '@ world xmin ', I5, /,
     &              '@ world xmax ', I5 )

            T = 1.0
            DO J = 1, NRD
               PHIR( J ) = REAL(  PHI( J ) ) / T
               PHII( J ) = AIMAG( PHI( J ) ) / T
               WRITE( *, * ) RD( J ), PHIR( J )
            END DO
            WRITE( *, * ) '&'

            DO J = 1, NRD
               WRITE( *, * ) RD( J ), PHII( J )
            END DO
            WRITE( *, * ) '&'
         END IF

      END DO   ! next plot

 4000 STOP
C            IF ( COCY(L:L) .EQ. 'H' ) CALL XNAME( 'Horiz. Displ.$', 100)
C            IF ( COCY(L:L) .EQ. 'V' ) CALL XNAME( 'Vert. Displ.$',  100)
C            IF ( COCY(L:L) .EQ. 'T' ) CALL XNAME( 'Shear Stress$',  100)
C            IF ( COCY(L:L) .EQ. 'N' ) CALL XNAME( 'Normal Stress$', 100)

      END
