<p><b>mpetersen@lanl.gov</b> 2011-02-07 13:17:24 -0700 (Mon, 07 Feb 2011)</p><p>BRANCH COMMIT: Add file src/operators/module_spline_interpolation.F for spline routines.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_mpas_core.F 2011-02-03 22:37:43 UTC (rev 712)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_mpas_core.F 2011-02-07 20:17:24 UTC (rev 713)
@@ -81,40 +81,10 @@
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
-! mrp 110131 temp commented out:
-! call rbfInterp_initialize(mesh)
-! call init_reconstruct(mesh)
-! call reconstruct(block % state % time_levs(1) % state, mesh)
-! mrp 110131 temp commented out end
+ call rbfInterp_initialize(mesh)
+ call init_reconstruct(mesh)
+ call reconstruct(block % state % time_levs(1) % state, mesh)
- ! mrp 110131 testing only
- block % mesh % u_src % array(:,:) = &
- block % mesh % u_src % array(:,:) *100
- print *, 'min,maxval(u_src) ',minval( block % mesh % u_src % array(:,:)),maxval( block % mesh % u_src % array(:,:))
-
- do iCell=1,block % mesh % nCells
- do k=1,mesh % nVertLevels
- block % state % time_levs(1) % state % tracers % array( &
- 1, k ,iCell) = 31-k
- block % state % time_levs(1) % state % tracers % array( &
- 2, k ,iCell) = 32 + k*0.32
- enddo
- block % state % time_levs(1) % state % tracers % array( &
- 3, :,iCell) = 1.
- block % state % time_levs(1) % state % tracers % array( &
- 4, :,iCell) = 1.
-
- block % state % time_levs(1) % state % tracers % array( &
- 3, 1:10,iCell) = 0.
- do k=1,mesh % nVertLevels,2
- block % state % time_levs(1) % state % tracers % array( &
- 4, k ,iCell) = 0.
- enddo
-
- end do
- ! mrp 110131 testing only end
-
-
! initialize velocities and tracers on land to be -1e34
! The reconstructed velocity on land will have values not exactly
! -1e34 due to the interpolation of reconstruction.
@@ -402,10 +372,6 @@
nVertLevels = block % mesh % nVertLevels
vertexDegree = block % mesh % vertexDegree
-! mrp 110131 for testing only, for basin domain where maxLevelCell is not initialized
- maxLevelCell(1:nCells) = nVertLevels
-! mrp 110131 for testing only end
-
! for z-grids, maxLevelCell should be in input state
! Isopycnal grid uses all vertical cells
if (config_vert_grid_type.eq.'isopycnal') then
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-03 22:37:43 UTC (rev 712)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-02-07 20:17:24 UTC (rev 713)
@@ -5,6 +5,7 @@
use constants
use dmpar
use vector_reconstruction
+ use spline_interpolation
contains
@@ -2051,349 +2052,4 @@
end subroutine equation_of_state_jm
-
- subroutine CubicSplineCoefficients(x,y,n,yp1,ypn,y2)
-
-! 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.
-
- implicit none
-
-! INPUT PARAMETERS:
-
- integer, intent(in) :: &
- n ! number of nodes
- real(kind=RKIND), intent(in), dimension(n) :: &
- x, &! depth 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(xa,ya,y2a,n,x,y)
-
-! 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
-! Mark Petersen, LANL, March 2009
-!
-! !REVISION HISTORY:
-! same as module
-
- implicit none
-
-! INPUT PARAMETERS:
-
- integer, intent(in) :: &
- n
- real(kind=RKIND), intent(in), dimension(n) :: &
- xa,ya,y2a
- real(kind=RKIND), intent(in) :: &
- x
-
-! OUTPUT PARAMETERS:
-
- real(kind=RKIND), intent(out) :: &
- y
-
-! 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(xa(k).gt.x)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
-
- end subroutine InterpolateCubicSpline
-
-
- subroutine InterpolateColumnCubicSpline( &
- phiIn, phiOut, phi2DerIn, &
- zIn, zOut, &
- kmaxIn, kmaxOut)
-
-! 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
-
- implicit none
-
-! 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(kmaxOut), intent(in) :: &
- zOut ! layer thickness, output grid
-
- integer, intent(in) :: &
- KmaxIn, &! number of layers, input grid
- KmaxOut ! number of layers, output grid
-
-! OUTPUT PARAMETERS:
-
- real (kind=RKIND), dimension(kmaxOut), intent(out) :: &
- phiOut ! interpolation variable, output grid
-
-! local variables:
-
- integer :: &
- kIn, kOut ! counters
-
- real (kind=RKIND) :: &
- a, b, h
-
- kOut = 1
-
- kInLoop: do kIn = 1,kmaxIn-1
-
- h = zIn(kIn+1)-zIn(kIn)
-
- do while(zOut(kOut) < zIn(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)) &
- *(h**2)/6.0
-
- kOut = kOut + 1
-
- if (kOut>kmaxOut) exit kInLoop
-
- enddo
-
- enddo kInLoop
-
- end subroutine InterpolateColumnCubicSpline
-
-
- 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.
-
- implicit none
-
-! INPUT PARAMETERS:
-
- integer, intent(in) :: &
- n ! number of nodes
- real(kind=RKIND), intent(in), dimension(n) :: &
- x, &! depth 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 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).
-! The integration formula was created by analytically integrating a
-! cubic spline between each node.
-
- implicit none
-
-! 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, &! depth 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
-
end module time_integration
Modified: branches/ocean_projects/vert_adv_mrp/src/operators/Makefile
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/operators/Makefile 2011-02-03 22:37:43 UTC (rev 712)
+++ branches/ocean_projects/vert_adv_mrp/src/operators/Makefile 2011-02-07 20:17:24 UTC (rev 713)
@@ -1,6 +1,6 @@
.SUFFIXES: .F .o
-OBJS = module_RBF_interpolation.o module_vector_reconstruction.o
+OBJS = module_RBF_interpolation.o module_vector_reconstruction.o module_spline_interpolation.o
all: operators
@@ -9,6 +9,7 @@
module_vector_reconstruction.o: module_RBF_interpolation.o
module_RBF_interpolation.o:
+module_spline_interpolation:
clean:
$(RM) *.o *.mod *.f90 libops.a
Added: branches/ocean_projects/vert_adv_mrp/src/operators/module_spline_interpolation.F
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/operators/module_spline_interpolation.F (rev 0)
+++ branches/ocean_projects/vert_adv_mrp/src/operators/module_spline_interpolation.F 2011-02-07 20:17:24 UTC (rev 713)
@@ -0,0 +1,344 @@
+module spline_interpolation
+
+ implicit none
+
+ private
+
+ public :: CubicSplineCoefficients, InterpolateCubicSpline, &
+ InterpolateColumnCubicSpline, IntegrateCubicSpline, &
+ IntegrateColumnCubicSpline
+
+ contains
+
+ subroutine CubicSplineCoefficients(x,y,n,yp1,ypn,y2)
+
+! 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.
+
+! INPUT PARAMETERS:
+
+ integer, intent(in) :: &
+ n ! number of nodes
+ real(kind=RKIND), intent(in), dimension(n) :: &
+ x, &! depth 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(xa,ya,y2a,n,x,y)
+
+! 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
+
+! INPUT PARAMETERS:
+
+ integer, intent(in) :: &
+ n
+ real(kind=RKIND), intent(in), dimension(n) :: &
+ xa,ya,y2a
+ real(kind=RKIND), intent(in) :: &
+ x
+
+! OUTPUT PARAMETERS:
+
+ real(kind=RKIND), intent(out) :: &
+ y
+
+! 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(xa(k).gt.x)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
+
+ end subroutine InterpolateCubicSpline
+
+
+ subroutine InterpolateColumnCubicSpline( &
+ phiIn, phiOut, phi2DerIn, &
+ zIn, zOut, &
+ kmaxIn, kmaxOut)
+
+! 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
+
+! 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(kmaxOut), intent(in) :: &
+ zOut ! layer thickness, output grid
+
+ integer, intent(in) :: &
+ KmaxIn, &! number of layers, input grid
+ KmaxOut ! number of layers, output grid
+
+! OUTPUT PARAMETERS:
+
+ real (kind=RKIND), dimension(kmaxOut), intent(out) :: &
+ phiOut ! interpolation variable, output grid
+
+! local variables:
+
+ integer :: &
+ kIn, kOut ! counters
+
+ real (kind=RKIND) :: &
+ a, b, h
+
+ kOut = 1
+
+ kInLoop: do kIn = 1,kmaxIn-1
+
+ h = zIn(kIn+1)-zIn(kIn)
+
+ do while(zOut(kOut) < zIn(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)) &
+ *(h**2)/6.0
+
+ kOut = kOut + 1
+
+ if (kOut>kmaxOut) exit kInLoop
+
+ enddo
+
+ enddo kInLoop
+
+end subroutine InterpolateColumnCubicSpline
+
+
+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.
+
+! INPUT PARAMETERS:
+
+ integer, intent(in) :: &
+ n ! number of nodes
+ real(kind=RKIND), intent(in), dimension(n) :: &
+ x, &! depth 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 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).
+! The integration formula was created by analytically integrating a
+! cubic spline between each node.
+
+! 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, &! depth 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
+
+end module spline_interpolation
+
</font>
</pre>