Hexahedral Subroutine

  SUBROUTINE Hexahedral( x, t, c, cimag, gradc, cxx, cyy, czz, cxy, cxz, cyz, rho, freq, Task )
    !! Trilinear hexahedral interpolation of SSP data

    ! assumes a rectilinear case (not the most general hexahedral)

    INTEGER,            PARAMETER     :: SSPFile = 40
    REAL      (KIND=8), INTENT( IN  ) :: freq
    REAL      (KIND=8), INTENT( IN  ) :: x( 3 )   ! x-y-z coordinate where sound speed is evaluated
    REAL      (KIND=8), INTENT( IN  ) :: t( 3 )   ! ray tangent; for edge cases of updating segments
    CHARACTER (LEN =3), INTENT( IN  ) :: Task
    REAL      (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 3 ), cxx, cyy, czz, cxy, cxz, cyz, rho ! sound speed and its derivatives
    INTEGER                           :: AllocateStatus, iSegxt, iSegyt, iy2, iz2
    REAL      (KIND=8)                :: c1, c2, c11, c12, c21, c22, cz11, cz12, cz21, cz22, cz1, cz2, &
                                         cx, cy, cz, s1, s2, s3

    IF ( Task == 'INI' ) THEN

       !  *** Section to read in SSP data ***

       ! Read dummy SSP information from the environmental file
       ! This is over-ridden by the info in the SSP file
       ! However, cz info is used in beam selection
       Depth = x( 3 )
       CALL ReadSSP( Depth, freq )

       ! Read the 3D SSP matrix

       WRITE( PRTFile, * )
       WRITE( PRTFile, * ) 'Reading sound speed profile from file'

       ! x coordinates
       READ( SSPFile,  * ) SSP%Nx
       WRITE( PRTFile, * )
       WRITE( PRTFile, * ) 'Number of points in x = ', SSP%Nx
       ALLOCATE( SSP%Seg%x( SSP%Nx ), STAT = AllocateStatus )
       IF ( AllocateStatus /= 0 ) CALL ERROUT( 'sspMod: Hexahedral', 'Insufficient memory to store SSP'  )
       READ( SSPFile,  * ) SSP%Seg%x
       !WRITE( PRTFile, * )
       !WRITE( PRTFile, * ) 'x-coordinates of SSP (km):'
       !WRITE( PRTFile, FMT="( F10.2 )"  ) SSP%Seg%x( 1 : SSP%Nx )
       IF ( .NOT. monotonic( SSP%Seg%x, SSP%Nx ) ) THEN
          CALL ERROUT( 'sspMod: Hexahedral', 'The x coordinates in the SSP must be monotone increasing' )
       END IF

       ! y coordinates
       READ( SSPFile,  * ) SSP%Ny
       WRITE( PRTFile, * )
       WRITE( PRTFile, * ) 'Number of points in y = ', SSP%Ny
       ALLOCATE( SSP%Seg%y( SSP%Ny ), STAT = AllocateStatus )
       IF ( AllocateStatus /= 0 ) CALL ERROUT( 'sspMod: Hexahedral', 'Insufficient memory to store SSP'  )
       READ( SSPFile,  * ) SSP%Seg%y
       !WRITE( PRTFile, * )
       !WRITE( PRTFile, * ) 'y-coordinates of SSP (km):'
       !WRITE( PRTFile, FMT="( F10.2 )"  ) SSP%Seg%y( 1 : SSP%Ny )
       IF ( .NOT. monotonic( SSP%Seg%y, SSP%Ny ) ) THEN
          CALL ERROUT( 'sspMod: Hexahedral', 'The y coordinates in the SSP must be monotone increasing' )
       END IF

       ! z coordinates
       READ( SSPFile,  * ) SSP%Nz
       WRITE( PRTFile, * )
       WRITE( PRTFile, * ) 'Number of points in z = ', SSP%Nz
       ALLOCATE( SSP%Seg%z( SSP%Nz ), STAT = AllocateStatus )
       IF ( AllocateStatus /= 0 ) CALL ERROUT( 'sspMod: Hexahedral', 'Insufficient memory to store SSP'  )
       READ( SSPFile,  * ) SSP%Seg%z
       !WRITE( PRTFile, * )
       !WRITE( PRTFile, * ) 'z-coordinates of SSP (km):'
       !WRITE( PRTFile, FMT="( F10.2 )"  ) SSP%Seg%z( 1 : SSP%Nz )
       IF ( .NOT. monotonic( SSP%Seg%z, SSP%Nz ) ) THEN
          CALL ERROUT( 'sspMod: Hexahedral', 'The z coordinates in the SSP must be monotone increasing' )
       END IF

       ! SSP matrix should be bigger than 2x2x2
       IF ( SSP%Nx < 2 .OR. SSP%Ny < 2 .OR. SSP%Nz < 2 ) THEN
          CALL ERROUT( 'sspMod: Hexahedral', &
               'You must have at least two points in x, y, z directions in your 3D SSP field'  )
       END IF

       IF ( SSP%Nz >= MaxSSP ) THEN
          ! LP: SSP%Nz / SSP%Seg%z will get assigned to SSP%NPts / SSP%z.
          CALL ERROUT( 'sspMod: Hexahedral', &
               'Number of SSP points in Z exceeds limit' )
       END IF

       ALLOCATE( SSP%cMat3( SSP%Nx, SSP%Ny, SSP%Nz ), SSP%czMat3( SSP%Nx, SSP%Ny, SSP%Nz - 1 ), STAT = AllocateStatus )
       IF ( AllocateStatus /= 0 ) CALL ERROUT( 'sspMod: Hexahedral', 'Insufficient memory to store SSP'  )

       WRITE( PRTFile, * )
       ! WRITE( PRTFile, * ) 'Sound speed matrix:'
       DO iz2 = 1, SSP%Nz
          DO iy2 = 1, SSP%Ny
             READ(  SSPFile, * ) SSP%cMat3( :, iy2, iz2 )
             ! WRITE( PRTFile, FMT="( 'z-section = ', F10.2, '  m' )"  ) SSP%Seg%z( iz2 )
             ! WRITE( PRTFile, FMT="( 'y-section = ', F10.2, ' km' )"  ) SSP%Seg%y( iy2 )
             ! WRITE( PRTFile, FMT="( 12F10.2 )"  ) SSP%cMat3( :, iy2, iz2 )
          END DO
       END DO

       CLOSE( SSPFile )

       SSP%Seg%x = 1000.0 * SSP%Seg%x   ! convert km to m
       SSP%Seg%y = 1000.0 * SSP%Seg%y   ! convert km to m

       ! calculate cz
       DO iSegxt = 1, SSP%Nx
          DO iSegyt = 1, SSP%Ny
             DO iz2 = 2, SSP%Nz
                SSP%czMat3( iSegxt, iSegyt, iz2 - 1 ) = &
                     ( SSP%cMat3( iSegxt, iSegyt, iz2 ) - SSP%cMat3( iSegxt, iSegyt, iz2 - 1 ) ) / &
                     ( SSP%Seg%z(                 iz2 ) - SSP%Seg%z(                 iz2 - 1 ) )
             END DO
          END DO
       END DO

       ! over-ride the SSP%z values read in from the environmental file with these new values
       SSP%NPts = SSP%Nz
       SSP%z( 1 : SSP%Nz ) = SSP%Seg%z

       RETURN
    ELSE

       ! *** Section to return SSP info ***

       ! Check that x is inside the box where the sound speed is defined
       IF ( x( 1 ) < SSP%Seg%x( 1 ) .OR. x( 1 ) > SSP%Seg%x( SSP%Nx ) .OR. &
            x( 2 ) < SSP%Seg%y( 1 ) .OR. x( 2 ) > SSP%Seg%y( SSP%Ny ) ) THEN ! .OR. &
!!$          x( 3 ) < SSP%Seg%z( 1 ) .OR. x( 3 ) > SSP%Seg%z( SSP%Nz ) ) THEN
          WRITE( PRTFile, * ) 'ray is outside the box where the ocean soundspeed is defined'
          WRITE( PRTFile, * ) ' x = ( x, y, z ) = ', x
          CALL ERROUT( 'Hexahedral', 'ray is outside the box where the soundspeed is defined' )
       END IF

       CALL Update3DXSegmentT( x, t )
       CALL Update3DYSegmentT( x, t )
       CALL Update3DZSegmentT( x, t )

       ! cz at the corners of the current rectangle
       cz11 = SSP%czMat3( iSegx,     iSegy    , iSegz )
       cz12 = SSP%czMat3( iSegx + 1, iSegy    , iSegz )
       cz21 = SSP%czMat3( iSegx,     iSegy + 1, iSegz )
       cz22 = SSP%czMat3( iSegx + 1, iSegy + 1, iSegz )

       ! for this depth, x( 3 ) get the sound speed at the corners of the current rectangle
       s3 = x( 3 ) - SSP%Seg%z( iSegz )
       c11 = SSP%cMat3( iSegx,     iSegy    , iSegz ) + s3 * cz11
       c21 = SSP%cMat3( iSegx + 1, iSegy    , iSegz ) + s3 * cz12
       c12 = SSP%cMat3( iSegx,     iSegy + 1, iSegz ) + s3 * cz21
       c22 = SSP%cMat3( iSegx + 1, iSegy + 1, iSegz ) + s3 * cz22

       ! s1 = proportional distance of x( 1 ) in x
       s1 = ( x( 1 ) - SSP%Seg%x( iSegx ) ) / ( SSP%Seg%x( iSegx + 1 ) - SSP%Seg%x( iSegx ) )
       s1 = MIN( s1, 1.0D0 )   ! force piecewise constant extrapolation for points outside the box
       s1 = MAX( s1, 0.0D0 )   ! "

       ! s2 = proportional distance of x( 2 ) in y
       s2 = ( x( 2 ) - SSP%Seg%y( iSegy ) ) / ( SSP%Seg%y( iSegy + 1 ) - SSP%Seg%y( iSegy ) )
       s2 = MIN( s2, 1.0D0 )   ! force piecewise constant extrapolation for points outside the box
       s2 = MAX( s2, 0.0D0 )   ! "

       ! interpolate the soundspeed in the x direction, at the two endpoints in y (top and bottom sides of rectangle)
       c1 = c11 + s1 * ( c21 - c11 )
       c2 = c12 + s1 * ( c22 - c12 )
       !c  = ( 1.0D0 - s2 ) * c1  + s2 * c2   ! interpolation in y
       cy = ( c2 - c1 ) / ( SSP%Seg%y( iSegy + 1 ) - SSP%Seg%y( iSegy ) )

       ! interpolate the soundspeed in the y direction, at the two endpoints in x (left and right sides of rectangle)
       c1 = c11 + s2 * ( c12 - c11 )
       c2 = c21 + s2 * ( c22 - c21 )

       c  = c1  + s1 * ( c2 - c1   )   ! interpolation in x

       ! interpolate the attenuation !!!! This will use the wrong segment if the ssp in the envfil is sampled at different depths
       s3 = s3 / ( SSP%z( iSegz + 1 ) - SSP%z( iSegz ) )   ! convert s3 to a proportional distance in the layer
       cimag = AIMAG( ( 1.0D0 - s3 ) * SSP%c( Isegz )  + s3 * SSP%c( Isegz + 1 ) )   ! volume attenuation is taken from the single c(z) profile

       cx = ( c2 - c1 ) / ( SSP%Seg%x( iSegx + 1 ) - SSP%Seg%x( iSegx ) )

       ! same thing on cz
       cz1 = cz11 + s2 * ( cz21 - cz11 )
       cz2 = cz12 + s2 * ( cz22 - cz12 )
       cz  = cz1  + s1 * ( cz2  - cz1  )   ! interpolation in z

       !gradc = [ cx, cy, cz ]
       gradc( 1 ) = cx
       gradc( 2 ) = cy
       gradc( 3 ) = cz

       cxx   = 0.0D0
       cyy   = 0.0D0
       czz   = 0.0D0
       cxy   = 0.0D0
       cxz   = 0.0D0
       cyz   = 0.0D0
       ! write( *, FMT="( 'SSP', I3, 2X, 4F10.2, F10.4 )" ) layer, x, c, cz

       ! linear interpolation for density
       W   = ( x( 3 ) - SSP%z( iSegz ) ) / ( SSP%z( iSegz + 1 ) - SSP%z( iSegz ) )
       rho = ( 1.0D0 - W ) * SSP%rho( iSegz ) + W * SSP%rho( iSegz + 1 )
    END IF

  END SUBROUTINE Hexahedral