<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
</head>
<body text="#000000" bgcolor="#ffffff">
Hi Michael and Mark,<br>
<br>
Here's a quote form the NR copyright page:<br>
<a class="moz-txt-link-freetext" href="http://www.nr.com/com/info-copyright.html">http://www.nr.com/com/info-copyright.html</a><br>
<br>
"Of course, the mathematical <i>algorithms</i> that underly the
programs are
not copyrightable. So, the test of whether a similar program is
actually infringing on the Numerical Recipes copyright is whether
the
<i>expression</i> of the underlying algorithm is so similar to that
in
Numerical Recipes as to make clear that the infringing program was
in
fact derived from the Numerical Recipes program. It is usually very
easy to tell if something was copied from Numerical Recipes, even if
variable names are changed (for example), since all programs are
full
of arbitrary stylistic choices as to the order in which things are
done, the way expressions are written, etc."<br>
<br>
Whenever I've used algorithms from NR, I've made sure that I wrote
code that was stylistically my own. That seems to be the important
distinction. One way this has worked for me is to work from NR in
C++ and rewrite the code in F90, though I guess that's not the
preferred option for any of you :) <br>
<br>
-Xylar<br>
<br>
On 2/14/11 6:31 PM, Michael Duda wrote:
<blockquote cite="mid:20110215013105.GA46267@ucar.edu" type="cite">
<pre wrap="">Hi, Mark.
I think the module looks great, and we'd definitely be able to make use
of this in the atmosphere cores, too. My only concern lies with the
copyright of the code in Numerical Recipes that this module is based on.
The authors of NR seem keen on selling licenses to use the code from
their book, whether the code is typed in or acquired electronically, and
they consider any modified versions of the code to be derivative works
falling under the same copyright. So, I wonder whether we could run into
problems by including code from NR in MPAS, even if it is for
non-commercial use.
Michael
On Wed, Feb 09, 2011 at 03:19:50PM -0700, Mark Petersen wrote:
</pre>
<blockquote type="cite">
<pre wrap="">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
</pre>
</blockquote>
<pre wrap="">
</pre>
<blockquote type="cite">
<pre wrap="">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
</pre>
</blockquote>
<pre wrap="">
</pre>
<blockquote type="cite">
<pre wrap="">_______________________________________________
mpas-developers mailing list
<a class="moz-txt-link-abbreviated" href="mailto:mpas-developers@mailman.ucar.edu">mpas-developers@mailman.ucar.edu</a>
<a class="moz-txt-link-freetext" href="http://mailman.ucar.edu/mailman/listinfo/mpas-developers">http://mailman.ucar.edu/mailman/listinfo/mpas-developers</a>
</pre>
</blockquote>
<pre wrap="">
_______________________________________________
mpas-developers mailing list
<a class="moz-txt-link-abbreviated" href="mailto:mpas-developers@mailman.ucar.edu">mpas-developers@mailman.ucar.edu</a>
<a class="moz-txt-link-freetext" href="http://mailman.ucar.edu/mailman/listinfo/mpas-developers">http://mailman.ucar.edu/mailman/listinfo/mpas-developers</a>
</pre>
</blockquote>
<br>
<br>
<pre class="moz-signature" cols="72">--
***********************
Xylar S. Asay-Davis
E-mail: <a class="moz-txt-link-abbreviated" href="mailto:xylar@lanl.gov">xylar@lanl.gov</a>
Phone: (505) 606-0025
Fax: (505) 665-2659
CNLS, MS B258
Los Alamos National Laboratory
Los Alamos, NM 87545
***********************
</pre>
</body>
</html>