!! Environment input file parsing and setup for BELLHOP

MODULE ReadEnvironmentBell
  !! Provides environment file reading and initialization

  USE BellhopMod
  USE MathConstants
  USE sspmod
  USE AttenMod
  USE FatalError
  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: ERROR_UNIT

  IMPLICIT NONE
  PUBLIC

CONTAINS

  SUBROUTINE ReadEnvironment( FileRoot, ThreeD )
    !! Reads and parses the main environment file

    ! Routine to read in and echo all the input data
    ! Note that default values of SSP, DENSITY, Attenuation will not work

    USE anglemod
    USE SourceReceiverPositions

    REAL      (KIND=8), PARAMETER   :: c0 = 1500.0
    LOGICAL,            INTENT(IN ) :: ThreeD
    CHARACTER (LEN=80), INTENT(IN ) :: FileRoot
    INTEGER            :: NPts, NMedia, iostat
    REAL               :: ZMin, ZMax
    REAL      (KIND=8) :: x( 2 ), t( 2 ), c, cimag, gradc( 2 ), crr, crz, czz, rho, sigma, Depth
    CHARACTER (LEN= 2) :: AttenUnit
    CHARACTER (LEN=10) :: PlotType

    WRITE( PRTFile, * ) 'BELLHOP/BELLHOP3D'
    WRITE( PRTFile, * )

    ! Open the environmental file
    OPEN( UNIT = ENVFile, FILE = TRIM( FileRoot ) // '.env', STATUS = 'OLD', IOSTAT = iostat, ACTION = 'READ' )
    IF ( IOSTAT /= 0 ) THEN   ! successful open?
       WRITE( PRTFile, * ) 'ENVFile = ', TRIM( FileRoot ) // '.env'
       CALL ERROUT( 'BELLHOP - READIN', 'Unable to open the environmental file' )
    END IF

    ! Prepend model name to title
    IF ( ThreeD ) THEN
       Title( 1 : 11 ) = 'BELLHOP3D- '
       READ(  ENVFile, * ) Title( 12 : 80 )
    ELSE
       Title( 1 :  9 ) = 'BELLHOP- '
       READ(  ENVFile, * ) Title( 10 : 80 )
    END IF

    WRITE( PRTFile, * ) Title

    READ(  ENVFile, *    ) freq
    WRITE( PRTFile, '('' frequency = '', G11.4, '' Hz'', / )' ) freq

    READ(  ENVFile, * ) NMedia
    WRITE( PRTFile, * ) 'Dummy parameter NMedia = ', NMedia
    IF ( NMedia /= 1 ) CALL ERROUT( 'READIN', &
         'Only one medium or layer is allowed in BELLHOP; sediment layers must be handled using a reflection coefficient' )

    CALL ReadTopOpt( Bdry%Top%HS%Opt, Bdry%Top%HS%BC, AttenUnit, FileRoot )

    ! *** Top BC ***

    IF ( Bdry%Top%HS%BC == 'A' ) THEN
       WRITE( PRTFile, "( //, '      z         alphaR      betaR     rho        alphaI     betaI', / )" )
       WRITE( PRTFile, "(     '     (m)         (m/s)      (m/s)   (g/cm^3)      (m/s)     (m/s)', / )" )
    END IF

    CALL TopBot( freq, AttenUnit, Bdry%Top%HS )

    ! ****** Read in ocean SSP data ******

    READ(  ENVFile, * ) NPts, Sigma, Bdry%Bot%HS%Depth
    WRITE( PRTFile, * )
    WRITE( PRTFile, FMT = "( ' Depth = ', F10.2, ' m' )" ) Bdry%Bot%HS%Depth

    IF ( Bdry%Top%HS%Opt( 1 : 1 ) == 'A' ) THEN
       WRITE( PRTFile, * ) 'Analytic SSP option'
       ! following is hokey, should be set in Analytic routine
       SSP%NPts = 2
       SSP%z( 1 ) = 0.0
       SSP%z( 2 ) = Bdry%Bot%HS%Depth
    ELSE
       x = [ 0.0D0, Bdry%Bot%HS%Depth ]   ! tells SSP Depth to read to
       t = [ 0.0, 1.0 ]
       CALL EvaluateSSP( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, 'INI' )
    END IF

    Bdry%Top%HS%Depth = SSP%z( 1 )   ! Depth of top boundary is taken from first SSP point
    ! bottom depth should perhaps be set the same way?

    ! *** Bottom BC ***
    Bdry%Bot%HS%Opt = '  '   ! initialize to blanks
    READ(  ENVFile, * ) Bdry%Bot%HS%Opt, Sigma
    WRITE( PRTFile, * )
    WRITE( PRTFile, FMT = "(33X, '( RMS roughness = ', G10.3, ' )' )" ) Sigma

    SELECT CASE ( Bdry%Bot%HS%Opt( 2 : 2 ) )
    CASE ( '~', '*' )
       WRITE( PRTFile, * ) '    Bathymetry file selected'
    CASE( '-', '_', ' ' )
    CASE DEFAULT
       CALL ERROUT( 'READIN', 'Unknown bottom option letter in second position' )
    END SELECT

    Bdry%Bot%HS%BC = Bdry%Bot%HS%Opt( 1 : 1 )
    CALL TopBot( freq, AttenUnit, Bdry%Bot%HS )

    ! *** source and receiver locations ***

    CALL ReadSxSy( ThreeD )     ! Read source/receiver x-y coordinates

    ZMin = SNGL( Bdry%Top%HS%Depth )
    ZMax = SNGL( Bdry%Bot%HS%Depth )
    ! CALL ReadSzRz( ZMin + 100 * SPACING( ZMin ), ZMax - 100 * SPACING( ZMax ) )   ! not sure why I had this
    CALL ReadSzRz( ZMin, ZMax )
    CALL ReadRcvrRanges
    IF ( ThreeD ) CALL ReadRcvrBearings
    CALL ReadfreqVec( freq,  Bdry%Top%HS%Opt( 6:6 ) )
    CALL ReadRunType( Beam%RunType, PlotType )

    Depth = Zmax - Zmin   ! water depth
    CALL ReadRayElevationAngles( freq, Depth, Bdry%Top%HS%Opt, Beam%RunType )
    IF ( ThreeD ) CALL ReadRayBearingAngles(   freq, Bdry%Top%HS%Opt, Beam%RunType )

    WRITE( PRTFile, * )
    WRITE( PRTFile, * ) '__________________________________________________________________________'
    WRITE( PRTFile, * )

    ! Limits for tracing beams

    IF ( ThreeD ) THEN
       READ(  ENVFile, * ) Beam%deltas, Beam%Box%x, Beam%Box%y, Beam%Box%z
       Beam%Box%x = 1000.0 * Beam%Box%x   ! convert km to m
       Beam%Box%y = 1000.0 * Beam%Box%y   ! convert km to m

       IF ( Beam%deltas == 0.0 ) Beam%deltas = ( Bdry%Bot%HS%Depth - Bdry%Top%HS%Depth ) / 10.0   ! Automatic step size selection
       ! WRITE( PRTFile, '('' frequency = '', G11.4, '' Hz'', / )' ) freq

       WRITE( PRTFile, * )
       WRITE( PRTFile, fmt = '(  '' Step length,       deltas = '', G11.4, '' m'' )' ) Beam%deltas
       WRITE( PRTFile, * )
       WRITE( PRTFile, fmt = '(  '' Maximum ray x-range, Box%x = '', G11.4, '' m'' )' ) Beam%Box%x
       WRITE( PRTFile, fmt = '(  '' Maximum ray y-range, Box%y = '', G11.4, '' m'' )' ) Beam%Box%y
       WRITE( PRTFile, fmt = '(  '' Maximum ray z-range, Box%z = '', G11.4, '' m'' )' ) Beam%Box%z
    ELSE
       READ(  ENVFile, * ) Beam%deltas, Beam%Box%z, Beam%Box%r

       WRITE( PRTFile, * )
       WRITE( PRTFile, fmt = '(  '' Step length,       deltas = '', G11.4, '' m'' )' ) Beam%deltas
       WRITE( PRTFile, * )
       WRITE( PRTFile, fmt = '(  '' Maximum ray depth, Box%z  = '', G11.4, '' m'' )' ) Beam%Box%z
       WRITE( PRTFile, fmt = '(  '' Maximum ray range, Box%r  = '', G11.4, ''km'' )' ) Beam%Box%r

       Beam%Box%r = 1000.0 * Beam%Box%r   ! convert km to m
    END IF

    ! *** Beam characteristics ***
    Beam%Type( 4 : 4 ) = Beam%RunType( 7 : 7 )   ! selects beam shift option

    SELECT CASE ( Beam%Type( 4 : 4 ) )
    CASE ( 'S' )
          WRITE( PRTFile, * ) 'Beam shift in effect'
    CASE DEFAULT
          WRITE( PRTFile, * ) 'No beam shift in effect'
    END SELECT

    IF ( Beam%RunType( 1 : 1 ) /= 'R' ) THEN   ! no worry about the beam type if this is a ray trace run

       ! Beam%Type( 1 : 1 ) is
       !   'G' or '^' Geometric hat beams in Cartesian coordinates
       !   'g' Geometric hat beams in ray-centered coordinates
       !   'B' Geometric Gaussian beams in Cartesian coordinates
       !   'b' Geometric Gaussian beams in ray-centered coordinates
       !   'S' Simple Gaussian beams
       !   'C' Cerveny Gaussian beams in Cartesian coordinates
       !   'R' Cerveny Gaussian beams in Ray-centered coordinates
       ! Beam%Type( 2 : 2 ) controls the setting of the beam width
       !   'F' space Filling
       !   'M' minimum width
       !   'W' WKB beams
       ! Beam%Type( 3 : 3 ) controls curvature changes on boundary reflections
       !   'D' Double
       !   'S' Single
       !   'Z' Zero
       ! Beam%Type( 4 : 4 ) selects whether beam shifts are implemented on boundary reflection
       !   'S' yes
       !   'N' no

       ! Curvature change can cause overflow in grazing case
       ! Suppress by setting BeamType( 3 : 3 ) = 'Z'

       Beam%Type( 1 : 1 ) = Beam%RunType( 2 : 2 )
       SELECT CASE ( Beam%Type( 1 : 1 ) )
       CASE ( 'G', 'g' , '^', 'B', 'b', 'S' )   ! geometric hat beams, geometric Gaussian beams, or simple Gaussian beams
       CASE ( 'R', 'C' )   ! Cerveny Gaussian Beams; read extra lines to specify the beam options
          READ(  ENVFile, * ) Beam%Type( 2 : 3 ), Beam%epsMultiplier, Beam%rLoop
          WRITE( PRTFile, * )
          WRITE( PRTFile, * )
          WRITE( PRTFile, * ) 'Type of beam = ', Beam%Type( 1 : 1 )

          SELECT CASE ( Beam%Type( 3 : 3 ) )
          CASE ( 'D' )
             WRITE( PRTFile, * ) 'Curvature doubling invoked'
          CASE ( 'Z' )
             WRITE( PRTFile, * ) 'Curvature zeroing invoked'
          CASE ( 'S' )
             WRITE( PRTFile, * ) 'Standard curvature condition'
          CASE DEFAULT
             CALL ERROUT( 'READIN', 'Unknown curvature condition' )
          END SELECT

          WRITE( PRTFile, * ) 'Epsilon multiplier', Beam%epsMultiplier
          WRITE( PRTFile, * ) 'Range for choosing beam width', Beam%rLoop

          ! Images, windows
          ! LP: These values are not initialized if not written in the file,
          ! and Component is not always written in the test env files.
          Beam%Nimage = 1
          Beam%iBeamWindow = 4
          Beam%Component = 'P'
          READ(  ENVFile, * ) Beam%Nimage, Beam%iBeamWindow, Beam%Component
          WRITE( PRTFile, * )
          WRITE( PRTFile, * ) 'Number of images, Nimage  = ', Beam%Nimage
          WRITE( PRTFile, * ) 'Beam windowing parameter  = ', Beam%iBeamWindow
          WRITE( PRTFile, * ) 'Component                 = ', Beam%Component
       CASE DEFAULT
          CALL ERROUT( 'READIN', 'Unknown beam type (second letter of run type)' )
       END SELECT
    END IF

    WRITE( PRTFile, * )
    CLOSE( ENVFile )

  END SUBROUTINE ReadEnvironment

  ! **********************************************************************!

  SUBROUTINE ReadTopOpt( TopOpt, BC, AttenUnit, FileRoot )

    CHARACTER (LEN= 6), INTENT( OUT ) :: TopOpt
    CHARACTER (LEN= 1), INTENT( OUT ) :: BC                     ! Boundary condition type
    CHARACTER (LEN= 2), INTENT( OUT ) :: AttenUnit
    CHARACTER (LEN=80), INTENT( IN  ) :: FileRoot
    INTEGER            :: iostat

    TopOpt = '      '   ! initialize to blanks
    READ(  ENVFile, * ) TopOpt
    WRITE( PRTFile, * )

    SSP%Type  = TopOpt( 1 : 1 )
    BC        = TopOpt( 2 : 2 )
    AttenUnit = TopOpt( 3 : 4 )
    SSP%AttenUnit = AttenUnit

    ! SSP approximation options

    SELECT CASE ( SSP%Type )
    CASE ( 'N' )
       WRITE( PRTFile, * ) '    N2-linear approximation to SSP'
    CASE ( 'C' )
       WRITE( PRTFile, * ) '    C-linear approximation to SSP'
    CASE ( 'P' )
       WRITE( PRTFile, * ) '    PCHIP approximation to SSP'
    CASE ( 'S' )
       WRITE( PRTFile, * ) '    Spline approximation to SSP'
    CASE ( 'Q' )
       WRITE( PRTFile, * ) '    Quad approximation to SSP'
       OPEN ( FILE = TRIM( FileRoot ) // '.ssp', UNIT = SSPFile, FORM = 'FORMATTED', STATUS = 'OLD', IOSTAT = iostat )
       IF ( IOSTAT /= 0 ) THEN   ! successful open?
          WRITE( PRTFile, * ) 'SSPFile = ', TRIM( FileRoot ) // '.ssp'
          CALL ERROUT( 'BELLHOP - READIN', 'Unable to open the SSP file' )
       END IF
    CASE ( 'H' )
       WRITE( PRTFile, * ) '    Hexahedral approximation to SSP'
       OPEN ( FILE = TRIM( FileRoot ) // '.ssp', UNIT = SSPFile, FORM = 'FORMATTED', STATUS = 'OLD', IOSTAT = iostat )
       IF ( IOSTAT /= 0 ) THEN   ! successful open?
          WRITE( PRTFile, * ) 'SSPFile = ', TRIM( FileRoot ) // '.ssp'
          CALL ERROUT( 'BELLHOP - READIN', 'Unable to open the SSP file' )
       END IF
    CASE ( 'A' )
       WRITE( PRTFile, * ) '    Analytic SSP option'
    CASE DEFAULT
       CALL ERROUT( 'READIN', 'Unknown option for SSP approximation' )
    END SELECT

    ! Attenuation options

    SELECT CASE ( AttenUnit( 1 : 1 ) )
    CASE ( 'N' )
       WRITE( PRTFile, * ) '    Attenuation units: nepers/m'
    CASE ( 'F' )
       WRITE( PRTFile, * ) '    Attenuation units: dB/mkHz'
    CASE ( 'M' )
       WRITE( PRTFile, * ) '    Attenuation units: dB/m'
    CASE ( 'W' )
       WRITE( PRTFile, * ) '    Attenuation units: dB/wavelength'
    CASE ( 'Q' )
       WRITE( PRTFile, * ) '    Attenuation units: Q'
    CASE ( 'L' )
       WRITE( PRTFile, * ) '    Attenuation units: Loss parameter'
    CASE DEFAULT
       CALL ERROUT( 'READIN', 'Unknown attenuation units' )
    END SELECT

    ! optional addition of volume attenuation using standard formulas

    SELECT CASE ( AttenUnit( 2 : 2 ) )
    CASE ( 'T' )
       WRITE( PRTFile, * ) '    THORP volume attenuation added'
    CASE ( 'F' )
       WRITE( PRTFile, * ) '    Francois-Garrison volume attenuation added'
       READ(  ENVFile, * ) T, Salinity, pH, z_bar
       WRITE( PRTFile, "( ' T = ', G11.4, 'degrees   S = ', G11.4, ' psu   pH = ', G11.4, ' z_bar = ', G11.4, ' m' )" ) &
            T, Salinity, pH, z_bar
    CASE ( 'B' )
       WRITE( PRTFile, * ) '    Biological attenaution'
       READ( ENVFile, *  ) NBioLayers
       WRITE( PRTFile, * ) '      Number of Bio Layers = ', NBioLayers

       DO iBio = 1, NBioLayers
          READ( ENVFile, *  ) bio( iBio )%Z1, bio( iBio )%Z2, bio( iBio )%f0, bio( iBio )%Q, bio( iBio )%a0
          WRITE( PRTFile, * ) '      Top    of layer = ', bio( iBio )%Z1, ' m'
          WRITE( PRTFile, * ) '      Bottom of layer = ', bio( iBio )%Z2, ' m'
          WRITE( PRTFile, * ) '      Resonance frequency = ', bio( iBio )%f0, ' Hz'
          WRITE( PRTFile, * ) '      Q  = ', bio( iBio )%Q
          WRITE( PRTFile, * ) '      a0 = ', bio( iBio )%a0
       END DO
    CASE ( ' ' )
    CASE DEFAULT
       CALL ERROUT( 'READIN', 'Unknown top option letter in fourth position' )
    END SELECT

    SELECT CASE ( TopOpt( 5 : 5 ) )
    CASE ( '~', '*' )
       WRITE( PRTFile, * ) '    Altimetry file selected'
    CASE ( '-', '_', ' ' )
    CASE DEFAULT
       CALL ERROUT( 'READIN', 'Unknown top option letter in fifth position' )
    END SELECT

    SELECT CASE ( TopOpt( 6 : 6 ) )
    CASE ( 'I' )
       WRITE( PRTFile, * ) '    Development options enabled'
    CASE ( ' ' )
    CASE DEFAULT
       CALL ERROUT( 'READIN', 'Unknown top option letter in sixth position' )
    END SELECT

  END SUBROUTINE ReadTopOpt

  ! **********************************************************************!

  SUBROUTINE ReadRunType( RunType, PlotType )
    !! Reads and validates the run type parameters

    ! Read the RunType variable and echo with explanatory information to the print file

    USE SourceReceiverPositions

    CHARACTER (LEN= 7), INTENT( OUT ) :: RunType
    CHARACTER (LEN=10), INTENT( OUT ) :: PlotType

    READ(  ENVFile, * ) RunType
    WRITE( PRTFile, * )

    SELECT CASE ( RunType( 1 : 1 ) )
    CASE ( 'R' )
       WRITE( PRTFile, * ) 'Ray trace run'
    CASE ( 'E' )
       WRITE( PRTFile, * ) 'Eigenray trace run'
    CASE ( 'I' )
       WRITE( PRTFile, * ) 'Incoherent TL calculation'
    CASE ( 'S' )
       WRITE( PRTFile, * ) 'Semi-coherent TL calculation'
    CASE ( 'C' )
       WRITE( PRTFile, * ) 'Coherent TL calculation'
    CASE ( 'A' )
       WRITE( PRTFile, * ) 'Arrivals calculation, ASCII  file output'
    CASE ( 'a' )
       WRITE( PRTFile, * ) 'Arrivals calculation, binary file output'
    CASE DEFAULT
       CALL ERROUT( 'READIN', 'Unknown RunType selected' )
    END SELECT

    SELECT CASE ( RunType( 2 : 2 ) )
    CASE ( 'C' )
       WRITE( PRTFile, * ) 'Cartesian beams'
    CASE ( 'R' )
       WRITE( PRTFile, * ) 'Ray centered beams'
    CASE ( 'S' )
       WRITE( PRTFile, * ) 'Simple gaussian beams'
    CASE ( 'b' )
       WRITE( PRTFile, * ) 'Geometric gaussian beams in ray-centered coordinates'
    CASE ( 'B' )
       WRITE( PRTFile, * ) 'Geometric gaussian beams in Cartesian coordinates'
    CASE ( 'g' )
       WRITE( PRTFile, * ) 'Geometric hat beams in ray-centered coordinates'
    CASE DEFAULT
       RunType( 2 : 2 ) = 'G'
       WRITE( PRTFile, * ) 'Geometric hat beams in Cartesian coordinates'
    END SELECT

    SELECT CASE ( RunType( 4 : 4 ) )
    CASE ( 'R' )
       WRITE( PRTFile, * ) 'Point source (cylindrical coordinates)'
    CASE ( 'X' )
       WRITE( PRTFile, * ) 'Line source (Cartesian coordinates)'
    CASE DEFAULT
       RunType( 4 : 4 ) = 'R'
       WRITE( PRTFile, * ) 'Point source (cylindrical coordinates)'
    END SELECT

    SELECT CASE ( RunType( 5 : 5 ) )
    CASE ( 'R' )
       WRITE( PRTFile, * ) 'Rectilinear receiver grid: Receivers at ( Rr( ir ), Rz( ir ) ) )'
       PlotType = 'rectilin  '
    CASE ( 'I' )
       WRITE( PRTFile, * ) 'Irregular grid: Receivers at Rr( : ) x Rz( : )'
       IF ( Pos%NRz /= Pos%NRr ) CALL ERROUT( 'READIN', 'Irregular grid option selected with NRz not equal to Nr' )
       PlotType = 'irregular '
    CASE DEFAULT
       WRITE( PRTFile, * ) 'Rectilinear receiver grid: Receivers at Rr( : ) x Rz( : )'
       RunType( 5 : 5 ) = 'R'
       PlotType = 'rectilin  '
    END SELECT

    SELECT CASE ( RunType( 6 : 6 ) )
    CASE ( '2' )
       WRITE( PRTFile, * ) 'N x 2D calculation (neglects horizontal refraction)'
    CASE ( '3' )
       WRITE( PRTFile, * ) '3D calculation'
    CASE DEFAULT
       RunType( 6 : 6 ) = '2'
    END SELECT

  END SUBROUTINE ReadRunType

  ! **********************************************************************!

  SUBROUTINE TopBot( freq, AttenUnit, HS )
    !! Handles top and bottom boundary conditions

    REAL     (KIND=8), INTENT( IN    ) :: freq     ! center / nominal frequency (wideband not supported)
    CHARACTER (LEN=2), INTENT( IN    ) :: AttenUnit
    TYPE( HSInfo ),    INTENT( INOUT ) :: HS
    REAL     (KIND=8) :: Mz, vr, alpha2_f          ! values related to grain size

    ! Echo to PRTFile user's choice of boundary condition

    SELECT CASE ( HS%BC )
    CASE ( 'V' )
       WRITE( PRTFile, * ) '    VACUUM'
    CASE ( 'R' )
       WRITE( PRTFile, * ) '    Perfectly RIGID'
    CASE ( 'A' )
       WRITE( PRTFile, * ) '    ACOUSTO-ELASTIC half-space'
    CASE ( 'G' )
       WRITE( PRTFile, * ) '    Grain size to define half-space'
    CASE ( 'F' )
       WRITE( PRTFile, * ) '    FILE used for reflection loss'
    CASE ( 'W' )
       WRITE( PRTFile, * ) '    Writing an IRC file'
    CASE ( 'P' )
       WRITE( PRTFile, * ) '    reading PRECALCULATED IRC'
    CASE DEFAULT
       CALL ERROUT( 'TopBot', 'Unknown boundary condition type' )
    END SELECT

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

    HS%cp  = 0.0
    HS%cs  = 0.0
    HS%rho = 0.0

    SELECT CASE ( HS%BC )
    CASE ( 'A' )                  ! *** Half-space properties ***
       zTemp = 0.0
       READ(  ENVFile, *    ) zTemp, alphaR, betaR, rhoR, alphaI, betaI
       WRITE( PRTFile, FMT = "( F10.2, 3X, 2F10.2, 3X, F6.2, 3X, 2F10.4 )" ) &
            zTemp, alphaR, betaR, rhoR, alphaI, betaI

       ! dummy parameters for a layer with a general power law for attenuation
       ! these are not in play because the AttenUnit for this is not allowed yet
       !freq0         = freq
       betaPowerLaw  = 1.0
       ft            = 1000.0

       HS%cp  = CRCI( zTemp, alphaR, alphaI, freq, freq, AttenUnit, betaPowerLaw, ft )
       HS%cs  = CRCI( zTemp, betaR,  betaI,  freq, freq, AttenUnit, betaPowerLaw, ft )

       HS%rho = rhoR
    CASE ( 'G' )                  ! *** Grain size (formulas from UW-APL HF Handbook)

       ! These formulas are from the UW-APL Handbook
       ! The code is taken from older Matlab and is unnecesarily verbose
       ! vr   is the sound speed ratio
       ! rhor is the density ratio
       READ(  ENVFile, *    ) zTemp, Mz
       WRITE( PRTFile, FMT = "( F10.2, 3X, F10.2 )" ) zTemp, Mz

       IF ( Mz >= -1 .AND. Mz < 1 ) THEN
          vr   = 0.002709 * Mz ** 2 - 0.056452 * Mz + 1.2778
          rhor = 0.007797 * Mz ** 2 - 0.17057  * Mz + 2.3139
       ELSE IF ( Mz >= 1 .AND. Mz < 5.3 ) THEN
          vr   = -0.0014881 * Mz ** 3 + 0.0213937 * Mz ** 2 - 0.1382798 * Mz + 1.3425
          rhor = -0.0165406 * Mz ** 3 + 0.2290201 * Mz ** 2 - 1.1069031 * Mz + 3.0455
       ELSE
          vr   = -0.0024324 * Mz + 1.0019
          rhor = -0.0012973 * Mz + 1.1565
       END IF

       IF ( Mz >= -1 .AND. Mz < 0 ) THEN
          alpha2_f = 0.4556
       ELSE IF ( Mz >= 0 .AND. Mz < 2.6 ) THEN
          alpha2_f = 0.4556 + 0.0245 * Mz
       ELSE IF( Mz >= 2.6 .AND. Mz < 4.5 ) THEN
          alpha2_f = 0.1978 + 0.1245 * Mz
       ELSE IF( Mz >= 4.5 .AND. Mz < 6.0 ) THEN
          alpha2_f = 8.0399 - 2.5228 * Mz + 0.20098 * Mz ** 2
       ELSE IF( Mz >= 6.0 .AND. Mz < 9.5 ) THEN
          alpha2_f = 0.9431 - 0.2041 * Mz + 0.0117 * Mz ** 2
       ELSE
          alpha2_f =  0.0601
       END IF

       ! AttenUnit = 'L'   ! loss parameter
!!! following uses a reference sound speed of 1500 ???
!!! should be sound speed in the water, just above the sediment
       ! the term vr / 1000 converts vr to units of m per ms
       alphaR = vr * 1500.0
       alphaI = alpha2_f * ( vr / 1000 ) * 1500.0 * log( 10.0 ) / ( 40.0 * pi )   ! loss parameter Sect. IV., Eq. (4) of handbook

       HS%cp  = CRCI( zTemp, alphaR, alphaI, freq, freq, 'L ', betaPowerLaw, ft )
       HS%cs  = 0.0
       HS%rho = rhoR
       WRITE( PRTFile, FMT = "( 'Converted sound speed =', 2F10.2, 3X, 'density = ', F10.2, 3X, 'loss parm = ', F10.4 )" ) &
            HS%cp, rhor, alphaI

    CASE ( 'V', 'R', 'F', 'W', 'P' )
       CONTINUE
    CASE DEFAULT
       WRITE( ERROR_UNIT, * ) 'TopBot: Unknown boundary condition type: ', HS%BC
       CALL ERROUT( 'TopBot', 'Unknown boundary condition type' )
    END SELECT

  END SUBROUTINE TopBot

  ! **********************************************************************!

  SUBROUTINE OpenOutputFiles( FileRoot, ThreeD )
    !! Opens output files based on run type

    ! Write appropriate header information

    USE SourceReceiverPositions
    USE angleMod
    USE bdryMod
    USE RWSHDFile

    LOGICAL,            INTENT( IN ) :: ThreeD
    CHARACTER (LEN=80), INTENT( IN ) :: FileRoot
    REAL               :: atten
    CHARACTER (LEN=10) :: PlotType

    SELECT CASE ( Beam%RunType( 1 : 1 ) )
    CASE ( 'R', 'E' )   ! Ray trace or Eigenrays
       OPEN ( FILE = TRIM( FileRoot ) // '.ray', UNIT = RAYFile, FORM = 'FORMATTED' )
       WRITE( RAYFile, * ) '''', Title( 1 : 50 ), ''''
       WRITE( RAYFile, * ) freq
       WRITE( RAYFile, * ) Pos%NSx, Pos%NSy, Pos%NSz
       WRITE( RAYFile, * ) Angles%Nalpha, Angles%Nbeta
       WRITE( RAYFile, * ) Bdry%Top%HS%Depth
       WRITE( RAYFile, * ) Bdry%Bot%HS%Depth

       IF ( ThreeD ) THEN
          WRITE( RAYFile, * ) '''xyz'''
       ELSE
          WRITE( RAYFile, * ) '''rz'''
       END IF

    CASE ( 'A' )        ! arrival file in ascii format
       OPEN ( FILE = TRIM( FileRoot ) // '.arr', UNIT = ARRFile, FORM = 'FORMATTED' )

       IF ( ThreeD ) THEN
          WRITE( ARRFile, * ) '''3D'''
       ELSE
          WRITE( ARRFile, * ) '''2D'''
       END IF

       WRITE( ARRFile, * ) freq

       ! write source locations

       IF ( ThreeD ) THEN
          WRITE( ARRFile, * ) Pos%NSx,    Pos%Sx(    1 : Pos%NSx )
          WRITE( ARRFile, * ) Pos%NSy,    Pos%Sy(    1 : Pos%NSy )
          WRITE( ARRFile, * ) Pos%NSz,    Pos%Sz(    1 : Pos%NSz )
       ELSE
          WRITE( ARRFile, * ) Pos%NSz,    Pos%Sz(    1 : Pos%NSz )
       END IF

       ! write receiver locations

       WRITE( ARRFile, *    ) Pos%NRz,    Pos%Rz(    1 : Pos%NRz )
       WRITE( ARRFile, *    ) Pos%NRr,    Pos%Rr(    1 : Pos%NRr )
       IF ( ThreeD ) THEN
          WRITE( ARRFile, * ) Pos%Ntheta, Pos%theta( 1 : Pos%Ntheta )
       END IF

    CASE ( 'a' )        ! arrival file in binary format
       OPEN ( FILE = TRIM( FileRoot ) // '.arr', UNIT = ARRFile, FORM = 'UNFORMATTED' )

       IF ( ThreeD ) THEN
          WRITE( ARRFile ) '''3D'''
       ELSE
          WRITE( ARRFile ) '''2D'''
       END IF

       WRITE( ARRFile    ) SNGL( freq )

       ! write source locations

       IF ( ThreeD ) THEN
          WRITE( ARRFile    ) Pos%NSx,    Pos%Sx(    1 : Pos%NSx )
          WRITE( ARRFile    ) Pos%NSy,    Pos%Sy(    1 : Pos%NSy )
          WRITE( ARRFile    ) Pos%NSz,    Pos%Sz(    1 : Pos%NSz )
       ELSE
          WRITE( ARRFile    ) Pos%NSz,    Pos%Sz(    1 : Pos%NSz )
       END IF

       ! write receiver locations

       WRITE( ARRFile       ) Pos%NRz,    Pos%Rz(    1 : Pos%NRz )
       WRITE( ARRFile       ) Pos%NRr,    Pos%Rr(    1 : Pos%NRr )
       IF ( ThreeD ) THEN
          WRITE( ARRFile    ) Pos%Ntheta, Pos%theta( 1 : Pos%Ntheta )
       END IF

    CASE DEFAULT
       atten = 0.0

       ! following to set PlotType has already been done in READIN if that was used for input
       SELECT CASE ( Beam%RunType( 5 : 5 ) )
       CASE ( 'R' )
          PlotType = 'rectilin  '
       CASE ( 'I' )
          PlotType = 'irregular '
       CASE DEFAULT
          PlotType = 'rectilin  '
       END SELECT

       CALL WriteHeader( TRIM( FileRoot ) // '.shd', Title, REAL( freq ), atten, PlotType )
    END SELECT

  END SUBROUTINE OpenOutputFiles

END MODULE ReadEnvironmentBell