      PROGRAM PLOTFIELD

C     Produces a shade plot
C     Use the option 'RF' for a FOM plot

#include "fgl.h"

      INTEGER SHDFIL, PRTFIL
      PARAMETER ( PLPFIL = 1, PRTFIL = 6, SHDFIL = 35,
     &   MAXNSD = 501, MAXNRD = 4096, MAXNR = 4094, MAXNRZ = 500000 )

      LOGICAL SEG, SMOOTHFL, GRIDFL, COLFL, FINAL, PEAKFL, SNAPFL
      INTEGER WINIDENT
      REAL    R( MAXNR ), SD( MAXNSD), RD( MAXNRD ), TL( MAXNRZ ),
     &        CVEC( 2 * MAXNRZ)
      CHARACTER TITLE*80, TITLEX*80, TITLEY*80, TITLEZ*80,
     &        OPTION*80, COORD*10, FILNAM*80,
     &        DEVCOL*9, PLTTYP*10

      DATA GRIDFL, COLFL, SMOOTHFL, FINAL, PEAKFL, BARFL, SNAPFL
     &     / 7 * .FALSE. /

C     *** Read in user input ***

      OPEN( FILE = 'PLPFIL', UNIT = PLPFIL, STATUS = 'OLD' )

      READ( PLPFIL, * ) OPTION
      READ( PLPFIL, * ) COORD
      READ( PLPFIL, * ) DEVCOL
      READ( PLPFIL, * ) ISEG
      READ( PLPFIL, * ) DR3DB
      READ( PLPFIL, * ) NPLOTS

      IF (  COORD(3:3) .EQ. 'G' )   GRIDFL = .TRUE. ! Plot grid
      IF ( DEVCOL(4:4) .EQ. 'C' )    COLFL = .TRUE. ! Color flag
      IF ( DEVCOL(5:5) .EQ. 'S' ) SMOOTHFL = .TRUE. ! Smooth color cell
      IF ( DEVCOL(6:6) .EQ. 'F' )    FINAL = .TRUE. ! Final labels
      IF ( DEVCOL(7:7) .EQ. 'P' )   PEAKFL = .TRUE. ! Box around peak
      IF ( DEVCOL(8:8) .NE. 'N' )    BARFL = .TRUE. ! Draw color bar
      IF ( DEVCOL(9:9) .EQ. 'S' )   SNAPFL = .TRUE. ! Take snapshot

C     *** Window size ***

      IF      ( NPLOTS .EQ. 1 ) THEN
         CALL PREFPO( 0,  850, 5,  900 )
C        CALL PREFPO( 0,  650, 5,  350 )
      ELSE IF ( NPLOTS .EQ. 2 ) THEN
         CALL PREFPO( 0,  850, 5,  900 )
      ELSE IF ( NPLOTS .EQ. 3 ) THEN
         CALL PREFPO( 0,  750, 5,  988 )
      ELSE IF ( NPLOTS .EQ. 4 ) THEN
         CALL PREFPO( 0, 1260, 5,  988 )
      ELSE IF ( NPLOTS .EQ. 5 ) THEN
         CALL PREFPO( 0, 1260, 5,  988 )
      ELSE
         CALL PREFPO( 0, 1260, 5,  988 )
      END IF

      CALL FOREGR()
      WINIDENT = WINOPE( 'PLOTFIELD', 9 )
C     CALL RGBMOD()      ! Use RGB mode instead of default (map)

      CALL ONEMAP()
      CALL CMODE()
      CALL GCONFI()
      CALL WINSET( WINIDENT )
      CALL GFLUSH()

C     --- set foreground/background
C      READ( 78, * ) ICOL, IRED, IGREEN, IBLUE
C      CALL MAPCOL(  ICOL, IRED, IGREEN, IBLUE )
C      READ( 78, * ) ICOL, IRED, IGREEN, IBLUE
C      CALL MAPCOL(  ICOL, IRED, IGREEN, IBLUE )
C      CALL SYSTEM( "makemap" )  ! ensure color table loaded

      CALL COLOR( BLACK )
      CALL COLOR( WHITE )

      CALL CLEAR()

C     *** Axis labels ***

      IF ( OPTION(2:2) .EQ. 'F' ) THEN
         TITLEX = 'Figure of Merit (dB)$'
         TITLEY = 'Depth (m)$'
         TITLEZ = 'Radius (km)$'
      ELSE
         TITLEX = 'r (km)$'
         TITLEY = 'z (m)$'
         IF ( OPTION(2:2) .EQ. 'L' ) THEN
            TITLEZ = ' $'
         ELSE
            TITLEZ = 'dB$'
         END IF
      END IF

C     *** Loop over plot number ***

      DO IPLOT = 1, NPLOTS

C     ------ Get source depth and plot window info
      READ( PLPFIL, *, END = 9500 ) FILNAM, SRCD
C     --- Append the extension '.shd' to the filename
      JJ = INDEX( FILNAM, ' ' )
      FILNAM( JJ:JJ+3 ) = '.shd'

      XAXL = 0.0
      YAXL = 0.0
      READ( PLPFIL, * ) XMIN, XMAX, XINC, XAXL, TITLEX
      READ( PLPFIL, * ) YMIN, YMAX, YINC, YAXL, TITLEY

C     *** Read header ***

      CALL RDHEAD( SHDFIL, FILNAM, MAXNR, MAXNSD, MAXNRD,
     &   FREQ, R, NR, DELTAR, PLTTYP, XS, YS, THETA,
     &   SD, NSD, RD, NRD, TITLE )

      IF ( NRD .LE. 1 ) CALL ERROUT( PRTFIL, 'F', 'PLOTFIELD',
     &                  'The number of receiver depths must be > 1' )

C     ------ FOM plots don't need conversion to km
      IF ( OPTION(2:2) .NE. 'F' ) THEN
         DO J = 1, NR
            R( J ) = R( J ) / 1000.0
         END DO
      ENDIF

C     *** Calculate source depth index ***

      ISD = NEIGH( SRCD, SD, NSD )

C     --- Check to see if within tolerance

      IF ( FREQ .EQ. 0.0 ) THEN
         TOL = 0.001
      ELSE
         TOL = 1500.0 / ( 6.0 * FREQ )  ! one wavelength
      ENDIF

      IF ( ABS( SD( ISD ) - SRCD ) .GT. TOL ) THEN
         WRITE( PRTFIL, * ) 'Sources exist at the following depths:'
         WRITE( PRTFIL, * ) ( SD( I ), I = 1, NSD )
         CALL ERROUT( PRTFIL, 'F', 'PLOTFIELD',
     &                'No source at specified depth' )
      ENDIF

C     Check program has enough space for the field

      NRZ = NRD * NR
      IF ( NRZ .GT. MAXNRZ ) THEN
         WRITE( PRTFIL, * ) 'NRZ = ', NRZ, '  MAXNRZ = ', MAXNRZ
         CALL ERROUT( PRTFIL, 'F', 'PLOTFIELD',
     &                'Too many field points' )
      ENDIF

C     *** Read and smooth the surface ********************************

      CALL READSURF( XMIN, XMAX, YMIN, YMAX,
     &   RD, NRD, R, NR, DELTAR, ISD, OPTION, DR3DB,
     &   TL, RMINT, ZMINT, TLMINT, RMAXT, ZMAXT, TLMAXT )

C     For snapshot, scale the field to have unit maximum

      IF ( SNAPFL ) THEN
        TLT = MAX( ABS( TLMINT ), ABS( TLMAXT ) )

        SCALE = 1.0
        IF ( TLT .GT. 0.0 ) SCALE = 1.0 / TLT

        DO JJ = 1, NRD * NR
           TL( JJ ) = SCALE * TL( JJ )
        END DO

        TLMINT = -0.9
        TLMAXT =  0.9
      ENDIF

      TLMIN = 0.0
      TLMAX = 0.0
      READ( PLPFIL, * ) TLMIN, TLMAX, DELTL, TITLEZ

C     Autoscale?

      IF ( TLMIN .EQ. 0.0 .AND. TLMAX .EQ. 0.0 ) THEN
         TLMIN = TLMAXT
         TLMAX = TLMINT
        tlmax = tlmin - 5.0
      ENDIF

C     *** Plot one of the subplots ***

      CALL PLOTONE( R, RD, TL, TITLE,
     &   COORD, IPLOT, NPLOTS,
     &   CVEC, NR, NRD,
     &   XMIN, XMAX, XINC, XAXL, XOR, TITLEX,
     &   YMIN, YMAX, YINC, YAXL, YOR, TITLEY,
     &   TLMIN, TLMAX, TLINC, TITLEZ,
     &   SMOOTHFL, GRIDFL, COLFL, FINAL, PEAKFL, BARFL, SNAPFL,
     &   SRCD, TLT, FREQ, PLTTYP, XS, YS, THETA,
     &   RMAXT, ZMAXT )

      END DO    ! Next plot

 9500 CONTINUE

      CALL GFLUSH()   ! flush graphics buffer

      IF ( DEVCOL(1:1) .EQ. 'S' ) THEN
C        CALL SYSTEM( "scrsave plotfield.rgb 40, 680,  50,  530" )
C        CALL SYSTEM( "scrsave plotfield.rgb 40, 360,  50,  290" )
C        CALL SYSTEM( "scrsave plotfield.rgb 0, 650,  5,  350" )
C        CALL SYSTEM( "scrsave plotfield.rgb 30, 620, 15,  795" )
         CALL SYSTEM( "scrsave plotfield.rgb 10, 1060, 55,  950" )
C        CALL SYSTEM( "scrsave plotfield.rgb 10, 600,  65,  800" )

      ELSE
         WRITE( *, * ) 'Type return to continue'
         READ( *, 9999 ) JUNK
 9999    FORMAT( A1 )
      END IF

      END

C-----------------------------------------------------------------------

      SUBROUTINE READSURF( XMIN, XMAX, YMIN, YMAX,
     &   RD, NRD, R, NR, DELTAR, ISD, OPTION, DR3DB,
     &   TL, RMINT, ZMINT, TLMINT, RMAXT, ZMAXT, TLMAXT )

C     *** Read pressure surfaces and convert to smoothed TL ***

      INTEGER SHDFIL, PRTFIL
      PARAMETER ( PRTFIL = 6, SHDFIL = 35, MAXNR = 4094 )

      REAL    RD( * ), R( * ), TL( * ), TLRAW( MAXNR ), TLSM( MAXNR )
      COMPLEX P( MAXNR )
      CHARACTER OPTION*80

      TLMINT =  1.0E9
      TLMAXT = -1.0E9

      IREC =  6 + NRD * ( ISD - 1 )

      DO I = 1, NRD
         IREC = IREC + 1
         READ( SHDFIL, REC = IREC ) ( P( J ), J = 1, NR )

         DO J = 1, NR

            IF ( OPTION(2:2) .EQ. 'L' .OR. OPTION(2:2) .EQ. 'F' ) THEN
C              ------ Linear scale
               PSQ = REAL( P( J ) )
               IF ( OPTION(1:1) .EQ. 'S' ) THEN
                  TLRAW( J ) = SQRT( 1000.0 * R( J ) ) * PSQ
               ELSE
                  TLRAW( J ) = PSQ
               ENDIF
            ELSE
C              ------ dB scale
               PSQ = P( J ) * CONJG( P( J ) )
               IF ( OPTION(1:1) .EQ. 'S' ) PSQ = 1000.0 * R( J ) * PSQ

               TLRAW( J ) = 200.0
               IF ( PSQ .GT. 1.0E-20 ) TLRAW( J ) = -10.0 *ALOG10( PSQ )
            ENDIF

         END DO   ! Next range

C        *** Smooth ***

         NROWS = 1
         CALL SMOOTH( TLRAW, TLSM, NROWS, 1, NR, DELTAR, DR3DB )

C        ------ Copy smoothed data in matrix
         DO J = 1, NR
            JJ = ( I - 1 ) * NR + J
            TL( JJ ) = TLSM( J )

C           --- check point is in plot window
            IF ( ABS( R(  J ) - XMIN ) + ABS( R(  J ) - XMAX )
     &       .LE. 1.001 * ABS( XMAX - XMIN ) .AND.
     &           ABS( RD( I ) - YMIN ) + ABS( RD( I ) - YMAX )
     &       .LE. 1.001 * ABS( YMAX - YMIN ) ) THEN

            IF ( TL( JJ ) .GT. TLMAXT ) THEN   ! log maximum
               TLMAXT = TL( JJ )
               RMAXT  = R(   J )
               ZMAXT  = RD(  I )
            END IF

            IF ( TL( JJ ) .LT. TLMINT ) THEN   ! log minimum
               TLMINT = TL( JJ )
               RMINT  = R(   J )
               ZMINT  = RD(  I )
            END IF
            END IF

         END DO   ! next range J

      END DO   ! Next depth IRD

      CLOSE( SHDFIL )

      WRITE( *, * ) 'R, Z, MAX = ', RMAXT, ZMAXT, TLMAXT
      WRITE( *, * ) 'R, Z, MIN = ', RMINT, ZMINT, TLMINT
      WRITE( *, * )

      RETURN
      END

C-----------------------------------------------------------------------

      SUBROUTINE PLOTONE( R, RD, TL, TITLE,
     &   COORD, IPLOT, NPLOTS,
     &   CVEC, NR, NRD,
     &   XMIN, XMAX, XINC, XAXL, XOR, TITLEX,
     &   YMIN, YMAX, YINC, YAXL, YOR, TITLEY,
     &   ZMIN, ZMAX, ZINC, TITLEZ,
     &   SMOOTHFL, GRIDFL, COLFL, FINAL, PEAKFL, BARFL, SNAPFL,
     &   SRCD, TLT, FREQ, PLTTYP, XS, YS, THETA,
     &   XMAXPOS, YMAXPOS )

C     Plots one of the subplots on a page

#include "fgl.h"

      PARAMETER ( NPLT = 10 )

      LOGICAL SMOOTHFL, GRIDFL, COLFL, FIRST, FINAL, PEAKFL, BARFL,
     &   SNAPFL
      REAL R( * ), RD( * ), TL( * )
      REAL XORM( NPLT ), YORM( NPLT ), XAXLM( NPLT ), YAXLM( NPLT )
      CHARACTER TITLE*80, TITLEX*80, TITLEY*80, TITLEZ*80, TITLE2*80,
     &          COORD*10, PLTTYP*10
      CHARACTER ABC*10

      SAVE FIRST
      DATA FIRST /.TRUE./
      DATA ABC /'abcdefghij'/

C     Default size in cm

C      DATA XORM  /  1.0,  1.0,  1.0,  1.0, 0.2,  0.2 /  ! X origin
C      DATA YORM  /  1.0,  1.0,  1.0,  0.0, 1.0,  1.0 /  ! Y origin
C      DATA XAXLM / 20.0, 20.0, 14.0, 12.0,15.0, 15.0 /  ! x-size
C      DATA YAXLM / 10.0, 10.0,  6.0,  6.0, 7.5,  7.5 /  ! y-size

      DATA XORM  /  1.0,  1.0,  1.0,  1.0,  1.0, 5*0.2 /       ! X origin
      DATA YORM  /  1.0,  1.0,  1.0,  0.0,  1.0,  1.0, 4*0.2 / ! Y origin
      DATA XAXLM / 20.0, 20.0, 14.0, 12.0, 15.0, 15.0, 4*8.0 / ! x-size
      DATA YAXLM / 10.0, 10.0,  6.0,  6.0,  7.5,  7.5, 4*5.0 / ! y-size

C     Calculate size of plot

      IF ( XAXL .EQ. 0.0 ) XAXL = XAXLM( NPLOTS )
      IF ( YAXL .EQ. 0.0 ) YAXL = YAXLM( NPLOTS )

C     Convert to pixels

      PIXCM = 40.0       ! Roughly 40 pixels/cm (98 pixels/inch)

      NFULL = NPLOTS - MOD( NPLOTS, 3 )   ! # plots in full columns

      XOR = PIXCM * ( XORM( NPLOTS ) +
     &               ( IPLOT - 1 ) / 3 * XAXL )

      IF ( IPLOT .LE. NFULL ) THEN
         YOR = PIXCM * ( YORM( NPLOTS ) +
     &          (2 - MOD( IPLOT - 1, 3 ) ) * YAXL )
      ELSE
         YOR = PIXCM * ( YORM( NPLOTS ) +
     &          ( NPLOTS - IPLOT ) * YAXL )
      ENDIF

      XAXL = PIXCM * XAXL
      YAXL = PIXCM * YAXL

      IF ( NPLOTS .GE. 5 .AND. IPLOT .GE. 4 ) THEN
C        XOR = 2.0 * XOR + XAXL
         XOR = XOR - PIXCM * 3.0   ! shift left
      ENDIF

C      IF ( NPLOTS .GE. 5 .AND. IPLOT .LE. 3 ) THEN
C         YOR = PIXCM * YORM( NPLOTS ) + ( 3 - IPLOT ) * YAXL
C      ENDIF

C     *** Titles ***

         IF ( SNAPFL ) THEN
C            WRITE( TITLE2, 8000 ) SRCD, TLT
C 8000       FORMAT( 'Time = ', F6.4, '  Max = ', F4.2, '$' )
            WRITE( TITLE2, 8000 ) SRCD
 8000       FORMAT( 'Time = ', F5.3, '$' )

         ELSE IF ( PLTTYP(1:5) .EQ. 'POLAR' ) THEN
            WRITE( TITLE2, 9001 ) FREQ, SRCD, THETA
 9001       FORMAT( 'FREQ = ', F4.0, '    Sd = ', F4.0,
     &              '   Rd = ', F4.0, '$' )

         ELSE IF ( PLTTYP(1:5) .EQ. 'CARTE' ) THEN
            WRITE( TITLE2, 9002 ) FREQ, SRCD, THETA
 9002       FORMAT( 'FREQ = ', F4.0, '    Sd = ', F4.0,
     &              '   Bearing = ', F4.0, '$' )
         ELSE IF ( PLTTYP(1:5) .EQ. 'OCEAN' ) THEN
            WRITE( TITLE2, 9003 ) SRCD
 9003       FORMAT( 'Depth = ', F4.0, '$' )
         ELSE
            WRITE( TITLE2, 9004 ) FREQ, SRCD
 9004       FORMAT( 'FREQ = ', G9.4, '    Sd = ', G9.4, '$' )
            TITLE2 = ' $'   ! suppress label
         ENDIF

      IF ( FINAL ) THEN
         IF ( NPLOTS .GT. 1 ) THEN
            TITLE = '  '
C            WRITE( TITLE, 7000 ) ABC( IPLOT:IPLOT )
 7000       FORMAT( '(', A1, ')' )
         ENDIF
      END IF

      FIRST = .FALSE.

C     *** Do the plot ***

      IF ( COORD(1:2) .EQ. 'PO' ) THEN
C        --- polar coordinates

         TITLEX = 'X (km)$'
         TITLEY = 'Y (km)$'

         CALL PLOTPO( CVEC, TL, R, NR, RD, NRD,
     &      XMIN, XMAX, XINC, XAXL, XOR, TITLEX,
     &      YMIN, YMAX, YINC, YAXL, YOR, TITLEY,
     &      XS, YS, ZMIN, ZMAX, TITLE, TITLE2 )
      ELSE
C        --- cartesian coordinates
         CALL PLOTCA( CVEC, TL, R, NR, RD, NRD,
     &      XMIN, XMAX, XINC, XAXL, XOR, TITLEX,
     &      YMIN, YMAX, YINC, YAXL, YOR, TITLEY,
     &      ZMIN, ZMAX, TITLE, TITLE2, TITLEZ,
     &      XMAXPOS, YMAXPOS, PEAKFL, BARFL )
      ENDIF

      RETURN
      END

C-----------------------------------------------------------------------

      SUBROUTINE PLOTCA( CVEC, F, X, NX, Y, NY,
     &   XMIN, XMAX, XINC, XAXL, XOR, TITLEX,
     &   YMIN, YMAX, YINC, YAXL, YOR, TITLEY,
     &   FMIN, FMAX, TITLE1, TITLE2, TITLEZ,
     &   XMAXPOS, YMAXPOS, PEAKFL, BARFL )

C     Routine to plot a rectangular grid of numbers
C     Note F is loaded as F( NX, NY )

#include "fgl.h"

      LOGICAL PEAKFL, BARFL
      INTEGER FIRST, LAST
      REAL F( NX, * ), CVEC( 2, NX, * ), X( * ), Y( * ), XY( 2 )
      CHARACTER ALABEL*10, TITLEX*80, TITLEY*80, TITLE1*80, TITLE2*80,
     &          TITLEZ*80

C     Viewport size

      IX1P = XOR
      IX2P = XOR + XAXL
      IY1P = YOR
      IY2P = YOR + YAXL

      CALL VIEWPO( IX1P, IX2P, IY1P, IY2P )

      BORDERXL = 0.15 * ( XMAX - XMIN ) * ( 600.0 / XAXL )
      BORDERXR = 0.40 * ( XMAX - XMIN ) * ( 600.0 / XAXL )
      BORDERY  = 0.20 * ( YMAX - YMIN ) * ( 320.0 / YAXL )

      CALL ORTHO2( XMIN - BORDERXL, XMAX + BORDERXR,
     &             YMAX + BORDERY , YMIN - BORDERY ) ! flip vertically

      CALL READMAP( FIRST, LAST, NCOLS )   ! *** Read the color map ***
C     CALL SHADEM( FLAT )   ! not Gouraud shading

C     *** Construct vectors with xy-coordinates ***

      DO IY = 1, NY
         DO IX = 1, NX
            CVEC( 1, IX, IY ) = X( IX  )
            CVEC( 2, IX, IY ) = Y( IY  )
         END DO
      END DO

      DO IX = 1, NX - 1

C        --- check X( IX ) in [XMIN, XMAX]
         IF ( ABS( X( IX ) - XMIN ) + ABS( X( IX ) - XMAX )
     &       .LE. 1.001 * ABS( XMAX - XMIN ) ) THEN

         CALL BGNTME()    ! Vertices of a triangular mesh
         DO IY = 1, NY

C        --- check Y( IY ) in [YMIN, YMAX]
         IF ( ABS( Y( IY ) - YMIN ) + ABS( Y( IY ) - YMAX )
     &       .LE. 1.001 * ABS( YMAX - YMIN ) ) THEN

            ICOL = FIRST + 0.5 +
     &        ( F( IX, IY )   - FMIN ) * ( NCOLS - 1 ) / ( FMAX - FMIN )
            ICOL = MAX( FIRST, ICOL )
            ICOL = MIN( LAST,  ICOL )

            CALL COLOR( ICOL )
            CALL V2F( CVEC( 1, IX,     IY    ) )

            ICOL = FIRST + 0.5 +
     &        ( F( IX+1, IY ) - FMIN ) * ( NCOLS - 1 ) / ( FMAX - FMIN )
            ICOL = MAX( FIRST, ICOL )
            ICOL = MIN( LAST,  ICOL )

            CALL COLOR( ICOL )
            CALL V2F( CVEC( 1, IX + 1, IY    ) )
         END IF

         END DO
         CALL ENDTME()
         END IF
      END DO

C     *** Draw rectangle around peak ***

      IF ( PEAKFL ) THEN
         SIZEX = ABS( XMAX - XMIN ) / 20.0
         SIZEY = ABS( YMAX - YMIN ) / 20.0

         CALL LINEWI( 3 )
         CALL COLOR( BLACK )
         CALL RECT( XMAXPOS - SIZEX, YMAXPOS - SIZEY,
     &              XMAXPOS + SIZEX, YMAXPOS + SIZEY )
         CALL LINEWI( 1 )
      END IF

C     *** Frame ***

      CALL COLOR( BLACK )

      CALL BGNLIN()
      CALL V2FXY( XMIN, YMIN )
      CALL V2FXY( XMAX, YMIN )
      CALL V2FXY( XMAX, YMAX )
      CALL V2FXY( XMIN, YMAX )
      CALL V2FXY( XMIN, YMIN )
      CALL ENDLIN()
 
C     *** X-axis ***

      WRITE( ALABEL, 9001 ) XMIN   ! NINT( XMIN )
      CALL CMOV2( XMIN, YMAX + 0.5 * BORDERY )
      CALL CHARST( ALABEL, 9 )

      WRITE( ALABEL, 9001 ) XMAX   ! NINT( XMAX )
      CALL CMOV2( XMAX - 0.5 * BORDERXL, YMAX + 0.5 * BORDERY )
      CALL CHARST( ALABEL, 9 )

      ILEN = INDEX( TITLEX, '$' ) - 1
      CALL CMOV2( XMIN + 0.4 * ( XMAX - XMIN ),
     &            YMAX + 0.5 * BORDERY )
      CALL CHARST( TITLEX, ILEN )

C     *** Y-axis ***

      WRITE( ALABEL, 9002 ) NINT( YMIN )
      CALL CMOV2( XMIN - 0.7 * BORDERXL, YMIN )
      CALL CHARST( ALABEL, 9 )

      WRITE( ALABEL, 9002 ) NINT( YMAX )
      CALL CMOV2( XMIN - 0.7 * BORDERXL, YMAX )
      CALL CHARST( ALABEL, 9 )

      ILEN = INDEX( TITLEY, '$' ) - 1
      CALL CMOV2( XMIN - BORDERXL, 0.5 * ( YMIN + YMAX ) )
      CALL CHARST( TITLEY, ILEN )

 9001 FORMAT( F4.2 )
 9002 FORMAT( I4 )

C     *** Plot title ***

      ILEN = 80
      CALL CMOV2( XMIN, YMIN - 0.5 * BORDERY )
      CALL CHARST( TITLE1, ILEN )

      ILEN = INDEX( TITLE2, '$' ) - 1
      CALL CMOV2( XMIN, YMIN - 0.1 * BORDERY )
      CALL CHARST( TITLE2, ILEN )

C     *** Draw color bar ***

      IF ( BARFL )
     &   CALL COLORBAR( XMIN, YMIN, XMAX, YMAX, NCOLS, FIRST, LAST,
     &              BORDERXR, FMIN, FMAX, TITLEZ )

      RETURN
      END

C**********************************************************************C

      SUBROUTINE PLOTPO( CVEC, F, X, NX, Y, NY,
     &   XMIN, XMAX, XINC, XAXL, XOR, TITLEX,
     &   YMIN, YMAX, YINC, YAXL, YOR, TITLEY,
     &   XS, YS, FMIN, FMAX, TITLE1, TITLE2 )

C     Routine to plot a rectangular grid of numbers
C     Note F is loaded as F( NX, NY )

#include "fgl.h"

      PARAMETER( DEGRAD = 3.141592 / 180.0 )

      INTEGER FIRST, LAST
      REAL F( NX, * ), CVEC( 2, NX, * ), X( * ), Y( * ), XY( 2 )
      CHARACTER ALABEL*10, TITLEX*80, TITLEY*80, TITLE1*80, TITLE2*80

C     Viewport size

      IX1P = XOR
      IX2P = XOR + XAXL
      IY1P = YOR
      IY2P = YOR + YAXL

      CALL VIEWPO( IX1P, IX2P, IY1P, IY2P )

      BORDERX = 0.30 * ABS( XMAX - XMIN )
      BORDERY = 0.25 * ABS( YMAX - YMIN )

      CALL ORTHO2( XMIN - BORDERX, XMAX + BORDERX,
     &             YMIN - BORDERY, YMAX + BORDERY ) ! flip vertically

      CALL READMAP( FIRST, LAST, NCOLS )   ! *** Read the color map ***
      CALL COLOR( WHITE )

C     *** Construct vectors with xy-coordinates ***
C     ( X is a range and Y is an angle )

      DO IY = 1, NY
         DO IX = 1, NX
            XT = XS + X( IX ) * COS( DEGRAD * Y( IY ) )
            YT = YS + X( IX ) * SIN( DEGRAD * Y( IY ) )

            XT = MIN( XT, XMAX )
            CVEC( 1, IX, IY ) = MAX( XT, XMIN )

            YT = MIN( YT, YMAX )
            CVEC( 2, IX, IY ) = MAX( YT, YMIN )
         END DO
      END DO

      DO IX = 1, NX - 1

         CALL BGNTME()    ! Vertices of a triangular mesh

         DO IY = 1, NY


            ICOL = FIRST + 0.5 +
     &        ( F( IX, IY )   - FMIN ) * ( NCOLS - 1 ) / ( FMAX - FMIN )
            ICOL = MAX( FIRST, ICOL )
            ICOL = MIN( LAST,  ICOL )

            CALL COLOR( ICOL )
            CALL V2F( CVEC( 1, IX    , IY ) )

            ICOL = FIRST + 0.5 +
     &        ( F( IX+1, IY ) - FMIN ) * ( NCOLS - 1 ) / ( FMAX - FMIN )
            ICOL = MAX( FIRST, ICOL )
            ICOL = MIN( LAST,  ICOL )

            CALL COLOR( ICOL )
            CALL V2F( CVEC( 1, IX + 1, IY ) )

         END DO
         CALL ENDTME()
      END DO

C     *** Frame ***

C     CALL COLOR( WHITE )
      CALL COLOR( BLACK )

      CALL BGNLIN()
         CALL V2FXY( XMIN, YMIN )
         CALL V2FXY( XMAX, YMIN )
         CALL V2FXY( XMAX, YMAX )
         CALL V2FXY( XMIN, YMAX )
         CALL V2FXY( XMIN, YMIN )
      CALL ENDLIN()
 
C     *** X-axis ***

      WRITE( ALABEL, 9001 ) XMIN
      CALL CMOV2( XMIN, YMIN - 0.5 * BORDERY )
      CALL CHARST( ALABEL, 9 )

      WRITE( ALABEL, 9001 ) XMAX
      CALL CMOV2( XMAX - 0.5 * BORDERX, YMIN - 0.5 * BORDERY )
      CALL CHARST( ALABEL, 9 )

      ILEN = INDEX( TITLEX, '$' ) - 1
      CALL CMOV2( 0.5 * ( XMIN + XMAX ), YMIN - 0.5 * BORDERY )
      CALL CHARST( TITLEX, ILEN )

C     *** Y-axis ***

      WRITE( ALABEL, 9001 ) YMIN
      CALL CMOV2( XMIN - 0.7 * BORDERX, YMIN )
      CALL CHARST( ALABEL, 9 )

      WRITE( ALABEL, 9001 ) YMAX
      CALL CMOV2( XMIN - 0.7 * BORDERX, YMAX )
      CALL CHARST( ALABEL, 9 )

      ILEN = INDEX( TITLEY, '$' ) - 1
      CALL CMOV2( XMIN - BORDERX, 0.5 * ( YMIN + YMAX ) )
      CALL CHARST( TITLEY, ILEN )

 9001 FORMAT( F5.0 )

C     *** Plot title ***

      ILEN = 80
      CALL CMOV2( XMIN, YMAX + 0.5 * BORDERY )
      CALL CHARST( TITLE1, ILEN )

      ILEN = INDEX( TITLE2, '$' ) - 1
      CALL CMOV2( XMIN, YMAX + 0.1 * BORDERY )
      CALL CHARST( TITLE2, ILEN )

C     *** Color bar ***

      CALL BGNTME()    ! Vertices of a triangular mesh

      DO ICOL = FIRST, LAST
         CALL COLOR( ICOL )

         XY( 1 ) = XMAX + 0.3 * BORDERX
         XY( 2 ) = YMAX - ( ICOL - FIRST + 1 ) * ( YMAX - YMIN ) / NCOLS
         CALL V2F( XY )

         XY( 1 ) = XMAX + 0.6 * BORDERX
         XY( 2 ) = YMAX - ( ICOL - FIRST + 1 ) * ( YMAX - YMIN ) / NCOLS
         CALL V2F( XY )
      END DO

      CALL ENDTME()

C     *** Labels for color bar ***

C     CALL COLOR( WHITE )
      CALL COLOR( BLACK )

      WRITE( ALABEL, 9002 ) FMIN
      CALL CMOV2( XMAX + 0.7 * BORDERX, YMAX  )
      CALL CHARST( ALABEL, 9 )

      WRITE( ALABEL, 9002 ) FMAX
      CALL CMOV2( XMAX + 0.7 * BORDERX, YMIN  )
      CALL CHARST( ALABEL, 9 )

 9002 FORMAT( F4.0 )

      CALL CMOV2( XMAX + 0.7 * BORDERX, 0.5 * ( YMIN + YMAX ) )
      CALL CHARST( ' dB', 3 )

      RETURN
      END


C**********************************************************************C

      SUBROUTINE COLORBAR( XMIN, YMIN, XMAX, YMAX, NCOLS, FIRST, LAST,
     &   BORDERXR, FMIN, FMAX, TITLEZ )

C     Draws the color bar

#include "fgl.h"

      INTEGER FIRST, LAST
      REAL XY( 2 )
      CHARACTER ALABEL*10, TITLEZ*80

      DATA ICOUNT / 0 /

      ICOUNT = ICOUNT + 1
C     IF ( ICOUNT .LT. 4 ) RETURN   ! color bar for last figure only

      CALL BGNTME()    ! Vertices of a triangular mesh

      DO ICOL = FIRST, LAST
         CALL COLOR( ICOL )

         XY( 1 ) = XMAX + 0.2 * BORDERXR
         XY( 2 ) = YMIN + ( ICOL - FIRST + 1 ) * ( YMAX - YMIN ) / NCOLS
         CALL V2F( XY )

         XY( 1 ) = XMAX + 0.4 * BORDERXR
         XY( 2 ) = YMIN + ( ICOL - FIRST + 1 ) * ( YMAX - YMIN ) / NCOLS
         CALL V2F( XY )
      END DO

      CALL ENDTME()

C     *** Labels for color bar ***

      CALL COLOR( WHITE )
      CALL COLOR( BLACK )

      WRITE( ALABEL, 9002 ) NINT( FMIN )
      IF ( FMIN .EQ. 0.0 ) ALABEL = ' 0'
      CALL CMOV2( XMAX + 0.5 * BORDERXR, YMIN  )
      CALL CHARST( ALABEL, 9 )

      WRITE( ALABEL, 9002 ) NINT( FMAX )
      IF ( FMAX .EQ. 0.0 ) ALABEL = ' 0'
      CALL CMOV2( XMAX + 0.5 * BORDERXR, YMAX  )
      CALL CHARST( ALABEL, 9 )

 9001 FORMAT( G8.3 )  ! ***
 9002 FORMAT( I4 ) ! G9.3  ! ***

      ILEN = INDEX( TITLEZ, '$' ) - 1
      CALL CMOV2( XMAX + 0.5 * BORDERXR, 0.5 * ( YMIN + YMAX ) )
      CALL CHARST( TITLEZ, ILEN )

      RETURN
      END

C**********************************************************************C

     SUBROUTINE V2FXY( X, Y )

C     Calls V2F for a point (x,y)

      REAL CVEC( 2)

      CVEC( 1 ) = X
      CVEC( 2 ) = Y

      CALL V2F( CVEC )

      RETURN
      END

C**********************************************************************C

      SUBROUTINE READMAP( FIRST, LAST, NCOLS )

C     *** Read and load the color map ***

      PARAMETER( MAPFIL = 77 )

      INTEGER FIRST, LAST

      FIRST = 9999
      LAST  = 1

      OPEN( FILE = 'MAPFIL', UNIT = MAPFIL, STATUS = 'OLD' )
      DO I = 1, 9999
         READ( MAPFIL, *, END = 700 ) ICOL, IRED, IGREEN, IBLUE
         CALL MAPCOL( ICOL, IRED, IGREEN, IBLUE )
         FIRST = MIN( ICOL, FIRST )
         LAST  = MAX( ICOL, LAST  )
      END DO

  700 CLOSE( MAPFIL )
      NCOLS = LAST - FIRST + 1

      CALL SYSTEM( "makemap" )  ! ensure color table loaded

      RETURN
      END
