<p><b>mpetersen@lanl.gov</b> 2011-02-09 09:07:59 -0700 (Wed, 09 Feb 2011)</p><p>BRANCH COMMIT: Improved comments and variable names in the spline module to be general and uniform.<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-09 00:24:15 UTC (rev 721)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-02-09 16:07:59 UTC (rev 722)
@@ -903,6 +903,13 @@
!TDR: still, this is a valid approximation at those interfaces and better than
!TDR: what will be grabbed from uninitialized memory.
+!mrp: I initialize tracerTop at the top and bottom as:
+!mrp: tracerTop(:,1,iCell)=0.0
+!mrp: tracerTop(:,maxLevelCell(iCell)+1,iCell)=0.0
+!mrp: after the following if statement, because there should be no contribution
+!mrp: from those layers. As you say, w should be aero anyway.
+
+
if (config_vert_tracer_adv.eq.'stencil'.and. &
config_vert_tracer_adv_order.eq.2) then
@@ -1029,9 +1036,8 @@
maxLevelCell(iCell), 2.0e30, 2.0e30, tracer2ndDer)
call InterpolateColumnCubicSpline( &
- tracersIn, tracersOut, tracer2ndDer, &
- posZMidZLevel, posZTopZLevel, &
- maxLevelCell(iCell), maxLevelCell(iCell)-1 )
+ posZMidZLevel, tracersIn, tracer2ndDer, maxLevelCell(iCell), &
+ posZTopZLevel, tracersOut, maxLevelCell(iCell)-1 )
tracerTop(iTracer,2:maxLevelCell(iCell),iCell) = tracersOut(1:maxLevelCell(iCell)-1)
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-09 00:24:15 UTC (rev 721)
+++ branches/ocean_projects/vert_adv_mrp/src/operators/module_spline_interpolation.F 2011-02-09 16:07:59 UTC (rev 722)
@@ -1,14 +1,3 @@
-!TDR: grep for "TDR" to find specific comments
-!global suggestions (see specific comments for reasoning)
-! zIn goes to xIn
-! zOut goes to xOut
-! better explanation of subroutines
-! remove "ocean specific" nomenclature (like "layer thickness")
-
-! TDR: I am happy to implement these changes, but it would be better
-! for you to do it, then I can review it.
-
-
module spline_interpolation
implicit none
@@ -20,35 +9,37 @@
IntegrateColumnCubicSpline, InterpolateColumnLinear, &
TestInterpolateColumn
-! TDR: is there a general reference for this entire module?
-! TDR: provide a brief explanation of what each public interface does
-! TDR: i.e. what is a typical use of each subroutine?
-! TDR: e.g. : given data (y) at locations (x), determine spline coefficients that can be passed into SSS.
+! Short Descriptions:
-! CubicSplineCoefficients:
+! CubicSplineCoefficients: Compute second derivatives at nodes.
+! This must be run before any of the other cubic spine functions.
-! InterpolateCubicSpline:
+! InterpolateCubicSpline: Compute cubic spline interpolation at a single point.
-! InterpolateColumnCubicSpline:
+! InterpolateColumnCubicSpline: Compute cubic spline interpolation at multiple points.
-! IntegrateCubicSpline:
+! IntegrateCubicSpline: Compute a single integral from spline data.
-! IntegrateColumnCubicSpline:
+! IntegrateColumnCubicSpline: Compute multiple integrals from spline data.
-! InterpolateColumnLinear:
+! InterpolateColumnLinear: Compute linear interpolation at multiple points.
-! TestInterpolateColumn:
+! TestInterpolateColumn: Test spline interpolation subroutines.
contains
subroutine CubicSplineCoefficients(x,y,n,yp1,ypn,y2)
-! TDR: The below description is pretty confusing
+! 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 n pairs (x_i, y_i), compute values of y2, which are
-! dy^2/dx^2 at each node used to create a cubic spline on each
-! interval. yp1 and ypn are the first derivatives at each end,
-! and a natural spline boundary condition is used when these are >1e30.
! This routine adapted from Numerical Recipes in Fortran, 2nd Ed, p. 109
! and was originally named spline.
@@ -57,9 +48,7 @@
integer, intent(in) :: &
n ! number of nodes
real(kind=RKIND), intent(in), dimension(n) :: &
-!TDR: since this will be in frameworks, there should be no
-!TDR "core specific" annotation. I suggest "location of nodes" throughout
- x, &! depth of nodes
+ x, &! location of nodes
y ! value at nodes
real(kind=RKIND), intent(in) :: &
yp1, &! first derivative of y at first endpoint
@@ -112,14 +101,16 @@
end subroutine CubicSplineCoefficients
- subroutine InterpolateCubicSpline(xa,ya,y2a,n,x,y)
+ subroutine InterpolateCubicSpline( &
+ x,y,y2,n, &
+ xOut,yOut)
-! TDR: this description is pretty confusing.
+! 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.
-! Given arrays xa and ya of length n, which tabulate a function (with
-! xa in order), and given the array y2a, which is the output from spline,
-! and given a value of x this routine returns a cubic-spline interpolated
-! value y.
! This routine adapted from Numerical Recipes in Fortran, 2nd Ed, p 110
! and was originally named splint, for SPLine INTerpolation
@@ -128,14 +119,14 @@
integer, intent(in) :: &
n
real(kind=RKIND), intent(in), dimension(n) :: &
- xa,ya,y2a
+ x,y,y2
real(kind=RKIND), intent(in) :: &
- x
+ xOut
! OUTPUT PARAMETERS:
real(kind=RKIND), intent(out) :: &
- y
+ yOut
! local variables:
@@ -148,60 +139,54 @@
khi=n
do while (khi-klo.gt.1)
k=(khi+klo)/2
- if(xa(k).gt.x)then
+ if(x(k).gt.xOut)then
khi=k
else
klo=k
endif
enddo
- h=xa(khi)-xa(klo)
- ! if (h.eq.0.) pause 'bad xa input in splint'
- a=(xa(khi)-x)/h
- b=(x-xa(klo))/h
- y=a*ya(klo) + b*ya(khi) &
- + ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.0
+ 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( &
- phiIn, phiOut, phi2DerIn, &
- zIn, zOut, &
- kmaxIn, kmaxOut)
+ x,y,y2,n, &
+ xOut,yOut,nOut)
-! TDR: this description is better, but we need to refrain from
-! TDR using core specific language. why switch to zIn/zOut. I think
-! TDR: that we should use xIn/xOut to be consistent with the other
-! TDR: subroutines. In that, zIn is replaced by xIn and is described
-! TDR: as "x locations of nodes, assumed to be monotonically increasing"
+! 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.
-! Given kmaxIn pairs (zIn,phiIn), compute cubic spline interpolation
-! phiOut at each zOut.
-! This subroutine assumes that both zIn and zOut are monotonically
-! increasing. This should occur naturally if the z coordinates are
-! computed from layer thicknesses in the calling routine.
-! It also assumes that all values of zOut are within the first
-! and last values of zIn.
-! This routine adapted from Numerical Recipes in Fortran, 2nd Ed, p 110
+! This routine adapted from Numerical Recipes in Fortran, 2nd Ed, p 110
! INPUT PARAMETERS:
- real (kind=RKIND), dimension(kmaxIn), intent(in) :: &
- zIn, &! layer thickness, input grid
- phiIn, &! interpolation variable, input grid
- phi2DerIn ! 2nd derivative of phi at nodes
+ 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(kmaxOut), intent(in) :: &
- zOut ! layer thickness, output grid
+ real (kind=RKIND), dimension(nOut), intent(in) :: &
+ xOut ! node location, output grid
integer, intent(in) :: &
- KmaxIn, &! number of layers, input grid
- KmaxOut ! number of layers, output grid
+ n, &! number of nodes, input grid
+ nOut ! number of nodes, output grid
! OUTPUT PARAMETERS:
- real (kind=RKIND), dimension(kmaxOut), intent(out) :: &
- phiOut ! interpolation variable, output grid
+ real (kind=RKIND), dimension(nOut), intent(out) :: &
+ yOut ! interpolation variable, output grid
! local variables:
@@ -213,21 +198,21 @@
kOut = 1
- kInLoop: do kIn = 1,kmaxIn-1
+ kInLoop: do kIn = 1,n-1
- h = zIn(kIn+1)-zIn(kIn)
+ h = x(kIn+1)-x(kIn)
- do while(zOut(kOut) < zIn(kIn+1))
+ do while(xOut(kOut) < x(kIn+1))
- a = (zIn(kIn+1)-zOut(kOut))/h
- b = (zOut(kOut)-zIn (kIn) )/h
- phiOut(kOut) = a*phiIn(kIn) + b*phiIn(kIn+1) &
- + ((a**3-a)*phi2DerIn(kIn) + (b**3-b)*phi2DerIn(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>kmaxOut) exit kInLoop
+ if (kOut>nOut) exit kInLoop
enddo
@@ -238,17 +223,20 @@
subroutine IntegrateCubicSpline(x,y,y2,n,x1,x2,y_integral)
-! Given n pairs (x_i, y_i), and dy^2/dx^2 at each node, compute
-! the integral of y from x1 to x2.
-! The integration formula was created by analytically integrating a
-! cubic spline between each node.
+! 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, &! depth of nodes
+ x, &! location of nodes
y, &! value at nodes
y2 ! dy^2/dx^2 at each node
real(kind=RKIND), intent(in) :: &
@@ -311,30 +299,35 @@
subroutine IntegrateColumnCubicSpline( &
x,y,y2,n, &
- xout,y_integral, nout)
+ xOut,y_integral, nOut)
-! Given n pairs (x_i, y_i), and dy^2/dx^2 at each node, compute
-! the integral of y. This is a cumulative integration, so that
-! y_integral(j) holds the integral of y from 0 to xout(j).
+! 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
+ nOut ! number of output locations to compute integral
real(kind=RKIND), intent(in), dimension(n) :: &
- x, &! depth of nodes
+ 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
+ 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
+ real(kind=RKIND), dimension(nOut), intent(out) :: &
+ y_integral ! integral from 0 to xOut
! local variables:
@@ -346,13 +339,13 @@
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)
+ eps1 = 0.0 ! note: could use 1e-12*xOut(nOut)
- k_loop: do k = 1,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)
+ 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
@@ -360,20 +353,20 @@
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
+ 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
+ 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
+ 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 )
@@ -385,38 +378,33 @@
subroutine InterpolateColumnLinear( &
- phiIn,phiOut, &
- zIn,zOut, &
- kmaxIn,kmaxOut)
+ x,y,n, &
+ xOut,yOut,nOut)
-! TDR: Desciptions need to be "core independent"
+! 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.
-! !DESCRIPTION:
-! Given kmaxIn pairs (zIn,phiIn), compute the linear interpolation
-! phiOut at each zOut.
-! This subroutine assumes that both zIn and zOut are monotonically
-! increasing. This should occur naturally if the z coordinates are
-! computed from layer thicknesses in the calling routine.
-! It also assumes that all values of zOut are within the first
-! and last values of zIn.
-
! !INPUT PARAMETERS:
- real (kind=RKIND), dimension(kmaxIn), intent(in) :: &
- zIn, &! layer thickness, input grid
- phiIn ! interpolation variable, input grid
+ real (kind=RKIND), dimension(n), intent(in) :: &
+ x, &! node location, input grid
+ y ! interpolation variable, input grid
- real (kind=RKIND), dimension(kmaxOut), intent(in) :: &
- zOut ! layer thickness, output grid
+ real (kind=RKIND), dimension(nOut), intent(in) :: &
+ xOut ! node location, output grid
integer, intent(in) :: &
- KmaxIn, &! number of layers, input grid
- KmaxOut ! number of layers, output grid
+ N, &! number of nodes, input grid
+ NOut ! number of nodes, output grid
! !OUTPUT PARAMETERS:
- real (kind=RKIND), dimension(kmaxOut), intent(out) :: &
- phiOut ! interpolation variable, output grid
+ real (kind=RKIND), dimension(nOut), intent(out) :: &
+ yOut ! interpolation variable, output grid
!-----------------------------------------------------------------------
!
@@ -429,18 +417,18 @@
kOut = 1
- kInLoop: do kIn = 1,kmaxIn-1
+ kInLoop: do kIn = 1,n-1
- do while(zOut(kOut) < zIn(kIn+1))
+ do while(xOut(kOut) < x(kIn+1))
- phiOut(kOut) = phiIn(kIn) &
- + (phiIn(kIn+1)-phiIn(kIn)) &
- /(zIn(kIn+1) -zIn(kIn) ) &
- *(zOut(kOut) -zIn(kIn) )
+ yOut(kOut) = y(kIn) &
+ + (y(kIn+1)-y(kIn)) &
+ /(x(kIn+1) -x(kIn) ) &
+ *(xOut(kOut) -x(kIn) )
kOut = kOut + 1
- if (kOut>kmaxOut) exit kInLoop
+ if (kOut>nOut) exit kInLoop
enddo
@@ -454,34 +442,34 @@
! Test function to show how to operate the cubic spline subroutines
integer, parameter :: &
- kmaxIn = 10
- real (kind=RKIND), dimension(kmaxIn) :: &
- phiIn, zIn, phi2DerIn
+ n = 10
+ real (kind=RKIND), dimension(n) :: &
+ y, x, y2
integer, parameter :: &
- kmaxOut = 100
- real (kind=RKIND), dimension(kmaxOut) :: &
- phiOut, zOut
+ nOut = 100
+ real (kind=RKIND), dimension(nOut) :: &
+ yOut, xOut
integer :: &
k
!-----------------------------------------------------------------------
!
-! Create zIn, PhiIn, zOut
+! Create x, y, xOut
!
!-----------------------------------------------------------------------
- do k=1,kmaxIn
- zIn(k) = k-4
+ do k=1,n
+ x(k) = k-4
! power function:
- !phiIn(k) = zIn(k)**2
+ !y(k) = x(k)**2
! trig function:
- phiIn(k) = sin(zIn(k)/2)
+ y(k) = sin(x(k)/2)
enddo
- do k=1,kmaxOut
- zOut(k) = zIn(1) + k/(kmaxOut+1.0)*(zIn(kmaxIn)-zIn(1))
+ do k=1,nOut
+ xOut(k) = x(1) + k/(nOut+1.0)*(x(n)-x(1))
enddo
!-----------------------------------------------------------------------
@@ -490,13 +478,14 @@
!
!-----------------------------------------------------------------------
- ! First, compute second derivative values at each node, phi2DerIn.
- call CubicSplineCoefficients(zIn,phiIn,kmaxIn, &
- 2.0e30, 2.0e30, phi2DerIn)
+ ! First, compute second derivative values at each node, y2.
+ call CubicSplineCoefficients(x,y,n, &
+ 2.0e30, 2.0e30, y2)
- ! Compute interpolated values phiOut.
- call InterpolateColumnCubicSpline(phiIn, phiOut, &
- phi2DerIn, zIn, zOut, kmaxIn, kmaxOut)
+ ! Compute interpolated values yOut.
+ call InterpolateColumnCubicSpline( &
+ x,y,y2,n, &
+ xOut,yOut,nOut)
!-----------------------------------------------------------------------
!
@@ -505,11 +494,11 @@
!-----------------------------------------------------------------------
! The following output can be copied directly into Matlab
- print '(a,10f8.4,a)', 'zIn = [',zIn,'];'
- print '(a,10f8.4,a)', 'phiIn = [',phiIn,'];'
- print '(a,100f8.4,a)', 'zOut = [',zOut,'];'
- print '(a,100f8.4,a)', 'phiOut = [',phiOut,'];'
- print *, "plot(zIn,phiIn,'-*r',zOut,phiOut,'x')"
+ 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
</font>
</pre>