PCHIP for interpolation of sound speed
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=8), | intent(in) | :: | x(2) | |||
real(kind=8), | intent(in) | :: | t(2) | |||
real(kind=8), | intent(out) | :: | c | |||
real(kind=8), | intent(out) | :: | cimag | |||
real(kind=8), | intent(out) | :: | gradc(2) | |||
real(kind=8), | intent(out) | :: | crr | |||
real(kind=8), | intent(out) | :: | crz | |||
real(kind=8), | intent(out) | :: | czz | |||
real(kind=8), | intent(out) | :: | rho | |||
real(kind=8), | intent(in) | :: | freq | |||
character(len=3), | intent(in) | :: | Task |
SUBROUTINE cPCHIP( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task ) !! PCHIP for interpolation of sound speed ! This implements the new monotone piecewise cubic Hermite interpolating ! polynomial (PCHIP) algorithm for the interpolation of the sound speed c. USE pchipMod 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 REAL (KIND=8) :: xt COMPLEX (KIND=8) :: c_cmplx IF ( Task == 'INI' ) THEN ! read in SSP data Depth = x( 2 ) CALL ReadSSP( Depth, freq ) ! 2 3 ! compute coefficients of std cubic polynomial: c0 + c1*x + c2*x + c3*x ! CALL PCHIP( SSP%z, SSP%c, SSP%NPts, SSP%cCoef, SSP%CSWork ) ELSE ! return SSP info CALL UpdateDepthSegmentT( x, t ) xt = x( 2 ) - SSP%z( iSegz ) c_cmplx = SSP%cCoef( 1, iSegz ) & + ( SSP%cCoef( 2, iSegz ) & + ( SSP%cCoef( 3, iSegz ) & + SSP%cCoef( 4, iSegz ) * xt ) * xt ) * xt c = REAL( c_cmplx ) cimag = AIMAG( c_cmplx ) gradc = [ 0.0D0, REAL( SSP%cCoef( 2, iSegz ) & + ( 2.0D0 * SSP%cCoef( 3, iSegz ) & + 3.0D0 * SSP%cCoef( 4, iSegz ) * xt ) * xt ) ] crr = 0.0D0 crz = 0.0D0 czz = REAL( 2.0D0 * SSP%cCoef( 3, iSegz ) + 6.0D0 * SSP%cCoef( 4, iSegz ) * xt ) 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 ) ! linear interp of density END IF END SUBROUTINE cPCHIP