[mpas-developers] review of new operators/module_spline_interpolation.F
Mark Petersen
mpetersen at lanl.gov
Wed Feb 9 15:19:50 MST 2011
Hi All,
We have been working on implementing various high-order transport
algorithms in the vertical direction (i.e. 1D). During this process we
have developed a new operator module called
"module_spline_interpolation.F" that is attached. This module contains as
set of general 1D cubic spline interpolation functions. While it can be
used as a general interpolator, we are using it to interpolate layer
center data to layer edges in order to compute vertical transport within
the continuity and tracer equations. This functionality is immediately
applicable to the hydrostatic and non-hydrostatic dynamical core, if
desired.
We propose to add this module to the operators directory in the trunk. It
is a completely stand-alone module (i.e. contains no "use" statements) and
includes a test subroutine at the end to show how it is used.
If possible, we would like to commit this to the trunk early next week,
and appreciate any feedback you may have.
Thanks,
Mark and Todd
-------------- next part --------------
module spline_interpolation
implicit none
private
public :: CubicSplineCoefficients, InterpolateCubicSpline, &
InterpolateColumnCubicSpline, IntegrateCubicSpline, &
IntegrateColumnCubicSpline, InterpolateColumnLinear, &
TestInterpolateColumn
! Short Descriptions:
! CubicSplineCoefficients: Compute second derivatives at nodes.
! This must be run before any of the other cubic spine functions.
! InterpolateCubicSpline: Compute cubic spline interpolation at a single point.
! InterpolateColumnCubicSpline: Compute cubic spline interpolation at multiple points.
! IntegrateCubicSpline: Compute a single integral from spline data.
! IntegrateColumnCubicSpline: Compute multiple integrals from spline data.
! InterpolateColumnLinear: Compute linear interpolation at multiple points.
! TestInterpolateColumn: Test spline interpolation subroutines.
contains
subroutine CubicSplineCoefficients(x,y,n,yp1,ypn,y2)
! Given arrays x(1:n) and y(1:n) containing a tabulated function,
! i.e., y(i) = f(x(i)), with x(1) < x(2) < ... < x(n), and given values
! yp1 and ypn for the first derivative of the interpolating function
! at points 1 and n, respectively, this routine returns an array
! y2(1:n) that contains the second derivatives of the interpolating
! function at the tabulated points x(1:n). If yp1 and/or ypn are equal
! to 1e30 or larger, the routine is signaled to set the corresponding
! boundary condition for a natural spline, with zero second derivative
! on that boundary.
! This routine adapted from Numerical Recipes in Fortran, 2nd Ed, p. 109
! and was originally named spline.
! INPUT PARAMETERS:
integer, intent(in) :: &
n ! number of nodes
real(kind=RKIND), intent(in), dimension(n) :: &
x, &! location of nodes
y ! value at nodes
real(kind=RKIND), intent(in) :: &
yp1, &! first derivative of y at first endpoint
ypn ! first derivative of y at last endpoint
! OUTPUT PARAMETERS:
real(kind=RKIND), intent(out), dimension(n) :: &
y2 ! dy^2/dx^2 at each node
! local variables:
integer, parameter :: &
NMAX=500
integer :: &
i,k
real(kind=RKIND) :: &
p,qn,sig,un,u(NMAX)
if (yp1.gt..99e30) then
y2(1)=0.0
u(1)=0.0
else
y2(1)=-0.5
u(1)=(3.0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
endif
do i=2,n-1
sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
p=sig*y2(i-1)+2.0
y2(i)=(sig-1.0)/p
u(i)=(6.0*((y(i+1)-y(i))/(x(i+1)-x(i)) &
-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1)) &
-sig*u(i-1))/p
enddo
if (ypn.gt..99e30) then
qn=0.0
un=0.0
else
qn=0.5
un=(3.0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
endif
y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.0)
do k=n-1,1,-1
y2(k)=y2(k)*y2(k+1)+u(k)
enddo
end subroutine CubicSplineCoefficients
subroutine InterpolateCubicSpline( &
x,y,y2,n, &
xOut,yOut)
! Given the arrays x(1:n) and y(1:n), which tabulate a function
! (with x(1) < x(2) < ... < x(n)), and given the array y2(1:n),
! which is the output from CubicSplineCoefficients above, and given
! a single value of xOut, this routine returns a cubic-spline interpolated
! value yOut at xOut.
! This routine adapted from Numerical Recipes in Fortran, 2nd Ed, p 110
! and was originally named splint, for SPLine INTerpolation
! INPUT PARAMETERS:
integer, intent(in) :: &
n
real(kind=RKIND), intent(in), dimension(n) :: &
x,y,y2
real(kind=RKIND), intent(in) :: &
xOut
! OUTPUT PARAMETERS:
real(kind=RKIND), intent(out) :: &
yOut
! local variables:
integer :: &
k,khi,klo
real(kind=RKIND) :: &
a,b,h
klo=1
khi=n
do while (khi-klo.gt.1)
k=(khi+klo)/2
if(x(k).gt.xOut)then
khi=k
else
klo=k
endif
enddo
h=x(khi)-x(klo)
! if (h.eq.0.) pause 'bad x input in splint'
a=(x(khi)-xOut)/h
b=(xOut-x(klo))/h
yOut=a*y(klo) + b*y(khi) &
+ ((a**3-a)*y2(klo)+(b**3-b)*y2(khi))*(h**2)/6.0
end subroutine InterpolateCubicSpline
subroutine InterpolateColumnCubicSpline( &
x,y,y2,n, &
xOut,yOut,nOut)
! Given the arrays x(1:n) and y(1:n), which tabulate a function,
! and given the array y2(1:n), which is the output from
! CubicSplineCoefficients above, this routine returns the
! cubic-spline interpolated values of yOut(1:nOut) at xOut(1:nOut).
! This subroutine assumes that both x and xOut are monotonically
! increasing, and that all values of xOut are within the first and
! last values of x.
! This routine adapted from Numerical Recipes in Fortran, 2nd Ed, p 110
! INPUT PARAMETERS:
real (kind=RKIND), dimension(n), intent(in) :: &
x, &! node location, input grid
y, &! interpolation variable, input grid
y2 ! 2nd derivative of y at nodes
real (kind=RKIND), dimension(nOut), intent(in) :: &
xOut ! node location, output grid
integer, intent(in) :: &
n, &! number of nodes, input grid
nOut ! number of nodes, output grid
! OUTPUT PARAMETERS:
real (kind=RKIND), dimension(nOut), intent(out) :: &
yOut ! interpolation variable, output grid
! local variables:
integer :: &
kIn, kOut ! counters
real (kind=RKIND) :: &
a, b, h
kOut = 1
kInLoop: do kIn = 1,n-1
h = x(kIn+1)-x(kIn)
do while(xOut(kOut) < x(kIn+1))
a = (x(kIn+1)-xOut(kOut))/h
b = (xOut(kOut)-x (kIn) )/h
yOut(kOut) = a*y(kIn) + b*y(kIn+1) &
+ ((a**3-a)*y2(kIn) + (b**3-b)*y2(kIn+1)) &
*(h**2)/6.0
kOut = kOut + 1
if (kOut>nOut) exit kInLoop
enddo
enddo kInLoop
end subroutine InterpolateColumnCubicSpline
subroutine IntegrateCubicSpline(x,y,y2,n,x1,x2,y_integral)
! Given the arrays x(1:n) and y(1:n), which tabulate a function,
! and given the array y2(1:n), which is the output from
! CubicSplineCoefficients above, this routine returns y_integral,
! the integral of y from x1 to x2. The integration formula was
! created by analytically integrating a cubic spline between each node.
! This subroutine assumes that x is monotonically increasing, and
! that x1 < x2.
! INPUT PARAMETERS:
integer, intent(in) :: &
n ! number of nodes
real(kind=RKIND), intent(in), dimension(n) :: &
x, &! location of nodes
y, &! value at nodes
y2 ! dy^2/dx^2 at each node
real(kind=RKIND), intent(in) :: &
x1,x2 ! limits of integration
! OUTPUT PARAMETERS:
real(kind=RKIND), intent(out) :: &
y_integral ! integral of y
! local variables:
integer :: i,j,k
real(kind=RKIND) :: h,h2, A2,B2, F1,F2, eps1
if (x1<x(1).or.x2>x(n).or.x1>x2) then
print *, 'error on integration bounds'
endif
y_integral = 0.0
eps1 = 1e-14*x2
do j=1,n-1 ! loop through sections
! section x(j) ... x(j+1)
if (x2<=x(j) +eps1) exit
if (x1>=x(j+1)-eps1) cycle
h = x(j+1) - x(j)
h2 = h**2
! left side:
if (x1<x(j)) then
F1 = -y(j)*h*0.5 + y2(j)*h**3/24.0
else
A2 = (x(j+1)-x1 )**2/h2
B2 = (x1 -x(j))**2/h2
F1 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ y2(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ y2(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
endif
! right side:
if (x2>x(j+1)) then
F2 = y(j+1)*h*0.5 - y2(j+1)*h**3/24.0
else
A2 = (x(j+1)-x2 )**2/h2
B2 = (x2 -x(j))**2/h2
F2 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ y2(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ y2(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
endif
y_integral = y_integral + F2 - F1
enddo ! j
end subroutine IntegrateCubicSpline
subroutine IntegrateColumnCubicSpline( &
x,y,y2,n, &
xOut,y_integral, nOut)
! Given the arrays x(1:n) and y(1:n), which tabulate a function,
! and given the array y2(1:n), which is the output from
! CubicSplineCoefficients above, this routine returns
! y_integral(1:nOut), the integral of y.
! This is a cumulative integration, so that
! y_integral(j) holds the integral of y from x(1) to xOut(j).
! The integration formula was created by analytically integrating a
! cubic spline between each node.
! This subroutine assumes that both x and xOut are monotonically
! increasing, and that all values of xOut are within the first and
! INPUT PARAMETERS:
integer, intent(in) :: &
n, &! number of nodes
nOut ! number of output locations to compute integral
real(kind=RKIND), intent(in), dimension(n) :: &
x, &! location of nodes
y, &! value at nodes
y2 ! dy^2/dx^2 at each node
real(kind=RKIND), dimension(nOut), intent(in) :: &
xOut ! output locations to compute integral
! OUTPUT PARAMETERS:
real(kind=RKIND), dimension(nOut), intent(out) :: &
y_integral ! integral from 0 to xOut
! local variables:
integer :: i,j,k
real(kind=RKIND) :: h,h2, A2,B2, F1,F2, eps1
y_integral = 0.0
j = 1
h = x(j+1) - x(j)
h2 = h**2
F1 = -y(j)*h*0.5 + y2(j)*h**3/24.0
eps1 = 0.0 ! note: could use 1e-12*xOut(nOut)
k_loop: do k = 1,nOut
if (k>1) y_integral(k) = y_integral(k-1)
do while(xOut(k) > x(j+1)-eps1)
F2 = y(j+1)*h*0.5 - y2(j+1)*h**3/24.0
y_integral(k) = y_integral(k) + F2 - F1
j = j+1
h = x(j+1) - x(j)
h2 = h**2
F1 = -y(j)*h*0.5 + y2(j)*h**3/24.0
if (abs(xOut(k) - x(j+1))<eps1) cycle k_loop
enddo
A2 = (x(j+1) - xOut(k))**2/h2
B2 = (xOut(k) - x(j) )**2/h2
F2 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ y2(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ y2(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
y_integral(k) = y_integral(k) + F2 - F1
if (k < nOut) then
A2 = (x(j+1) -xOut(k))**2/h2
B2 = (xOut(k) -x(j) )**2/h2
F1 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ y2(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ y2(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
endif
enddo k_loop
end subroutine IntegrateColumnCubicSpline
subroutine InterpolateColumnLinear( &
x,y,n, &
xOut,yOut,nOut)
! Given the arrays x(1:n) and y(1:n), which tabulate a function,
! this routine returns the linear interpolated values of yOut(1:nOut)
! at xOut(1:nOut).
! This subroutine assumes that both x and xOut are monotonically
! increasing, and that all values of xOut are within the first and
! last values of x.
! !INPUT PARAMETERS:
real (kind=RKIND), dimension(n), intent(in) :: &
x, &! node location, input grid
y ! interpolation variable, input grid
real (kind=RKIND), dimension(nOut), intent(in) :: &
xOut ! node location, output grid
integer, intent(in) :: &
N, &! number of nodes, input grid
NOut ! number of nodes, output grid
! !OUTPUT PARAMETERS:
real (kind=RKIND), dimension(nOut), intent(out) :: &
yOut ! interpolation variable, output grid
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer :: &
kIn, kOut ! counters
kOut = 1
kInLoop: do kIn = 1,n-1
do while(xOut(kOut) < x(kIn+1))
yOut(kOut) = y(kIn) &
+ (y(kIn+1)-y(kIn)) &
/(x(kIn+1) -x(kIn) ) &
*(xOut(kOut) -x(kIn) )
kOut = kOut + 1
if (kOut>nOut) exit kInLoop
enddo
enddo kInLoop
end subroutine InterpolateColumnLinear
subroutine TestInterpolateColumn
! Test function to show how to operate the cubic spline subroutines
integer, parameter :: &
n = 10
real (kind=RKIND), dimension(n) :: &
y, x, y2
integer, parameter :: &
nOut = 100
real (kind=RKIND), dimension(nOut) :: &
yOut, xOut
integer :: &
k
!-----------------------------------------------------------------------
!
! Create x, y, xOut
!
!-----------------------------------------------------------------------
do k=1,n
x(k) = k-4
! power function:
!y(k) = x(k)**2
! trig function:
y(k) = sin(x(k)/2)
enddo
do k=1,nOut
xOut(k) = x(1) + k/(nOut+1.0)*(x(n)-x(1))
enddo
!-----------------------------------------------------------------------
!
! Interpolate
!
!-----------------------------------------------------------------------
! First, compute second derivative values at each node, y2.
call CubicSplineCoefficients(x,y,n, &
2.0e30, 2.0e30, y2)
! Compute interpolated values yOut.
call InterpolateColumnCubicSpline( &
x,y,y2,n, &
xOut,yOut,nOut)
!-----------------------------------------------------------------------
!
! Output
!
!-----------------------------------------------------------------------
! The following output can be copied directly into Matlab
print '(a,10f8.4,a)', 'x = [',x,'];'
print '(a,10f8.4,a)', 'y = [',y,'];'
print '(a,100f8.4,a)', 'xOut = [',xOut,'];'
print '(a,100f8.4,a)', 'yOut = [',yOut,'];'
print *, "plot(x,y,'-*r',xOut,yOut,'x')"
end subroutine TestInterpolateColumn
end module spline_interpolation
More information about the mpas-developers
mailing list