      SUBROUTINE EVALAD( RPROF, NPROF, PHIS, PHIR, RD, NRD, 
     &   R, NR, M, OPT, P )

C     Computes pressure field using adiabatic mode theory
C     Normalized to pressure at 1 meter

C     OPT
C
C     X     Cartesian   (x, z) coordinates
C     T     Translationally invariant ocean
C     R     Cylindrical (r, z) coordinates
C     S     Scaled cylindrical coordinates ( 1/r fall-off removed )

      COMPLEX CI
      PARAMETER ( MAXM = 1500, MAXNRD = 101,
     &    CI = ( 0.0, 1.0 ), PI = 3.1415926 )

      REAL RD( * ), RPROF( * ), R( * )
      COMPLEX PHIS( * ), PHIL( MAXM, MAXNRD ), PHIR( MAXM, * ), 
     &   CONST( MAXM ), HANK( MAXM ),
     &   SUM, CKINT, PHIINT, 
     &   P( NRD, * ), CKL( MAXM ), CKR( MAXM ), CKMID

      COMPLEX*16 SUMK( MAXM ), SUMKINV( MAXM )
      CHARACTER OPT*( * ), TITLE*80

      IF ( NRD. GT. MAXNRD )
     &   STOP 'FATAL ERROR: NRD > MAXNRD in routine EVALAD'

C     *** Initialization ***

      RPROF( NPROF + 1 ) = 1.0E9
      IPROF = 1

C     ------ Receiver depths at left  of segment
      CALL GETMOD( IPROF, ' ', MAXM, RD, NRD, 'N',
     &                       CKL, PHIL, M1, FREQ, TITLE )
      M = MIN( M, M1 )
C     ------ Receiver depths at right of segment
      CALL GETMOD( IPROF+1, ' ', MAXM, RD, NRD, 'N',
     &                       CKR, PHIR, M1, FREQ, TITLE )
      M = MIN( M, M1 )

      DO MODE = 1, M
         CONST( MODE ) = CI * SQRT( 2.0 * PI ) * EXP( CI * PI / 4.0 )
     &                   * PHIS( MODE )
         SUMK(    MODE ) = 0.0
         SUMKINV( MODE ) = 0.0
         IF ( OPT(1:1) .EQ. 'T' )
     &      CONST( MODE ) = CONST( MODE ) / SQRT( CKL( MODE ) )
      END DO

C     *** March forward in range ***

      DO IR = 1, NR

C        ------ Crossing into new range segment?

   50    CONTINUE
         IF ( R( IR ) .GT. 1000.0 * RPROF( IPROF + 1 ) ) THEN

            CALL NEWSEG( R, IR, RPROF, IPROF, NPROF,
     &         CKL, CKR, PHIL, PHIR, M, RD, NRD, SUMK, SUMKINV )

            GO TO 50   ! Check again
         ENDIF

C        ------ Compute proportional distance W and interpolate

         IF ( IR .GT. 1 ) THEN
            RLEFT = MAX( R( IR - 1 ), 1000.0 * RPROF( IPROF ) )
         ELSE
            RLEFT = 1000.0 * RPROF( IPROF )
         ENDIF
         RMID = 0.5 * ( R( IR ) + RLEFT )

         W    = ( R( IR ) / 1000.0 - RPROF( IPROF ) ) /
     &                  ( RPROF( IPROF + 1 ) - RPROF( IPROF ) )
         WMID = ( RMID    / 1000.0 - RPROF( IPROF ) ) /
     &                  ( RPROF( IPROF + 1 ) - RPROF( IPROF ) )

         DO L = 1, M

            CKINT = CKL( L ) + W    * ( CKR( L ) - CKL( L ) )
            CKMID = CKL( L ) + WMID * ( CKR( L ) - CKL( L ) )
            SUMK(    L ) = SUMK(    L ) + CKMID * ( R( IR ) - RLEFT )
            SUMKINV( L ) = SUMKINV( L ) + ( R( IR ) - RLEFT ) / CKMID

            IF ( OPT(4:4) .NE. 'I' ) THEN   ! coherent case
               HANK( L ) = CONST( L ) *       EXP( -CI * SUMK( L ) )
            ELSE
               HANK( L ) = CONST( L ) * EXP( REAL( -CI * SUMK( L ) ) ) ! inc. case
            ENDIF

            IF ( OPT(1:1) .EQ. 'R' ) THEN             ! Cylindrical coords.
               HANK( L ) = HANK( L ) / SQRT( CKINT * R( IR ) )
            ELSE IF ( OPT(1:1) .EQ. 'X' ) THEN        ! Cartesian coords.
               HANK( L ) = HANK( L ) / CKINT
            ELSE IF ( OPT(1:1) .EQ. 'T' ) THEN        ! Translationally invariant
               HANK( L ) = HANK( L ) / SQRT( CKINT * SUMKINV( L ) )
            ELSE                                      ! Scaled cylindrical coords.
               HANK( L ) = HANK( L ) / SQRT( CKINT )
            ENDIF

         END DO

C        *** For each rcvr, add up modal contributions ***

         DO I = 1, NRD
            SUM = 0.0

            IF ( OPT(4:4) .NE. 'I' )  THEN   ! coherent case
               DO MODE = 1, M
                  PHIINT = PHIL( MODE, I ) +
     &                      W * ( PHIR( MODE, I ) - PHIL( MODE, I ) )
                  SUM = SUM + PHIINT * HANK( MODE )
               END DO
            ELSE                                ! incoherent case
               DO MODE = 1, M
                  PHIINT = PHIL( MODE, I ) +
     &                      W * ( PHIR( MODE, I ) - PHIL( MODE, I ) )
                  SUM = SUM + ABS(  PHIINT * HANK( MODE ) ) ** 2
               END DO
               SUM = SQRT( SUM )
            ENDIF

            P( I, IR ) = SUM
         END DO

      END DO   ! Next range step

      RETURN
      END
C**********************************************************************C
      SUBROUTINE NEWSEG( R, IR, RPROF, IPROF, NPROF,
     &   CKL, CKR, PHIL, PHIR, M, RD, NRD, SUMK, SUMKINV )

C     Treats the crossing into a new range segment

      PARAMETER ( MAXM = 1500 )

      REAL RD( * ), RPROF( * ), R( * )
      COMPLEX PHIL( MAXM, * ), PHIR( MAXM, * ), 
     &   CKL( * ), CKR( * ), CKMID
      COMPLEX*16 SUMK( * ), SUMKINV( * )
      CHARACTER TITLE*80

C     *** Do phase integral up to the new range ***

      IF ( IR .GT. 1 ) THEN
         RLEFT = MAX( R( IR - 1 ), 1000.0 * RPROF( IPROF ) )
      ELSE
         RLEFT = 1000.0 * RPROF( IPROF )
      ENDIF

      RMID = 0.5 * ( 1000.0 * RPROF( IPROF + 1 ) + RLEFT )
      WMID = ( RMID / 1000.0 - RPROF( IPROF ) ) /
     &                  ( RPROF( IPROF + 1 ) - RPROF( IPROF ) )

      DO L = 1, M
         CKMID     = CKL(  L ) + WMID  * ( CKR( L ) - CKL( L ) )
         SUMK( L ) = SUMK( L ) +
     &           CKMID * ( 1000.0 * RPROF( IPROF + 1 ) - RLEFT )
         SUMKINV( L ) = SUMKINV( L ) +
     &       ( 1000.0 * RPROF( IPROF + 1 ) - RLEFT ) / ckmid
      END DO

C     ***  Copy right modes to left ***

      DO MODE = 1, M
         CKL( MODE ) = CKR( MODE )
         DO J = 1, NRD
            PHIL( MODE, J ) = PHIR( MODE, J )
         END DO
      END DO

C     *** Read in the new right mode set ***

      IPROF = IPROF + 1

      IF ( IPROF + 1  .LE. NPROF ) THEN
         CALL GETMOD( IPROF+1, ' ', MAXM, RD, NRD, 'N',
     &                           CKR, PHIR, M1, FREQ, TITLE )
         M = MIN( M, M1 )
         WRITE( *, * ) 'New profile read',
     &       R(IR), IPROF+1, RPROF(IPROF+1), ' #modes=', M
      ENDIF


      RETURN
      END
