      SUBROUTINE GETMOD( IPROF, FILNAM, MAXM, RD, NRD, COMP,
     &   CK, PHIR, M, FREQ, TITLE )

C     Read in modes and extract values at rcvr depths

C     INPUT:
C        IPROF  unit number  for the mode set
C        FILNAM name of file for the mode set
C        MAXM   row dimension of PHIR in calling program
C        RD     vector of receiver depths where modes are to be evaluated
C        NRD number of such receivers
C        COMP   component (vertical, horizontal,...) for elastic problems
C               ignored for purely acoustic problems

C     OUTPUT:
C        CK     vector of eigenvalues
C        PHIR   matrix of tabulated modes
C        M      number of modes
C        FREQ   frequency
C        TITLE  title in mode file

C     If IPROF <> 0
C        we assume the mode file is assigned to unit
C        29 + IPROF and no open is necessary
C     else
C        we open a file with name FILNAM to read the modes

      INTEGER PRTFIL
      PARAMETER ( PRTFIL = 6, MAXN = 16001, MAXMED = 51, MAXNRD = 501 )

      INTEGER
     &   N( MAXMED ), IRD( MAXNRD )

      REAL
     &   RD( * ), Z( MAXN ), W( MAXNRD ), DEPTH( MAXMED ), RHO( MAXMED )

      COMPLEX
     &   PHIR( MAXM, * ), PHIT( MAXN ), CK( * ), KTOP2, KBOT2

      CHARACTER
     &   TITLE*( * ), FILNAM*( * ), MATER( MAXMED )*8,
     &   COMP*( *), BCTOP*1, BCBOT*1

C     ------ OK, Let's go. ------

      IF ( NRD .GT. MAXNRD ) THEN
         WRITE( PRTFIL, * ) 'NRD, MAXNRD = ', NRD, MAXNRD
         CALL ERROUT( PRTFIL, 'F', 'GETMOD', 'Too many receiver depths')
      ENDIF

C     *** Read the header data from the mode file ***

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

      IF ( M .LE. 0 ) RETURN

C     *** Locate indices of receiver points ***

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

C     *** Loop over receiver depths to check for safe interpolation ***

C     Receivers must be within a fraction of a wavelength
C     of tabulated pts. We accept one wavelength a 1500 m/s

      TOL = 1500.0 / FREQ

      DO IR = 1, NRD

         IZ = IRD( IR )
         WT = ABS( MIN( W( IR ), 1.0 - W( IR ) ) )

         IF ( RD( IR ) .LT. DEPTHT ) THEN
C           ------ Rcvr in upper halfspace
            CST = 0.0 ! should be passed by MODHDR
            IF ( CST .NE. 0.0 .OR. BCTOP(1:1) .NE. 'A' ) THEN
               WRITE( PRTFIL, * ) 'Receiver depth: ', RD( IR )
               WRITE( PRTFIL, * ) 'Highest valid depth: ', DEPTHT
               CALL ERROUT( PRTFIL, 'F', 'GETMOD',
     &                                'Rcvr above highest valid depth' )
            ENDIF

         ELSE IF ( RD( IR ) .GT. DEPTHB ) THEN
C           ------ Rcvr in lower halfspace
            CSB = 0.0 ! should be passed by MODHDR
            IF ( CSB .NE. 0.0 .OR. BCBOT(1:1) .NE. 'A' ) THEN
               WRITE( PRTFIL, * ) 'Receiver depth: ', RD( IR )
               WRITE( PRTFIL, * ) 'Lowest valid depth: ', DEPTHB
               CALL ERROUT( PRTFIL, 'F', 'GETMOD',
     &                                 'Rcvr below lowest valid depth' )
            ENDIF

         ELSE IF ( NTOT .GT. 1 ) THEN
C           ------ Rcvr between two grid points
            IF ( WT * ( Z( IZ + 1 ) - Z( IZ ) ) .GT. TOL ) THEN
            WRITE( PRTFIL, * ) 'Nearest depths: ', Z( IZ ), Z( IZ + 1 )
            WRITE( PRTFIL, * ) 'Tolerance: ', TOL
            CALL ERROUT( PRTFIL, 'F', 'GETMOD',
     &                        'Modes not tabulated near requested pt.' )
            ENDIF
         ELSE
C           ------ Rcvr near a single grid point
            IF ( ABS( RD( IR ) - Z( IZ ) ) .GT. TOL ) THEN
            WRITE( PRTFIL, * ) 'Rd, Tabulation depth ', RD( IR), Z( IZ )
            WRITE( PRTFIL, * ) 'Tolerance: ', TOL
            CALL ERROUT( PRTFIL, 'F', 'GETMOD',
     &                        'Modes not tabulated near requested pt.' )
            ENDIF
         ENDIF

      END DO   ! next receiver depth

C     *** Read in the modes ***

      DO MODE = 1, M

         CALL GETONE( MODE, MODFIL, NTOT, NMAT, W, IRD,
     &      N, MATER, NMEDIA, COMP,
     &      KTOP2, DEPTHT, BCTOP,
     &      KBOT2, DEPTHB, BCBOT,
     &      RD, NRD, CK, PHIT )

         DO IZ = 1, NRD
            PHIR( MODE, IZ ) = PHIT( IZ )
         END DO

      END DO   ! next mode

      CLOSE( MODFIL )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE MODHDR( IPROF, FILNAM,
     &   TITLE, FREQ, NMEDIA, NTOT, NMAT,
     &   N, MATER, DEPTH, RHO,
     &   BCTOP, RHOT, DEPTHT,
     &   BCBOT, RHOB, DEPTHB,
     &   M, MAXM, Z, CK, KTOP2, KBOT2, MODFIL )

C     Reads the header information from MODFIL
C     Note T suffix means top
C          B suffix means bottom

C     The mode file may be specified either by a profile number (IPROF)
C        or a name FILNAM
C     If IPROF =  0 then we assume you provide FILNAM
C     If IPROF <> 0 then a file name of the form MODFILxxx is built

C     FILNAM is the user-provided file name
C     FILNAMT is the temporary name we build
C     These have to be two separate variables do ensure there
C     is space allocated to construct the file name even when
C     the user calls us with FILNAM = ' '


      INTEGER PRTFIL
      PARAMETER ( PRTFIL = 6, PI = 3.141592 )

      INTEGER   N( * )
      REAL      DEPTH( * ), Z( * ), RHO( * )
      COMPLEX   CK( * ), CPT, CST, CPB, CSB, KTOP2, KBOT2
      CHARACTER TITLE*( * ), FILNAM*( * ), MATER( * )*8,
     &          BCTOP*1, BCBOT*1, FILNAMT*80

C     *** open MODFIL ***

      IF ( IPROF .EQ. 0 ) THEN
         MODFIL = 30
         FILNAMT = FILNAM   ! copy filename
      ELSE
         MODFIL = 29 + IPROF
         IF ( IPROF .LE. 9 ) THEN
            WRITE( FILNAMT, 1000 ) IPROF
 1000       FORMAT( 'MODFIL', I1.1 )
         ELSE IF ( IPROF .LE. 99 ) THEN
            WRITE( FILNAMT, 2000 ) IPROF
 2000       FORMAT( 'MODFIL', I2.2 )
         ELSE IF ( IPROF .LE. 999 ) THEN
            WRITE( FILNAMT, 3000 ) IPROF
 3000       FORMAT( 'MODFIL', I3.3 )
         ELSE
            STOP "FATAL ERROR in GETMOD"
         END IF
      ENDIF

C     ------ get record length
      OPEN ( UNIT = MODFIL, FILE = FILNAMT, STATUS = 'OLD', 
     &   ACCESS = 'DIRECT', FORM = 'UNFORMATTED', RECL = 20 )
      READ( MODFIL, REC = 1 ) LRECL
      CLOSE( MODFIL )

C     ------ re-open file with the correct record length
      OPEN ( UNIT = MODFIL, FILE = FILNAMT, STATUS = 'OLD', 
     &   ACCESS = 'DIRECT', FORM = 'UNFORMATTED', RECL = 4 * LRECL )

C     *** Read header info ***

      READ( MODFIL, REC = 1 ) LRECL, TITLE(1:80),
     &                        FREQ, NMEDIA, NTOT, NMAT
      READ( MODFIL, REC = 2 ) ( N( MED ), MATER( MED ), MED = 1, NMEDIA)
      READ( MODFIL, REC = 3 ) BCTOP(1:1), CPT, CST, RHOT, DEPTHT,
     &                        BCBOT(1:1), CPB, CSB, RHOB, DEPTHB
      READ( MODFIL, REC = 4 )
     &      ( DEPTH( MED ), RHO( MED ), MED = 1, NMEDIA )
      READ( MODFIL, REC = 5 ) M, LRECL

      IF ( M .GT. MAXM ) THEN
         WRITE( PRTFIL, * ) 'M = ', M, '   MAXM = ', MAXM
         CALL ERROUT( PRTFIL, 'F', 'MODHDR',
     &   'Insufficient storage to read all the modes: increase MAXM' )
      ENDIF

      IF ( M .EQ. 0 ) RETURN
      READ( MODFIL, REC = 6 ) ( Z(   IZ ), IZ = 1, NTOT )

C     *** Read in eigenvalues, CK( I ) ***

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

         READ( MODFIL, REC = 6 + M + IREC )
     &                        ( CK( MODE ), MODE = IFIRST, ILAST )

         IFIRST = ILAST + 1
      END DO

      IF ( BCTOP(1:1) .EQ. 'A' ) KTOP2 = ( 2.0 * PI * FREQ / CPT ) **2
      IF ( BCBOT(1:1) .EQ. 'A' ) KBOT2 = ( 2.0 * PI * FREQ / CPB ) **2

      RETURN
      END

C**********************************************************************C

      SUBROUTINE GETONE( MODE, MODFIL, NTOT, NMAT, W, IRD,
     &   N, MATER, NMEDIA, COMP,
     &   KTOP2, DEPTHT, BCTOP,
     &   KBOT2, DEPTHB, BCBOT,
     &   RD, NRD, CK, PHIR )

C     Read in a single eigenfunction and extract receiver values
C     Results are returned in PHIR

      PARAMETER ( MAXN = 16001 )

      LOGICAL    TUFLUK
      INTEGER    N( * ), IRD( * )
      REAL       RD( * ), W( * )
      COMPLEX    PHIR( * ), CK( * ), PHI( MAXN ),
     &           GAMT, GAMB, KTOP2, KBOT2
      COMPLEX*16 PEKRT, GAMMA2
      CHARACTER  MATER( * )*8, COMP*( *), BCTOP*1, BCBOT*1

      READ( MODFIL, REC = 6 + MODE ) ( PHI( J ), J = 1, NMAT )

C     *** Is there an elastic medium in the problem? ***

      TUFLUK = .FALSE.

      DO MED = 1, NMEDIA
         IF ( MATER( MED ) .EQ. 'ELASTIC' ) TUFLUK = .TRUE.
      END DO

C     ------ Extract the component specified by 'COMP'

      IF ( TUFLUK ) CALL EXTRACT( PHI, N, MATER, NMEDIA, COMP )

C     *** Extract values at receiver depths ***

      GAMT = 0.0
      GAMB = 0.0

C     --- n.b. should be using real( ck(mode)) for KRAKEN
      IF ( BCTOP(1:1) .EQ. 'A' ) THEN
         GAMMA2 = CK( MODE ) ** 2 - KTOP2
         GAMT = PEKRT( GAMMA2 )
      END IF

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

      DO IR = 1, NRD

         IF ( RD( IR ) .LT. DEPTHT ) THEN
C           ------ Rcvr in upper halfspace
            PHIR( IR ) = PHI( 1    ) *
     &                     EXP( -GAMT * ( DEPTHT - RD( IR  ) ) )
         ELSE IF ( RD( IR ) .GT. DEPTHB ) THEN
C           ------ Rcvr in lower halfspace
            PHIR( IR ) = PHI( NTOT ) *
     &                     EXP( -GAMB * ( RD( IR ) - DEPTHB ) )
         ELSE IF ( NTOT .GT. 1 ) THEN
            IZ = IRD( IR )
            PHIR( IR ) = PHI( IZ ) +
     &                     W( IR ) * ( PHI( IZ + 1 ) - PHI( IZ ) )
         ELSE   ! case where mode is tabulated at only one depth
            IZ = IRD( IR )
            PHIR( IR ) = PHI( IZ )
         ENDIF

      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE EXTRACT( PHI, N, MATER, NMEDIA, COMP )

C     For elastic problems where a stress-displacement vector is output,
C     extracts the desired component

      INTEGER   N( * )
      COMPLEX   PHI( * )
      CHARACTER MATER( * )*8, COMP*( *)

      J = 1
      K = 1

      DO MED = 1, NMEDIA

         DO I = 1, N( MED ) + 1
            IF ( MATER( MED ) .EQ. 'ACOUSTIC' ) THEN
               PHI( J ) = PHI( K )
               K = K + 1
            ELSE
               IF ( COMP .EQ. 'H' ) PHI( J ) = PHI( K     )
               IF ( COMP .EQ. 'V' ) PHI( J ) = PHI( K + 1 )
               IF ( COMP .EQ. 'T' ) PHI( J ) = PHI( K + 2 )
               IF ( COMP .EQ. 'N' ) PHI( J ) = PHI( K + 3 )
               K = K + 4
            ENDIF
            J = J + 1
         END DO

      END DO

      RETURN
      END
