      PROGRAM PLTTLR

C     Plot the propagation loss along a specified track

      INTEGER PLPFIL, PRTFIL
      PARAMETER ( PLPFIL = 1, PRTFIL = 6, MAXM = 1500, MAXN = 5000,
     &   MAXNR = 5000, MAXNPROF = 51 )

      REAL      R( MAXNR ), RM( MAXNR), TLRAW( MAXNR ),
     &          TL( MAXNR ), RPROF( MAXNPROF ), RR( 1 )
      COMPLEX   CK( MAXM ), PHIS( MAXM ), PHIR( MAXM ), P( MAXNR )
      CHARACTER OPT*50, TITLE*80, COMP*1

      DATA RMIN, RMAX, TLMIN, TLMAX / 0.0, 100.0, 40.0, 120.0 /

C     Read in control data

      OPEN( FILE = 'PLPFIL', UNIT = PLPFIL, STATUS = 'OLD' )
      READ( PLPFIL, * ) OPT
      COMP = OPT(3:3)

      M = 99999
      READ( PLPFIL, * ) M

C     ------ Read profile ranges
      RPROF( 3 ) = -999.9
      READ( PLPFIL, * ) NPROF, ( RPROF( I ), I = 1, NPROF )

      IF ( NPROF .GT. MAXNPROF ) THEN
         WRITE( PRTFIL, * ) 'NPROF, MAXNPROF = ', NPROF, MAXNPROF
         CALL ERROUT( PRTFIL, 'F', 'PLOTTLR', 'Too many profiles' )
      ENDIF

      CALL SUBTAB( RPROF, NPROF )

      READ( PLPFIL, * ) SRCD, RCVRD
      READ( PLPFIL, * ) RMIN, RMAX, NR

      NSD = 1
      NRD = 1
      RR( 1 ) = 0.0

C     Compute range vector

      DELTAR = ( RMAX - RMIN ) / ( NR - 1 )
      DO J = 1, NR
         R(  J ) = RMIN + ( J - 1 ) * DELTAR
         RM( J ) = 1000.0 * R( J )
      END DO

      IF ( RM( 1 ) .LE. 0.0 ) RM( 1 ) = MIN( 1.0, RM( 2 ) )
 
C     Read in modes

      IPROF = 1

      CALL GETMOD( IPROF, ' ', MAXM, SRCD,  NSD,  'N',
     &   CK, PHIS, MSRC, FREQ, TITLE )

      CALL GETMOD( IPROF, ' ', MAXM, RCVRD, NRD, COMP,
     &   CK, PHIR, MSRC, FREQ, TITLE )

      M = MIN( M, MSRC )

C     *** Evaluate the field ***

      IF ( NPROF .EQ. 1 ) THEN                  ! Range independent
         CALL EVAL( PHIS, PHIR, NRD,
     &                 RM, NR, RR, CK, M, OPT, P )
      ELSE IF ( OPT(2:2) .EQ. 'C' ) THEN        ! Coupled mode theory
         CALL EVALCM( RPROF, NPROF, PHIS, PHIR, RCVRD, NRD,
     &                     RM, NR, CK, M, OPT, P )
      ELSE                                      !Adiabatic mode theory
         CALL EVALAD( RPROF, NPROF, PHIS, PHIR, RCVRD, NRD, 
     &                     RM, NR, M, OPT, P )
      ENDIF

C     *** Compute scales for TL axis ***

      TLMINT = +1.0E9
      TLMAXT = -1.0E9

      DO J = 1, NR
         PSQ = P( J ) * CONJG( P( J ) )
         TLRAW( J ) = 200.0
         IF ( PSQ .GT. 1.0E-20 ) TLRAW( J ) = -10.0 * ALOG10( PSQ )
         TLMINT = MIN( TLMINT, TLRAW( J ) )
         TLMAXT = MAX( TLMAXT, TLRAW( J ) )
      END DO

      READ( PLPFIL, * ) TLMIN, TLMAX

C     --- Compute max, min automatically?
      IF ( TLMIN .EQ. TLMAX ) THEN
        TLMIN = 10.0 * INT( TLMINT / 10.0 )
        TLMAX = 10.0 * INT( TLMAXT / 10.0 ) + 10.0
      END IF

C     Smooth

      READ( PLPFIL, * ) DR3DB
      DELTAR = 1000.0 * DELTAR
      NROWS = 1
      CALL SMOOTH( TLRAW, TL, NROWS, NRD, NR, DELTAR, DR3DB )

C     *** Write data ***

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

      WRITE( PRTFIL, 2000 ) NINT( FREQ ), NINT( SRCD ), NINT( RCVRD )
 2000 FORMAT(  '@ subtitle "',
     &         'F  =', I4, ' Hz, ',
     &         'Sd =', I4, ' m,  ',
     &         'Rd =', I4, ' m"' )
 
      WRITE( PRTFIL, 3000 ) INT( RMIN  ), INT( RMAX  ),
     &                      INT( TLMIN ), INT( TLMAX )
 3000 FORMAT( '@ world xmin ', I5, /,
     &        '@ world xmax ', I5, /,
     &        '@ world ymin ', I5, /,
     &        '@ world ymax ', I5 )

      DO IR = 1, NR
         WRITE( PRTFIL, * ) R( IR ), TL( IR )
      END DO

      STOP

      END
