<p><b>mpetersen@lanl.gov</b> 2011-02-22 10:48:59 -0700 (Tue, 22 Feb 2011)</p><p>Altered NR-like sections on spline interpolations. Changed some spacing in module_time_integration.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-02-18 18:32:25 UTC (rev 740)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-02-22 17:48:59 UTC (rev 741)
@@ -893,186 +893,189 @@
!
! tracer tendency: vertical advection term -d/dz( h \phi w)
!
+!mrp 110222
+ call TestInterpolate
+ stop
+!mrp 110222
if (config_vert_grid_type.eq.'zlevel') then
- allocate(tracerTop(num_tracers,nVertLevels+1,nCells))
+ allocate(tracerTop(num_tracers,nVertLevels+1,nCells))
- ! Tracers at the top and bottom boundary are assigned nearest
- ! cell-centered value, regardless of tracer interpolation method.
- ! wTop=0 at top and bottom sets the boundary condition.
- do iCell=1,nCellsSolve
- do iTracer=1,num_tracers
- tracerTop(iTracer,1,iCell) = tracers(iTracer,1,iCell)
- tracerTop(iTracer,maxLevelCell(iCell)+1,iCell) = &
+ ! Tracers at the top and bottom boundary are assigned nearest
+ ! cell-centered value, regardless of tracer interpolation method.
+ ! wTop=0 at top and bottom sets the boundary condition.
+ do iCell=1,nCellsSolve
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,1,iCell) = tracers(iTracer,1,iCell)
+ tracerTop(iTracer,maxLevelCell(iCell)+1,iCell) = &
tracers(iTracer,maxLevelCell(iCell),iCell)
+ end do
end do
- end do
- if (config_vert_tracer_adv.eq.'stencil'.and. &
- config_vert_tracer_adv_order.eq.2) then
+ if (config_vert_tracer_adv.eq.'stencil'.and. &
+ config_vert_tracer_adv_order.eq.2) then
- ! Compute tracerTop using centered stencil, a simple average.
+ ! Compute tracerTop using centered stencil, a simple average.
- do iCell=1,nCellsSolve
- do k=2,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tracerTop(iTracer,k,iCell) = &
- ( tracers(iTracer,k-1,iCell) &
- +tracers(iTracer,k ,iCell))/2.0
+ do iCell=1,nCellsSolve
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ ( tracers(iTracer,k-1,iCell) &
+ +tracers(iTracer,k ,iCell))/2.0
+ end do
end do
end do
- end do
- elseif (config_vert_tracer_adv.eq.'stencil'.and. &
- config_vert_tracer_adv_order.eq.3) then
+ elseif (config_vert_tracer_adv.eq.'stencil'.and. &
+ config_vert_tracer_adv_order.eq.3) then
- ! Compute tracerTop using 3rd order stencil. This is the same
- ! as 4th order, but includes upwinding.
+ ! Compute tracerTop using 3rd order stencil. This is the same
+ ! as 4th order, but includes upwinding.
- ! Hardwire flux3Coeff at 1.0 for now. Could add this to the
- ! namelist, if desired.
- flux3Coef = 1.0
- do iCell=1,nCellsSolve
- k=2
+ ! Hardwire flux3Coeff at 1.0 for now. Could add this to the
+ ! namelist, if desired.
+ flux3Coef = 1.0
+ do iCell=1,nCellsSolve
+ k=2
do iTracer=1,num_tracers
tracerTop(iTracer,k,iCell) = &
hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
end do
- do k=3,maxLevelCell(iCell)-1
- cSignWTop = sign(flux3Coef,wTop(k,iCell))
- do iTracer=1,num_tracers
- tracerTop(iTracer,k,iCell) = &
- ( (-1.+ cSignWTop)*tracers(iTracer,k-2,iCell) &
- +( 7.-3.*cSignWTop)*tracers(iTracer,k-1,iCell) &
- +( 7.+3.*cSignWTop)*tracers(iTracer,k ,iCell) &
- +(-1.- cSignWTop)*tracers(iTracer,k+1,iCell) &
- )/12.
+ do k=3,maxLevelCell(iCell)-1
+ cSignWTop = sign(flux3Coef,wTop(k,iCell))
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ ( (-1.+ cSignWTop)*tracers(iTracer,k-2,iCell) &
+ +( 7.-3.*cSignWTop)*tracers(iTracer,k-1,iCell) &
+ +( 7.+3.*cSignWTop)*tracers(iTracer,k ,iCell) &
+ +(-1.- cSignWTop)*tracers(iTracer,k+1,iCell) &
+ )/12.
+ end do
end do
+ k=maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
end do
- k=maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tracerTop(iTracer,k,iCell) = &
- hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
- + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
- end do
- end do
- elseif (config_vert_tracer_adv.eq.'stencil'.and. &
- config_vert_tracer_adv_order.eq.4) then
+ elseif (config_vert_tracer_adv.eq.'stencil'.and. &
+ config_vert_tracer_adv_order.eq.4) then
- ! Compute tracerTop using 4rd order stencil [-1 7 7 -1]
+ ! Compute tracerTop using 4rd order stencil [-1 7 7 -1]
- do iCell=1,nCellsSolve
- k=2
- do iTracer=1,num_tracers
- tracerTop(iTracer,k,iCell) = &
- hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
- + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ do iCell=1,nCellsSolve
+ k=2
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ do k=3,maxLevelCell(iCell)-1
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ (- tracers(iTracer,k-2,iCell) &
+ +7.*tracers(iTracer,k-1,iCell) &
+ +7.*tracers(iTracer,k ,iCell) &
+ - tracers(iTracer,k+1,iCell) &
+ )/12.
+ end do
end do
- do k=3,maxLevelCell(iCell)-1
- do iTracer=1,num_tracers
- tracerTop(iTracer,k,iCell) = &
- (- tracers(iTracer,k-2,iCell) &
- +7.*tracers(iTracer,k-1,iCell) &
- +7.*tracers(iTracer,k ,iCell) &
- - tracers(iTracer,k+1,iCell) &
- )/12.
- end do
+ k=maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k,iCell) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
end do
- k=maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tracerTop(iTracer,k,iCell) = &
- hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
- + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
- end do
- end do
- elseif (config_vert_tracer_adv.eq.'spline'.and. &
- config_vert_tracer_adv_order.eq.2) then
+ elseif (config_vert_tracer_adv.eq.'spline'.and. &
+ config_vert_tracer_adv_order.eq.2) then
- ! Compute tracerTop using linear interpolation.
+ ! Compute tracerTop using linear interpolation.
- do iCell=1,nCellsSolve
- do k=2,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! Note hRatio on the k side is multiplied by tracer at k-1
- ! and hRatio on the Km1 (k-1) side is mult. by tracer at k.
- tracerTop(iTracer,k,iCell) = &
- hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
- + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ do iCell=1,nCellsSolve
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ ! Note hRatio on the k side is multiplied by tracer at k-1
+ ! and hRatio on the Km1 (k-1) side is mult. by tracer at k.
+ tracerTop(iTracer,k,iCell) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
end do
end do
- end do
- elseif (config_vert_tracer_adv.eq.'spline'.and. &
- config_vert_tracer_adv_order.eq.3) then
+ elseif (config_vert_tracer_adv.eq.'spline'.and. &
+ config_vert_tracer_adv_order.eq.3) then
- ! Compute tracerTop using cubic spline interpolation.
+ ! Compute tracerTop using cubic spline interpolation.
- allocate(tracer2ndDer(nVertLevels))
- allocate(tracersIn(nVertLevels),tracersOut(nVertLevels), &
- posZMidZLevel(nVertLevels), posZTopZLevel(nVertLevels-1))
+ allocate(tracer2ndDer(nVertLevels))
+ allocate(tracersIn(nVertLevels),tracersOut(nVertLevels), &
+ posZMidZLevel(nVertLevels), posZTopZLevel(nVertLevels-1))
- ! For the ocean, zlevel coordinates are negative and decreasing,
- ! but spline functions assume increasing, so flip to positive.
+ ! For the ocean, zlevel coordinates are negative and decreasing,
+ ! but spline functions assume increasing, so flip to positive.
- posZMidZLevel = -zMidZLevel(1:nVertLevels)
- posZTopZLevel = -zTopZLevel(2:nVertLevels)
+ posZMidZLevel = -zMidZLevel(1:nVertLevels)
+ posZTopZLevel = -zTopZLevel(2:nVertLevels)
- do iCell=1,nCellsSolve
- ! mrp 110201 efficiency note: push tracer loop down
- ! into spline subroutines to improve efficiency
- do iTracer=1,num_tracers
+ do iCell=1,nCellsSolve
+ ! mrp 110201 efficiency note: push tracer loop down
+ ! into spline subroutines to improve efficiency
+ do iTracer=1,num_tracers
- ! Place data in arrays to avoid creating new temporary arrays for every
- ! subroutine call.
- tracersIn(1:maxLevelCell(iCell))=tracers(iTracer,1:maxLevelCell(iCell),iCell)
+ ! Place data in arrays to avoid creating new temporary arrays for every
+ ! subroutine call.
+ tracersIn(1:maxLevelCell(iCell))=tracers(iTracer,1:maxLevelCell(iCell),iCell)
- call CubicSplineCoefficients(posZMidZLevel, &
- tracersIn, &
- maxLevelCell(iCell), 2.0e30, 2.0e30, tracer2ndDer)
+ call CubicSplineCoefficients(posZMidZLevel, &
+ tracersIn, maxLevelCell(iCell), tracer2ndDer)
- call InterpolateColumnCubicSpline( &
- posZMidZLevel, tracersIn, tracer2ndDer, maxLevelCell(iCell), &
- posZTopZLevel, tracersOut, maxLevelCell(iCell)-1 )
+ call InterpolateCubicSpline( &
+ posZMidZLevel, tracersIn, tracer2ndDer, maxLevelCell(iCell), &
+ posZTopZLevel, tracersOut, maxLevelCell(iCell)-1 )
- tracerTop(iTracer,2:maxLevelCell(iCell),iCell) = tracersOut(1:maxLevelCell(iCell)-1)
+ tracerTop(iTracer,2:maxLevelCell(iCell),iCell) = tracersOut(1:maxLevelCell(iCell)-1)
+ end do
end do
- end do
- deallocate(tracer2ndDer)
- deallocate(tracersIn,tracersOut, posZMidZLevel, posZTopZLevel)
+ deallocate(tracer2ndDer)
+ deallocate(tracersIn,tracersOut, posZMidZLevel, posZTopZLevel)
- else
+ else
- print *, 'Abort: Incorrect combination of ', &
- 'config_vert_tracer_adv and config_vert_tracer_adv_order.'
- print *, 'Use:'
- print *, 'config_vert_tracer_adv=''stencil'' and_vert_tracer_adv_order=2,3,4 or'
- print *, 'config_vert_tracer_adv=''spline'' and_vert_tracer_adv_order=2,3'
- call dmpar_abort(dminfo)
+ print *, 'Abort: Incorrect combination of ', &
+ 'config_vert_tracer_adv and config_vert_tracer_adv_order.'
+ print *, 'Use:'
+ print *, 'config_vert_tracer_adv=''stencil'' and_vert_tracer_adv_order=2,3,4 or'
+ print *, 'config_vert_tracer_adv=''spline'' and_vert_tracer_adv_order=2,3'
+ call dmpar_abort(dminfo)
- endif ! vertical tracer advection method
+ endif ! vertical tracer advection method
- do iCell=1,nCellsSolve
- ! For top layer, contribution is from bottom only.
- k=1
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ do iCell=1,nCellsSolve
+ ! For top layer, contribution is from bottom only.
+ k=1
+ do iTracer=1,num_tracers
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ wTop(k+1,iCell)*tracerTop(iTracer,k+1,iCell)
+ end do
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ - ( wTop(k ,iCell)*tracerTop(iTracer,k ,iCell) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1,iCell))
+ end do
end do
- do k=2,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( wTop(k ,iCell)*tracerTop(iTracer,k ,iCell) &
- - wTop(k+1,iCell)*tracerTop(iTracer,k+1,iCell))
- end do
end do
- end do
- deallocate(tracerTop)
+ deallocate(tracerTop)
endif ! ZLevel
Modified: branches/ocean_projects/vert_adv_mrp/src/operators/module_spline_interpolation.F
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/operators/module_spline_interpolation.F 2011-02-18 18:32:25 UTC (rev 740)
+++ branches/ocean_projects/vert_adv_mrp/src/operators/module_spline_interpolation.F 2011-02-22 17:48:59 UTC (rev 741)
@@ -5,44 +5,35 @@
private
public :: CubicSplineCoefficients, InterpolateCubicSpline, &
- InterpolateColumnCubicSpline, IntegrateCubicSpline, &
- IntegrateColumnCubicSpline, InterpolateColumnLinear, &
- TestInterpolateColumn
+ IntegrateCubicSpline, IntegrateColumnCubicSpline, InterpolateLinear, &
+ TestInterpolate
! 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.
+! InterpolateCubicSpline: Compute cubic spline interpolation.
-! 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.
+! InterpolateLinear: Compute linear interpolation.
-! TestInterpolateColumn: Test spline interpolation subroutines.
+! TestInterpolate: Test spline interpolation subroutines.
contains
- subroutine CubicSplineCoefficients(x,y,n,yp1,ypn,y2)
+ subroutine CubicSplineCoefficients(x,y,n,y2ndDer)
-! 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.
+! Given arrays x(1:n) and y(1:n) containing a function,
+! i.e., y(i) = f(x(i)), with x monotonically increasing
+! this routine returns an array y2ndDer(1:n) that contains
+! the second derivatives of the interpolating function at x(1:n).
+! This routine uses boundary conditions 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) :: &
@@ -50,131 +41,56 @@
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
+ y2ndDer ! dy^2/dx^2 at each node
! local variables:
- integer, parameter :: &
- NMAX=500
- integer :: &
- i,k
+ integer :: i
real(kind=RKIND) :: &
- p,qn,sig,un,u(NMAX)
+ temp,xRatio,a(n)
- 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
+ y2ndDer(1)=0.0
+ y2ndDer(n)=0.0
+ a(1)=0.0
- 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
+ do i=2,n-1
+ xRatio=(x(i)-x(i-1))/(x(i+1)-x(i-1))
+ temp=1.0/(2.0+xRatio*y2ndDer(i-1))
+ y2ndDer(i)=temp*(xRatio-1.0)
+ a(i) = temp*(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)) &
+ -xRatio*a(i-1))
+ 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 i=n-1,1,-1
+ y2ndDer(i)=y2ndDer(i)*y2ndDer(i+1)+a(i)
+ enddo
- 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, &
+ x,y,y2ndDer,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
+! and given the array y2ndDer(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
+ y2ndDer ! 2nd derivative of y at nodes
real (kind=RKIND), dimension(nOut), intent(in) :: &
xOut ! node location, output grid
@@ -207,7 +123,7 @@
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)) &
+ + ((a**3-a)*y2ndDer(kIn) + (b**3-b)*y2ndDer(kIn+1)) &
*(h**2)/6.0
kOut = kOut + 1
@@ -218,13 +134,13 @@
enddo kInLoop
-end subroutine InterpolateColumnCubicSpline
+end subroutine InterpolateCubicSpline
-subroutine IntegrateCubicSpline(x,y,y2,n,x1,x2,y_integral)
+subroutine IntegrateCubicSpline(x,y,y2ndDer,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
+! and given the array y2ndDer(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.
@@ -238,7 +154,7 @@
real(kind=RKIND), intent(in), dimension(n) :: &
x, &! location of nodes
y, &! value at nodes
- y2 ! dy^2/dx^2 at each node
+ y2ndDer ! dy^2/dx^2 at each node
real(kind=RKIND), intent(in) :: &
x1,x2 ! limits of integration
@@ -270,24 +186,24 @@
! left side:
if (x1<x(j)) then
- F1 = -y(j)*h*0.5 + y2(j)*h**3/24.0
+ F1 = -y(j)*h*0.5 + y2ndDer(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 )
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(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
+ F2 = y(j+1)*h*0.5 - y2ndDer(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 )
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
endif
y_integral = y_integral + F2 - F1
@@ -298,11 +214,11 @@
subroutine IntegrateColumnCubicSpline( &
- x,y,y2,n, &
+ x,y,y2ndDer,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
+! and given the array y2ndDer(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
@@ -320,7 +236,7 @@
real(kind=RKIND), intent(in), dimension(n) :: &
x, &! location of nodes
y, &! value at nodes
- y2 ! dy^2/dx^2 at each node
+ y2ndDer ! dy^2/dx^2 at each node
real(kind=RKIND), dimension(nOut), intent(in) :: &
xOut ! output locations to compute integral
@@ -338,7 +254,7 @@
j = 1
h = x(j+1) - x(j)
h2 = h**2
- F1 = -y(j)*h*0.5 + y2(j)*h**3/24.0
+ F1 = -y(j)*h*0.5 + y2ndDer(j)*h**3/24.0
eps1 = 0.0 ! note: could use 1e-12*xOut(nOut)
k_loop: do k = 1,nOut
@@ -346,21 +262,21 @@
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
+ F2 = y(j+1)*h*0.5 - y2ndDer(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
+ F1 = -y(j)*h*0.5 + y2ndDer(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 )
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
y_integral(k) = y_integral(k) + F2 - F1
@@ -368,8 +284,8 @@
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 )
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
endif
enddo k_loop
@@ -377,7 +293,7 @@
end subroutine IntegrateColumnCubicSpline
- subroutine InterpolateColumnLinear( &
+ subroutine InterpolateLinear( &
x,y,n, &
xOut,yOut,nOut)
@@ -434,17 +350,17 @@
enddo kInLoop
- end subroutine InterpolateColumnLinear
+ end subroutine InterpolateLinear
-subroutine TestInterpolateColumn
+ subroutine TestInterpolate
! Test function to show how to operate the cubic spline subroutines
integer, parameter :: &
n = 10
real (kind=RKIND), dimension(n) :: &
- y, x, y2
+ y, x, y2ndDer
integer, parameter :: &
nOut = 100
@@ -462,8 +378,6 @@
do k=1,n
x(k) = k-4
- ! power function:
- !y(k) = x(k)**2
! trig function:
y(k) = sin(x(k)/2)
enddo
@@ -478,29 +392,36 @@
!
!-----------------------------------------------------------------------
- ! First, compute second derivative values at each node, y2.
- call CubicSplineCoefficients(x,y,n, &
- 2.0e30, 2.0e30, y2)
+ ! First, compute second derivative values at each node, y2ndDer.
+ call CubicSplineCoefficients(x,y,n,y2ndDer)
! Compute interpolated values yOut.
- call InterpolateColumnCubicSpline( &
- x,y,y2,n, &
- xOut,yOut,nOut)
+ call InterpolateCubicSpline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
-!-----------------------------------------------------------------------
-!
-! Output
-!
-!-----------------------------------------------------------------------
+ ! The following output can be copied directly into Matlab
+ print *, 'subplot(2,1,1)'
+ 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')"
-! 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')"
+ ! Compute interpolated values yOut.
+ call IntegrateColumnCubicSpline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
-end subroutine TestInterpolateColumn
+ ! The following output can be copied directly into Matlab
+ print *, 'subplot(2,1,2)'
+ print '(a,10f8.4,a)', 'x = [',x,'];'
+ print '(a,10f8.4,a)', 'y = 2*cos(-3/2) -2*cos(x/2);'
+ print '(a,100f8.4,a)', 'xOut = [',xOut,'];'
+ print '(a,100f8.4,a)', 'yOut = [',yOut,'];'
+ print *, "plot(x,y,'-*r',xOut,yOut,'x')"
+ end subroutine TestInterpolate
+
end module spline_interpolation
</font>
</pre>