[mpas-developers] review of new operators/module_spline_interpolation.F
Mark Petersen
mpetersen at lanl.gov
Tue Feb 15 09:01:32 MST 2011
Xylar and Michael,
Thanks for the copyright statement. Only two of the seven subroutines are
from NR, while the other five (column and integral forms) are fully
derived and coded by me.
One from NR, InterpolateCubicSpline, is a single-point version of my
InterpolateColumnCubicSpline, so is redundant and can be deleted. I will
alter the other, CubicSplineCoefficients, so it is not similar to NR.
Mark
On Tue, 15 Feb 2011, Xylar Asay-Davis wrote:
> Hi Michael and Mark,
>
> Here's a quote form the NR copyright page:
> http://www.nr.com/com/info-copyright.html
>
> "Of course, the mathematical algorithms 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 expression 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."
>
> 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 :)
>
> -Xylar
>
> On 2/14/11 6:31 PM, Michael Duda wrote:
>
> 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:
>
> 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
>
> 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
>
>
> _______________________________________________
> mpas-developers mailing list
> mpas-developers at mailman.ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/mpas-developers
>
> _______________________________________________
> mpas-developers mailing list
> mpas-developers at mailman.ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/mpas-developers
>
>
>
> --
>
> ***********************
> Xylar S. Asay-Davis
> E-mail: xylar at lanl.gov
> Phone: (505) 606-0025
> Fax: (505) 665-2659
> CNLS, MS B258
> Los Alamos National Laboratory
> Los Alamos, NM 87545
> ***********************
>
>
>
More information about the mpas-developers
mailing list