ReadEnvironment Subroutine

  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' )
    ENDIF

    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