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

C     Reads in the info in ENVFIL

      IMPLICIT NONE

      INTEGER MAXSSP

      PARAMETER ( MAXSSP = 2001 )

      INTEGER ENVFIL, PRTFIL, NMEDIA, NG( * ), MAXMED,
     &   NELTS, MED
      REAL*8 ALPHAR, BETAR, ALPHAI, BETAI, RHOR, FREQ, RHOT, RHOB,
     &   SIGMA( * ), DEPTH( * ),
     &   BUMDEN, ETA, XI, RHO( 1 ), C, DELTAZ
      COMPLEX*16 CPT, CST, CPB, CSB, CP( MAXSSP ), CS( MAXSSP )
      CHARACTER TOPOPT*( * ), BOTOPT*( * ), TITLE*( * ),
     &   SSPTYP*1, ATUNIT*1, BCTYPE*1

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

C     OPEN ( FILE = 'ENVFIL', UNIT = ENVFIL, STATUS = 'OLD',
C    &       FORM = 'FORMATTED' )

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

C     *** Title, Freq, Nmedia ***

      READ(  ENVFIL, * ) TITLE(9:80)
      WRITE( PRTFIL, * ) TITLE

      READ( ENVFIL, * ) FREQ
      READ( ENVFIL, * ) NMEDIA

      WRITE( PRTFIL, 1000 ) FREQ, NMEDIA
 1000 FORMAT( ' Frequency = ', G11.4, ' NMEDIA = ', I3, // )

      IF ( NMEDIA .GT. MAXMED ) THEN
         WRITE( PRTFIL, * ) 'MAXMED = ', MAXMED
         CALL ERROUT( PRTFIL, 'F', 'READIN', 'Too many media' )
      ENDIF

C     ******** TOP OPTIONS ********

      READ( ENVFIL, * ) TOPOPT(1:8)
C      CALL STR$UPCASE( TOPOPT, TOPOPT )

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

C     *** SSP approximation options ***

      IF ( TOPOPT(1:1) .EQ. 'N' ) THEN 
         WRITE( PRTFIL, * ) '    N2-LINEAR approximation to SSP'

      ELSE IF ( TOPOPT(1:1) .EQ. 'C' ) THEN 
         WRITE( PRTFIL, * ) '    C-LINEAR approximation to SSP'

      ELSE IF ( TOPOPT(1:1) .EQ. 'S' ) THEN 
         WRITE( PRTFIL, * ) '    SPLINE approximation to SSP'

      ELSE IF ( TOPOPT(1:1) .EQ. 'A' ) THEN 
         WRITE( PRTFIL, * ) '    ANALYTIC SSP option'

      ELSE
         CALL ERROUT( PRTFIL, 'F', 'READIN',
     &                'Unknown option for SSP approximation' )
      ENDIF

C     *** Attenuation options ***

      IF ( TOPOPT(3:3) .EQ. 'N' ) THEN 
         WRITE( PRTFIL, * ) '    Attenuation units: nepers/m'

      ELSE IF ( TOPOPT(3:3) .EQ. 'F' ) THEN 
         WRITE( PRTFIL, * ) '    Attenuation units: dB/mkHz'

      ELSE IF ( TOPOPT(3:3) .EQ. 'M' ) THEN 
         WRITE( PRTFIL, * ) '    Attenuation units: dB/m'

      ELSE IF ( TOPOPT(3:3) .EQ. 'W' ) THEN 
         WRITE( PRTFIL, * ) '    Attenuation units: dB/wavelength'

      ELSE IF ( TOPOPT(3:3) .EQ. 'Q' ) THEN 
         WRITE( PRTFIL, * ) '    Attenuation units: Q'

      ELSE IF ( TOPOPT(3:3) .EQ. 'T' ) THEN 
         WRITE( PRTFIL, * ) '    THORP attenuation added'

      ELSE
         CALL ERROUT( PRTFIL, 'F', 'READIN',
     &                'Unknown attenuation units' )
      ENDIF

C     *** CALL TOPBOT  to read top BC ***

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

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

C     ****** Internal media ******

      IF ( BCTYPE(1:1) .NE. 'A' ) WRITE( PRTFIL, 2000 )

      DO MED = 1, NMEDIA
         READ(  ENVFIL, *    ) NG( MED ), SIGMA( MED ), DEPTH( MED + 1 )
         WRITE( PRTFIL, 3000 ) NG( MED ), SIGMA( MED )
 3000    FORMAT( /, '          ( Number of pts = ', I5, 
     &                  '  RMS roughness = ', G10.3, ' )')

C        *** Call PROFIL to read in SSP ***

         CALL PROFIL( DEPTH, CP, CS, RHO, MED, NELTS,
     &      FREQ, SSPTYP, ATUNIT, 'INIT', ENVFIL, PRTFIL  )

C        *** automatic calculation of f.d. mesh (10 pts/wavelength) ***
         IF ( NG( MED ) .EQ. 0 ) THEN
            C = alphar
            IF ( betar .GT. 0.0 ) C = betar ! shear?

            DELTAZ = 0.02 * C / FREQ
            NG( MED ) = ( DEPTH( MED + 1 ) - DEPTH( MED ) ) / DELTAZ
            WRITE( PRTFIL, * ) 'Number of pts = ', NG( MED )
         ENDIF

      END DO   ! next MED

C     ****** Bottom properties ******

      READ( ENVFIL, * ) BOTOPT(1:8), SIGMA( NMEDIA + 1 )
C      CALL STR$UPCASE( BOTOPT, BOTOPT )
      BCTYPE(1:1) = BOTOPT(1:1)

      WRITE( PRTFIL, * )
      WRITE( PRTFIL, 5000 ) SIGMA( NMEDIA + 1 )
 5000 FORMAT( 33X, '( RMS roughness = ', G10.3, ' )' )

C     *** CALL TOPBOT  to read bottom BC ***

      CALL TOPBOT( ENVFIL, PRTFIL, FREQ, BCTYPE, ATUNIT,
     &    CPB, CSB, RHOB, BUMDEN, ETA, XI )

      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
      SUBROUTINE PROFIL( DEPTH, CP, CS, RHOT, MED, N1,
     &   FREQ, SSPTYP, ATUNIT, TASK, ENVFIL, PRTFIL  )

C     Call the particular SSP routine specified by SSPTYP

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

C     Note that FREQ is only need if TASK = 'INIT'

      IMPLICIT REAL*8 (A-H, O-Z)
      INTEGER ENVFIL, PRTFIL

      REAL*8     RHOT( * ), DEPTH( * )
      COMPLEX*8  CPT, CST
      COMPLEX*16 CP( * ), CS( * )
      CHARACTER  SSPTYP*1, ATUNIT*1, TASK*8

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

C        *** Analytic profile option ***
         IF ( TASK(1:4) .EQ. 'INIT' ) THEN
            N1 = 21
            CALL ANALYT( DEPTH, CP, CS, RHOT,
     &                                     MED, N1, FREQ, ATUNIT, TASK )
            H = ( DEPTH( MED+1 ) - DEPTH( MED ) ) / ( N1 - 1 )

            DO I = 1, N1
               Z = DEPTH( MED ) + ( I - 1 ) * H
C              --- for f77/f90 compatibility in AIMAG/DIMAG:
               CPT =  CP( I )
               CST =  CS( I )
               WRITE( PRTFIL, 2000 ) Z,  REAL( CPT ),  REAL( CST ),
     &                         RHOT(I), AIMAG( CPT ), AIMAG( CST )
 2000          FORMAT( F10.2, 3X, 2F10.2, 3X, F6.2, 3X, 2F10.4 )
            END DO
         ELSE
            CALL ANALYT( DEPTH, CP, CS, RHOT,
     &                                     MED, N1, FREQ, ATUNIT, TASK )
         ENDIF

      ELSEIF ( SSPTYP(1:1) .EQ. 'N' ) THEN
C        *** N2-linear profile option ***
         CALL N2LIN(  DEPTH, CP, CS, RHOT, MED, N1,
     &                FREQ, ATUNIT, TASK, ENVFIL, PRTFIL )

      ELSEIF ( SSPTYP(1:1) .EQ. 'C' ) THEN
C        *** C-linear profile option ***
         CALL CLIN(   DEPTH, CP, CS, RHOT, MED, N1,
     &                FREQ, ATUNIT, TASK, ENVFIL, PRTFIL )

      ELSEIF ( SSPTYP(1:1) .EQ. 'S' ) THEN
C        *** Cubic spline profile option ***
         CALL CCUBIC( DEPTH, CP, CS, RHOT, MED, N1,
     &                FREQ, ATUNIT, TASK, ENVFIL, PRTFIL )

      ELSE
C        *** Non-existent profile option ***
         WRITE( PRTFIL, * ) 'Profile option: ', SSPTYP(1:1)
         CALL ERROUT( PRTFIL, 'F', 'PROFIL',
     &                'Unknown profile option' )

      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 (with positive imaginary part)

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.1415926535897932D0 * 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.6858896D0

      ELSE IF ( ATUNIT(1:1) .EQ. 'F' ) THEN
         ALPHAT = ALPHA * FREQ / 8685.8896D0

      ELSE IF ( ATUNIT(1:1) .EQ. 'W' ) THEN
         IF ( C .NE. 0.0 ) ALPHAT = ALPHA * FREQ / ( 8.6858896D0 * C )

C        --- the following lines give f^1.25 frequency dependence
C        FAC = SQRT( SQRT( FREQ / 50.0 ) )
C        IF ( C .NE. 0.0 )
C     &      ALPHAT = FAC * ALPHA * FREQ / ( 8.6858896D0 * 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.4D0                 ! dB / m
         ALPHAT = ALPHAT / 8.6858896D0             ! Nepers / m

      ENDIF

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

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

      RETURN
      END
C**********************************************************************C
      SUBROUTINE N2LIN( DEPTH, CP, CS, RHOT,
     &                  MED, N1, FREQ, ATUNIT, TASK, ENVFIL, PRTFIL )

C     Tabulate CP, CS, RHO for specified MED

C     Uses N2-linear segments for P and S-wave speeds
C     Uses rho-linear segments for density

      IMPLICIT REAL*8 (A-H, O-Z)
      INTEGER ENVFIL, PRTFIL
      PARAMETER ( MAXMED = 51, MAXSSP = 2001, PREC = 1.0E-6 )

      INTEGER LOC( MAXMED ), NPTS( MAXMED )
      REAL*8 DEPTH( * ), RHOT( * ), Z( MAXSSP ), RHO( MAXSSP )
      COMPLEX*16 CP( * ), CS( * ), ALPHA( MAXSSP ), BETA( MAXSSP ),
     &   N2BOT, N2TOP, CRCI
      CHARACTER ATUNIT*1, TASK*8

      SAVE Z, ALPHA, BETA, RHO, LOC, NPTS
      COMMON /CPREV/ ALPHAR, BETAR, RHOR, ALPHAI, BETAI

C     If TASK = 'INIT' then this is the first call and SSP is read.
C     Any other call is a request for SSP subtabulation.

      IF ( TASK(1:4) .EQ. 'INIT' ) THEN
C        ------ TASK 'INIT' for initialization

C        The variable LOC( MED ) points to the starting point for the
C        data in the arrays Z, ALPHA, BETA and RHO

         IF ( MED .EQ. 1 ) THEN
            LOC( MED ) = 0
         ELSE
            LOC( MED ) = LOC( MED - 1 ) + NPTS( MED - 1 )
         ENDIF
         ILOC = LOC( MED )

C        *** Read in data and convert attenuation to Nepers/m ***

         N1 = 1
         DO I = 1, MAXSSP
            READ( ENVFIL, *)      Z( ILOC + I ), ALPHAR, BETAR,
     &         RHOR, ALPHAI, BETAI
            WRITE( PRTFIL, 2000 ) Z( ILOC + I ), ALPHAR, BETAR,
     &         RHOR, ALPHAI, BETAI
 2000       FORMAT(F10.2, 3X, 2F10.2, 3X, F6.2, 3X, 2F10.4)

            ALPHA(ILOC + I) = CRCI( ALPHAR, ALPHAI, FREQ, ATUNIT )
            BETA( ILOC + I) = CRCI( BETAR,  BETAI,  FREQ, ATUNIT )
            RHO(  ILOC + I) = RHOR

C           ------ Did we read the last point?
            IF ( ABS( Z( ILOC + I ) - DEPTH( MED+1 ) ) .LT. PREC ) THEN
               NPTS( MED ) = N1
               IF ( MED .EQ. 1 ) DEPTH( 1 ) = Z( 1 )
               RETURN
            ENDIF

            N1 = N1 + 1
         END DO

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

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

      ELSE
C        ------ TASK = 'TABULATE'
         ILOC = LOC( MED )
         N    = N1 - 1
         H    = ( Z( ILOC + NPTS( MED ) ) - Z( ILOC + 1 ) ) / N
         LAY  = 1

         DO I = 1, N1
            ZT = Z( ILOC + 1 ) + ( I - 1 ) * H
C           ------ Make sure no overshoot
            IF ( I .EQ. N1 ) ZT = Z( ILOC + NPTS( MED ) )

 7000       IF ( ZT .GT. Z( ILOC + LAY + 1 ) ) THEN
               LAY = LAY + 1
               GO TO 7000
            ENDIF

            R = ( ZT - Z( ILOC + LAY ) ) / 
     &         ( Z( ILOC + LAY+1 ) - Z( ILOC + LAY ) )

C           ------ P-wave
            N2TOP   = 1.0 / ALPHA( ILOC + LAY     )**2
            N2BOT   = 1.0 / ALPHA( ILOC + LAY + 1 )**2
            CP( I ) = 1.0 / SQRT( ( 1.0 - R ) * N2TOP + R * N2BOT )

C           ------ S-wave
            IF ( BETA(ILOC + LAY) .NE. 0.0 ) THEN
               N2TOP   = 1.0 / BETA( ILOC + LAY     )**2
               N2BOT   = 1.0 / BETA( ILOC + LAY + 1 )**2
               CS( I ) = 1.0 / SQRT( ( 1.0 - R ) * N2TOP + R * N2BOT )
            ELSE
               CS( I ) = 0.0
            ENDIF

            RHOT( I ) = ( 1.0 - R ) * RHO( ILOC + LAY ) +
     &                            R * RHO( ILOC + LAY + 1 )
         END DO

      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE CLIN( DEPTH, CP, CS, RHOT,
     &                 MED, N1, FREQ, ATUNIT, TASK, ENVFIL, PRTFIL  )

C     Tabulate CP, CS, RHO for specified MED

C     Uses c-linear segments for P and S-wave speeds
C     Uses rho-linear segments for density

      IMPLICIT REAL*8 (A-H, O-Z)
      INTEGER ENVFIL, PRTFIL
      PARAMETER ( MAXMED = 51, MAXSSP = 2001, PREC = 1.0E-6 )

      INTEGER LOC( MAXMED ), NPTS( MAXMED )
      REAL*8 DEPTH( * ), RHOT( * ), Z( MAXSSP ), RHO( MAXSSP )
      COMPLEX*16 CP( * ), CS( * ), ALPHA( MAXSSP ), BETA( MAXSSP ), CRCI
      CHARACTER ATUNIT*1, TASK*8

      SAVE Z, ALPHA, BETA, RHO, LOC, NPTS
      COMMON /CPREV/ ALPHAR, BETAR, RHOR, ALPHAI, BETAI

C     If TASK = 'INIT' then this is the first call and SSP is read.
C     Any other call is a request for SSP subtabulation.

      IF ( TASK(1:4) .EQ. 'INIT' ) THEN
C        ------ TASK 'INIT' FOR INITIALIZATION
         NPTS( MED ) = N1

C        The variable LOC(MED) points to the starting point for the
C        data in the arrays Z, ALPHA, BETA and RHO

         IF ( MED .EQ. 1 ) THEN
            LOC( MED ) = 0
         ELSE
            LOC( MED ) = LOC( MED - 1 ) + NPTS( MED - 1 )
         ENDIF
         ILOC = LOC( MED )

C        *** Read in data and convert attenuation to Nepers/m ***

         N1 = 1
         DO I = 1, MAXSSP
            READ( ENVFIL, * )     Z( ILOC + I ), ALPHAR, BETAR,
     &         RHOR, ALPHAI, BETAI
            WRITE( PRTFIL, 2000 ) Z( ILOC + I ), ALPHAR, BETAR,
     &         RHOR, ALPHAI, BETAI
 2000       FORMAT( F10.2, 3X, 2F10.2, 3X, F6.2, 3X, 2F10.4 )

            ALPHA(ILOC + I) = CRCI( ALPHAR, ALPHAI, FREQ, ATUNIT )
            BETA( ILOC + I) = CRCI( BETAR,  BETAI,  FREQ, ATUNIT )
            RHO(  ILOC + I) = RHOR

C           ------ Did we read the last point?
            IF ( ABS( Z( ILOC + I ) - DEPTH( MED+1 ) ) .LT. PREC ) THEN
               NPTS( MED ) = N1
               IF ( MED .EQ. 1 ) DEPTH( 1 ) = Z( 1 )
               RETURN
            ENDIF

            N1 = N1 + 1
         END DO

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

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

      ELSE
C        ------ TASK = 'TABULATE'
         ILOC = LOC( MED )
         N    = N1 - 1
         H    = ( Z( ILOC + NPTS( MED ) ) - Z( ILOC + 1 ) ) / N
         LAY  = 1

         DO I = 1, N1
            ZT = Z( ILOC + 1 ) + ( I - 1 ) * H
C           ------ Make sure no overshoot
            IF ( I .EQ. N1 ) ZT = Z( ILOC + NPTS( MED ) )

 7000       IF ( ZT .GT. Z( ILOC + LAY + 1 ) ) THEN
               LAY = LAY + 1
               GO TO 7000
            ENDIF

            R = ( ZT - Z( ILOC + LAY ) ) /
     &         ( Z( ILOC + LAY + 1 ) - Z( ILOC + LAY ) )
            CP(I)   = (1.0-R)*ALPHA(ILOC + LAY) + R*ALPHA(ILOC + LAY+1)
            CS(I)   = (1.0-R)* BETA(ILOC + LAY) + R* BETA(ILOC + LAY+1)
            RHOT(I) = (1.0-R)*  RHO(ILOC + LAY) + R*  RHO(ILOC + LAY+1)
         END DO
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE CCUBIC( DEPTH, CP, CS, RHOT,
     &                   MED, N1, FREQ, ATUNIT, TASK, ENVFIL, PRTFIL  )

C     Tabulate CP, CS, RHO for specified MED
C     using cubic spline interpolation

      IMPLICIT REAL*8 (A-H, O-Z)
      INTEGER ENVFIL, PRTFIL
      PARAMETER ( MAXMED = 51, MAXSSP = 2001, PREC = 1.0E-6 )

      INTEGER LOC( MAXMED ), NPTS( MAXMED )
      REAL*8 DEPTH( * ), Z( MAXSSP ), RHOT( * )
      COMPLEX*16 CP( * ), CS( * ), ESPLINE, CRCI,
     &    ALPHA( 4, MAXSSP ), BETA( 4, MAXSSP ), RHO( 4, MAXSSP )
      CHARACTER ATUNIT*1, TASK*8

      SAVE Z, ALPHA, BETA, RHO, LOC, NPTS
      COMMON /CPREV/ ALPHAR, BETAR, RHOR, ALPHAI, BETAI

C     If TASK = 'INIT' then this is the first call and SSP is read.
C     Any other call is a request for SSP subtabulation.

      IF ( TASK(1:4) .EQ. 'INIT' ) THEN
C        ------ Task 'INIT' for initialization
         NPTS( MED ) = N1

C        The variable LOC(MED) points to the starting point for the
C        data in the arrays Z, ALPHA, BETA and RHO

         IF ( MED .EQ. 1 ) THEN
            LOC( MED ) = 0
         ELSE
            LOC( MED ) = LOC( MED - 1 ) + NPTS( MED - 1 )
         ENDIF
         ILOC = LOC( MED )

C        *** Read in data and convert attenuation to Nepers/m ***

         N1 = 1

         DO I = 1, MAXSSP
            READ( ENVFIL, *)      Z( ILOC + I ), ALPHAR, BETAR,
     &         RHOR, ALPHAI, BETAI
            WRITE( PRTFIL, 2000 ) Z( ILOC + I ), ALPHAR, BETAR,
     &         RHOR, ALPHAI, BETAI
 2000       FORMAT(F10.2, 3X, 2F10.2, 3X, F6.2, 3X, 2F10.4)

            ALPHA(1, ILOC + I) = CRCI( ALPHAR, ALPHAI, FREQ, ATUNIT )
            BETA( 1, ILOC + I) = CRCI( BETAR,  BETAI,  FREQ, ATUNIT )
            RHO(  1, ILOC + I) = RHOR

C           ------ Did we read the last point?
            IF ( ABS( Z( ILOC + I ) - DEPTH(MED+1) ) .LT. PREC ) THEN
               NPTS( MED ) = N1
               IF ( MED .EQ. 1 ) DEPTH( 1 ) = Z( 1 )

C              *** Compute spline coefs ***

               IBCBEG = 0
               IBCEND = 0
               CALL CSPLINE( Z( ILOC + 1 ), ALPHA(1, ILOC + 1),
     &                       NPTS( MED ), IBCBEG, IBCEND, NPTS( MED ) )
               CALL CSPLINE( Z( ILOC + 1 ),  BETA( 1, ILOC + 1 ),
     &                       NPTS( MED ), IBCBEG, IBCEND, NPTS( MED ) )
               CALL CSPLINE( Z( ILOC + 1 ),   RHO( 1, ILOC + 1 ),
     &                       NPTS( MED ), IBCBEG, IBCEND, NPTS( MED ) )

               RETURN
            ENDIF

            N1 = N1 + 1
         END DO

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

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

      ELSE
C        ------ Task = 'TABULATE'
         ILOC = LOC( MED )
         N    = N1 - 1
         H    = ( Z( ILOC + NPTS( MED ) ) - Z( ILOC + 1 ) ) / N
         LAY  = 1
         DO I = 1, N1
            ZT = Z( ILOC + 1 ) + ( I - 1 ) * H
C           ------ Make sure no overshoot
            IF ( I .EQ. N1 ) ZT = Z( ILOC + NPTS( MED ) )
 7000       IF ( ZT .GT. Z( ILOC + LAY + 1 ) ) THEN
               LAY = LAY + 1
               GO TO 7000
            ENDIF

            HSPLNE = ZT - Z( ILOC + LAY )

            CP(   I ) =        ESPLINE( ALPHA(1, ILOC + LAY), HSPLNE )
            CS(   I ) =        ESPLINE(  BETA(1, ILOC + LAY), HSPLNE )
            RHOT( I ) = DBLE( ESPLINE(   RHO(1, ILOC + LAY), HSPLNE ) )

         END DO
      ENDIF

      RETURN
      END
