      PROGRAM KRAKEN

C     Program for solving for ocean acoustic normal modes

C     Michael B. Porter

      INCLUDE 'common.f'

      REAL*4 ZMIN, ZMAX

      NV( 1 ) = 1
      NV( 2 ) = 2
      NV( 3 ) = 4
      NV( 4 ) = 8
      NV( 5 ) = 16

C     *** Read in environmental info ***

      TITLE = 'KRAKEN- '

      CALL READIN( TITLE, FREQ, MAXMED, NMEDIA,
     &   TOPOPT, CPT, CST, RHOT, BUMDEN, ETA, XI, NG, SIGMA, DEPTH, 
     &   BOTOPT, CPB, CSB, RHOB,
     &   ENVFIL, PRTFIL )

C     *** Spectral limits ***

      READ(  ENVFIL, *    ) CLOW, CHIGH
      WRITE( PRTFIL, 500  ) CLOW, CHIGH
  500 FORMAT( /, ' CLOW = ', G12.5, '  CHIGH = ', G12.5 )

C     *** Maximum range for calculations ***

      READ(  ENVFIL, * ) RMAX
      WRITE( PRTFIL, * ) 'RMAX = ', RMAX

C     *** Read source/receiver depths ***

      ZMIN = DEPTH( 1 )
      ZMAX = DEPTH( NMEDIA + 1 )

      CALL SDRD( ENVFIL, PRTFIL, ZMIN, ZMAX,
     &   MAXNSD, NSD, SD,
     &   MAXNRD, NRD, RD )

      OMEGA2 = ( 2.0 * PI * FREQ ) ** 2
      CLOSE( ENVFIL )

      M = MAXM

C     *** Main loop: solve the problem for a sequence of meshes ***

      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * ) 'Mesh multiplier   CPU seconds'

      DO ISET = 1, NSETS

         DO MED = 1, NMEDIA
            N( MED ) = NG( MED ) * NV( ISET )
            H( MED ) = ( DEPTH( MED + 1 ) - DEPTH( MED ) ) / N( MED )
         END DO

         HV( ISET ) = H( 1 )
         CALL SOLVE( ERROR )
         IF ( RMAX .EQ. 0.0 ) GOTO 3000
         IF ( ERROR .LT. 1.0 / ( 1000.0 * RMAX ) ) GOTO 3000
      END DO

C     ------ Fall through indicates failure to converge
      CALL ERROUT( PRTFIL, 'W', 'KRAKEN',
     &             'Too many meshes needed: check convergence' )

C     Solution complete: discard modes with phase velocity above CHIGH

 3000 OMEGA = SQRT( OMEGA2 )

      DO MODE = M, 1, -1
         IF ( EXTRAP( 1, MODE ) .GT. OMEGA2 / CHIGH ** 2 ) GOTO 5000
      END DO

 5000 M = MODE

C     *** Write eigenvalues to PRTFIL and MODFIL ***

      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * )
     &   '   I           K                ALPHA         PHASE SPEED'

      DO MODE = 1, M
         CK( MODE ) = SQRT( EXTRAP( 1, MODE ) + CK( MODE ) )
         WRITE( PRTFIL, 6000 ) MODE,
     &                         CK( MODE ), OMEGA / DBLE( CK( MODE ) )
 6000    FORMAT( I5, 3G19.10 )
      END DO

      WRITE( MODFIL, REC = 5 ) M, LRECL

      IFIRST = 1
      DO IREC = 1, 1 + ( 2 * M - 1 ) / LRECL

         ILAST  = MIN( M, IFIRST + LRECL / 2 - 1 )
         WRITE( MODFIL, REC = 6 + M + IREC )
     &      ( CMPLX( CK( MODE ) ), MODE = IFIRST, ILAST )

         IFIRST = ILAST + 1
      END DO

      STOP
      END
C**********************************************************************C
      SUBROUTINE INIT

C     Initializes arrays defining difference equations

      INCLUDE 'common.f'

      LOGICAL ELFLAG
      REAL*8 CP2, CS2
      COMPLEX*16 CP( MAXN ), CS( MAXN )

      COMMON /SPACE/ CP, CS

      ELFLAG   = .FALSE.
      CMIN     = 1.0E6
      NFACT    = 0
      LOC( 1 ) = 0

C     *** Loop over media ***

      DO MED = 1, NMEDIA
         IF ( MED .NE. 1 ) LOC( MED ) = LOC( MED-1 ) + N( MED-1 ) + 1
         N1 = N(   MED ) + 1
         I  = LOC( MED ) + 1

         IF ( LOC( MED ) + N1 .GT. MAXN )
     &      CALL ERROUT( PRTFIL, 'F', 'KRAKEN',
     &             'Insufficient storage for mesh: '//
     &             'Reduce mesh size or change program dimensions' )

C        *** PROFIL reads in the data for a medium ***

         CALL PROFIL( DEPTH, CP( I ), CS( I ), RHO( I ),
     &      MED, N1, FREQ, TOPOPT(1:1), TOPOPT(3:3), 'TAB',
     &      ENVFIL, PRTFIL )

C        *** Load diagonals of the finite-difference equations ***

         IF ( REAL( CS( I ) ) .EQ. 0.0 ) THEN

C           --- Case of an acoustic medium ---

            MATER( MED ) = 'ACOUSTIC'
            IF ( NFACT .EQ. 0 ) NFACT = MED
            NLACT = MED

            DO J = I, I + N( MED )
               CMIN = MIN( DBLE( CP( J ) ), CMIN )
               B1(  J ) = -2.0 + H( MED ) ** 2
     &                   * DBLE( OMEGA2 / CP( J ) ** 2 )
               B1C( J ) = DIMAG( OMEGA2 / CP( J ) ** 2 )
C              --- use DIMAG in f77, AIMAG in f90
            END DO

         ELSE

C           --- Case of an elastic medium ---

            IF ( SIGMA( MED ) .NE. 0.0 )
     &         CALL ERROUT( PRTFIL, 'F', 'KRAKEN',
     &               'Rough elastic interfaces are not allowed' )

            MATER( MED ) = 'ELASTIC'
            ELFLAG = .TRUE.
            TWOH   = 2.0 * H( MED )

            DO J = I, I + N( MED )
               CMIN = MIN( DBLE( CS( J ) ), CMIN )

               CP2 = CP( J ) ** 2
               CS2 = CS( J ) ** 2

               B1( J ) = TWOH / ( RHO( J ) * CS2 )
               B2( J ) = TWOH / ( RHO( J ) * CP2 )
               B3( J ) = 4.0 * TWOH * RHO(J) * CS2 * ( CP2 - CS2 ) / CP2
               B4( J ) = TWOH * ( CP2 - 2.0*CS2 ) / CP2

               RHO( J ) = TWOH * OMEGA2 * RHO( J )
            END DO

         ENDIF
C        ------ Next MED:
      END DO

C     (CLOW, CHIGH) = phase speed interval for the mode search
C     user specified interval is reduced if it exceeds domain
C     of possible mode phase speeds

C     *** Bottom properties ***

      IF ( BOTOPT(1:1) .EQ. 'A' ) THEN
         IF ( REAL( CSB ) .GT. 0.0 ) THEN
C           --- Elastic bottom:
            ELFLAG = .TRUE.
            CMIN  = MIN( CMIN,  DBLE( CSB ) )
            CHIGH = MIN( CHIGH, DBLE( CSB ) )
         ELSE
C           --- Acoustic bottom:
            CMIN  = MIN( CMIN,  DBLE( CPB ) )
            CHIGH = MIN( CHIGH, DBLE( CPB ) )
         ENDIF
      ENDIF

C     *** Top properties ***

      IF ( TOPOPT(2:2) .EQ. 'A' ) THEN
         IF ( REAL( CST ) .GT. 0.0 ) THEN
C           --- Elastic top half-space:
            ELFLAG = .TRUE.
            CMIN =  MIN( CMIN,  DBLE( CST ) )
            CHIGH = MIN( CHIGH, DBLE( CST ) )
         ELSE
C           --- Acoustic top half-space:
            CMIN  = MIN( CMIN,  DBLE( CPT ) )
            CHIGH = MIN( CHIGH, DBLE( CPT ) )
         ENDIF
      ENDIF

C     --- If elastic medium then reduce CMIN for Scholte wave
      IF ( ELFLAG ) CMIN = 0.85 * CMIN
      CLOW = MAX( CLOW, CMIN )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE SOLVE( ERROR )

C     Solves the eigenvalue problem at the current mesh
C     and produces a new extrapolation

      INCLUDE 'common.f'

      CALL CLTIME
      CALL INIT

C     Choose a solver ...
C     Use bisection if this is the 1st or 2nd set and
C     there are no internal elastic media.

      IF ( ( ISET .LE. 2 ) .AND. ( NMEDIA .LE. NLACT-NFACT+1 ) ) THEN
         CALL SOLVE1
      ELSE
         CALL SOLVE2
      ENDIF

      DO MODE = 1, M
         EXTRAP( ISET, MODE ) = EVMAT( ISET, MODE )
      END DO

C     *** If this is the first mesh, compute the eigenvectors ***

      IF ( ISET .EQ. 1 ) CALL VECTOR

C     *** Now do the Richardson extrapolation ***

      ERROR = 10            ! initialize error to a large number
      KEY = 2 * M / 3 + 1   ! KEY value used to check convergence

      IF ( ISET .GT. 1 ) THEN
         T1 = EXTRAP( 1, KEY )

         DO J = ISET - 1, 1, -1
            DO MODE = 1, M
               X1 = NV( J    ) ** 2
               X2 = NV( ISET ) ** 2
               F1 = EXTRAP( J,     MODE )
               F2 = EXTRAP( J + 1, MODE )
               EXTRAP( J, MODE ) = F2 - ( F1 - F2 ) / ( X2 / X1 - 1.0 )
            END DO
         END DO

         T2 = EXTRAP( 1, KEY )
         ERROR = ABS( T2 - T1 )
      ENDIF

      CALL RDTIME( TEND )
      ET( ISET ) = TEND
      WRITE( PRTFIL, 5000 ) NV( ISET ), ET( ISET )
 5000 FORMAT( 1X, I8, 6X, G15.3 )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE SOLVE2

C     Provides initial guess to root finder for each EVMAT(I)

      INCLUDE 'common.f'

      REAL*8 P( 10 )
      CHARACTER*80 ERRMSG

      X = OMEGA2 / CLOW ** 2
      MAXIT = 500

      DO MODE = 1, M

C        For first or second meshes, use a high guess
C        Otherwise use extrapolation to produce an initial guess

         X = 1.00001 * X
 
         IF ( ISET .GE. 2 ) THEN

            DO I = 1, ISET - 1
               P( I ) = EVMAT( I, MODE )
            END DO

            IF ( ISET .GE. 3 ) THEN
               DO I = 1, ISET - 2
                  DO J = 1, ISET - I - 1
                     X1 = HV( J     ) ** 2
                     X2 = HV( J + I ) ** 2

                     P( J ) = ( ( HV( ISET ) ** 2 - X2 ) * P( J     ) -
     &                          ( HV( ISET ) ** 2 - X1 ) * P( J + 1 ) )
     &                                  / ( X1 - X2 )
                  END DO
               END DO
               X = P( 1 )
            ENDIF

         ENDIF

C        *** Use the secant method to refine the eigenvalue ***

         TOL = 100.0 * MACHEP * ABS( X )
         CALL ZSECX( X, TOL, IT, MAXIT, ERRMSG )

         IF ( ERRMSG .NE. ' ' ) THEN
            WRITE( PRTFIL, * ) 'ISET, MODE = ', ISET, MODE
            CALL ERROUT( PRTFIL, 'W', 'KRAKEN-ZSECX', ERRMSG )
            X = 1.0E-9   ! make sure value is discarded
         ENDIF

         EVMAT( ISET, MODE ) = X

C        --- Toss out modes outside user specified spectrum
         IF ( OMEGA2 / X .GT. CHIGH ** 2 ) THEN
            M = MODE - 1
            RETURN
         ENDIF

      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE FUNCT( X, DELTA, IPOW )

C     FUNCT( X ) = 0 is the dispersion relation

      INCLUDE 'common.f'

      PARAMETER ( ROOF = 1.0E5, FLOOR = 1.0E-5, IPOWR = 5, IPOWF = -5 )
      CHARACTER BCTYPE*1 

C     For a k below the cts spectrum limit, force a zero

      IF ( X .LE. OMEGA2 / CHIGH ** 2 ) THEN
         DELTA = 0.0
         IPOW  = 0
         RETURN
      ENDIF

      MODECT = 0

C     --- Bottom impedance
      BCTYPE(1:1) = BOTOPT(1:1)
      CALL BCIMP( X, BCTYPE, 'BOT', CPB, CSB, RHOB, F, G, IPOW )

C     --- Shoot through acoustic layers
      CALL ACOUST( X, F,  G,  IPOW  )

C     --- Top impedance
      BCTYPE(1:1) = TOPOPT(2:2)
      CALL BCIMP( X, BCTYPE, 'TOP', CPT, CST, RHOT, F1, G1, IPOW1 )

      DELTA = F * G1 - G * F1
      IF ( G * DELTA .GT. 0.0 ) MODECT = MODECT + 1
      IPOW = IPOW + IPOW1

C     *** Deflate previous roots ***

      IF ( ( MODE .GT. 1 ) .AND. ( NMEDIA .GT. NLACT - NFACT+1 ) ) THEN
         DELTA = DELTA / ( X - EVMAT( ISET, MODE-1 ) )

         IF ( MODE .GT. 2 ) THEN
            DO J = 1, MODE - 2
               DELTA = DELTA / ( X - EVMAT(ISET, J) )

C              ------ Scale if necessary
   50          IF ( ABS( DELTA ) .LT. FLOOR .AND.
     &                     ABS( DELTA )   .GT. 0.0) THEN
                  DELTA = ROOF * DELTA
                  IPOW  = IPOW - IPOWR
                  GOTO 50
               ENDIF

   60          IF ( ABS( DELTA ) .GE. ROOF ) THEN
                  DELTA = FLOOR * DELTA
                  IPOW  = IPOW - IPOWF
                  GOTO 60
               ENDIF

            END DO
         ENDIF
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE ACOUST( X, F, G, IPOW )

C     Shoot through acoustic layers

      INCLUDE 'common.f'

      PARAMETER ( ROOF = 1.0E5, FLOOR = 1.0E-5, IPOWR = 5, IPOWF = -5 )

      IF ( NFACT .EQ. 0 ) RETURN

C     *** Loop over successive acoustic media ***

      DO MED = NLACT, NFACT, -1
         H2K2 = H( MED ) ** 2 * X
         I = LOC( MED ) + N( MED ) + 1
         RHOM = RHO( I )
         P1 = -2.0 * G
         P2 = ( B1( I ) - H2K2 ) * G - 2.0 * H( MED ) * F * RHOM

C        *** Shoot through a single medium ***

         DO I = LOC( MED ) + N( MED ), LOC( MED ) + 1, -1

            P0 = P1
            P1 = P2
            P2 = ( H2K2 - B1( I ) ) * P1 - P0

            IF( P0 * P1 .LE. 0.0 ) MODECT = MODECT + 1
C            IF ( ( P0 .GT. 0.0 ) .AND. (P1 .LT. 0.0) ) MODECT = MODECT+1
C            IF ( ( P0 .LT. 0.0 ) .AND. (P1 .GT. 0.0) ) MODECT = MODECT+1
C            IF (   P1 .EQ. 0.0 )                       MODECT = MODECT+1

C           ------ Scale if necessary
   60       IF ( ABS( P2 ) .GE. ROOF ) THEN
               P0 = FLOOR * P0
               P1 = FLOOR * P1
               P2 = FLOOR * P2
               IPOW = IPOW - IPOWF
               GOTO 60
            ENDIF

         END DO

C        ------ F = P'/RHO and G = -P since FP+GP'/RHO = 0
         RHOM = RHO( LOC( MED ) + 1 )
         F = -( P2 - P0 ) / ( 2.0 * H( MED ) ) / RHOM
         G = -P1

      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE VECTOR

C     Do inverse iteration to compute each of the eigenvectors
C     and write these to the disk file

      INCLUDE 'common.f'

      INTEGER IZTAB( MAXNZT )
      REAL Z( MAXN ), ZTAB( MAXNZT ), PHITAB( MAXNZT ), WTS( MAXNZT )
      REAL*8 PHI( MAXN ), D( MAXN ), E( MAXN ), W
      CHARACTER BCTOP*1, BCBOT*1

      COMMON /SPACE/ RV1( MAXN ), RV2( MAXN ), RV3( MAXN ), RV4( MAXN )

      EQUIVALENCE ( Z, D, PHI )

      BCTOP(1:1) = TOPOPT(2:2)
      BCBOT(1:1) = BOTOPT(1:1)

C     *** Tabulate z-coordinates and off-diagonals of matrix ***

      NTOT = 0
      J    = 1
      Z( 1 ) = DEPTH( NFACT )

      DO MED = NFACT, NLACT

         HRHO = H( MED ) * RHO( LOC( MED ) + 1 )

         DO I = 1, N( MED )
            J = J + 1
            Z( J ) = Z( J - 1 ) + H( MED )
            E( J ) = 1.0 / HRHO
         END DO

         NTOT = NTOT + N( MED )
      END DO

      NTOT1 = NTOT + 1

C     Calculate the indices, weights, ... for mode interpolation

      CALL MERGEV( SD, NSD, RD, NRD, ZTAB, NZTAB )
      CALL WEIGHT( Z, NTOT1, ZTAB, NZTAB, WTS, IZTAB )

C     Open MODFIL and write header

C     ------ Logical record length in `longwords' (4 bytes)
      LRECL = MAX( 2 * NZTAB, 32, 3 * ( NLACT - NFACT + 1 ) )

      OPEN ( FILE = 'MODFIL', UNIT = MODFIL, STATUS = 'NEW',
     &   ACCESS = 'DIRECT', RECL = 4 * LRECL, FORM = 'UNFORMATTED' )

      WRITE( MODFIL, REC = 1 ) LRECL, TITLE, REAL( FREQ ),
     &   NLACT - NFACT + 1, NZTAB, NZTAB
      WRITE( MODFIL, REC = 2 )
     &   ( N( MED ), MATER( MED ), MED = NFACT, NLACT )
      WRITE( MODFIL, REC = 3 )
     &   BCTOP(1:1), CMPLX( CPT ), CMPLX( CST ),
     &     REAL( RHOT ), REAL( DEPTH( 1          ) ),
     &   BCBOT(1:1), CMPLX( CPB ), CMPLX( CSB ),
     &     REAL( RHOB ), REAL( DEPTH( NMEDIA + 1 ) )
      WRITE( MODFIL, REC = 4 ) ( REAL( DEPTH( MED ) ),
     &   REAL( RHO( LOC( MED ) + 1 ) ), MED = NFACT, NLACT )
      WRITE( MODFIL, REC = 6 ) ( ZTAB( IZ ), IZ = 1, NZTAB )

C     *** Main loop: for each eigenvalue call SINVIT to get eigenvector

      DO MODE = 1, M
         W = 0.0
         X = EVMAT( 1, MODE )

C        *** Corner elt requires top impedance ***

         CALL BCIMP( X, BCTOP, 'TOP', CPT, CST, RHOT, F, G, IPOW )

         IF ( G .EQ. 0.0 ) THEN
            D( 1 ) = 1.0
            E( 2 ) = 0.0
         ELSE
            L = LOC( NFACT ) + 1
            XH2    = X * H( NFACT ) * H( NFACT )
            HRHO   = H( NFACT ) * RHO( L )
            D( 1 ) = ( B1( L ) - XH2 ) / HRHO / 2.0 + F / G
         ENDIF

C        *** Set up the diagonal ***

         ITP = NTOT
         J   = 1
         L   = LOC( NFACT ) + 1

         DO MED = NFACT, NLACT
            XH2  = X * H( MED ) ** 2
            HRHO = H( MED ) * RHO( LOC( MED ) + 1 )

            IF ( MED .GE. NFACT + 1 ) THEN
               L = L + 1
               D( J ) = ( D( J ) + ( B1( L ) - XH2 ) / HRHO ) / 2.0
            ENDIF

            DO I = 1, N( MED )
               J = J + 1
               L = L + 1
               D( J ) = ( B1( L ) - XH2 ) / HRHO
C              ------ Locate index of turning point nearest top
               IF ( B1( L ) - XH2 + 2.0 .GT. 0.0 ) THEN
                  ITP = MIN( J, ITP )
               ENDIF
            END DO

         END DO

C        *** Corner elt requires bottom impedance ***

         CALL BCIMP( X, BCBOT, 'BOT', CPB, CSB, RHOB, F, G, IPOW )

         IF ( G .EQ. 0.0 ) THEN
            D( NTOT1 ) = 1.0
            E( NTOT1 ) = 0.0
         ELSE
            D( NTOT1 ) =  D( NTOT1 ) / 2.0 - F / G
         ENDIF

C        *** Use inverse iteration to compute eigenvector ***

         CALL SINVIT( NTOT1, D, E, W, PHI, IERR,
     &      RV1, RV2, RV3, RV4, PHI )

         IF ( IERR .NE. 0 ) THEN
            WRITE( PRTFIL, * ) 'MODE = ', MODE
            CALL ERROUT( PRTFIL, 'W', 'KRAKEN-SINVIT',
     &                   'Inverse iteration failed to converge' )
         ENDIF

C        *** Normalize the eigenvector ***

         CALL NORMIZ( PHI, ITP, NTOT1, X )

C        Tabulate the modes at the source depths and write to disk

         DO IZ = 1, NZTAB
            I = IZTAB( IZ )
            PHITAB( IZ ) = PHI( I ) +
     &         WTS(IZ) * ( PHI( I + 1 ) - PHI( I ) )
         END DO

         WRITE( MODFIL, REC = 6 + MODE )
     &      ( PHITAB( IZ ), 0.0, IZ = 1, NZTAB )

      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE NORMIZ( PHI, ITP, NTOT1, X )

C     Normalize the eigenvector:
C        SQNRM = Integral(PHI ** 2) by the trapezoidal rule:
C        Integral(F) = H*( F(1)+...+F(N-1) + 0.5*(F(0)+F(N)) )

C     Compute perturbation due to material absorption
C     Call SCAT to figure interfacial scatter loss

      INCLUDE 'common.f'

      REAL*8 PHI( * )
      COMPLEX*16 PERK, DEL
      CHARACTER BCTYPE*1

      SQNRM = 0.0
      PERK  = 0.0

C     *** Compute perturbation due to loss in top half-space ***

      IF ( TOPOPT(2:2) .EQ. 'A' ) THEN
         DEL = -0.5*( OMEGA2 / CPT ** 2 - DBLE( OMEGA2 / CPT ** 2 ) ) /
     &          SQRT( X                 - DBLE( OMEGA2 / CPT ** 2 ) )
         PERK = PERK - DEL * PHI( 1 ) ** 2 / RHOT
      ENDIF

C     *** Compute norm and pertubation due to material absorption ***

      L = LOC( NFACT )
      J = 1

      DO MED = NFACT, NLACT
         L = L + 1
         RHOM = RHO( L )

         SQNRM = SQNRM + 0.5*H( MED ) *             PHI( J ) ** 2 / RHOM
         PERK  = PERK  + 0.5*H( MED ) * CI*B1C(L) * PHI( J ) ** 2 / RHOM

         DO I = 1, N( MED ) - 1
            L = L + 1
            J = J + 1

            SQNRM = SQNRM +  H( MED ) *             PHI( J ) ** 2 / RHOM
            PERK  = PERK  +  H( MED ) * CI*B1C(L) * PHI( J ) ** 2 / RHOM

         END DO

         L = L + 1
         J = J + 1

         SQNRM = SQNRM + 0.5*H( MED ) *             PHI( J ) ** 2 / RHOM
         PERK  = PERK  + 0.5*H( MED ) * CI*B1C(L) * PHI( J ) ** 2 / RHOM

      END DO

C     *** Compute perturbation due to loss in bottom half-space ***

      IF ( BOTOPT(1:1) .EQ. 'A' ) THEN
         DEL = -0.5*( OMEGA2 / CPB ** 2 - DBLE( OMEGA2 / CPB ** 2 ) ) /
     &          SQRT( X                 - DBLE( OMEGA2 / CPB ** 2 ) )
         PERK = PERK - DEL * PHI( J ) ** 2 / RHOB
      ENDIF

C     *** Compute derivative of top admitance ***

      X1 = 0.9999999D0 * X
      X2 = 1.0000001D0 * X

      BCTYPE(1:1) = TOPOPT(2:2)

      CALL BCIMP( X1, BCTYPE, 'TOP', CPT, CST, RHOT, F1, G1, IPOW )
      CALL BCIMP( X2, BCTYPE, 'TOP', CPT, CST, RHOT, F2, G2, IPOW )

      DRHODX = 0.0
      IF ( G1 .NE. 0.0 ) THEN
         DRHODX = -( F2 / G2 - F1 / G1 ) / ( X2 - X1 )
      ENDIF

C     *** Compute derivative of bottom admitance ***

      BCTYPE(1:1) = BOTOPT(1:1)

      CALL BCIMP( X1, BCTYPE, 'BOT', CPB, CSB, RHOB, F1, G1, IPOW )
      CALL BCIMP( X2, BCTYPE, 'BOT', CPB, CSB, RHOB, F2, G2, IPOW )

      DETADX = 0.0
      IF ( G1 .NE. 0.0 ) THEN
         DETADX = -( F2 / G2 - F1 / G1 ) / ( X2 - X1 )
      ENDIF

C     *** Scale the mode ***

      RN = SQNRM + DRHODX * PHI( 1 ) ** 2 - DETADX * PHI( NTOT1 ) ** 2

      IF ( RN .LE. 0.0 ) THEN
         RN = -RN
         WRITE( PRTFIL, * ) 'MODE = ', MODE
         CALL ERROUT( PRTFIL, 'W', 'KRAKEN',
     &      'Normalization constant non-positive' )
      ENDIF

      SCALE = 1.0 / SQRT( RN )
      IF ( PHI( ITP ) .LT. 0.0 ) SCALE = -SCALE

      DO J = 1, NTOT1
         PHI( J ) = SCALE * PHI( J )
      END DO

      PERK = SCALE ** 2 * PERK

C     *** Compute interfacial scatter loss ***

      CALL SCAT( PERK, PHI, X )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE SCAT( PERK, PHI, X )

C     Figure scatter loss

      INCLUDE 'common.f'

      REAL*8 PHI( * )
      COMPLEX*16 CIMPED, KX, TWERSK, PERK,
     &   ETA1SQ, ETA2SQ, KUPING, U, PHIC
      CHARACTER BCTYPE*1

      OMEGA = SQRT( OMEGA2 )
      KX    = SQRT( X )

C     *** Top Twersky roughness ***

      BCTYPE(1:1) = TOPOPT(2:2)
      IF ( BCTYPE(1:1) .EQ. 'S' .OR. BCTYPE(1:1) .EQ. 'H' .OR.
     &     BCTYPE(1:1) .EQ. 'T' .OR. BCTYPE(1:1) .EQ. 'I' ) THEN

         I = LOC( NFACT ) + N( NFACT ) + 1
         RHOINS = RHO( I )
         CINS = SQRT( OMEGA2 * H( NFACT ) ** 2 / ( 2.0 + B1( NFACT ) ) )

         CIMPED = TWERSK( BCTYPE(1:1), OMEGA,
     &                       BUMDEN, XI, ETA, KX, RHOINS, CINS )

         CIMPED = CIMPED / ( -CI * OMEGA * RHOINS )
         DPHIDZ = PHI( 2 ) / H( NFACT )
         PERK   = PERK - CIMPED * DPHIDZ ** 2
      ENDIF

C     *** Bottom Twersky roughness ***

      BCTYPE(1:1) = BOTOPT(1:1)
      IF ( BCTYPE(1:1) .EQ. 'S' .OR. BCTYPE(1:1) .EQ. 'H' .OR.
     &     BCTYPE(1:1) .EQ. 'T' .OR. BCTYPE(1:1) .EQ. 'I' ) THEN

         I = LOC( NLACT ) + N( NLACT ) + 1
         RHOINS = RHO( I )
         CINS   = SQRT( OMEGA2 * H(NLACT) ** 2 / ( 2.0 + B1(NLACT) ) )

         CIMPED = TWERSK( BCTYPE(1:1), OMEGA,
     &                       BUMDEN, XI, ETA, KX, RHOINS, CINS )

         CIMPED = CIMPED / ( -CI * OMEGA * RHOINS )
         DPHIDZ = PHI( 2 ) / H( NFACT )
         PERK   = PERK - CIMPED * DPHIDZ ** 2
      ENDIF

      J = 1
      L = LOC( NFACT )

C     *** Loop over media ***

      DO MED = NFACT - 1, NLACT

C        *** Calculate RHO1, ETA1SQ, PHI, U ***

         IF ( MED .EQ. NFACT - 1 ) THEN
C           *** Top properties ***

            BCTYPE(1:1) = TOPOPT(2:2)

            IF ( BCTYPE(1:1) .EQ. 'A' ) THEN
C              ------ Acousto-elastic
               RHO1   = RHOT
               ETA1SQ = X - OMEGA2 / CPT ** 2
               U = SQRT( ETA1SQ ) * PHI( 1 ) / RHOT

            ELSE IF ( BCTYPE(1:1) .EQ. 'V' ) THEN
C              ------ Vacuum
               RHO1   = 1.0E-9
               ETA1SQ = 1.0
               RHOINS = RHO( LOC( NFACT ) + 1 )
               U      = PHI( 2 ) / H( NFACT ) / RHOINS

            ELSE IF ( BCTYPE(1:1) .EQ. 'R' ) THEN
C              ------ Rigid
               RHO1   = 1.0E+9
               ETA1SQ = 1.0
               U      = 0.0
            ENDIF
         ELSE
            H2 = H( MED ) ** 2
            J  = J + N( MED )
            L  = LOC( MED ) + N( MED ) + 1

            RHO1   = RHO( L )
            ETA1SQ = ( 2.0 + B1( L ) ) / H2 - X
            U = ( -PHI( J-1 ) - 0.5 * ( B1(L) - H2*X ) * PHI( J ) )
     &               / ( H( MED ) * RHO1 )
         ENDIF

C        *** Calculate RHO2, ETA2 ***

         IF ( MED .EQ. NLACT ) THEN
C           *** Bottom properties ***

            BCTYPE(1:1) = BOTOPT(1:1)

            IF ( BCTYPE(1:1) .EQ. 'A' ) THEN
C              ------ Acousto-elastic
               RHO2   = RHOB
               ETA2SQ = OMEGA2 / CPB ** 2 - X

            ELSE IF ( BCTYPE(1:1) .EQ. 'V' ) THEN
C              ------ Vacuum
               RHO2   = 1.0E-9
               ETA2SQ = 1.0

            ELSE IF ( BCTYPE(1:1) .EQ. 'R' ) THEN
C              ------ Rigid
               RHO2   = 1.0E+9
               ETA2SQ = 1.0
            ENDIF
         ELSE
            RHO2   = RHO( L + 1 )
            ETA2SQ = ( 2.0 + B1( L + 1 ) ) / H( MED + 1 ) ** 2 - X
         ENDIF

         PHIC = PHI( J )   ! convert to complex*16
         PERK = PERK + KUPING( SIGMA( MED + 1 ), ETA1SQ, RHO1,
     &                         ETA2SQ, RHO2, PHIC, U )

      END DO

      CK( MODE ) = PERK

      RETURN
      END
C**********************************************************************C
      SUBROUTINE SOLVE1

C     Uses Sturm sequences to isolate the eigenvalues
C     and Brent's method to refine them

      INCLUDE 'common.f'

      REAL*8 XL( MAXM + 1 ), XR( MAXM + 1 )
      CHARACTER*80 ERRMSG

C     *** Determine number of modes ***

      XMIN = 1.00001D0 * OMEGA2 / CHIGH ** 2

      CALL FUNCT( XMIN, DELTA, IPOW )
      M = MODECT

      XMAX = OMEGA2 / CLOW ** 2
      CALL FUNCT( XMAX, DELTA, IPOW )
      M = M - MODECT

      IF ( M .EQ. 0 ) THEN

C        Open a dummy MODFIL for possible use by FIELD3D

         LRECL = 32
         NZTAB = 0
         OPEN ( FILE = 'MODFIL', UNIT = MODFIL, STATUS = 'NEW',
     &          ACCESS = 'DIRECT', RECL = 4 * LRECL,
     &          FORM = 'UNFORMATTED' )
         WRITE( MODFIL, REC = 1 )
     &          LRECL, TITLE, REAL( FREQ ), 1, NZTAB, NZTAB
         WRITE( MODFIL, REC = 5 ) M, LRECL

         CALL ERROUT( PRTFIL, 'F', 'KRAKEN',
     &                'No modes for given phase speed interval' )

      ENDIF

      IF( M .GT. MAXM ) THEN
         WRITE( PRTFIL, * ) 'Number of modes = ', M
         WRITE( PRTFIL, * ) 'Number of modes allowable = ', MAXM
         CALL ERROUT( PRTFIL, 'F', 'KRAKEN',
     &     'Too many modes'//
     &     'Reduce phase speed interval or increase dimension MAXM' )
      ENDIF

      NTOT = 0
      DO MED = NFACT, NLACT
         NTOT = NTOT + N( MED )
      END DO

      IF ( M .GT. NTOT / 3 ) THEN
         WRITE( PRTFIL, * ) 'Approximate number of modes =', M
         CALL ERROUT( PRTFIL, 'W', 'KRAKEN',
     &     'Mesh too coarse to sample the modes adequately' )
      ENDIF

C     *** Initialize upper and lower bounds ***

      CALL BISECT( XMIN, XMAX, XL, XR )

C     *** Call ZBRENT to refine each eigenvalue in turn ***

      DO MODE = 1, M
         X1  = XL( MODE )
         X2  = XR( MODE )
         EPS = 100.0 * MACHEP * X2

         CALL ZBRENTX( X, X1, X2, EPS, ERRMSG )

         IF ( ERRMSG .NE. ' ' ) THEN
            WRITE( PRTFIL, * ) 'ISET, MODE = ', ISET, MODE
            CALL ERROUT( PRTFIL, 'W', 'KRAKEN-ZBRENTX', ERRMSG )
         ENDIF

         EVMAT( ISET, MODE ) = X
      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE BISECT( XMIN, XMAX, XL, XR )

C     Returns an isolating interval (XL, XR) for each eigenvalue

      PARAMETER( MXBIS = 50 )
      INCLUDE 'common.f'

      REAL*8 XL( * ), XR( * )

      DO I = 1, M + 1
         XL( I ) = XMIN
         XR( I ) = XMAX
      END DO

      CALL FUNCT( XMAX, DELTA, IPOW )
      NZER1 = MODECT

C     *** For each eigenvalue, do: ***

      IF ( M .EQ. 1 ) RETURN

      DO 5000 MODE = 1, M - 1

C        Obtain initial guesses for X1 and X2

         IF ( XL( MODE ) .EQ. XMIN ) THEN
            X2 = XR( MODE )
            DO J = MODE + 1, M
               IF ( XL( J ) .GT. XMIN ) GOTO 3000
            END DO
 3000       X1 = XL( J )

C           Begin bisection (assuming we will never
C           need more than MXBIS bisections per mode)

            DO J = 1, MXBIS
               X = X1 + ( X2 - X1 ) / 2
               CALL FUNCT( X, DELTA, IPOW )
               NZEROS = MODECT - NZER1

               IF ( NZEROS .LT. MODE ) THEN
                  X2 = X
                  XR( MODE ) = X
               ELSE
                  X1 = X
                  IF ( XR( NZEROS + 1 ) .GE. X ) XR( NZEROS + 1 ) = X
                  IF ( XL( NZEROS     ) .LE. X ) XL( NZEROS     ) = X
               ENDIF

               IF ( XL( MODE ) .NE. XMIN ) GOTO 5000
            END DO
         ENDIF
 5000 CONTINUE

      RETURN
      END
