      SUBROUTINE EVALPDQ( CK, PHIR, PHIS, M, IELTSC, ADJELT, ISET,
     &   X, Y, NNODES, NODE, NELTS, XS, YS, THETA, NTHETA,
     &   RMINM, RMAXM, NR, MACT, OPT, P )

C     Computes 3-D pressure field using adiabatic mode theory
C     Normalized to pressure of point source at 1 meter

C     Note RMINM must be zero

      COMPLEX CI
      PARAMETER ( PI = 3.141592, DEGRAD = PI/180.0, CI = ( 0.0, 1.0 ),
     &   MAXM = 100 )

      INTEGER M( * ), NODE( 3, * ), ADJELT( 3, * ), ISET( * ), OUTSIDE
      REAL    X( * ), Y( * ), THETA( * )
      COMPLEX PHIR( MAXM, * ), PHIIN( MAXM ), PHIOUT( MAXM ),
     &        PHIS( MAXM, * ), COEF( MAXM ),
     &        CK(   MAXM, * ), CKIN( MAXM ), CKOUT( MAXM ),
     &        P( NTHETA, * ), PHASEINC( MAXM ), T, CKAVG
      CHARACTER OPT*( * )

C     *** Loop over angle ***

      DELTAR = ( RMAXM - RMINM ) / ( NR - 1 )

      DO ITHETA = 1, NTHETA
         TSX = COS( DEGRAD * THETA( ITHETA ) )
         TSY = SIN( DEGRAD * THETA( ITHETA ) )

C        --- Get modal values
         IELT = IELTSC
         WRITE( *, * )
         WRITE( *, * ) 'ITHETA, TSX, TSY', ITHETA, TSX, TSY

         CALL FIRST( IELT, OUTSIDE, RIN, ROUT, XS, YS, TSX, TSY,
     &      X, Y, MPROP, M, MAXM, CK, PHIR, PHIS, COEF, NODE,
     &      ISET, CKIN, PHIIN, CKOUT, PHIOUT )

         MPROP = MIN( MACT, MPROP )

         DO MODE = 1, MPROP
            COEF( MODE ) = SQRT( 2.0 * PI ) * EXP( CI * PI / 4.0 )
     &                     * COEF( MODE )
            CKAVG = 0.5 * ( CKIN( MODE ) + CKOUT( MODE ) )
            PHASEINC( MODE ) = EXP( -CI * CKAVG * DELTAR )
         END DO

C        *** March forward in range ***

         DO IR = 1, NR
            RM = RMINM + ( IR - 1 ) * DELTAR
            IF ( RM .EQ. 0.0 ) RM = DELTAR

C           Crossing into new element?

  100       CONTINUE
            IF ( RM .GT. ROUT ) THEN

C              --- Copy outside info to inside
               NEWELT = ADJELT( OUTSIDE, IELT )
               RIN = ROUT

               DO MODE = 1, MPROP
                  CKIN(  MODE ) =  CKOUT( MODE )
                  PHIIN( MODE ) = PHIOUT( MODE )
                  CKAVG = 0.5 * ( CKIN( MODE ) + CKOUT( MODE ) )
                  PHASEINC( MODE ) = EXP( -CI * CKAVG * DELTAR )
               END DO

C              --- Get new outside info
               CALL OUT( IELT, NEWELT, OUTSIDE, ROUT, ADJELT, ISET,
     &            XS, YS, TSX, TSY, X, Y, NODE,
     &            MPROP, M, MAXM, CK, PHIR, CKOUT, PHIOUT )
               IELT = NEWELT

               GOTO 100   ! Check again

            END IF

C           *** Compute modal contribution at this range ***

            T = 0.0
            DO L = 1, MPROP
               COEF( L ) = COEF( L ) * PHASEINC( L )
               T = T + COEF( L ) * PHIIN( L )
            END DO

            P( ITHETA, IR ) = T / SQRT( RM * CKIN( 1 ) )

         END DO   ! next range

      END DO   ! next bearing

      RETURN
      END
