MODULE AttenMod

  ! Attenuation module
  ! Routines to convert a sound speed and attenuation in user units to a complex sound speed
  ! Includes a formula for volume attenuation

  USE FatalError
  IMPLICIT NONE
  INTEGER, PRIVATE, PARAMETER      :: PRTFile = 6
  INTEGER, PARAMETER               :: MaxBioLayers = 200
  INTEGER                          :: iBio, NBioLayers
  REAL    (KIND=8) :: T = 20, Salinity = 35, pH = 8, z_bar = 0, FG   ! Francois-Garrison volume attenuation; temperature, salinity, ...

  TYPE bioStructure
     REAL (KIND=8) :: Z1, Z2, f0, Q, a0
  END TYPE bioStructure

  TYPE( bioStructure ) :: bio( MaxBioLayers )

CONTAINS
  FUNCTION CRCI( z, c, alpha, freq, freq0, AttenUnit, beta, fT )

    ! Converts real wave speed and attenuation to a single
    !  complex wave speed (with positive imaginary part)

    ! AttenUnit
    ! 6 Cases:    N Nepers/meter
    !             M dB/meter      (M for Meters)
    !             m dB/meter with a power law
    !             F dB/m-kHz      (F for frequency dependent)
    !             W dB/wavelength (W for Wavelength)
    !             Q Q
    !             L Loss parameter
    !
    ! second letter adds volume attenuation according to standard laws:
    !             T for Thorp
    !             F for Francois Garrison
    !             B for biological
    !
    ! freq is the current frequency
    ! freq0 is the reference frequency for which the dB/meter was specified (used only for 'm')

    ! Returns
    ! c     real      part of sound speed
    ! alpha imaginary part of sound speed

    USE MathConstants

    REAL     (KIND=8), INTENT( IN )  :: freq, freq0, alpha, c, z, beta, fT
    CHARACTER (LEN=2), INTENT( IN )  :: AttenUnit
    REAL     (KIND=8)                :: f2, omega, alphaT, Thorp, a, FG
    COMPLEX  (KIND=8)                :: CRCI

    omega = 2.0 * pi * freq

    !  Convert to Nepers/m 
    alphaT = 0.0
    SELECT CASE ( AttenUnit( 1 : 1 ) )
    CASE ( 'N' )
       alphaT = alpha
    CASE ( 'M' )   ! dB/m
       alphaT = alpha / 8.6858896D0
    CASE ( 'm' )   ! dB/m with power law
       alphaT = alpha / 8.6858896D0
       IF ( freq < fT ) THEN   ! frequency raised to the power beta
          alphaT = alphaT * ( freq / freq0 ) ** beta
       ELSE                    ! linear in frequency
          alphaT = alphaT * ( freq / freq0 ) * ( fT / freq0 ) ** ( beta - 1 )
       END IF
    CASE ( 'F' )   ! dB/(m kHz)
       alphaT = alpha * freq / 8685.8896D0
    CASE ( 'W' )   ! dB/wavelength
       IF ( c /= 0.0         ) alphaT = alpha * freq / ( 8.6858896D0 * c )
       !        The following lines give f^1.25 frequency dependence
       !        FAC = SQRT( SQRT( freq / 50.0 ) )
       !        IF ( c /= 0.0 ) alphaT = FAC * alpha * freq / ( 8.6858896D0 * c )
    CASE ( 'Q' )   ! Quality factor
       IF ( c * alpha /= 0.0 ) alphaT = omega / ( 2.0 * c * alpha )
    CASE ( 'L' )   ! loss parameter
       IF ( c /= 0.0         ) alphaT = alpha * omega / c
    END SELECT

    ! added volume attenuation
    SELECT CASE ( AttenUnit( 2 : 2 ) )
    CASE ( 'T' )
       f2 = ( freq / 1000.0 ) ** 2

       ! Original formula from Thorp 1967
       ! Thorp = 40.0 * f2 / ( 4100.0 + f2 ) + 0.1 * f2 / ( 1.0 + f2 )   ! dB/kyard
       ! Thorp = Thorp / 914.4D0                 ! dB / m
       ! Thorp = Thorp / 8.6858896D0             ! Nepers / m

       ! Updated formula from JKPS Eq. 1.34
       Thorp = 3.3d-3 + 0.11 * f2 / ( 1.0 + f2 ) + 44.0 * f2 / ( 4100.0 + f2 ) + 3d-4 * f2   ! dB/km
       Thorp = Thorp / 8685.8896 ! Nepers / m

       alphaT = alphaT + Thorp
    CASE ( 'F' )   ! Francois-Garrison
       FG     = Franc_Garr( freq / 1000 );   ! dB/km
       FG     = FG / 8685.8896;              ! Nepers / m
       alphaT = alphaT + FG
    CASE ( 'B' )   ! biological attenuation per Orest Diachok
       DO iBio = 1, NBioLayers
          IF ( z >= bio( iBio )%Z1 .AND. z <= bio( iBio )%Z2 ) THEN
             a = bio( iBio )%a0 / ( ( 1.0 - bio( iBio )%f0 ** 2 / freq ** 2  ) ** 2 + 1.0 / bio( iBio )%Q ** 2 )   ! dB/km
             a = a / 8685.8896   ! Nepers / m
             alphaT = alphaT + a
          END IF
       END DO
    END SELECT

    ! Convert Nepers/m to equivalent imaginary sound speed 
    alphaT = alphaT * c * c / omega
    CRCI   = CMPLX( c, alphaT, KIND=8 )

    IF ( alphaT > c ) THEN
       WRITE( PRTFile, * ) 'Complex sound speed: ', CRCI
       WRITE( PRTFile, * ) 'Usually this means you have an attenuation that is way too high'
       CALL ERROUT( 'AttenMod : CRCI ', 'The complex sound speed has an imaginary part > real part' )
    END IF

    RETURN
  END FUNCTION CRCI

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

  FUNCTION Franc_Garr( f )

    ! Francois Garrison formulas for attenuation
    ! Based on a Matlab version by D. Jackson APL-UW

    ! mbp Feb. 2019
    ! Verified using F-G Table IV

    ! alpha = attenuation   (dB/km)
    ! f     = frequency     (kHz)
    ! T     = temperature   (deg C)
    ! S     = salinity      (psu)
    ! pH    = 7 for neutral water
    ! z_bar = depth         (m)

    !     Returns
    !        alpha = volume attenuation in dB/km

    REAL (KIND=8) :: f, Franc_Garr
    REAL (KIND=8) :: c, A1, A2, A3, P1, P2, P3, f1, f2

    c = 1412 + 3.21 * T + 1.19 * Salinity + 0.0167 * z_bar

    ! Boric acid contribution
    A1 = 8.86 / c * 10 ** ( 0.78 * pH - 5 )
    P1 = 1
    f1 = 2.8 * sqrt( Salinity / 35 ) * 10 ** ( 4 - 1245 / ( T + 273 ) )

    ! Magnesium sulfate contribution
    A2 = 21.44 * Salinity / c * ( 1 + 0.025 * T )
    P2 = 1 - 1.37D-4 * z_bar + 6.2D-9 * z_bar ** 2
    f2 = 8.17 * 10 ** ( 8 - 1990 / ( T + 273 ) ) / ( 1 + 0.0018 * ( Salinity - 35 ) )

    ! Viscosity
    P3 = 1 - 3.83D-5 * z_bar + 4.9D-10 * z_bar ** 2
    if ( T < 20 ) THEN
       A3 = 4.937D-4 - 2.59D-5 * T + 9.11D-7 * T ** 2 - 1.5D-8  * T ** 3
    else
       A3 = 3.964D-4 -1.146D-5 * T + 1.45D-7 * T ** 2 - 6.5D-10 * T ** 3
    end if

    Franc_Garr = A1 * P1 * ( f1 * f ** 2 ) / ( f1 ** 2 + f ** 2 ) + A2 * P2 * ( f2 * f ** 2 ) / ( f2 ** 2 + f ** 2 ) + &
         A3 * P3 * f ** 2

  END FUNCTION Franc_Garr


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

  SUBROUTINE GrainSize_to_Geoacoustic( Mz, vr, rhor, alpha2_f )
    ! These formulas are from the UW-APL Handbook
 
    ! vr   is the sound speed ratio
    ! rhor is the density ratio

    REAL (KIND=8), INTENT( IN  ) :: Mz                            ! Grain size (if selected)
    REAL (KIND=8), INTENT( OUT ) :: vr, rhor, alpha2_f            ! grain size info (sound speed ratio, density, attenuation)
 
    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
  END SUBROUTINE GrainSize_to_Geoacoustic

END MODULE AttenMod
