Cubic spline interpolation for 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 cCubic( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task ) !! Cubic spline interpolation for sound speed 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 :: iBCBeg, iBCEnd REAL (KIND=8) :: hSpline COMPLEX (KIND=8) :: c_cmplx, cz_cmplx, czz_cmplx IF ( Task == 'INI' ) THEN ! *** Task 'INIT' for initialization *** Depth = x( 2 ) CALL ReadSSP( Depth, freq ) SSP%cSpline( 1, 1 : SSP%NPts ) = SSP%c( 1 : SSP%NPts ) ! Compute spline coefs iBCBeg = 0 iBCEnd = 0 CALL CSpline( SSP%z, SSP%cSpline( 1, 1 ), SSP%NPts, iBCBeg, iBCEnd, SSP%NPts ) ELSE ! *** Section to return SSP info *** CALL UpdateDepthSegmentT( x, t ) hSpline = x( 2 ) - SSP%z( iSegz ) ! c = Spline( SSP%cSpline( 1, iSegz ), hSpline ) ! cz = SplineX( SSP%cSpline( 1, iSegz ), hSpline ) ! czz = SplineXX( SSP%cSpline( 1, iSegz ), hSpline ) CALL SplineALL( SSP%cSpline( 1, iSegz ), hSpline, c_cmplx, cz_cmplx, czz_cmplx ) c = DBLE( c_cmplx ) cimag = AIMAG( c_cmplx ) gradc = [ 0.0D0, DBLE( cz_cmplx ) ] czz = DBLE( czz_cmplx ) crr = 0.0d0 crz = 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 cCubic