<p><b>mpetersen@lanl.gov</b> 2011-02-24 17:05:27 -0700 (Thu, 24 Feb 2011)</p><p>OCEAN CORE TRUNK COMMIT: <br>
Merge high order vertical advection branch to trunk. <br>
<br>
This revision allows one to choose the method and order of the<br>
vertical tracer advection using namelist flags<br>
config_vert_tracer_adv <br>
config_vert_tracer_adv_order <br>
This changes the method of approximating the tracer value at the top<br>
of each cell, where the tracer value is multiplied by wTop to obtain<br>
the upward tracer flux.<br>
<br>
Revisions by Mark, reviewed by Todd. <br>
Tested on global x3 and double gyre quad and hex domains. <br>
</p><hr noshade><pre><font color="gray">
Property changes on: trunk/mpas
___________________________________________________________________
Added: svn:mergeinfo
+ /branches/ocean_projects/vert_adv_mrp:704-745
Modified: trunk/mpas/namelist.input.ocean
===================================================================
--- trunk/mpas/namelist.input.ocean 2011-02-23 21:50:59 UTC (rev 745)
+++ trunk/mpas/namelist.input.ocean 2011-02-25 00:05:27 UTC (rev 746)
@@ -45,7 +45,8 @@
config_eos_type = 'linear'
/
&advection
- config_vert_tracer_adv = 'upwind'
+ config_vert_tracer_adv = 'spline'
+ config_vert_tracer_adv_order = 3
config_tracer_adv_order = 2
config_thickness_adv_order = 2
config_positive_definite = .false.
Modified: trunk/mpas/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2011-02-23 21:50:59 UTC (rev 745)
+++ trunk/mpas/src/core_ocean/Registry 2011-02-25 00:05:27 UTC (rev 746)
@@ -32,7 +32,8 @@
namelist real vmix config_vmixTanhZMid -100
namelist real vmix config_vmixTanhZWidth 100
namelist character eos config_eos_type linear
-namelist character advection config_vert_tracer_adv 'centered'
+namelist character advection config_vert_tracer_adv stencil
+namelist integer advection config_vert_tracer_adv_order 4
namelist integer advection config_tracer_adv_order 2
namelist integer advection config_thickness_adv_order 2
namelist logical advection config_positive_definite false
@@ -130,6 +131,9 @@
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
var persistent real zMidZLevel ( nVertLevels ) 0 - zMidZLevel mesh - -
var persistent real zTopZLevel ( nVertLevelsP1 ) 0 - zTopZLevel mesh - -
+var persistent real hMeanTopZLevel ( nVertLevels ) 0 - hMeanTopZLevel mesh - -
+var persistent real hRatioZLevelK ( nVertLevels ) 0 - hRatioZLevelK mesh - -
+var persistent real hRatioZLevelKm1 ( nVertLevels ) 0 - hRatioZLevelKm1 mesh - -
# Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
Modified: trunk/mpas/src/core_ocean/module_mpas_core.F
===================================================================
--- trunk/mpas/src/core_ocean/module_mpas_core.F 2011-02-23 21:50:59 UTC (rev 745)
+++ trunk/mpas/src/core_ocean/module_mpas_core.F 2011-02-25 00:05:27 UTC (rev 746)
@@ -76,14 +76,36 @@
type (block_type), intent(inout) :: block
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
- integer :: i
+ integer :: i, iEdge, iCell, k
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
-
+
call rbfInterp_initialize(mesh)
call init_reconstruct(mesh)
call reconstruct(block % state % time_levs(1) % state, mesh)
+
+ ! 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.
+
+ do iEdge=1,block % mesh % nEdges
+ ! mrp 101115 note: in order to include flux boundary conditions, the following
+ ! line will need to change. Right now, set boundary edges between land and
+ ! water to have zero velocity.
+ block % state % time_levs(1) % state % u % array( &
+ block % mesh % maxLevelEdgeTop % array(iEdge)+1 &
+ :block % mesh % maxLevelEdgeBot % array(iEdge), iEdge) = 0.0
+
+ block % state % time_levs(1) % state % u % array( &
+ block % mesh % maxLevelEdgeBot % array(iEdge)+1: &
+ block % mesh % nVertLevels,iEdge) = -1e34
+ end do
+ do iCell=1,block % mesh % nCells
+ block % state % time_levs(1) % state % tracers % array( &
+ :, block % mesh % maxLevelCell % array(iCell)+1 &
+ :block % mesh % nVertLevels,iCell) = -1e34
+ end do
if (.not. config_do_restart) then
block % state % time_levs(1) % state % xtime % scalar = 0.0
@@ -224,7 +246,7 @@
block_ptr => domain % blocklist
if (associated(block_ptr % next)) then
write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
- 'that there is only one block per processor.'
+ 'that there is only one block per processor.'
end if
call timer_start("global diagnostics")
@@ -253,7 +275,8 @@
integer :: iTracer, cell
real (kind=RKIND), dimension(:), pointer :: &
- hZLevel, zMidZLevel, zTopZLevel
+ hZLevel, zMidZLevel, zTopZLevel, &
+ hMeanTopZLevel, hRatioZLevelK, hRatioZLevelKm1
real (kind=RKIND), dimension(:,:), pointer :: h
integer :: nVertLevels
@@ -266,6 +289,9 @@
zMidZLevel => block % mesh % zMidZLevel % array
zTopZLevel => block % mesh % zTopZLevel % array
nVertLevels = block % mesh % nVertLevels
+ hMeanTopZLevel => block % mesh % hMeanTopZLevel % array
+ hRatioZLevelK => block % mesh % hRatioZLevelK % array
+ hRatioZLevelKm1 => block % mesh % hRatioZLevelKm1 % array
! These should eventually be in an input file. For now
! I just read them in from h(:,1).
@@ -281,6 +307,15 @@
zTopZLevel(k+1) = zTopZLevel(k)- hZLevel(k)
end do
+ hMeanTopZLevel(1) = 0.0
+ hRatioZLevelK(1) = 0.0
+ hRatioZLevelKm1(1) = 0.0
+ do k = 2,nVertLevels
+ hMeanTopZLevel(k) = 0.5*(hZLevel(k-1) + hZLevel(k))
+ hRatioZLevelK(k) = 0.5*hZLevel(k)/hMeanTopZLevel(k)
+ hRatioZLevelKm1(k) = 0.5*hZLevel(k-1)/hMeanTopZLevel(k)
+ end do
+
block => block % next
end do
Modified: trunk/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_ocean/module_time_integration.F 2011-02-23 21:50:59 UTC (rev 745)
+++ trunk/mpas/src/core_ocean/module_time_integration.F 2011-02-25 00:05:27 UTC (rev 746)
@@ -5,6 +5,7 @@
use constants
use dmpar
use vector_reconstruction
+ use spline_interpolation
contains
@@ -71,7 +72,7 @@
type (block_type), pointer :: block
type (state_type) :: provis
- integer :: rk_step
+ integer :: rk_step, iEdge
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
@@ -265,7 +266,7 @@
type (mesh_type), intent(in) :: grid
integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
- vertex1, vertex2, eoe, i, j, kmax
+ vertex1, vertex2, eoe, i, j
integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
@@ -357,27 +358,39 @@
! for z-level, only compute height tendency for top layer.
if (config_vert_grid_type.eq.'isopycnal') then
- kmax = nVertLevels
- elseif (config_vert_grid_type.eq.'zlevel') then
- kmax = 1
- endif
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,kmax
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend_h(k,cell1) = tend_h(k,cell1) - flux
- tend_h(k,cell2) = tend_h(k,cell2) + flux
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell1) = tend_h(k,cell1) - flux
+ tend_h(k,cell2) = tend_h(k,cell2) + flux
+ end do
end do
- end do
+ do iCell=1,nCells
+ do k=1,nVertLevels
+ tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
+ end do
+ end do
- do iCell=1,nCells
- do k=1,kmax
- tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
+ elseif (config_vert_grid_type.eq.'zlevel') then
+
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,min(1,maxLevelEdgeTop(iEdge))
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell1) = tend_h(k,cell1) - flux
+ tend_h(k,cell2) = tend_h(k,cell2) + flux
+ end do
end do
- end do
+ do iCell=1,nCells
+ tend_h(1,iCell) = tend_h(1,iCell) / areaCell(iCell)
+ end do
+ endif ! config_vert_grid_type
+
!
! height tendency: vertical advection term -d/dz(hw)
!
@@ -599,10 +612,6 @@
do iEdge=1,grid % nEdgesSolve
- ! forcing in top layer only
- tend_u(1,iEdge) = tend_u(1,iEdge) &
- + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
-
k = maxLevelEdgeTop(iEdge)
! efficiency note: it would be nice to avoid this
@@ -612,18 +621,17 @@
if (k>0) then
- ! bottom drag is the same as POP:
- ! -c |u| u where c is unitless and 1.0e-3.
- ! see POP Reference guide, section 3.4.4.
+ ! forcing in top layer only
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- - 1.0e-3*u(k,iEdge) &
- *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
+ ! bottom drag is the same as POP:
+ ! -c |u| u where c is unitless and 1.0e-3.
+ ! see POP Reference guide, section 3.4.4.
- ! old bottom drag, just linear friction
- ! du/dt = u/tau where tau=100 days.
- !tend_u(k,iEdge) = tend_u(k,iEdge) &
- ! - u(k,iEdge)/(100.0*86400.0)
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ - 1.0e-3*u(k,iEdge) &
+ *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
endif
@@ -692,7 +700,7 @@
type (mesh_type), intent(in) :: grid
integer :: i, k, iCell, iEdge, iTracer, cell1, cell2, upwindCell,&
- nEdges, nCells, nVertLevels, num_tracers
+ nEdges, nCells, nCellsSolve, nVertLevels, num_tracers
real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
real (kind=RKIND) :: flux, tracer_edge, r
real (kind=RKIND), dimension(:), pointer :: &
@@ -707,17 +715,18 @@
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
- real (kind=RKIND), dimension(:), pointer :: zTopZLevel
- real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, tracerTop, boundaryMask
- real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer
+ 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, boundaryMask
+ real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer, tracerTop
- real (kind=RKIND), dimension(:), allocatable:: vertDiffTop
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND) :: coef_3rd_order
+ real (kind=RKIND) :: coef_3rd_order, flux3Coef, cSignWTop
-
u => s % u % array
h => s % h % array
boundaryCell=> grid % boundaryCell % array
@@ -732,6 +741,9 @@
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
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
@@ -739,6 +751,7 @@
nEdges = grid % nEdges
nCells = grid % nCells
+ nCellsSolve = grid % nCellsSolve
nVertLevels = grid % nVertLevels
num_tracers = s % num_tracers
@@ -880,44 +893,183 @@
!
! tracer tendency: vertical advection term -d/dz( h \phi w)
!
- allocate(tracerTop(num_tracers,nVertLevels+1))
- tracerTop(:,1)=0.0
- do iCell=1,grid % nCellsSolve
- if (config_vert_tracer_adv.eq.'centered') then
- do k=2,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tracerTop(iTracer,k) = ( tracers(iTracer,k-1,iCell) &
- +tracers(iTracer,k ,iCell))/2.0
- end do
- end do
-
- elseif (config_vert_tracer_adv.eq.'upwind') then
- do k=2,maxLevelCell(iCell)
- if (wTop(k,iCell)>0.0) then
- upwindCell = k
- else
- upwindCell = k-1
- endif
- do iTracer=1,num_tracers
- tracerTop(iTracer,k) = tracers(iTracer,upwindCell,iCell)
- end do
- end do
+ if (config_vert_grid_type.eq.'zlevel') then
- endif
- tracerTop(:,maxLevelCell(iCell)+1)=0.0
+ allocate(tracerTop(num_tracers,nVertLevels+1,nCells))
- do k=1,maxLevelCell(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
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
- - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
+ tracerTop(iTracer,1,iCell) = tracers(iTracer,1,iCell)
+ tracerTop(iTracer,maxLevelCell(iCell)+1,iCell) = &
+ tracers(iTracer,maxLevelCell(iCell),iCell)
end do
end do
- enddo ! iCell
- deallocate(tracerTop)
+ if (config_vert_tracer_adv.eq.'stencil'.and. &
+ config_vert_tracer_adv_order.eq.2) then
+ ! 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
+ end do
+ end do
+ end do
+
+ 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.
+
+ ! 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.
+ 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
+
+ ! 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)
+ 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
+ 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
+
+ ! 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)
+ end do
+ end do
+ end do
+
+ elseif (config_vert_tracer_adv.eq.'spline'.and. &
+ config_vert_tracer_adv_order.eq.3) then
+
+ ! Compute tracerTop using cubic spline interpolation.
+
+ 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.
+
+ 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
+
+ ! 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), tracer2ndDer)
+
+ call InterpolateCubicSpline( &
+ posZMidZLevel, tracersIn, tracer2ndDer, maxLevelCell(iCell), &
+ posZTopZLevel, tracersOut, maxLevelCell(iCell)-1 )
+
+ tracerTop(iTracer,2:maxLevelCell(iCell),iCell) = tracersOut(1:maxLevelCell(iCell)-1)
+
+ end do
+ end do
+
+ deallocate(tracer2ndDer)
+ deallocate(tracersIn,tracersOut, posZMidZLevel, posZTopZLevel)
+
+ 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 config_vert_tracer_adv_order=2,3,4 or'
+ print *, 'config_vert_tracer_adv=''spline'' and config_vert_tracer_adv_order=2,3'
+ call dmpar_abort(dminfo)
+
+ endif ! vertical tracer advection method
+
+ do iCell=1,nCellsSolve
+ 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 ,iCell) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1,iCell))
+ end do
+ end do
+ end do
+
+ deallocate(tracerTop)
+
+ endif ! ZLevel
+
!
! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="gray">abla \phi)
!
@@ -1056,7 +1208,7 @@
allocate(fluxVertTop(num_tracers,nVertLevels+1))
fluxVertTop(:,1) = 0.0
- do iCell=1,grid % nCellsSolve
+ do iCell=1,nCellsSolve
do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
@@ -1354,17 +1506,19 @@
! Compute kinetic energy in each cell
!
ke(:,:) = 0.0
- do iCell=1,nCells
- do i=1,nEdgesOnCell(iCell)
- iEdge = edgesOnCell(i,iCell)
- do k=1,maxLevelCell(iCell)
- ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
- end do
- end do
- do k=1,maxLevelCell(iCell)
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeBot(iEdge)
+ ke(k,cell1) = ke(k,cell1) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ ke(k,cell2) = ke(k,cell2) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ enddo
+ end do
+ do iCell = 1,nCells
+ do k = 1,maxLevelCell(iCell)
ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
- end do
- end do
+ enddo
+ enddo
!
! Compute v (tangential) velocities
@@ -1373,7 +1527,9 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- do k = 1,maxLevelEdgeBot(iEdge)
+ ! mrp 101115 note: in order to include flux boundary conditions,
+ ! the following loop may need to change to maxLevelEdgeBot
+ do k = 1,maxLevelEdgeTop(iEdge)
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
end do
@@ -1539,7 +1695,7 @@
endif
!
- ! vertical velocity
+ ! vertical velocity through layer interface
!
if (config_vert_grid_type.eq.'isopycnal') then
! set vertical velocity to zero in isopycnal case
@@ -1556,7 +1712,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(iEdge)
+ do k=2,maxLevelEdgeBot(iEdge)
flux = u(k,iEdge) * dvEdge(iEdge)
div_u(k,cell1) = div_u(k,cell1) + flux
div_u(k,cell2) = div_u(k,cell2) - flux
@@ -1564,16 +1720,14 @@
end do
do iCell=1,nCells
- do k=1,maxLevelCell(iCell)
- div_u(k,iCell) = div_u(k,iCell) / areaCell(iCell)
- end do
-
- ! Vertical velocity at bottom is zero.
+ ! Vertical velocity through layer interface at top and
+ ! bottom is zero.
+ wTop(1,iCell) = 0.0
wTop(maxLevelCell(iCell)+1,iCell) = 0.0
- do k=maxLevelCell(iCell),1,-1
- wTop(k,iCell) = wTop(k+1,iCell) - div_u(k,iCell)*h(k,iCell)
+ do k=maxLevelCell(iCell),2,-1
+ wTop(k,iCell) = wTop(k+1,iCell) &
+ - div_u(k,iCell)/areaCell(iCell)*h(k,iCell)
end do
-
end do
deallocate(div_u)
Modified: trunk/mpas/src/operators/Makefile
===================================================================
--- trunk/mpas/src/operators/Makefile 2011-02-23 21:50:59 UTC (rev 745)
+++ trunk/mpas/src/operators/Makefile 2011-02-25 00:05:27 UTC (rev 746)
@@ -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
Copied: trunk/mpas/src/operators/module_spline_interpolation.F (from rev 745, branches/ocean_projects/vert_adv_mrp/src/operators/module_spline_interpolation.F)
===================================================================
--- trunk/mpas/src/operators/module_spline_interpolation.F (rev 0)
+++ trunk/mpas/src/operators/module_spline_interpolation.F 2011-02-25 00:05:27 UTC (rev 746)
@@ -0,0 +1,427 @@
+module spline_interpolation
+
+ implicit none
+
+ private
+
+ public :: CubicSplineCoefficients, InterpolateCubicSpline, &
+ 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.
+
+! IntegrateCubicSpline: Compute a single integral from spline data.
+
+! IntegrateColumnCubicSpline: Compute multiple integrals from spline data.
+
+! InterpolateLinear: Compute linear interpolation.
+
+! TestInterpolate: Test spline interpolation subroutines.
+
+ contains
+
+ subroutine CubicSplineCoefficients(x,y,n,y2ndDer)
+
+! 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.
+
+! INPUT PARAMETERS:
+
+ integer, intent(in) :: &
+ n ! number of nodes
+ real(kind=RKIND), intent(in), dimension(n) :: &
+ x, &! location of nodes
+ y ! value at nodes
+
+! OUTPUT PARAMETERS:
+
+ real(kind=RKIND), intent(out), dimension(n) :: &
+ y2ndDer ! dy^2/dx^2 at each node
+
+! local variables:
+
+ integer :: i
+ real(kind=RKIND) :: &
+ temp,xRatio,a(n)
+
+ y2ndDer(1)=0.0
+ y2ndDer(n)=0.0
+ a(1)=0.0
+
+ 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
+
+ do i=n-1,1,-1
+ y2ndDer(i)=y2ndDer(i)*y2ndDer(i+1)+a(i)
+ enddo
+
+ end subroutine CubicSplineCoefficients
+
+
+ subroutine InterpolateCubicSpline( &
+ 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 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.
+
+! INPUT PARAMETERS:
+
+ real (kind=RKIND), dimension(n), intent(in) :: &
+ x, &! node location, input grid
+ y, &! interpolation variable, input grid
+ y2ndDer ! 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)*y2ndDer(kIn) + (b**3-b)*y2ndDer(kIn+1)) &
+ *(h**2)/6.0
+
+ kOut = kOut + 1
+
+ if (kOut>nOut) exit kInLoop
+
+ enddo
+
+ enddo kInLoop
+
+end subroutine InterpolateCubicSpline
+
+
+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 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.
+! 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
+ y2ndDer ! 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 + 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 &
+ + 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 - 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 &
+ + 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
+
+ enddo ! j
+
+ end subroutine IntegrateCubicSpline
+
+
+ subroutine IntegrateColumnCubicSpline( &
+ 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 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
+! 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
+ y2ndDer ! 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 + y2ndDer(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 - 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 + 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 &
+ + 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
+
+ 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 &
+ + 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
+
+ end subroutine IntegrateColumnCubicSpline
+
+
+ subroutine InterpolateLinear( &
+ 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 InterpolateLinear
+
+
+ subroutine TestInterpolate
+
+! Test function to show how to operate the cubic spline subroutines
+
+ integer, parameter :: &
+ n = 10
+ real (kind=RKIND), dimension(n) :: &
+ y, x, y2ndDer
+
+ 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
+ ! 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, y2ndDer.
+ call CubicSplineCoefficients(x,y,n,y2ndDer)
+
+ ! Compute interpolated values yOut.
+ call InterpolateCubicSpline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
+
+ ! 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')"
+
+ ! Compute interpolated values yOut.
+ call IntegrateColumnCubicSpline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
+
+ ! 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>