      SUBROUTINE READIN( TITLE, FREQ, ALPHA, NBEAMS, MXNBEAMS, ISINGL,
     &   NIMAGE, IBWIN, DELTAS, MXN, ZBOX, RBOX, EPMULT, RLOOP,
     &   SD, NSD, MXNSD, RD, NRD, MXNRD, R, NR, MXNR,
     &   RNG, DEP, TBR, TBZ, NBTYPTS, MXBTY,
     &   THETATOP, RTOP, PHITOP, NTOPPTS,
     &   THETABOT, RBOT, PHIBOT, NBOTPTS, MXREF,
     &   TOPOPT, DEPTHT, CPT, RHOT,
     &   BOTOPT,  DEPTHB, CPB, RHOB, RUNTYP, BEAMTYP )

C     Routine to read in and echo all the input data

C     Note that default values of SSP, DENSITY, Attenuation
C     will not work

      IMPLICIT REAL*4 ( A-H, O-Z )
      INTEGER ENVFIL, PRTFIL, BRCFIL, TRCFIL, BTYFIL, RAYFIL, ARRFIL
      PARAMETER ( ENVFIL =  5, PRTFIL =  6, BRCFIL = 31, TRCFIL = 32,
     &            BTYFIL = 41, RAYFIL = 21, ARRFIL = 36,
     &            PI = 3.14159265, DEGRAD = 2.0 * PI / 360.0 )

      REAL      SD( * ), RD( * ), R( * ), ALPHA( * ), LAMX10, LAMX66
      REAL*8    FREQD, RHOT, RHOB, BUMDEN, ETA, XI
      COMPLEX*16 CPT, CPB, CST, CSB
      CHARACTER TITLE*80, TOPOPT*4, BOTOPT*2, RUNTYP*2, BEAMTYP*2,
     &          SSPTYP*1, BCTYPE*1, ATUNIT*1

      REAL*8 RNG( * ), DEP( * ), TBR( * ), TBZ( * )
      REAL*8 THETABOT( * ), RBOT( * ), PHIBOT( * ),
     &       THETATOP( * ), RTOP( * ), PHITOP( * )

      REAL*8 ALPHAR, BETAR, ALPHAI, BETAI, RHOR
      COMMON /CPREV/ ALPHAR, BETAR, RHOR, ALPHAI, BETAI

      ALPHAR = 1500.0
      BETAR  = 0.0
      RHOR   = 1.0
      ALPHAI = 0.0
      BETAI  = 0.0
      ISINGL = 0

      TITLE(1:9) = 'BELLHOP- '
      READ(  ENVFIL, * ) TITLE(10:80)
      WRITE( PRTFIL, * ) TITLE

      READ(  ENVFIL, *    ) FREQ
      WRITE( PRTFIL, 1000 ) FREQ
 1000 FORMAT(' Frequency = ', G11.4, / )

      FREQD = FREQ

      READ(  ENVFIL, * ) NMEDIA
      WRITE( PRTFIL, * ) 'Dummy parameter NMEDIA = ', NMEDIA

C     ****** Read in SSP data ******

      READ(  ENVFIL, * ) TOPOPT
      WRITE( PRTFIL, * )

      SSPTYP(1:1) = TOPOPT(1:1)
      BCTYPE(1:1) = TOPOPT(2:2)
      ATUNIT(1:1) = TOPOPT(3:3)

C     *** Top BC ***

      IF ( BCTYPE(1:1) .EQ. 'A' ) WRITE( PRTFIL, 2000 )
 2000 FORMAT( //, '      Z          ',
     &      'ALPHAR     BETAR      RHO       ALPHAI     BETAI', / )

      CALL TOPBOT( ENVFIL, PRTFIL, FREQD, BCTYPE, ATUNIT,
     &   CPT, CST, RHOT, BUMDEN, ETA, XI )

      DEPTHT = 0.0   ! should be read in somewhere

C     *** Ocean SSP ***

      READ(  ENVFIL, * ) NPTS, SIGMA, DEPTHB
      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * ) 'Depth = ', DEPTHB

      IF ( TOPOPT(1:1) .EQ. 'A' ) THEN
         WRITE( PRTFIL, * ) 'Analytic SSP option'
      ELSE
         CALL SSP( RDUMMY, DEPTHB, C, CR, CZ, CRR, CRZ, CZZ,
     &             TOPOPT, 'INIT' )
      ENDIF

C     *** Bottom BC ***

      READ( ENVFIL, * ) BOTOPT, SIGMA

      WRITE( PRTFIL, * )
      WRITE( PRTFIL, 3000 ) SIGMA
 3000 FORMAT( 33X, '( RMS roughness = ', G10.3, ' )' )

      CALL TOPBOT( ENVFIL, PRTFIL, FREQD, BOTOPT, ATUNIT,
     &    CPB, CSB, RHOB, BUMDEN, ETA, XI )

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

      CALL SDRD( ENVFIL, PRTFIL, DEPTHT, DEPTHB,
     &   MXNSD, NSD, SD,
     &   MXNRD, NRD, RD )

C     *** Read receiver ranges ***

      CALL RANGES( ENVFIL, PRTFIL, MXNR, NR, R )

C     *** Read run type ***

      READ(  ENVFIL, * ) RUNTYP(1:2)

      WRITE( PRTFIL, * )
      IF      ( RUNTYP(1:1) .EQ. 'R' ) THEN
         WRITE( PRTFIL, * ) 'Ray trace run'
      ELSE IF ( RUNTYP(1:1) .EQ. 'E' ) THEN
         WRITE( PRTFIL, * ) 'Eigenray trace run'
      ELSE IF ( RUNTYP(1:1) .EQ. 'I' ) THEN
         WRITE( PRTFIL, * ) 'Incoherent TL calculation'
      ELSE IF ( RUNTYP(1:1) .EQ. 'S' ) THEN
         WRITE( PRTFIL, * ) 'Semi-coherent TL calculation'
      ELSE IF ( RUNTYP(1:1) .EQ. 'C' ) THEN
         WRITE( PRTFIL, * ) 'Coherent TL calculation'
      ELSE IF ( RUNTYP(1:1) .EQ. 'A' ) THEN
         WRITE( PRTFIL, * ) 'Arrivals calculation'
      ELSE
         CALL ERROUT( PRTFIL, 'F', 'READIN',
     &                'Unknown RUNTYP selected' )
      ENDIF

      IF ( RUNTYP(2:2) .EQ. 'C' ) THEN
         WRITE( PRTFIL, * ) 'Cartesian beams'
      ELSE IF ( RUNTYP(2:2) .EQ. 'R' ) THEN
         WRITE( PRTFIL, * ) 'Ray centered beams'
      ELSE
         RUNTYP(2:2) = 'G'
         WRITE( PRTFIL, * ) 'Geometric beams'
      ENDIF

C     *** Read beam angles ***

      NBEAMS     = 0
      ALPHA( 3 ) = -999.9
      READ( ENVFIL, * ) NBEAMS, ( ALPHA( IBEAM ), IBEAM = 1, NBEAMS )

      IF ( NBEAMS .GT. MXNBEAMS ) THEN
         WRITE( PRTFIL, * ) 'NBEAMS, MXNBEAMS = ', NBEAMS, MXNBEAMS
         CALL ERROUT( PRTFIL, 'F', 'READIN', 'Too many beams'  )
      ENDIF

      CALL SUBTAB( ALPHA, NBEAMS )
      CALL SORT(   ALPHA, NBEAMS )

      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * ) 'Number of beams   = ', NBEAMS
      IF ( NBEAMS .GE. 1 ) WRITE( PRTFIL, * )
     &                ( ALPHA( IBEAM ), IBEAM = 1, MIN( NBEAMS, 51 ) )

C     --- convert to radians
      DO IBEAM = 1, NBEAMS
         ALPHA( IBEAM ) = DEGRAD * ALPHA( IBEAM )
      END DO

C     *** Limits for tracing beams ***

      READ(  ENVFIL, * ) DELTAS, ZBOX, RBOX

C     --- Automatic step size selection < 10 wavelengths

      IF ( DELTAS .EQ. 0.0 ) THEN

         LAMX10 =  15000.0 / FREQ               ! 10 wavelengths
         LAMX66 = 100000.0 / FREQ               ! 66 wavelengths
         DBY20  = ( DEPTHB - DEPTHT ) / 20.0
         IF      ( RUNTYP(1:1) .EQ. 'R' ) THEN
           DELTAS = MIN( LAMX66, DBY20 )
         ELSE IF ( RUNTYP(1:1) .EQ. 'E' ) THEN
           DELTAS = MIN( LAMX66, DBY20 )
         ELSE IF ( RUNTYP(1:1) .EQ. 'I' ) THEN
           DELTAS = MIN( LAMX66, DBY20 )
         ELSE IF ( RUNTYP(1:1) .EQ. 'S' ) THEN
           DELTAS = MIN( LAMX66, DBY20 )
         ELSE IF ( RUNTYP(1:1) .EQ. 'C' ) THEN
           DELTAS = MIN( LAMX10, DBY20 )
         ELSE IF ( RUNTYP(1:1) .EQ. 'A' ) THEN
           DELTAS = MIN( LAMX10, DBY20 )
         ENDIF

         DELTAS = MAX( DELTAS, RBOX / MXN )      ! at most do MXN steps
         DELTAS = MIN( DELTAS, DEPTHB - DEPTHT ) ! DELTAS < channel height 
      ENDIF

      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * ) 'Step length,     DELTAS = ', DELTAS
      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * ) 'Maximum ray depth, ZBOX = ', ZBOX
      WRITE( PRTFIL, * ) 'Maximum ray range, RBOX = ', RBOX

      RBOX = 1000.0 * RBOX

C     ****** Beam characteristics ******

C     note: curvature change can cause overflow in grazing case
C     suppress by setting BEAMTYP(2:2) = 'Z'
      IF ( RUNTYP(2:2) .EQ. 'G' ) THEN
         BEAMTYP(1:2) = 'MS'
         RLOOP        = 1.0
         EPMULT       = 1.0
      ELSE

         READ(  ENVFIL, * ) BEAMTYP(1:2), EPMULT, RLOOP, ISINGL
         WRITE( PRTFIL, * )
         WRITE( PRTFIL, * )
         WRITE( PRTFIL, * ) 'Type of beam = ', BEAMTYP(1:1)

         IF ( BEAMTYP(2:2) .EQ. 'D' ) THEN
            WRITE( PRTFIL, * ) 'Curvature doubling invoked'
         ELSE IF ( BEAMTYP(2:2) .EQ. 'Z' ) THEN
            WRITE( PRTFIL, * ) 'Curvature zeroing invoked'
         ELSE IF ( BEAMTYP(2:2) .EQ. 'S' ) THEN
            WRITE( PRTFIL, * ) 'Standard curvature condition'
         ELSE
            CALL ERROUT( PRTFIL, 'F', 'READIN',
     &                'Unknown curvature condition' )
         ENDIF

         WRITE( PRTFIL, * ) 'Epsilon multiplier', EPMULT
         WRITE( PRTFIL, * ) 'Range for choosing beam width', RLOOP
         WRITE( PRTFIL, * ) 'Index of beam to trace = ', ISINGL

C        ****** Images, windows ******

         READ(  ENVFIL, * ) NIMAGE, IBWIN
         WRITE( PRTFIL, * )
         WRITE( PRTFIL, * ) 'Number of images, NIMAGE = ', NIMAGE
         WRITE( PRTFIL, * ) 'Beam windowing parameter,  IBWIN = ', IBWIN

      ENDIF

      WRITE( PRTFIL, * )

C     *** Write appropriate header information ***

      IF ( RUNTYP(1:1) .EQ. 'R' .OR. RUNTYP(1:1) .EQ. 'E' ) THEN
         OPEN ( FILE = 'RAYFIL', UNIT = RAYFIL, FORM = 'FORMATTED' )
         WRITE( RAYFIL, * ) '''', TITLE(1:50), ''''
         WRITE( RAYFIL, * ) FREQ
         WRITE( RAYFIL, * ) NBEAMS
         WRITE( RAYFIL, * ) DEPTHT
         WRITE( RAYFIL, * ) DEPTHB
      ELSE IF ( RUNTYP(1:1) .EQ. 'A' ) THEN
         OPEN ( FILE = 'ARRFIL', UNIT = ARRFIL, FORM = 'FORMATTED' )
         WRITE( ARRFIL, * ) FREQ, NSD, NRD, NR
         WRITE( ARRFIL, * ) ( SD( I ), I = 1, NSD )
         WRITE( ARRFIL, * ) ( RD( I ), I = 1, NRD )
         WRITE( ARRFIL, * ) ( R( IR ), IR = 1, NR )
      ELSE
         CALL HEADER( TITLE, SD, NSD, RD, NRD,
     &               R, NR, FREQ, 0.0, '          ', 0.0, 0.0, 0.0 )
      END IF

      CLOSE( ENVFIL )

C     *** Optionally read in bottom bathymetry ***

      IF ( BOTOPT(2:2) .EQ. 'B' ) THEN
         WRITE( PRTFIL, * ) '*********************************'
         WRITE( PRTFIL, * )
     &      'Using bottom-bathymetry file'

         OPEN ( FILE = 'BTYFIL', UNIT = BTYFIL, STATUS = 'OLD' )
         READ(  BTYFIL, * ) NBTYPTS
         WRITE( PRTFIL, * ) 'Number of bathymetry points', NBTYPTS

         IF ( NBTYPTS .GT. MXBTY ) THEN
           WRITE( PRTFIL, * ) 'NBTYPTS, MXBTY = ', NBTYPTS, MXBTY
           CALL ERROUT( PRTFIL, 'F', 'READIN',
     &        'Too many bathymetry points'  )
         ENDIF

         WRITE( PRTFIL, * )
         WRITE( PRTFIL, * ) ' Range (km)  Depth (m)'
         DO I = 1, NBTYPTS
            READ(  BTYFIL, * ) RNG( I ), DEP( I )
            WRITE( PRTFIL, 4000 ) RNG( I ), DEP( I )
 4000       FORMAT( 2G11.3 )
            RNG( I ) = 1000.0 * RNG( I )   ! Convert km to m
         ENDDO

C        --- compute tangent to bottom
         DO I = 1, NBTYPTS - 1
            TBR( I ) = RNG( I + 1 ) - RNG( I )
            TBZ( I ) = DEP( I + 1 ) - DEP( I )
            RLEN = SQRT( TBR( I ) * TBR( I ) + TBZ( I ) * TBZ( I ) )
            TBR( I ) = TBR( I ) / RLEN
            TBZ( I ) = TBZ( I ) / RLEN
         END DO
      ELSE   ! no bathymetry give, use SSP depth for flat bottom
         NBTYPTS = 2
         RNG( 1 ) = 0.0
         DEP( 1 ) = DEPTHB
         RNG( 2 ) = RBOX
         DEP( 2 ) = DEPTHB
         TBR( 1 ) = 1.0
         TBZ( 1 ) = 0.0
      ENDIF

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

      IF ( BOTOPT(1:1) .EQ. 'F' ) THEN
         WRITE( PRTFIL, * )
     &      'Using tabulated bottom reflection coefficient'

         OPEN ( FILE = 'BRCFIL', UNIT = BRCFIL )
         READ( BRCFIL, * ) NBOTPTS,
     &      ( THETABOT( I ), RBOT( I ), PHIBOT( I ), I = 1, NBOTPTS )

         IF ( NBOTPTS .GT. MXREF ) THEN
           WRITE( PRTFIL, * ) 'NBOTPTS, MXREF = ', NBOTPTS, MXREF
           CALL ERROUT( PRTFIL, 'F', 'READIN',
     &        'Too many bottom reflection-coefficient points'  )
         ENDIF

         DO I = 1, NBOTPTS
            PHIBOT( I ) = DEGRAD * PHIBOT( I )
         END DO
      ENDIF

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

      IF ( TOPOPT(2:2) .EQ. 'F' ) THEN
         WRITE( PRTFIL, * ) 'Using tabulated top reflection coefficient'
         OPEN ( FILE = 'TRCFIL', UNIT = TRCFIL )
         READ( TRCFIL, * ) NTOPPTS,
     &      ( THETATOP( I ), RTOP( I ), PHITOP( I ), I = 1, NTOPPTS )

         IF ( NTOPPTS .GT. MXREF ) THEN
           WRITE( PRTFIL, * ) 'NTOPPTS, MXREF = ', NTOPPTS, MXREF
           CALL ERROUT( PRTFIL, 'F', 'READIN',
     &        'Too many top reflection-coefficient points'  )
         ENDIF

         DO I = 1, NTOPPTS
            PHITOP( I ) = DEGRAD * PHITOP( I )
         END DO
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE SSP( R, Z, C, CR, CZ, CRR, CRZ, CZZ, OPT, TASK )

C     Call the particular profil routine based specified by TASK

C     SSP is expected to perform two tasks:
C        TASK = 'TAB'  then tabulate CP, CS, RHOT 
C        TASK = 'INIT' then initialize

      INTEGER PRTFIL
      PARAMETER ( PRTFIL = 6 )

      CHARACTER OPT*4, TASK*4

      IF      ( OPT(1:1) .EQ. 'N' ) THEN
         CALL N2LIN(  R, Z, C, CR, CZ, CRR, CRZ, CZZ, TASK )
      ELSE IF ( OPT(1:1) .EQ. 'C' ) THEN
         CALL CLIN(   R, Z, C, CR, CZ, CRR, CRZ, CZZ, TASK )
      ELSE IF ( OPT(1:1) .EQ. 'S' ) THEN
         CALL CCUBIC( R, Z, C, CR, CZ, CRR, CRZ, CZZ, TASK )
      ELSE IF ( OPT(1:1) .EQ. 'A' ) THEN
         CALL ANALYT( R, Z, C, CR, CZ, CRR, CRZ, CZZ, TASK )
      ELSE
         WRITE( PRTFIL, * ) 'Non-existent profile option: ', OPT(1:1)
         STOP
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE N2LIN( R, Z, C, CR, CZ, CRR, CRZ, CZZ, TASK )

C     N2-linear interpolation of SSP data

      INTEGER ENVFIL, PRTFIL
      PARAMETER ( ENVFIL = 5, PRTFIL = 6, MXSSP = 2001, PREC = 1.0E-6 )

      REAL ZV( MXSSP ), CV( MXSSP ), N2V( MXSSP ), N2ZV( MXSSP )
      CHARACTER TASK*4
      SAVE ZTMP, CTMP, ZV, CV, N2V, N2ZV, NSSP, LAYER

      IF ( TASK(1:4) .EQ. 'INIT' ) THEN

C        *** Section to read in SSP data ***

         WRITE( PRTFIL, * )
         WRITE( PRTFIL, * ) 'N2-Linear SSP option'
         WRITE( PRTFIL, * )
         WRITE( PRTFIL, * ) 'Sound speed profile:'

         DEPTH = Z
         NSSP  = 1

         DO IZ = 1, MXSSP

            READ(  ENVFIL, * ) ZTMP, CTMP
            ZV( IZ ) = ZTMP
            CV( IZ ) = CTMP
            WRITE( PRTFIL, * ) ZV( IZ ), CV( IZ )
            N2V( IZ ) = 1.0 / CV( IZ ) ** 2

            IF ( IZ .GT. 1 ) N2ZV( IZ - 1 )  = 
     &                ( N2V( IZ ) - N2V( IZ - 1 ) ) /
     &                (  ZV( IZ ) -  ZV( IZ - 1 ) )

C           ------ Have we read the last point?
            IF ( ABS( ZV( IZ ) - DEPTH ) .LT. PREC ) THEN
               LAYER = 1
               RETURN
            ENDIF

            NSSP = NSSP + 1

         END DO

C        ------ Fall through means too many points in the profile

         WRITE( PRTFIL, * ) 'Max. #SSP points: ', MXSSP
         CALL ERROUT( PRTFIL, 'F', 'N2LIN',
     &                'Number of SSP points exceeds limit' )

      ELSE

C        *** Section to return SSP info ***

         IF ( Z .LT. ZV( LAYER ) .OR. Z .GT. ZV( LAYER + 1 ) ) THEN

            DO IZ = 2, NSSP   ! Search for bracketting depths
               IF ( Z .LT. ZV( IZ ) ) THEN
                  LAYER = IZ - 1
                  GOTO 3000
               ENDIF
            END DO

         ENDIF


3000     CONTINUE

         W = ( Z - ZV( LAYER ) ) / ( ZV( LAYER + 1 ) - ZV( LAYER ) )

         C     = 1.0 / SQRT( ( 1.0 - W ) * N2V( LAYER )
     &                               + W * N2V( LAYER + 1 ) )
         CZ  = -0.5 * C * C * C * N2ZV( LAYER )
         CR  = 0.0
         CRR = 0.0
         CRZ = 0.0
         CRR = 0.0
         CZZ = 3.0 * CZ * CZ / C
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE CLIN( R, Z, C, CR, CZ, CRR, CRZ, CZZ, TASK )

C     C-linear interpolation of SSP data

      INTEGER ENVFIL, PRTFIL
      PARAMETER ( ENVFIL = 5, PRTFIL = 6, MXSSP = 2001, PREC = 1.0E-6 )

      REAL ZV( MXSSP ), CV( MXSSP ), CZV( MXSSP )
      CHARACTER TASK*4
      SAVE ZTMP, CTMP, ZV, CV, CZV, NSSP, LAYER

      IF ( TASK(1:4) .EQ. 'INIT' ) THEN

C        *** Section to read in SSP data ***

         WRITE( PRTFIL, * )
         WRITE( PRTFIL, * ) 'C-Linear SSP option'
         WRITE( PRTFIL, * )
         WRITE( PRTFIL, * ) 'Sound speed profile:'

         DEPTH = Z
         NSSP = 1

         DO IZ = 1, MXSSP

            READ(  ENVFIL, * ) ZTMP, CTMP
            ZV( IZ ) = ZTMP
            CV( IZ ) = CTMP
            WRITE( PRTFIL, * ) ZV( IZ ), CV( IZ )

            IF ( IZ .GT. 1 ) CZV( IZ - 1 )  = 
     &                ( CV( IZ ) - CV( IZ - 1 ) ) /
     &                ( ZV( IZ ) - ZV( IZ - 1 ) )

C           ------ Have we read the last point?
            IF ( ABS( ZV( IZ ) - DEPTH ) .LT. PREC ) THEN
               LAYER = 1
               RETURN
            ENDIF

            NSSP = NSSP + 1

         END DO

C        ------ Fall through means too many points in the profile

         WRITE( PRTFIL, * ) 'Max. #SSP points: ', MXSSP
         CALL ERROUT( PRTFIL, 'F', 'N2LIN',
     &                'Number of SSP points exceeds limit' )

      ELSE

C        *** Section to return SSP info ***

         IF ( Z .LT. ZV( LAYER ) .OR. Z .GT. ZV( LAYER + 1 ) ) THEN
            DO IZ = 2, NSSP   ! Search for bracketting depths
               IF ( Z .LT. ZV( IZ ) ) THEN
                  LAYER = IZ - 1
                  GOTO 3000
               ENDIF
            END DO

         ENDIF

 3000    CZ  = CZV( LAYER )
         C   = CV( LAYER ) + ( Z - ZV( LAYER ) ) * CZ
         CR  = 0.0
         CRR = 0.0
         CRZ = 0.0
         CRR = 0.0
         CZZ = 0.0
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE CCUBIC( R, Z, C, CR, CZ, CRR, CRZ, CZZ, TASK )

C     Cubic spline interpolation

      INTEGER ENVFIL, PRTFIL
      PARAMETER ( ENVFIL = 5, PRTFIL = 6, MXSSP = 2001, PREC = 1.E-6 )

      REAL ZV( MXSSP ), CV( MXSSP ), CVS( 4, MXSSP )
      CHARACTER TASK*4

      SAVE ZTMP, CTMP, ZV, CV, CVS, NSSP, LAYER

      IF ( TASK(1:4) .EQ. 'INIT' ) THEN

C        *** Task 'INIT' for initialization ***

         WRITE( PRTFIL, * )
         WRITE( PRTFIL, * ) 'Spline SSP option'
         WRITE( PRTFIL, * )
         WRITE( PRTFIL, * ) 'Sound speed profile:'

         DEPTH = Z
         NSSP  = 1

         DO IZ = 1, MXSSP

            READ(  ENVFIL, * ) ZTMP, CTMP
            ZV( IZ ) = ZTMP
            CV( IZ ) = CTMP
            WRITE( PRTFIL, * ) ZV( IZ ), CV( IZ )

            CVS( 1, IZ ) = CV( IZ )

            IF ( ABS( ZV( IZ ) - DEPTH ) .LT. PREC ) THEN
               LAYER = 1

C              ------ Compute spline coefs
               IBCBEG = 0
               IBCEND = 0
               CALL CSPLINE( ZV, CVS(1, 1), NSSP, IBCBEG, IBCEND, NSSP )

               RETURN
            ENDIF

            NSSP = NSSP + 1

         END DO

C        ------ Fall through means too many points in the profile

         WRITE( PRTFIL, * ) 'Max. #SSP points: ', MXSSP
         CALL ERROUT( PRTFIL, 'F', 'N2LIN',
     &                'Number of SSP points exceeds limit' )

      ELSE

C        *** Section to return SSP info ***

         IF ( Z .LT. ZV( LAYER ) .OR. Z .GT. ZV( LAYER+1 ) ) THEN
            DO IZ = 2, NSSP   ! Search for bracketting depths
               IF ( Z .LT. ZV( IZ ) ) THEN
                  LAYER = IZ - 1
                  GOTO 3000
               ENDIF
            END DO

         ENDIF

 3000    HSPLNE = Z - ZV( LAYER )

C         C   = SPLINE(   CVS( 1, LAYER ), HSPLNE )
C         CZ  = SPLINEX(  CVS( 1, LAYER ), HSPLNE )
C         CZZ = SPLINEXX( CVS( 1, LAYER ), HSPLNE )

         CALL SPLINEALL( CVS( 1, LAYER ), HSPLNE, C, CZ, CZZ )

         CR  = 0.0
         CRR = 0.0
         CRZ = 0.0
         CRR = 0.0
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE TOPBOT( ENVFIL, PRTFIL, FREQ, BCTYPE, ATUNIT,
     &   CPHS, CSHS, RHOHS, BUMDEN, ETA, XI )

C     Handles top and bottom boundary conditions
C
C     Input:
C         ENVFIL: Environmental file
C         PRTFIL: Print file
C         FREQ:   Frequency
C         BCTYPE: Boundary condition type
C
C     Output:
C        CPHS:    P-wave speed in halfspace
C        CSHS:    S-wave speed in halfspace
C        RHOHS:   density in halfspace
C
C        BUMDEN:  Bump density
C        ETA:     Principal radius 1
C        XI:      Principal radius 2

      IMPLICIT NONE

      INTEGER ENVFIL, PRTFIL
      REAL*8 ALPHAR, BETAR, ALPHAI, BETAI, RHOR, FREQ, RHOHS,
     &   BUMDEN, ETA, XI, ZTMP
      COMPLEX*16 CPHS, CSHS, CRCI
      CHARACTER BCTYPE*1, ATUNIT*1

      COMMON /CPREV/ ALPHAR, BETAR, RHOR, ALPHAI, BETAI

C     *** Echo to PRTFIL user's choice of boundary condition ***

      IF ( BCTYPE(1:1) .EQ. 'S' ) THEN
         WRITE( PRTFIL, * ) '    Twersky SOFT BOSS scatter model'

      ELSE IF ( BCTYPE(1:1) .EQ. 'H' ) THEN 
         WRITE( PRTFIL, * ) '    Twersky HARD BOSS scatter model'

      ELSE IF ( BCTYPE(1:1) .EQ. 'T' ) THEN 
         WRITE( PRTFIL, * ) '    Twersky (amplitude only) SOFT BOSS ', 
     &                                       'scatter model'

      ELSE IF ( BCTYPE(1:1) .EQ. 'I' ) THEN 
         WRITE( PRTFIL, * ) '    Twersky (amplitude only) HARD BOSS ', 
     &                                       'scatter model'

      ELSE IF ( BCTYPE(1:1) .EQ. 'V' ) THEN 
         WRITE( PRTFIL, * ) '    VACUUM'

      ELSE IF ( BCTYPE(1:1) .EQ. 'R' ) THEN 
         WRITE( PRTFIL, * ) '    Perfectly RIGID'

      ELSE IF ( BCTYPE(1:1) .EQ. 'A' ) THEN
         WRITE( PRTFIL, * ) '    ACOUSTO-ELASTIC half-space'

      ELSE IF ( BCTYPE(1:1) .EQ. 'F' ) THEN 
         WRITE( PRTFIL, * ) '    FILE used for reflection loss'

      ELSE IF ( BCTYPE(1:1) .EQ. 'W' ) THEN 
         WRITE( PRTFIL, * ) '    Writing an IRC file'

      ELSE IF ( BCTYPE(1:1) .EQ. 'P' ) THEN
         WRITE( PRTFIL, * ) '    reading PRECALCULATED IRC'

      ELSE
         CALL ERROUT( PRTFIL, 'F', 'TOPBOT',
     &                'Unknown boundary condition type' )
      ENDIF

C     ****** Read in BC parameters depending on particular choice ******

      CPHS  = 0.0
      CSHS  = 0.0
      RHOHS = 0.0

C     *** Twersky ice model parameters ***

      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

         READ(  ENVFIL, *    ) BUMDEN, ETA, XI
         WRITE( PRTFIL, 1000 ) BUMDEN, ETA, XI
 1000    FORMAT( /, ' Twersky ice model parameters:', /, 
     &   ' Bumden = ', G15.6, '  Eta = ', G11.3, '  Xi = ', G11.3, /)
      ENDIF

C     *** Half-space properties ***

      IF ( BCTYPE(1:1) .EQ. 'A' ) THEN

         READ(  ENVFIL, *    ) ZTMP, ALPHAR, BETAR, RHOR, ALPHAI, BETAI
         WRITE( PRTFIL, 2000 ) ZTMP, ALPHAR, BETAR, RHOR, ALPHAI, BETAI
 2000    FORMAT( F10.2, 3X, 2F10.2, 3X, F6.2, 3X, 2F10.4 )

         CPHS = CRCI( ALPHAR, ALPHAI, FREQ, ATUNIT )
         CSHS = CRCI( BETAR,  BETAI,  FREQ, ATUNIT )
         RHOHS = RHOR
      ENDIF

      RETURN
      END
C**********************************************************************C
      FUNCTION CRCI( C, ALPHA, FREQ, ATUNIT )

C     Converts real wave speed and attenuation to a single
C     complex wave speed

C     6 CASES:    N for Nepers/meter
C                 M for dB/meter      (M for Meters)
C                 F for dB/m-kHZ      (F for Frequency dependent)
C                 W for dB/wavelength (W for Wavelength)
C                 Q for Q
C                 T for Thorp

      IMPLICIT REAL*8 (A-H, O-Z)
      COMPLEX*16 CRCI
      CHARACTER ATUNIT*1

      OMEGA = 2.0 * 3.1415926535897932 * FREQ

C     *** Convert to Nepers/m ***

      ALPHAT = 0.0

      IF      ( ATUNIT(1:1) .EQ. 'N' ) THEN
         ALPHAT = ALPHA
      ELSE IF ( ATUNIT(1:1) .EQ. 'M' ) THEN
         ALPHAT = ALPHA / 8.6858896
      ELSE IF ( ATUNIT(1:1) .EQ. 'F' ) THEN
         ALPHAT = ALPHA * FREQ / 8685.8896
      ELSE IF ( ATUNIT(1:1) .EQ. 'W' ) THEN
         IF ( C .NE. 0.0 ) ALPHAT = ALPHA * FREQ / ( 8.6858896 * C )
      ELSE IF ( ATUNIT(1:1) .EQ. 'Q' ) THEN
         IF( C * ALPHA .NE. 0.0 ) ALPHAT = OMEGA / ( 2.0 * C * ALPHA )
      ELSE IF ( ATUNIT(1:1) .EQ. 'T' ) THEN
         F2     = ( FREQ / 1000.0 ) **2
         ALPHAT = 40.0 * F2 / ( 4100.0 + F2 ) + 0.1*F2 / ( 1.0 + F2 )
         ALPHAT = ALPHAT / 914.4                 ! dB / m
         ALPHAT = ALPHAT / 8.6858896             ! Nepers / m
      ENDIF

C     *** Convert Nepers/m to equivalent imaginary sound speed ***

      ALPHAT = ALPHAT * C**2 / OMEGA
      CRCI   = CMPLX( C, ALPHAT )

      RETURN
      END
