      PROGRAM KRAKENC

C     Program for solving for ocean acoustic normal modes

C     Michael B. Porter

      INCLUDE 'comc.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 = 'KRAKENC-'

      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

C     *** Optionally read in bottom reflection coefficient data ***

      DEGRAD = PI / 180.0

      IF ( BOTOPT(1:1) .EQ. 'F' ) THEN
         WRITE( PRTFIL, * )
     &      'Using tabulated bottom reflection coefficient'
         OPEN( FILE = 'BRCFIL', UNIT = BRCFIL )
         READ( BRCFIL, * )
     &      NBPTS, ( THETAB( I ), RB( I ), PHIB( I ), I = 1, NBPTS )

         DO I = 1, NBPTS
            PHIB( I ) = -DEGRAD * PHIB( I )
         END DO
      ENDIF

C     *** Optionally read in top reflection coefficient data ***

      IF ( TOPOPT(2:2) .EQ. 'F' ) THEN
         OPEN( FILE = 'TRCFIL', UNIT = TRCFIL )
         WRITE( PRTFIL, * ) 'Using tabulated top reflection coefficient'
         READ( TRCFIL, * )
     &      NTPTS, ( THETAT( I ), RT( I ), PHIT( I ), I = 1, NTPTS )

         DO I = 1, NTPTS
            PHIT( I ) = DEGRAD * PHIT( I )
         END DO
      ENDIF

C     *** Optionally read in internal reflection coefficient data ***

      IF ( BOTOPT(1:1) .EQ. 'P' ) THEN
         OPEN( FILE = 'IRCFIL', UNIT = IRCFIL )
         WRITE( PRTFIL, * ) 'Reading precalculated refl. coeff. table'
         READ( IRCFIL, * ) NKPTS,
     &      ( XTAB(IK), FTAB(IK), GTAB(IK), ITAB(IK), IK = 1, NKPTS )
         IF ( NKPTS .GT. MXREF ) THEN
            WRITE( PRTFIL, * ) 'NKPTS, MXREF = ', NKPTS, MXREF
            CALL ERROUT( PRTFIL, 'F', 'KRAKENC',
     &             'Too many k-points in IRC file' )
         ENDIF
      ENDIF

      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', 'KRAKENC',
     &             '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 ( OMEGA / DBLE( SQRT( EXTRAP( 1, MODE ) ) ) .LT. CHIGH )
     &                            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 )

C        Zero out positive imaginary part which would cause
C        growth in range. Should be a small number.

C        --- use DIMAG for f77, AIMAG for f90 below
         IF ( DIMAG( CK( MODE ) ) .GT. 0.0 )
     &               CK( MODE ) = REAL( CK( MODE ) )
      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 'comc.f'

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

      COMMON /SPACE/ CP, CS, DUMMY( 4 * MAXN )

      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', 'KRAKENC',
     &             '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' )

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

         IF ( CS( I ) .EQ. ( 0.0, 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 * OMEGA2 / CP( J ) ** 2
            END DO
         ELSE

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

            IF ( SIGMA( MED ) .NE. 0.0 )
     &         CALL ERROUT( PRTFIL, 'F', 'KRAKENC',
     &               '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     *** Bottom properties ***

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

C     *** Top properties ***

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

      IF ( ELFLAG ) CMIN = 0.85 * CMIN
      CLOW = MAX( CLOW, 0.99 * CMIN )

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

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

      INCLUDE 'comc.f'

      COMPLEX*16 T1, T2, F1, F2

      CALL CLTIME   ! clear the timer
      CALL INIT     ! set up the finite difference mesh
      CALL SOLVE2   ! solve for the eigenvalues

      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     *** Richardson extrapolation to improve the accuracy ***

      ERROR = 10
      KEY = 2 * M / 3 + 1   ! index of element 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 )   ! check elapsed time
      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)

C     A sequence of restarts is done to compute a complete set of modes
C     for the first two sets (ISET=1,2).
C     'Complete' in this case means that all the modes in the
C     user-specified interval are found. However, if, for ISET=2,
C     the same number of modes is found as for ISET=1 then it is assumed
C     that the set is complete and no restarts are performed.

      INCLUDE 'comc.f'

      INTEGER ISEED( 4 )
      REAL*4 RAN
      COMPLEX*16 X, P( 10 ), XTEMP( MAXM ), CTRY
      CHARACTER*80 ERRMSG

      SAVE ISEED
      DATA ISEED /2301, 4562, 3219, 1823/

      X = ( 1.0, 0.0 ) * OMEGA2 / CLOW ** 2
      MAXIT = 5000   ! maximum # of iterations in root-finder

C     --- MXTRIES = # of restarts for the root finder
C     --- ITRY = counter that keeps track of # of tries

      IF ( TOPOPT(4:4) .EQ. '.' .AND. ISET .LE. 2 ) THEN
         MXTRIES = MAXM
      ELSE
         MXTRIES = 1
      ENDIF
      ITRY = 1

C     --- compute size of matrix (used to check grid is fine enough)
      NTOT = 0
      DO MED = NFACT, NLACT
         NTOT = NTOT + N( MED )
      END DO

      DO MODE = 1, M

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

 5500    CONTINUE
         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 = 10000.0 * MACHEP * ABS( X )
         CALL ZSECCX( X, TOL, IT, MAXIT, ERRMSG )

C        --- any problems in root finder?

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

C        *** Is the mode inside the user specified spectral limits? ***

         IF ( SQRT( OMEGA2 ) / REAL( SQRT( X ) ) .GT. CHIGH ) THEN
C           --- Mode outside, restart at a random point
            CTRY = CLOW     + RAN( ISEED ) * ( CHIGH - CLOW )
     &      + ( 0.0, 0.01 ) * RAN( ISEED ) * CLOW
         
            X = ( 1.0, 0.0 ) * OMEGA2 / CTRY ** 2

            IF ( ITRY .LT. MXTRIES ) THEN
               MXTRIES = MAX( MXTRIES, 5 * MODE )
               ITRY = ITRY + 1
               WRITE( *, * ) 'Restart at phase speed = ', CTRY
               GOTO 5500
            ELSE
               M = MODE - 1
               GOTO 6000   ! done searching for modes
            ENDIF
         ENDIF

C        *** Are we getting too many modes for the grid? ***

         IF ( MODE .GT. NTOT / 5 ) THEN

            CALL ERROUT( PRTFIL, 'W', 'KRAKENC',
     &      'Mesh too coarse to compute all the modes in the spectrum: '
     &      //'Increase the mesh size' )
  
            M = MODE
            EVMAT( ISET, MODE ) = X
            XTEMP( MODE )       = X
            GOTO 6000   ! done searching for modes
         ENDIF

         EVMAT( ISET, MODE ) = X
         XTEMP( MODE )       = X

      END DO   ! next mode

C     --- If no modes, open a dummy MODFIL for use by FIELD3D

 6000 IF ( M .EQ. 0 ) THEN
         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', 'KRAKENC',
     &                'No modes for given phase speed interval' )

      ENDIF

C     --- Check whether we have hit the ceiling on max # modes

      IF ( M .EQ. MAXM ) THEN
         WRITE( PRTFIL, * ) 'Number of modes = ', M
         WRITE( PRTFIL, * ) 'Number of modes allowable = ', MAXM
         CALL ERROUT( PRTFIL, 'W', 'KRAKENC',
     &     'Too many modes'//
     &     'Program will compute as many as it can' )
      ENDIF

      CALL ORDER( XTEMP, M )   ! order eigenvalues by real part

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

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

C     FUNCT( X ) = 0 is the dispersion relation

      INCLUDE 'comc.f'

      PARAMETER ( ROOF = 1.0E5, FLOOR = 1.0E-5, IPOWR = 5, IPOWF = -5 )
      COMPLEX*16 X, DELTA, F, G, F1, G1
      CHARACTER BCTYPE*1

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
      IPOW  = IPOW + IPOW1
C     --- squaring the dispersion reln. seems to make the
C         root-finder more robust
      delta = delta ** 2
      ipow  = 2 * ipow

C     *** Deflate previous roots ***

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

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

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

   60          IF ( ABS( DBLE( 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 'comc.f'

      COMPLEX*16 X, F, G, P0, P1, P2, H2K2
      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

C           --- Scale if necessary
   60       IF ( ABS( DBLE( 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 'comc.f'

      INTEGER IZTAB( MAXNZT )
      REAL Z( MAXN ), ZTAB( MAXNZT ), WTS( MAXNZT )
      COMPLEX PHITAB( MAXNZT )
      COMPLEX*16 RV1( MAXN ), RV2( MAXN ), RV3( MAXN ), RV4( MAXN ),
     &   PHI( MAXN ), D( MAXN ), E( MAXN ), W, XH2, X, F, G
      CHARACTER BCTOP*1, BCBOT*1

      COMMON /SPACE/ RV1, RV2, RV3, RV4

      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 ( REAL( 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', 'KRAKENC-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 ), 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 'comc.f'

      COMPLEX*16 PHI(*), X, X1, F1, G1, X2, F2, G2,
     &   DRHODX, DETADX, SQNRM, RN, SCALE
      CHARACTER BCTYPE*1

      SQNRM = 0.0

C     *** Loop to compute norm of eigenvector ***

      J = 1

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

         SQNRM = SQNRM + 0.5 * H( MED ) * PHI( J ) ** 2 / RHOM

         DO I = 1, N( MED ) - 1
            J = J + 1
            SQNRM = SQNRM +    H( MED ) * PHI( J ) ** 2 / RHOM
         END DO

         J = J + 1
         SQNRM = SQNRM + 0.5 * H( MED ) * PHI( J ) ** 2 / RHOM

      END DO

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

      X1 = 0.9999999 * X
      X2 = 1.0000001 * 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
      SCALE = 1.0 / SQRT( RN )
      IF ( REAL( PHI( ITP ) ) .LT. 0.0 ) SCALE = -SCALE

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

C     *** Compute interfacial scatter loss ***

      CALL SCAT( PHI, X )

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

C     Figure scatter loss

      INCLUDE 'comc.f'

      COMPLEX*16 PHI( * ), X, PERK, ETA1SQ, ETA2SQ,
     &   KUPING, U, PEKRT
      CHARACTER BCTYPE*1

      PERK = 0.0
      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      = PEKRT( 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
            ELSE
C              --- Tabulated
               RHO1   = 0.0
               ETA1SQ = 0.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
            ELSE
C              --- Tabulated
               RHO2   = 0.0
               ETA2SQ = 0.0
            ENDIF
         ELSE
            RHO2   = RHO( L + 1 )
            ETA2SQ = ( 2.0 + B1( L + 1 ) ) / H( MED + 1 ) ** 2 - X
         ENDIF

         PERK = PERK + KUPING( SIGMA( MED + 1 ), ETA1SQ, RHO1,
     &                         ETA2SQ, RHO2, PHI( J ), U )

      END DO

      CK( MODE ) = PERK

      RETURN
      END
C**********************************************************************C
      SUBROUTINE ORDER( X, N )

C     Does an insertion sort of the complex vector X in order
C     of decreasing real part

C     At the Ith step, the first I-1 positions contain a sorted
C     vector.  We shall insert the Ith value into its place in that
C     vector shifting up to produce a new vector of length I.

      COMPLEX*16 X( * ), T

      IF ( N .EQ. 1 ) RETURN

      DO I = 2, N

         T = X( I )

         IF ( REAL( T ) .GT. REAL( X( 1 ) ) ) THEN
C           --- Goes in the first position
            CALL SHIFTINC( X, T, 1, I - 1 )

         ELSE IF ( REAL( T ) .GT. REAL( X( I - 1 ) ) ) THEN

C           *** Binary search for its place ***

            IRIGHT = I - 1
            ILEFT  = 1

 1000       IMID = ( ILEFT + IRIGHT ) / 2

            IF ( REAL( T ) .GT. REAL( X( IMID ) ) ) THEN
               IRIGHT = IMID
            ELSE
               ILEFT  = IMID
            ENDIF

            IF ( IRIGHT .GT. ILEFT + 1 ) GOTO 1000

C           --- Shift and insert
            CALL SHIFTINC( X, T, IRIGHT, I - 1 )

         ENDIF

      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE SHIFTINC( X, T, I1, I2 )

C     Shift values (I1, I2) one position to the right
C     and insert T at position I1

      COMPLEX*16 X( * ), T

      DO K = I2, I1, -1
         X( K + 1 ) = X( K )
      END DO

      X( I1 ) = T

      RETURN
      END
