[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