      SUBROUTINE EVAL( C, PHI, NZ, R, NR, RR,
     &   CK, M, OPTION, P )

C     Given modes and wavenumbers, compute pressure field
C     Normalized to pressure of point source at 1 meter
C     OPTION = X     Cartesian   (x, z) coordinates
C     OPTION = R     Cylindrical (r, z) coordinates

      COMPLEX CI
      PARAMETER ( MAXM = 1500, CI = (0.0, 1.0), PI = 3.1415926, 
     &   MAXNZ = 251, MINEXP = -37 )

      REAL RR( * ), R( * )
      COMPLEX C( * ), PHI( MAXM, * ), P( NZ, * ), T, HANK( MAXM ),
     &   CICK( MAXM ), CONST( MAXM ), CMAT( MAXM, MAXNZ ), CK( * ),
     &   FACTOR
      COMPLEX*16 PHASE
      CHARACTER OPTION*(*)

C     *** If no modes, return vanishing pressure ***

      IF ( M .LE. 0 ) THEN
         DO IR = 1, NR
            DO IZ = 1, NZ
               P( IZ, IR ) = 0.0
            END DO
         END DO
         RETURN
      END IF

C     Initialization

      FACTOR = CI * SQRT( 2.0 * PI ) * EXP( CI * PI / 4.0 )

      DO MODE = 1, M

         IF ( OPTION(1:1) .EQ. 'X' ) THEN
            CONST( MODE ) = FACTOR * C( MODE ) / CK( MODE )
         ELSE
            CONST( MODE ) = FACTOR * C( MODE ) / SQRT( CK( MODE ) )
         ENDIF

         CICK( MODE ) = -CI * CK( MODE )   ! use e{i(wt-kr)} form
         IF ( OPTION(4:4) .EQ. 'I' ) CICK( MODE ) = REAL( CICK( MODE ) )   ! Incoherent case

      END DO

      DO IZ = 1, NZ

         DO MODE = 1, M
            CMAT( MODE, IZ ) = CONST( MODE ) * 
     &         PHI( MODE, IZ ) * EXP( CICK( MODE ) * RR( IZ ) )
         END DO

      END DO

C     *** Loop over range ***

      DO IR = 1, NR

         DO MODE = 1, M

C           --- eliminate underflows (can raise CPU time)
            IF ( REAL( CICK( MODE ) * R( IR ) ) .LT. MINEXP ) THEN
               HANK( MODE ) = 0.0
            ELSE
C              HANK( MODE ) = EXP( CICK( MODE ) * R( IR ) )
C              handle an SGI bug in sin( large number )
               PHASE = CICK( MODE ) * R( IR )
               HANK( MODE ) = EXP( PHASE )
            ENDIF
         END DO

C        *** Loop over depth ***

         DO IZ = 1, NZ
            T = 0.0
            IF ( OPTION(4:4) .NE. 'I' )  THEN   ! coherent case
               DO MODE = 1, M
                  T = T + CMAT( MODE, IZ ) * HANK( MODE )
               END DO
            ELSE                                ! incoherent case
               DO MODE = 1, M
                  T = T + ABS(  CMAT( MODE, IZ ) * HANK( MODE ) ) ** 2
               END DO
               T = SQRT( T )
            ENDIF

C           ------ Cylindrical or cartesian coordinates?
            IF ( OPTION(1:1) .EQ. 'R' )
     &         T = T / SQRT( R( IR ) + RR( IZ ) )
            P( IZ, IR ) = T
         END DO

      END DO   ! next range step

      RETURN
      END
