TopBot Subroutine

public subroutine TopBot(freq, AttenUnit, HS)

Handles top and bottom boundary conditions

! following uses a reference sound speed of 1500 ??? ! should be sound speed in the water, just above the sediment

Arguments

Type IntentOptional Attributes Name
real(kind=8), intent(in) :: freq
character(len=2), intent(in) :: AttenUnit
type(HSInfo), intent(inout) :: HS

Calls

proc~~topbot~~CallsGraph proc~topbot TopBot proc~crci CRCI proc~topbot->proc~crci proc~errout ERROUT proc~topbot->proc~errout proc~crci->proc~errout proc~franc_garr Franc_Garr proc~crci->proc~franc_garr

Called by

proc~~topbot~~CalledByGraph proc~topbot TopBot proc~readenvironment ReadEnvironment proc~readenvironment->proc~topbot program~bellhop BELLHOP program~bellhop->proc~readenvironment program~bellhop3d BELLHOP3D program~bellhop3d->proc~readenvironment

Source Code

  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

    END SELECT

  END SUBROUTINE TopBot