[mpas-developers] review of new operators/module_spline_interpolation.F
Xylar Asay-Davis
xylar at lanl.gov
Tue Feb 15 08:48:26 MST 2011
Also, I would say something like "adapted from algorithms in Numerical
Recipes Section X.X" rather than referencing the code in NR.
On 2/15/11 8:46 AM, 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
> ***********************
>
--
***********************
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
***********************
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/mpas-developers/attachments/20110215/95642296/attachment-0001.html
More information about the mpas-developers
mailing list