SourceReceiverPositions.f90 Source File

!! Provides methods to read in depths, depths, bearing data

MODULE SourceReceiverPositions
!! Reads in source depths, receiver depths, receiver ranges, and receiver bearings

  USE monotonicMod
  USE SortMod
  USE SubTabulate
  USE FatalError

  IMPLICIT NONE
  PUBLIC

  INTEGER, PARAMETER          :: Number_to_Echo = 10
  INTEGER, PRIVATE            :: IAllocStat     ! used to capture status after allocation
  INTEGER, PRIVATE, PARAMETER :: ENVFile = 5, PRTFile = 6   ! unit 5 is usually (not always) the ENVFile
  INTEGER                     :: Nfreq = 1      ! number of frequencies
  REAL (KIND=8), ALLOCATABLE  :: freqVec( : )   ! frequency vector for broadband runs

  TYPE Position
     INTEGER              :: NSx = 1, NSy = 1, NSz = 1, NRz = 1, NRr = 1, Ntheta = 1    ! number of x, y, z, r, theta coordinates
     REAL                 :: Delta_r, Delta_theta
     INTEGER, ALLOCATABLE :: iSz( : ), iRz( : )
     REAL,    ALLOCATABLE :: Sx( : ), Sy( : ), Sz( : )          ! Source x, y, z coordinates
     REAL,    ALLOCATABLE :: Rr( : ), Rz( : ), ws( : ), wr( : ) ! Receiver r, z coordinates and weights for interpolation
     REAL,    ALLOCATABLE :: theta( : )                         ! Receiver bearings
  END TYPE Position

  TYPE ( Position ) :: Pos   ! structure containing source and receiver positions

CONTAINS

  FUNCTION RToIR( r )

    ! index of nearest rcvr before normal
    ! Compute upper index on rcvr line
    !!! Assumes Pos%r is a vector of equally spaced points

    REAL ( KIND=8 ), INTENT( IN ) :: r
    INTEGER                       :: RToIR, ti
    REAL ( KIND=8 )               :: temp

    temp = ( r - Pos%Rr( 1 ) ) / Pos%Delta_r
    ! LP: Added snapping to deal with floating-point error at the int boundaries.
    IF ( ABS( temp - REAL( NINT( temp ), 8 ) ) < 1D-6 ) THEN
       temp = REAL( NINT( temp ), 8 )
    END IF
    ! mbp: should be ", 0 )" ? [LP: on 2D Cerveny cart]
    RToIR = MAX( MIN( INT( temp ) + 1, Pos%NRr ), 1 )

  END FUNCTION RToIR

  SUBROUTINE ReadfreqVec( freq0, BroadbandOption )

    ! Optionally reads a vector of source frequencies for a broadband run
    ! If the broadband option is not selected, then the input freq (a scalar) is stored in the frequency vector

    REAL (KIND=8), INTENT( IN ) :: freq0             ! Nominal or carrier frequency
    CHARACTER(LEN=1), INTENT(IN) :: BroadbandOption
    INTEGER                     :: ifreq

    ! Broadband run?
    IF ( BroadbandOption == 'B' ) THEN
       READ( ENVFile, * ) Nfreq
       WRITE( PRTFile, * ) '__________________________________________________________________________'
       WRITE( PRTFile, * )
       WRITE( PRTFile, * )
       WRITE( PRTFile, * ) '   Number of frequencies =', Nfreq
       IF ( Nfreq <= 0 ) CALL ERROUT( 'ReadEnvironment', 'Number of frequencies must be positive'  )
    END IF

    IF ( ALLOCATED( freqVec ) ) DEALLOCATE( freqVec )
    ALLOCATE( freqVec( MAX( 3, Nfreq ) ), Stat = IAllocStat )
    IF ( IAllocStat /= 0 ) CALL ERROUT( 'ReadEnvironment', 'Too many frequencies'  )

    IF ( BroadbandOption == 'B' ) THEN
       WRITE( PRTFile, * ) '   Frequencies (Hz)'
       freqVec( 2 ) = -999.9
       freqVec( 3 ) = -999.9
       READ(  ENVFile, * ) freqVec( 1 : Nfreq )
       CALL SubTab( freqVec, Nfreq )

       WRITE( PRTFile, "( 5G14.6 )" ) ( freqVec( ifreq ), ifreq = 1, MIN( Nfreq, Number_to_Echo ) )
       IF ( Nfreq > Number_to_Echo ) WRITE( PRTFile,  "( G14.6 )" ) ' ... ', freqVec( Nfreq )
    ELSE
       freqVec( 1 ) = freq0
    END IF

    RETURN

  END SUBROUTINE ReadfreqVec

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

  SUBROUTINE ReadSxSy( ThreeD )

    ! Reads source x-y coordinates

    LOGICAL, INTENT( IN ) :: ThreeD   ! flag indicating whether this is a 3D run

    IF ( ThreeD ) THEN
       CALL ReadVector( Pos%NSx, Pos%Sx, 'Source   x-coordinates, Sx', 'km' )
       CALL ReadVector( Pos%NSy, Pos%Sy, 'Source   y-coordinates, Sy', 'km' )
    ELSE
       ALLOCATE( Pos%Sx( 1 ), Pos%Sy( 1 ) )
       Pos%Sx( 1 ) = 0.
       Pos%Sy( 1 ) = 0.
    END IF

    RETURN
  END SUBROUTINE ReadSxSy

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

  SUBROUTINE ReadSzRz( zMin, zMax )

    ! Reads source and receiver z-coordinates (depths)
    ! zMin and zMax are limits for those depths; sources and receivers are shifted to be within those limits

    REAL,    INTENT( IN ) :: zMin, zMax
    !LOGICAL               :: monotonic

    CALL ReadVector( Pos%NSz, Pos%Sz, 'Source   z-coordinates, Sz', 'm' )
    CALL ReadVector( Pos%NRz, Pos%Rz, 'Receiver z-coordinates, Rz', 'm' )

    IF ( ALLOCATED( Pos%ws ) ) DEALLOCATE( Pos%ws, Pos%iSz )
    ALLOCATE( Pos%ws( Pos%NSz ), Pos%iSz( Pos%NSz ), Stat = IAllocStat )
    IF ( IAllocStat /= 0 ) CALL ERROUT( 'ReadSzRz', 'Too many sources'  )

    IF ( ALLOCATED( Pos%wr ) ) DEALLOCATE( Pos%wr, Pos%iRz )
    ALLOCATE( Pos%wr( Pos%NRz ), Pos%iRz( Pos%NRz ), Stat = IAllocStat  )
    IF ( IAllocStat /= 0 ) CALL ERROUT( 'ReadSzRz', 'Too many receivers'  )

    ! *** Check for Sz/Rz in upper or lower halfspace ***

    IF ( ANY( Pos%Sz( 1 : Pos%NSz ) < zMin ) ) THEN
       WHERE ( Pos%Sz < zMin ) Pos%Sz = zMin
       WRITE( PRTFile, * ) 'Warning in ReadSzRz : Source above or too near the top bdry has been moved down'
    END IF

    IF ( ANY( Pos%Sz( 1 : Pos%NSz ) > zMax ) ) THEN
       WHERE( Pos%Sz > zMax ) Pos%Sz = zMax
       WRITE( PRTFile, * ) 'Warning in ReadSzRz : Source below or too near the bottom bdry has been moved up'
    END IF

    IF ( ANY( Pos%Rz( 1 : Pos%NRz ) < zMin ) ) THEN
       WHERE( Pos%Rz < zMin ) Pos%Rz = zMin
       WRITE( PRTFile, * ) 'Warning in ReadSzRz : Receiver above or too near the top bdry has been moved down'
    END IF

    IF ( ANY( Pos%Rz( 1 : Pos%NRz ) > zMax ) ) THEN
       WHERE( Pos%Rz > zMax ) Pos%Rz = zMax
       WRITE( PRTFile, * ) 'Warning in ReadSzRz : Receiver below or too near the bottom bdry has been moved up'
    END IF

!!$    IF ( .NOT. monotonic( Pos%sz, Pos%NSz ) ) THEN
!!$       CALL ERROUT( 'SzRzRMod', 'Source depths are not monotonically increasing' )
!!$    END IF
!!$
!!$    IF ( .NOT. monotonic( Pos%rz, Pos%NRz ) ) THEN
!!$       CALL ERROUT( 'SzRzRMod', 'Receiver depths are not monotonically increasing' )
!!$    END IF

    RETURN
  END SUBROUTINE ReadSzRz

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

  SUBROUTINE ReadRcvrRanges

    CALL ReadVector( Pos%NRr, Pos%Rr, 'Receiver r-coordinates, Rr', 'km' )

    ! calculate range spacing
    Pos%delta_r = 0.0
    IF ( Pos%NRr /= 1 ) Pos%delta_r = Pos%Rr( Pos%NRr ) - Pos%Rr( Pos%NRr - 1 )

    IF ( .NOT. monotonic( Pos%rr, Pos%NRr ) ) THEN
       CALL ERROUT( 'ReadRcvrRanges', 'Receiver ranges are not monotonically increasing' )
    END IF

    RETURN
  END SUBROUTINE ReadRcvrRanges

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

  SUBROUTINE ReadRcvrBearings

    CALL ReadVector( Pos%Ntheta, Pos%theta, 'Receiver bearings, theta', 'degrees' )

    ! full 360-degree sweep? remove duplicate angle
    IF ( Pos%Ntheta > 1 ) THEN
       IF ( ABS( MOD( Pos%theta( Pos%Ntheta ) - Pos%theta( 1 ), 360.0 ) ) < 10.0 * TINY( 1.0D0 ) ) &
          Pos%Ntheta = Pos%Ntheta - 1
    END IF

    ! calculate angular spacing
    Pos%Delta_theta = 0.0
    IF ( Pos%Ntheta /= 1 ) Pos%Delta_theta = Pos%theta( Pos%Ntheta ) - Pos%theta( Pos%Ntheta - 1 )

    IF ( .NOT. monotonic( Pos%theta, Pos%Ntheta ) ) THEN
       CALL ERROUT( 'ReadRcvrBearings', 'Receiver bearings are not monotonically increasing' )
    END IF

    RETURN
  END SUBROUTINE ReadRcvrBearings

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

  SUBROUTINE ReadVector( Nx, x, Description, Units )

    ! Read a vector x
    ! Description is something like 'receiver ranges'
    ! Units       is something like 'km'

    INTEGER,   INTENT( OUT ) :: Nx
    REAL,      ALLOCATABLE, INTENT( OUT ) :: x( : )
    CHARACTER, INTENT( IN  ) :: Description*( * ), Units*( * )
    INTEGER                  :: ix

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

    READ(  ENVFile, * ) Nx
    WRITE( PRTFile, * ) '   Number of ' // Description // ' = ', Nx

    IF ( Nx <= 0 ) CALL ERROUT( 'ReadVector', 'Number of ' // Description // 'must be positive'  )

    IF ( ALLOCATED( x ) ) DEALLOCATE( x )
    ALLOCATE( x( MAX( 3, Nx ) ), Stat = IAllocStat )
    IF ( IAllocStat /= 0 ) CALL ERROUT( 'ReadVector', 'Too many ' // Description )

    WRITE( PRTFile, * ) '   ', Description // ' (' // Units // ')'
    x( 3 ) = -999.9
    READ( ENVFile, * ) x( 1 : Nx )

    CALL SubTab( x, Nx )
    CALL Sort(   x, Nx )

    WRITE( PRTFile, "( 5G14.6 )" ) '   ', ( x( ix ), ix = 1, MIN( Nx, Number_to_Echo ) )
    IF ( Nx > Number_to_Echo ) WRITE( PRTFile, "( G14.6 )" ) ' ... ', x( Nx )

    WRITE( PRTFile, * )

    ! Vectors in km should be converted to m for internal use
    IF ( LEN_TRIM( Units ) >= 2 ) THEN
       IF ( Units( 1 : 2 ) == 'km' ) x = 1000.0 * x
    END IF

  END SUBROUTINE ReadVector

END MODULE SourceReceiverPositions