Quad Subroutine

  SUBROUTINE Quad( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
    !! Quadrilateral interpolation for sound speed profiles

    ! Bilinear quadrilateral interpolation of SSP data in 2D

    INTEGER,           PARAMETER      :: SSPFile = 40
    REAL      (KIND=8), INTENT( IN  ) :: freq
    REAL      (KIND=8), INTENT( IN  ) :: x( 2 )   ! r-z coordinate where sound speed is evaluated
    REAL      (KIND=8), INTENT( IN  ) :: t( 2 )   ! ray tangent; for edge cases of updating segments
    CHARACTER (LEN=3),  INTENT( IN  ) :: Task
    REAL      (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 2 ), crr, crz, czz, rho ! sound speed and its derivatives
    INTEGER                           :: AllocateStatus, iSegT, iz2
    REAL      (KIND=8)                :: c1, c2, cz1, cz2, cr, cz, s1, s2, delta_r, delta_z

    IF ( Task == 'INI' ) THEN

       !  *** read in SSP data ***

       Depth = x( 2 )
       CALL ReadSSP( Depth, freq )

      ! Read the 2D SSP matrix
       WRITE( PRTFile, * ) '__________________________________________________________________________'
       WRITE( PRTFile, * )
       WRITE( PRTFile, * ) 'Using range-dependent sound speed'

       READ( SSPFile,  * ) SSP%Nr
       WRITE( PRTFile, * ) 'Number of SSP ranges = ', SSP%Nr

       IF ( SSP%Nr < 2 ) THEN
          CALL ERROUT( 'sspMod: Quad', 'You must have at least two profiles in your 2D SSP field'  )
       END IF

       ALLOCATE( SSP%cMat( SSP%NPts, SSP%Nr ), SSP%czMat( SSP%NPts - 1, SSP%Nr ), SSP%Seg%r( SSP%Nr ), STAT = AllocateStatus )
       IF ( AllocateStatus /= 0 ) CALL ERROUT( 'sspMod: Quad', 'Insufficient memory to store SSP'  )

       READ( SSPFile,  * ) SSP%Seg%r( 1 : SSP%Nr )
       WRITE( PRTFile, * )
       WRITE( PRTFile, * ) 'Profile ranges (km):'
       WRITE( PRTFile, FMT="( F10.2 )"  ) SSP%Seg%r( 1 : SSP%Nr )

       IF ( .NOT. monotonic( SSP%Seg%r, SSP%Nr ) ) THEN
          CALL ERROUT( 'sspMod: Quad', 'The ranges in the SSP must be monotone increasing' )
       END IF

       SSP%Seg%r = 1000.0 * SSP%Seg%r   ! convert km to m

       WRITE( PRTFile, * )
       WRITE( PRTFile, * ) 'Sound speed matrix:'
       WRITE( PRTFile, * ) ' Depth (m )     Soundspeed (m/s)'
       DO iz2 = 1, SSP%NPts
          READ(  SSPFile, * ) SSP%cMat( iz2, : )
          ! WRITE( PRTFile, FMT="( 'iSegz depth = ', F10.2, ' m' )"  ) SSP%z( iz2 )
          ! WRITE( PRTFile, FMT="( 12F10.2 )"  ) SSP%cMat( iz2, : )
          WRITE( PRTFile, FMT="( 12F10.2 )"  ) SSP%z( iz2 ), SSP%cMat( iz2, : )
       END DO

       CLOSE( SSPFile )

       ! calculate cz
       DO iSegt = 1, SSP%Nr
          DO iz2 = 2, SSP%NPts
             delta_z = ( SSP%z( iz2 ) - SSP%z( iz2 - 1 ) )
             SSP%czMat( iz2 - 1, iSegt ) = ( SSP%cMat( iz2, iSegt ) - SSP%cMat( iz2 - 1, iSegt ) ) / delta_z
          END DO
       END DO

       SSP%Nz = SSP%NPts
       RETURN

    ELSE

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

       CALL UpdateDepthSegmentT( x, t )

       ! Check that x is inside the box where the sound speed is defined
       IF ( x( 1 ) < SSP%Seg%r( 1 ) .OR. x( 1 ) > SSP%Seg%r( SSP%Nr ) ) THEN ! .OR. &
          WRITE( PRTFile, * ) 'ray is outside the box where the ocean soundspeed is defined'
          WRITE( PRTFile, * ) ' x = ( r, z ) = ', x
          CALL ERROUT( 'sspMod: Quad', 'ray is outside the box where the soundspeed is defined' )
       END IF

       CALL UpdateRangeSegmentT( x, t )

       ! for this depth, x( 2 ), get the sound speed at both ends of the segment
       cz1 = SSP%czMat( iSegz, iSegr )
       cz2 = SSP%czMat( iSegz, iSegr + 1 )

       s2      = x( 2 ) - SSP%z( iSegz )
       delta_z = SSP%z( iSegz + 1 ) - SSP%z( iSegz )

       c1 = SSP%cMat( iSegz, iSegr     ) + s2 * cz1
       c2 = SSP%cMat( iSegz, iSegr + 1 ) + s2 * cz2

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

       c     = ( 1.0D0 - s1 ) * c1  + s1 * c2

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

       cz  = ( 1.0D0 - s1 ) * cz1 + s1 * cz2

       cr  = ( c2  - c1  ) / delta_r
       crz = ( cz2 - cz1 ) / delta_r

       gradc = [ cr, cz ]
       crr   = 0.0D0
       czz   = 0.0D0

       ! linear interpolation for density
       W   = ( x( 2 ) - 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 Quad