[mpas-developers] review of new operators/module_spline_interpolation.F
Michael Duda
duda at ucar.edu
Mon Feb 14 18:31:05 MST 2011
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
More information about the mpas-developers
mailing list