<p><b>mpetersen@lanl.gov</b> 2011-02-03 15:37:43 -0700 (Thu, 03 Feb 2011)</p><p>BRANCH COMMIT: Added Cubic Spline interpolation to vertical advection.<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-01 21:00:50 UTC (rev 711)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_mpas_core.F 2011-02-03 22:37:43 UTC (rev 712)
@@ -112,8 +112,6 @@
enddo
end do
-
-
! mrp 110131 testing only end
@@ -404,7 +402,7 @@
nVertLevels = block % mesh % nVertLevels
vertexDegree = block % mesh % vertexDegree
-! mrp 110131 for testing only
+! mrp 110131 for testing only, for basin domain where maxLevelCell is not initialized
maxLevelCell(1:nCells) = nVertLevels
! mrp 110131 for testing only end
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-01 21:00:50 UTC (rev 711)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-02-03 22:37:43 UTC (rev 712)
@@ -714,12 +714,13 @@
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
- real (kind=RKIND), dimension(:), pointer :: zTopZLevel, &
+ real (kind=RKIND), dimension(:), pointer :: zTopZLevel,zMidZLevel, &
hRatioZLevelK, hRatioZLevelKm1
+ real (kind=RKIND), dimension(:), allocatable:: vertDiffTop, tracer2ndDer, tracersIn, tracersOut, posZMidZLevel, &
+ posZTopZLevel
real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, tracerTop, boundaryMask
real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer
- real (kind=RKIND), dimension(:), allocatable:: vertDiffTop
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
@@ -739,6 +740,7 @@
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
zTopZLevel => grid % zTopZLevel % array
+ zMidZLevel => grid % zMidZLevel % array
hRatioZLevelK => grid % hRatioZLevelK % array
hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
boundaryEdge => grid % boundaryEdge % array
@@ -1034,6 +1036,48 @@
end do
end do
+ elseif (config_vert_tracer_adv.eq.'CubicSpline') then
+
+ allocate(tracer2ndDer(nVertLevels))
+ allocate(tracersIn(nVertLevels),tracersOut(nVertLevels), &
+ posZMidZLevel(nVertLevels), posZTopZLevel(nVertLevels))
+
+ do iCell=1,grid % 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. For the ocean, zlevel coordinates are negative and
+ ! decreasing, but spline functions assume increasing, so flip to positive
+ tracersIn(1:maxLevelCell(iCell))=tracers(iTracer,1:maxLevelCell(iCell),iCell)
+ posZMidZLevel(1:maxLevelCell(iCell)) = -zMidZLevel(1:maxLevelCell(iCell))
+ posZTopZLevel(1:maxLevelCell(iCell)-1) = -zTopZLevel(2:maxLevelCell(iCell))
+
+ call CubicSplineCoefficients(posZMidZLevel, &
+ tracersIn, &
+ maxLevelCell(iCell), 2.0e30, 2.0e30, tracer2ndDer)
+
+ call InterpolateColumnCubicSpline( &
+ tracersIn, tracersOut, tracer2ndDer, &
+ posZMidZLevel, posZTopZLevel, &
+ maxLevelCell(iCell), maxLevelCell(iCell)-1 )
+
+ tracerTop(iTracer,2:maxLevelCell(iCell)) = tracersOut(1:maxLevelCell(iCell)-1)
+
+ end do
+ tracerTop(:,maxLevelCell(iCell)+1)=0.0
+ do k=1,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
+ end do
+ end do
+ end do
+
+ deallocate(tracer2ndDer)
+ deallocate(tracersIn,tracersOut, posZMidZLevel, posZTopZLevel)
endif
deallocate(tracerTop)
@@ -2007,4 +2051,349 @@
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
</font>
</pre>