      SUBROUTINE EVALCM( RPROF, NPROF, PHIS, PHI, RD, NRD, 
     &   R, NR, CK, M, OPT, P )

C     Computes pressure field using coupled mode theory
C     Normalized to pressure of point source at 1 meter

C     OPT = X     Cartesian   (X, Z) coordinates
C     OPT = R     Cylindrical (R, Z) coordinates

C     Note number of propagating modes is reset after first segment.
C     Thus M restricts the number of modes in the source field but
C     thereafter energy is allowed to couple into higher-order modes.

C     Note flaws:
C        should half-space contribution include imaginary part?
C        Then we need to handle the growing part of leaky waves.

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

      LOGICAL FIRST
      REAL      RD( * ), RPROF( * ), R( * )
      COMPLEX   PHIS( * ), PHI( MAXM, * ), A( MAXM ), SUM,
     &          P( NRD, * ), CK( * )
      CHARACTER OPT*( * ), TITLE*80
      SAVE FIRST
      DATA FIRST /.TRUE./

C     *** Compute ranges (in meters) where new profiles are used ***

      IF ( FIRST ) THEN
         IF ( RPROF( 1 ) .NE. 0.0 )
     &      STOP 'FATAL ERROR: First profile must start at zero range'

         DO IPROF = NPROF, 2, -1
            RPROF( IPROF ) = 500.0*( RPROF( IPROF ) + RPROF( IPROF-1 ) )
         END DO

         RPROF(  NPROF + 1 ) = 1.0E9
         FIRST = .FALSE.
      ENDIF

C     *** Evaluate mode excitation coefficients, A(mode) ***

      IPROF = 1

C     THIS HAS ALREADY BEEN READ IN PLOTTLR!

      CALL GETMOD( IPROF, ' ', MAXM, RD, NRD, 'N',
     &                         CK, PHI,  M1, FREQ, TITLE )

      M = MIN( M, M1 )

      DO MODE = 1, M
         IF ( OPT(1:1) .EQ. 'X' ) THEN
            A( MODE ) = SQRT( 2.0*PI ) * EXP( CI * PI / 4.0 ) *
     &                  PHIS( MODE ) / CK( MODE )
         ELSE
C           ------ 'R' Cylindrical coordinates
            A( MODE ) = CI * SQRT( 2.0*PI ) * EXP( CI * PI / 4.0 ) *
     &                  PHIS( MODE ) / SQRT( CK(MODE) )
         ENDIF
      END DO

C     *** March forward in range ***

      DO IR = 1, NR
C        WRITE( *, * ) 'Range = ', R( IR )

C        *** Crossing into new range segment? ***

         IF ( R( IR ) .GT. RPROF( IPROF + 1 ) ) THEN

            IPROF = IPROF + 1

C           ------ Advance to interface
            IF ( IR .EQ. 1 ) THEN
               CALL ADVPHA( RPROF( IPROF )              , A, CK, M )
            ELSE
               CALL ADVPHA( RPROF( IPROF ) - R( IR - 1 ), A, CK, M )
            ENDIF

C           ------ Here's where we cross over
            IF ( IPROF  .LE. NPROF ) THEN
               CALL NEWPRO( IPROF, CK, PHI, M, RD, NRD, A )
               WRITE( *, * ) 'New profile read', R(I R ), IPROF,
     &                       ' #modes = ', M
            ENDIF

C           ------ Are there other segments to cross?
   50       CONTINUE
            IF ( R( IR ) .GT. RPROF( IPROF + 1 ) ) THEN
               IPROF = IPROF + 1
               CALL ADVPHA( RPROF( IPROF ) - RPROF( IPROF - 1 ),
     &                                                   A, CK, M )
               IF ( IPROF  .LE. NPROF )
     &            CALL NEWPRO( IPROF, CK, PHI, M, RD, NRD, A )

               GOTO 50       ! Check again
            ENDIF

C           ------ Advance the remaining distance
            CALL ADVPHA( R( IR ) - RPROF( IPROF ), A, CK, M )

         ELSE

            IF ( IR .EQ. 1 ) THEN
               CALL ADVPHA( R( IR )              , A, CK, M )
            ELSE
               CALL ADVPHA( R( IR ) - R( IR - 1 ), A, CK, M )
            ENDIF

         ENDIF

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

         DO IRD = 1, NRD

            SUM = 0.0
            DO MODE = 1, M
               SUM = SUM + A( MODE ) * PHI( MODE, IRD )
            END DO

            IF ( OPT(1:1) .EQ. 'R' .AND. R( IR ) .NE. 0.0 ) THEN
               P( IRD, IR ) = SUM / SQRT( R( IR ) )
            ELSE
               P( IRD, IR ) = SUM
            ENDIF

         END DO

C        ------ Next range step
      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE ADVPHA( DELTAR, A, CK, M )

C     *** ADVance the PHAse of the coefficients ***

      COMPLEX CI
      PARAMETER ( CI = ( 0.0, 1.0 ) )

      COMPLEX A( * ), CK( * )

      DO L = 1, M
         A( L ) = A( L ) * EXP( -CI * CK( L ) * DELTAR )
      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE NEWPRO( IPROF, CK, PHIR, MR, RD, NRD, A )

C     For a given profil number:
C         read in modes for current segment
C         project the pressure field onto the new modes
C         extract values of the modes at rcvr depths

      PARAMETER ( MAXN = 16001, MAXM = 1500, MAXNRD = 201 )

      INTEGER    IRD( MAXNRD )
      REAL       Z( MAXN ), RD( * ), W( MAXNRD )
      COMPLEX    CK( * ), PHI( MAXN ), PHIR( MAXM, * ),
     &           P( MAXN ), A( * ), SUM,
     &           GAMTL( MAXM ), GAMBL( MAXM ),
     &           PHITL( MAXM ), PHIBL( MAXM ),
     &           GAMTR, GAMBR, PHITR, PHIBR,
     &           KTOP2R, KBOT2R, TAIL, PHITMP
      COMPLEX*16 PEKRT, GAMMA2
      CHARACTER  BCBOTR * 1, BCTOPR * 1

      IF ( NRD. GT. MAXNRD )
     &   STOP 'NRD > MAXNRD in routine NEWPRO'

C     *** Compute pressure along the left of the interface ***

      CALL PLEFT( IPROF, A, CK, Z, PHI, MR, P, NR, NTOT,
     &   BCTOPR, RHOTR, KTOP2R, DEPTHTR,
     &   BCBOTR, RHOBR, KBOT2R, DEPTHBR,
     &   GAMTL, GAMBL, DEPTHTL, DEPTHBL, PHITL, PHIBL, ML, MAXM )

C     ***  Read in eigenfunctions and extract receiver values ***

C     Compute weights for mode interpolation at rcvr depths

      CALL WEIGHT( Z, NTOT, RD, NRD, W, IRD )

      sum2 = 0.0
      MODFIL = 29 + IPROF

      DO MODE = 1, MR
         READ( MODFIL, REC = 6 + MODE ) ( PHI( IZ ), IZ = 1, NR )

         IF ( BCTOPR .EQ. 'A' ) THEN
            PHITR = PHI( 1 )
            GAMMA2 = CK( MODE ) ** 2 - KTOP2R
            GAMTR = PEKRT( GAMMA2 )
         ENDIF

         IF ( BCBOTR .EQ. 'A' ) THEN
            PHIBR = PHI( NR )
            GAMMA2 = CK( MODE ) ** 2 - KBOT2R
            GAMBR = PEKRT( GAMMA2 )
         ENDIF

C        --- Compute new amplitudes:
C           A = Integral[ P( z ) * PHI( z ) dz ]
C           (Repeat for each mode PHI to produce excitation coef. A)
C           Integral is done using trapezoidal rule

         SUM = 0.0

         DO IZ = 1, NTOT

C           ------ Calculate PHI at that depth

            PHITMP = 0.0
            ZT = Z( IZ )

            IF ( ZT .GT. DEPTHBR ) THEN
               IF ( BCBOTR .EQ. 'A' ) PHITMP = PHIBR *
     &                          EXP( -GAMBR * ( ZT - DEPTHBR ) )
            ELSE IF ( ZT .LT. DEPTHTR ) THEN
               IF ( BCTOPR .EQ. 'A' ) PHITMP = PHITR *
     &                          EXP( -GAMTR * ( DEPTHTR - ZT ) )
            ELSE
               PHITMP = PHI( IZ )
            ENDIF

            SUM = SUM + P( IZ ) * PHITMP

         END DO

C        *** Compute the contribution from tails in the halfspaces ***

         IF ( BCTOPR .EQ. 'A' ) THEN
            SUM   = SUM + TAIL( Z( 1    ), PHITL, GAMTL, DEPTHTL, ML,
     &                             PHITR / RHOTR, GAMTR, DEPTHTR )
         ENDIF

         IF ( BCBOTR .EQ. 'A' ) THEN
            SUM   = SUM + TAIL( Z( NTOT ), PHIBL, GAMBL, DEPTHBL, ML,
     &                             PHIBR / RHOBR, GAMBR, DEPTHBR )
         ENDIF

         IF ( MODE .GT. ML ) A( MODE ) = 0.0
         WRITE( *, * ) MODE, ABS( A( MODE ) ), ABS( SUM )

         A( MODE ) = SUM
         sum2 = sum2 + abs( sum )**2

C        *** Subtabulate modes at receiver depths ***

         DO IR = 1, NRD

            PHIR( MODE, IR ) = 0.0
            IF      ( RD( IR ) .LT. DEPTHTR ) THEN
C              ------ Rcvr in upper halfspace
               IF ( BCTOPR .EQ. 'A' ) PHIR( MODE, IR ) = PHITR *
     &                     EXP( -GAMTR * ( DEPTHTR   - RD( IR  ) ) )
            ELSE IF ( RD( IR ) .GT. DEPTHBR ) THEN
C              ------ Rcvr in lower halfspace
               IF ( BCBOTR .EQ. 'A' ) PHIR( MODE, IR ) = PHIBR *
     &                     EXP( -GAMBR * ( RD( IR ) - DEPTHBR ) )
            ELSE
               IZ = IRD( IR )
               PHIR( MODE, IR ) = PHI( IZ ) +
     &                     W( IR ) * ( PHI( IZ + 1 ) - PHI( IZ ) )
            ENDIF

         END DO

      END DO   ! next mode

      CLOSE( MODFIL )

      WRITE( *, * ) 'sum2**', sum2

      RETURN
      END
C**********************************************************************C
      SUBROUTINE PLEFT( IPROF, A, CK, Z, PHI, M, P, NR, NTOT,
     &   BCTOP, RHOT, KTOP2, DEPTHT,
     &   BCBOT, RHOB, KBOT2, DEPTHB,
     &   GAMTL, GAMBL, DEPTHTL, DEPTHBL, PHITL, PHIBL, ML, MAXM )

C     Computes the pressure field along the interface.
C     Also returns information needed for the tails in the halfspaces

      PARAMETER ( MAXN = 16001, MAXMED = 50 )

      INTEGER    N( MAXMED )
      REAL       ZL( MAXN ), Z( * ), RHOL( MAXMED ), RHO( MAXMED ),
     &           DEPTHL( MAXMED ), DEPTH( MAXMED )
      COMPLEX    CK( * ), PHI( * ), PL( MAXN ), P( * ), A( * ),
     &           GAMTL( * ), GAMBL( * ), PHITL( * ), PHIBL( * ),
     &           KTOP2, KBOT2
      COMPLEX*16 PEKRT
      CHARACTER  TITLE * 80, MATER( MAXMED ) * 8, BCBOT * 1, BCTOP * 1,
     &           FILNAM * 8

C     *** Read modal info at end of last segment ***

      CALL MODHDR( IPROF - 1, FILNAM,
     &   TITLE, FREQ, NMEDIA, NL, NMAT,
     &   N, MATER, DEPTHL, RHOL,
     &   BCTOP, RHOT, DEPTHTL,
     &   BCBOT, RHOB, DEPTHBL,
     &   ML, MAXM, ZL, CK, KTOP2, KBOT2, MODFIL )

      ML = M   ! we only know M amplitudes

C     *** Compute pressure at the interface ***

      DO IZ = 1, NL
         PL( IZ ) = 0.0
      END DO

      DO MODE = 1, ML
         READ( MODFIL, REC = 6 + MODE ) ( PHI( IZ ), IZ = 1, NL )

         DO IZ = 1, NL
            PL( IZ ) = PL( IZ ) + A( MODE ) * PHI( IZ )
         END DO

C        ------ Halfspace information ------

         PHITL( MODE ) = A( MODE ) * PHI( 1  )
         PHIBL( MODE ) = A( MODE ) * PHI( NL )

         GAMTL( MODE ) = 0.0
         GAMBL( MODE ) = 0.0

C        we take the real part, GAMMA2, since leaky modes grow in depth
         IF ( BCTOP .EQ. 'A' ) THEN
            GAMMA2 = CK( MODE ) ** 2 - KTOP2
            GAMTL( MODE ) = PEKRT( GAMMA2 )
         END IF

         IF ( BCBOT .EQ. 'A' ) THEN
            GAMMA2 = CK( MODE ) ** 2 - KBOT2
            GAMBL( MODE ) = PEKRT( GAMMA2 )
         END IF

      END DO

      CLOSE( MODFIL )

C     *** Read modal data in new segment ***

      CALL MODHDR( IPROF, FILNAM,
     &   TITLE, FREQ, NMEDIA, NR, NMAT,
     &   N, MATER, DEPTH, RHO,
     &   BCTOP, RHOT, DEPTHT,
     &   BCBOT, RHOB, DEPTHB,
     &   M, MAXM, Z, CK, KTOP2, KBOT2, MODFIL )

C     *** Upslope? Extend the Z vector with data from ZL ***

      NTOT = NR

      DO IZL = 1, NL
         IF ( ZL( IZL ) .GT. Z( NTOT ) ) THEN
            NTOT  = NTOT + 1
            Z( NTOT ) = ZL( IZL )
         ENDIF
      END DO

C     *** Retabulate the pressure on the new grid ***

      IZL = 1
      MED = 1
      RHOMED = RHO( 1 )

C     --- Depth of next interface below this one
      IF ( MED .LT. NMEDIA ) THEN
         DBELOW = DEPTH( MED + 1 )
         RHOBEL = RHO( MED + 1 )
      ELSE
         DBELOW = DEPTHB
         RHOBEL = RHOB
      ENDIF

      DO IZ = 1, NTOT
         ZT = Z( IZ )

C        ------ Get medium density
         IF      ( ZT .LT. DEPTHT ) THEN
            RHOMED = RHOT
         ELSE IF ( ZT .GT. DEPTHB ) THEN
            RHOMED = RHOB
         ELSE IF ( MED .LT. NMEDIA ) THEN
            IF   ( ZT .GT. DBELOW ) THEN
               MED    = MED + 1
               RHOMED = RHO( MED )
            ENDIF
         ENDIF

C        --- DEPTH of next interface BELow this one
         IF ( MED .LT. NMEDIA ) THEN
            DBELOW = DEPTH( MED + 1 )
            RHOBEL = RHO( MED + 1 )
         ELSE
            DBELOW = DEPTHB
            RHOBEL = RHOB
         ENDIF

 4000    IF ( ZT .GT. ZL( IZL + 1 ) .AND. IZL .LT. NL - 1 ) THEN
            IZL = IZL + 1
            GOTO 4000
         ENDIF

C        ------ Calculate P at that depth

         IF      ( ZT .GT. DEPTHBL ) THEN   ! lower halfspace
            IF ( BCBOT .EQ. 'A' ) THEN
               P( IZ ) = 0.0
               DO MODE = 1, ML
                  P( IZ ) = P( IZ ) + PHIBL( MODE ) *
     &                 EXP( -GAMBL( MODE ) * ( ZT - DEPTHBL ) )
               END DO
            ENDIF
         ELSE IF ( ZT .LT. DEPTHTL ) THEN   ! upper halfpace
            IF ( BCTOP .EQ. 'A' ) THEN
               P( IZ ) = 0.0
               DO MODE = 1, ML
                  P( IZ ) = P( IZ ) + PHITL( MODE ) *
     &                 EXP( -GAMTL( MODE ) * ( DEPTHTL - ZT ) )
               END DO
            ENDIF
         ELSE
            P( IZ ) = PL( IZL ) +
     &         ( ZT - ZL( IZL ) ) / ( ZL( IZL + 1 ) - ZL( IZL ) ) *
     &                              ( PL( IZL + 1 ) - PL( IZL ) )
         ENDIF

C        *** compute mesh width, h ***

         IF ( IZ .EQ. 1 ) THEN
C           ------ First point
            H = 0.5 * ( Z(  2   ) - Z(   1      ) ) / RHOMED
         ELSE IF ( IZ .EQ. NTOT ) THEN
C           ------ Last point
            H = 0.5 * ( Z( NTOT ) - Z( NTOT - 1 ) ) / RHOMED
         ELSE
            IF ( Z( IZ - 1 ) .LT. DBELOW .AND.
     &           Z( IZ + 1 ) .GT. DBELOW ) THEN
C              ------ Point just above or below the interface
               H = 0.5 * ( Z( IZ + 1 ) / RHOBEL
     &                   - Z( IZ - 1 ) / RHOMED
     &                   - DBELOW / RHOBEL
     &                   + DBELOW / RHOMED )
            ELSE
               H = 0.5 * ( Z( IZ + 1 ) - Z( IZ   - 1 ) ) / RHOMED
            ENDIF
         ENDIF

         P( IZ ) = H * P( IZ )

      END DO   ! next depth IZ

      RETURN
      END
C**********************************************************************C
      COMPLEX FUNCTION TAIL( D, PHIL, GAML, DL, ML, PHIR, GAMR, DR )

      COMPLEX GAML( * ), PHIL( * ), GAMR, PHIR, FR

      TAIL = 0.0
      FR =  PHIR         * EXP( -GAMR         * ( D - DR ) )

      DO MODE = 1, ML

         IF ( D .EQ. DL ) THEN
            TAIL = TAIL + FR * PHIL( MODE ) / ( GAML( MODE ) + GAMR )
         ELSE
            TAIL = TAIL + FR * 
     &         PHIL( MODE ) * EXP( -GAML( MODE ) * ( D - DL ) )
     &               / ( GAML( MODE ) + GAMR )
         ENDIF

      END DO

      RETURN
      END
