<p><b>duda</b> 2011-06-29 14:02:04 -0600 (Wed, 29 Jun 2011)</p><p>BRANCH COMMIT<br>
<br>
Import spline interpolation module from trunk.<br>
<br>
<br>
A src/operators/module_spline_interpolation.F<br>
M src/operators/Makefile<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/operators/Makefile
===================================================================
--- branches/atmos_physics/src/operators/Makefile 2011-06-29 18:19:49 UTC (rev 907)
+++ branches/atmos_physics/src/operators/Makefile 2011-06-29 20:02:04 UTC (rev 908)
@@ -1,6 +1,6 @@
.SUFFIXES: .F .o
-OBJS = module_RBF_interpolation.o module_vector_reconstruction.o
+OBJS = module_RBF_interpolation.o module_vector_reconstruction.o module_spline_interpolation.o
all: operators
@@ -9,6 +9,7 @@
module_vector_reconstruction.o: module_RBF_interpolation.o
module_RBF_interpolation.o:
+module_spline_interpolation:
clean:
$(RM) *.o *.mod *.f90 libops.a
Copied: branches/atmos_physics/src/operators/module_spline_interpolation.F (from rev 847, trunk/mpas/src/operators/module_spline_interpolation.F)
===================================================================
--- branches/atmos_physics/src/operators/module_spline_interpolation.F (rev 0)
+++ branches/atmos_physics/src/operators/module_spline_interpolation.F 2011-06-29 20:02:04 UTC (rev 908)
@@ -0,0 +1,427 @@
+module spline_interpolation
+
+ implicit none
+
+ private
+
+ public :: CubicSplineCoefficients, InterpolateCubicSpline, &
+ IntegrateCubicSpline, IntegrateColumnCubicSpline, InterpolateLinear, &
+ TestInterpolate
+
+! 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.
+
+! IntegrateCubicSpline: Compute a single integral from spline data.
+
+! IntegrateColumnCubicSpline: Compute multiple integrals from spline data.
+
+! InterpolateLinear: Compute linear interpolation.
+
+! TestInterpolate: Test spline interpolation subroutines.
+
+ contains
+
+ subroutine CubicSplineCoefficients(x,y,n,y2ndDer)
+
+! Given arrays x(1:n) and y(1:n) containing a function,
+! i.e., y(i) = f(x(i)), with x monotonically increasing
+! this routine returns an array y2ndDer(1:n) that contains
+! the second derivatives of the interpolating function at x(1:n).
+! This routine uses boundary conditions for a natural spline,
+! with zero second derivative on that boundary.
+
+! INPUT PARAMETERS:
+
+ integer, intent(in) :: &
+ n ! number of nodes
+ real(kind=RKIND), intent(in), dimension(n) :: &
+ x, &! location of nodes
+ y ! value at nodes
+
+! OUTPUT PARAMETERS:
+
+ real(kind=RKIND), intent(out), dimension(n) :: &
+ y2ndDer ! dy^2/dx^2 at each node
+
+! local variables:
+
+ integer :: i
+ real(kind=RKIND) :: &
+ temp,xRatio,a(n)
+
+ y2ndDer(1)=0.0
+ y2ndDer(n)=0.0
+ a(1)=0.0
+
+ do i=2,n-1
+ xRatio=(x(i)-x(i-1))/(x(i+1)-x(i-1))
+ temp=1.0/(2.0+xRatio*y2ndDer(i-1))
+ y2ndDer(i)=temp*(xRatio-1.0)
+ a(i) = temp*(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)) &
+ -xRatio*a(i-1))
+ enddo
+
+ do i=n-1,1,-1
+ y2ndDer(i)=y2ndDer(i)*y2ndDer(i+1)+a(i)
+ enddo
+
+ end subroutine CubicSplineCoefficients
+
+
+ subroutine InterpolateCubicSpline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
+
+! Given the arrays x(1:n) and y(1:n), which tabulate a function,
+! and given the array y2ndDer(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.
+
+! INPUT PARAMETERS:
+
+ real (kind=RKIND), dimension(n), intent(in) :: &
+ x, &! node location, input grid
+ y, &! interpolation variable, input grid
+ y2ndDer ! 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)*y2ndDer(kIn) + (b**3-b)*y2ndDer(kIn+1)) &
+ *(h**2)/6.0
+
+ kOut = kOut + 1
+
+ if (kOut>nOut) exit kInLoop
+
+ enddo
+
+ enddo kInLoop
+
+end subroutine InterpolateCubicSpline
+
+
+subroutine IntegrateCubicSpline(x,y,y2ndDer,n,x1,x2,y_integral)
+
+! Given the arrays x(1:n) and y(1:n), which tabulate a function,
+! and given the array y2ndDer(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
+ y2ndDer ! 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 + y2ndDer(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 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(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 - y2ndDer(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 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(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,y2ndDer,n, &
+ xOut,y_integral, nOut)
+
+! Given the arrays x(1:n) and y(1:n), which tabulate a function,
+! and given the array y2ndDer(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
+ y2ndDer ! 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 + y2ndDer(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 - y2ndDer(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 + y2ndDer(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 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(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 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
+ endif
+
+ enddo k_loop
+
+ end subroutine IntegrateColumnCubicSpline
+
+
+ subroutine InterpolateLinear( &
+ 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 InterpolateLinear
+
+
+ subroutine TestInterpolate
+
+! Test function to show how to operate the cubic spline subroutines
+
+ integer, parameter :: &
+ n = 10
+ real (kind=RKIND), dimension(n) :: &
+ y, x, y2ndDer
+
+ 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
+ ! 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, y2ndDer.
+ call CubicSplineCoefficients(x,y,n,y2ndDer)
+
+ ! Compute interpolated values yOut.
+ call InterpolateCubicSpline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
+
+ ! The following output can be copied directly into Matlab
+ print *, 'subplot(2,1,1)'
+ 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')"
+
+ ! Compute interpolated values yOut.
+ call IntegrateColumnCubicSpline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
+
+ ! The following output can be copied directly into Matlab
+ print *, 'subplot(2,1,2)'
+ print '(a,10f8.4,a)', 'x = [',x,'];'
+ print '(a,10f8.4,a)', 'y = 2*cos(-3/2) -2*cos(x/2);'
+ print '(a,100f8.4,a)', 'xOut = [',xOut,'];'
+ print '(a,100f8.4,a)', 'yOut = [',yOut,'];'
+ print *, "plot(x,y,'-*r',xOut,yOut,'x')"
+
+ end subroutine TestInterpolate
+
+end module spline_interpolation
+
</font>
</pre>