<p><b>mpetersen@lanl.gov</b> 2010-05-14 12:09:25 -0600 (Fri, 14 May 2010)</p><p>Changed vertical edge variables to reference the top rather than bottom of cell. All are named *Top and are dimensioned nVertLevels+1. This lets me avoid any special cases for the surface and bottom in my loops.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry 2010-05-14 17:06:19 UTC (rev 274)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry 2010-05-14 18:09:25 UTC (rev 275)
@@ -35,6 +35,7 @@
dim R3 3
dim vertexDegree vertexDegree
dim nVertLevels nVertLevels
+dim nVertLevelsP1 nVertLevels+1
#
# var type name_in_file ( dims ) iro- name_in_code super-array array_class
@@ -93,7 +94,7 @@
var integer maxLevelsEdge ( nEdges ) iro kmaxEdge - -
var real hZLevel ( nVertLevels ) iro hZLevel - -
var real zMidZLevel ( nVertLevels ) iro zMidZLevel - -
-var real zBotZLevel ( nVertLevels ) iro zBotZLevel - -
+var real zTopZLevel ( nVertLevelsP1 ) iro zTopZLevel - -
# Boundary conditions: read from input, saved in restart and written to output
var integer boundaryEdge ( nVertLevels nEdges ) iro boundaryEdge - -
@@ -122,16 +123,14 @@
var real uReconstructY ( nVertLevels nCells Time ) o uReconstructY - -
var real uReconstructZ ( nVertLevels nCells Time ) o uReconstructZ - -
var real zMid ( nVertLevels nCells Time ) o zMid - -
-var real zBot ( nVertLevels nCells Time ) o zBot - -
-var real zSurface ( nCells Time ) o zSurface - -
+var real zTop ( nVertLevelsP1 nCells Time ) o zTop - -
var real zMidEdge ( nVertLevels nEdges Time ) o zMidEdge - -
-var real zBotEdge ( nVertLevels nEdges Time ) o zBotEdge - -
-var real zSurfaceEdge ( nEdges Time ) o zSurfaceEdge - -
-var real pMid ( nVertLevels nCells Time ) o pMid - -
-var real pBot ( nVertLevels nCells Time ) o pBot - -
-var real pMidZLevel ( nVertLevels nCells Time ) o pmidZLevel - -
+var real zTopEdge ( nVertLevelsP1 nEdges Time ) o zTopEdge - -
+var real p ( nVertLevels nCells Time ) o p - -
+var real pTop ( nVertLevelsP1 nCells Time ) o pTop - -
+var real pZLevel ( nVertLevels nCells Time ) o pZLevel - -
var real MontPot ( nVertLevels nCells Time ) o MontPot - -
-var real wBot ( nVertLevels nCells Time ) o wBot - -
+var real wTop ( nVertLevelsP1 nCells Time ) o wTop - -
# Other diagnostic variables: neither read nor written to any files
var real vh ( nVertLevels nEdges Time ) - vh - -
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_global_diagnostics.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_global_diagnostics.F 2010-05-14 17:06:19 UTC (rev 274)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_global_diagnostics.F 2010-05-14 18:09:25 UTC (rev 275)
@@ -36,7 +36,7 @@
real (kind=RKIND) :: areaCellGlobal, areaEdgeGlobal, areaTriangleGlobal
real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle, areaEdge
real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, circulation, vorticity, ke, pv_edge, pv_vertex, &
- pv_cell, gradPVn, gradPVt, zmid, zbot, pmid, pbot, MontPot, w, rho
+ pv_cell, gradPVn, gradPVt, zMid, zTop, p, pTop, MontPot, wTop, rho
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
real (kind=RKIND) :: localCFL, localSum
@@ -66,7 +66,7 @@
rho => state % rho % array
tracers => state % tracers % array
v => state % v % array
- w => state % w % array
+ wTop => state % wTop % array
h_edge => state % h_edge % array
circulation => state % circulation % array
vorticity => state % vorticity % array
@@ -76,10 +76,10 @@
pv_cell => state % pv_cell % array
gradPVn => state % gradPVn % array
gradPVt => state % gradPVt % array
- zmid => state % zmid % array
- zbot => state % zbot % array
- pmid => state % pmid % array
- pbot => state % pbot % array
+ zMid => state % zMid % array
+ zTop => state % zTop % array
+ p => state % p % array
+ pTop => state % pTop % array
MontPot => state % MontPot % array
variableIndex = 0
@@ -152,28 +152,28 @@
gradPVt(:,1:nEdgesSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
- ! zmid
+ ! zMid
variableIndex = variableIndex + 1
call computeFieldAreaWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), &
- zmid(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
+ zMid(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
- ! zbot
+ ! zTop
variableIndex = variableIndex + 1
call computeFieldAreaWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), &
- zbot(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
+ zTop(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
- ! pmid
+ ! p
variableIndex = variableIndex + 1
call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
- pmid(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
+ p(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
- ! pbot
+ ! pTop
variableIndex = variableIndex + 1
call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
- pbot(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
+ pTop(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
! MontPot
@@ -182,10 +182,10 @@
MontPot(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
- ! w vertical velocity
+ ! wTop vertical velocity
variableIndex = variableIndex + 1
- call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
- w(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
+ call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels+1, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
+ wTop(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
nVariables = variableIndex
@@ -294,19 +294,19 @@
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeEdgeGlobal
- ! zmid
+ ! zMid
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/(areaCellGlobal*nVertLevels)
- ! zbot
+ ! zTop
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/(areaCellGlobal*nVertLevels)
- ! pmid
+ ! p
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
- ! pbot
+ ! pTop
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
@@ -314,7 +314,7 @@
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
- ! w vertical velocity
+ ! wTop vertical velocity
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-05-14 17:06:19 UTC (rev 274)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-05-14 18:09:25 UTC (rev 275)
@@ -28,7 +28,7 @@
! mrp 100507: for diagnostic output
integer :: iTracer
real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
- hZLevel, zmidZLevel, zbotZLevel
+ hZLevel, zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: delta_rho
@@ -110,8 +110,8 @@
yCell => block_ptr % mesh % yCell % array
hZLevel => block_ptr % mesh % hZLevel % array
- zmidZLevel => block_ptr % mesh % zmidZLevel % array
- zbotZLevel => block_ptr % mesh % zbotZLevel % array
+ zMidZLevel => block_ptr % mesh % zMidZLevel % array
+ zTopZLevel => block_ptr % mesh % zTopZLevel % array
nCells = block_ptr % mesh % nCells
nEdges = block_ptr % mesh % nEdges
@@ -122,12 +122,11 @@
if (config_vert_grid_type.eq.'zlevel') then
! These should eventually be in an input file. For now
! I just read them in from h(:,1).
- hZLevel = h(:,1)
- zmidZLevel(1) = -0.5*hZLevel(1)
- zbotZLevel(1) = -hZLevel(1)
- do iLevel = 2,nVertLevels
- zmidZLevel(iLevel) = zbotZLevel(iLevel-1)-0.5*hZLevel(iLevel)
- zbotZLevel(iLevel) = zbotZLevel(iLevel-1)- hZLevel(iLevel)
+ hZLevel = h(:,1)
+ zTopZLevel(1) = 0.0
+ do iLevel = 1,nVertLevels
+ zMidZLevel(iLevel) = zTopZLevel(iLevel)-0.5*hZLevel(iLevel)
+ zTopZLevel(iLevel+1) = zTopZLevel(iLevel)- hZLevel(iLevel)
enddo
if (config_vert_grid_type.eq.'isopycnal') then
print *, ' Using isopycnal coordinates'
@@ -164,14 +163,14 @@
' min u max u ',&
' min u_src max u_src ', &
' min h max h ',&
- ' hZLevel zmidZlevel', &
- ' zbotZlevel'
+ ' hZLevel zMidZlevel', &
+ ' zTopZlevel'
do iLevel = 1,nVertLevels
print '(i5,20es12.4)', ilevel, rho(ilevel,1), &
minval(u(iLevel,1:nEdges)), maxval(u(iLevel,1:nEdges)), &
minval(u_src(iLevel,1:nEdges)), maxval(u_src(iLevel,1:nEdges)), &
minval(h(iLevel,1:nCells)), maxval(h(iLevel,1:nCells)), &
- hZLevel(iLevel),zmidZlevel(iLevel),zbotZlevel(iLevel)
+ hZLevel(iLevel),zMidZlevel(iLevel),zTopZlevel(iLevel)
enddo
! print some diagnostics
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-05-14 17:06:19 UTC (rev 274)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-05-14 18:09:25 UTC (rev 275)
@@ -229,7 +229,7 @@
! xsad 10-02-09:
! commenting this out until we incorporate the necessary lapack routines into mpas
- !call reconstruct(block % time_levs(2) % state, block % mesh)
+ call reconstruct(block % time_levs(2) % state, block % mesh)
! xsad 10-02-09 end
block => block % next
@@ -258,22 +258,21 @@
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
- upstream_bias, wBotEdge, rho0Inv
+ upstream_bias, wTopEdge, rho0Inv
real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zSurfaceEdge, &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
zMidZLevel
- real (kind=RKIND), dimension(grid % nVertLevels) :: w_dudzBotEdge
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, wBot, &
- tend_h, tend_u, hFluxBot, circulation, vorticity, ke, pv_edge, &
- divergence, MontPot, pmidZLevel, zMidEdge, zBotEdge
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, wTop, &
+ tend_h, tend_u, hFluxTop, circulation, vorticity, ke, pv_edge, &
+ divergence, MontPot, pZLevel, zMidEdge, zTopEdge
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
integer, dimension(:,:), pointer :: &
cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
edgesOnEdge, edgesOnVertex
- real (kind=RKIND) :: u_diffusion, visc, dist
- real (kind=RKIND), dimension(:), allocatable:: fluxVert
+ real (kind=RKIND) :: u_diffusion, visc
+ real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
@@ -283,7 +282,7 @@
h => s % h % array
u => s % u % array
v => s % v % array
- wBot => s % wBot % array
+ wTop => s % wTop % array
h_edge => s % h_edge % array
circulation => s % circulation % array
vorticity => s % vorticity % array
@@ -291,9 +290,8 @@
ke => s % ke % array
pv_edge => s % pv_edge % array
MontPot => s % MontPot % array
- pmidZLevel => s % pmidZLevel % array
- zBotEdge => s % zBotEdge % array
- zSurfaceEdge=> s % zSurfaceEdge % array
+ pZLevel => s % pZLevel % array
+ zTopEdge => s % zTopEdge % array
zMidEdge => s % zMidEdge % array
weightsOnEdge => grid % weightsOnEdge % array
@@ -313,11 +311,11 @@
h_s => grid % h_s % array
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
- zmidZLevel => grid % zmidZLevel % array
+ zMidZLevel => grid % zMidZLevel % array
tend_h => tend % h % array
tend_u => tend % u % array
- hFluxBot => tend % wBot % array
+ hFluxTop => tend % wTop % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -358,65 +356,61 @@
! height tendency: vertical advection term -d/dz(hw)
!
if (config_vert_grid_type.eq.'isopycnal') then
- hFluxBot=0.0 ! vertical fluxes are zero
+ hFluxTop=0.0 ! vertical fluxes are zero
elseif (config_vert_grid_type.eq.'zlevel') then
+
do iCell=1,nCells
+ hFluxTop(1,iCell) = 0.0
+ hFluxTop(nVertLevels+1,iCell) = 0.0
- ! change nvertlevels to kmt later
- ! this next line can be set permanently somewhere upon startup.
- hFluxBot(nVertLevels,iCell) = 0.0
-
+ ! compute hw, the vertical height flux, using
+ ! hw = -integral_zbot^ztop ( div(hu) )
do k=nVertLevels,2,-1
+ ! I use +tend_h because tend_h is -div.(hu)
+ hFluxTop(k,iCell) = hFluxTop(k+1,iCell) + tend_h(k,iCell)*h(k,iCell)
+ end do
- ! F^v_{k-1/2} =
- ! F^v_{k+1/2} - </font>
<font color="red">abla_h \cdot \left( F^h \right) h_k
- ! note +tend_h because tend_h is -div.(hu)
- hFluxBot(k-1,iCell) = hFluxBot(k,iCell) + tend_h(k,iCell)*h(k,iCell)
-
- ! now tend_h is div.(hu) + d/dz(hw):
- ! - </font>
<font color="gray">abla_h \cdot F^h_k
- ! - \frac{F^v_{k-1/2} - F^v_{k+1/2}}{ h^n_k}
- ! note if you substitute line above this is just tend_h=0.
- ! so this line can be changed to tend_h=0 in the future.
+ ! add -d/dz(hw) to heigh tendency
+ ! note if you substitute this into the line above this is
+ ! just tend_h=0 except for the top layer.
+ ! So this loop could be changed to just k=1 in the future.
+ do k=1,nVertLevels
tend_h(k,iCell) = tend_h(k,iCell) &
- - (hFluxBot(k-1,iCell) - hFluxBot(k,iCell))/h(k,iCell)
+ - (hFluxTop(k,iCell) - hFluxTop(k+1,iCell))/h(k,iCell)
end do
- k=1
- ! Surface tracer hFluxBot(k,iCell) = 0.0 is not used
- tend_h(k,iCell) = tend_h(k,iCell) &
- - ( - hFluxBot(k,iCell))/h(k,iCell)
-
end do
endif ! coordinate type
!
! velocity tendency: vertical advection term -w du/dz
!
+ allocate(w_dudzTopEdge(nVertLevels+1))
+ w_dudzTopEdge(1) = 0.0
+ w_dudzTopEdge(nVertLevels+1) = 0.0
tend_u(:,:) = 0.0
if (config_vert_grid_type.eq.'zlevel') then
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels-1
+ do k=2,nVertLevels
! Average w from cell center to edge
- wBotEdge = 0.5*(wBot(k,cell1)+wBot(k,cell2))
+ wTopEdge = 0.5*(wTop(k,cell1)+wTop(k,cell2))
! compute dudz at vertical interface with first order derivative.
- w_dudzBotEdge(k) = wBotEdge * (u(k,iEdge)-u(k+1,iEdge)) &
- / (zMidZLevel(k) - zMidZLevel(k+1))
+ w_dudzTopEdge(k) = wTopEdge * (u(k-1,iEdge)-u(k,iEdge)) &
+ / (zMidZLevel(k-1) - zMidZLevel(k))
end do
- w_dudzBotEdge(nVertLevels) = 0.0
! Average w*du/dz from vertical interface to vertical middle of cell
- tend_u(1,iEdge) = - 0.5 * w_dudzBotEdge(1)
- do k=2,nVertLevels
- tend_u(k,iEdge) = - 0.5 * (w_dudzBotEdge(k-1) + w_dudzBotEdge(k))
+ do k=1,nVertLevels
+ tend_u(k,iEdge) = - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
enddo
enddo
endif
+ deallocate(w_dudzTopEdge)
!
! velocity tendency: pressure gradient
@@ -437,8 +431,8 @@
cell2 = cellsOnEdge(2,iEdge)
do k=1,nVertLevels
tend_u(k,iEdge) = tend_u(k,iEdge) &
- - rho0Inv*( pmidZLevel(k,cell2) &
- - pmidZLevel(k,cell1) )/dcEdge(iEdge)
+ - rho0Inv*( pZLevel(k,cell2) &
+ - pZLevel(k,cell1) )/dcEdge(iEdge)
end do
enddo
endif
@@ -493,35 +487,30 @@
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
- allocate(fluxVert(0:nVertLevels))
+ allocate(fluxVertTop(1:nVertLevels+1))
+ fluxVertTop(1) = 0.0
+ fluxVertTop(nVertLevels+1) = 0.0
do iEdge=1,grid % nEdgesSolve
- fluxVert(0) = 0.0
- fluxVert(nVertLevels) = 0.0
- do k=nVertLevels-1,1,-1
- fluxVert(k) = config_vert_viscosity &
- * ( u(k,iEdge) - u(k+1,iEdge) ) &
- / (zMidEdge(k,iEdge) -zMidEdge(k+1,iEdge))
+ do k=2,nVertLevels
+ fluxVertTop(k) = config_vert_viscosity &
+ * ( u(k-1,iEdge) - u(k,iEdge) ) &
+ / (zMidEdge(k-1,iEdge) -zMidEdge(k,iEdge))
enddo
do k=1,nVertLevels
- if(k.eq.1) then
- dist = zSurfaceEdge(iEdge) - zBotEdge(k,iEdge)
- else
- dist = zBotEdge(k-1,iEdge) - zBotEdge(k,iEdge)
- endif
-
- tend_u(k,iEdge) = tend_u(k,iEdge) + (fluxVert(k-1) - fluxVert(k))/dist
-
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ + (fluxVertTop(k) - fluxVertTop(k+1)) &
+ /(zTopEdge(k,iEdge) - zTopEdge(k+1,iEdge))
enddo
end do
- deallocate(fluxVert)
+ deallocate(fluxVertTop)
-!print *, 'ncells,grid % nCellsSolve, size(hFluxBot)',ncells,grid % nCellsSolve, size(hFluxBot,1), size(hFluxBot,2)
+!print *, 'ncells,grid % nCellsSolve, size(hFluxTop)',ncells,grid % nCellsSolve, size(hFluxTop,1), size(hFluxTop,2)
! print '(a,i5,f10.5)', &
-! 'k, min(tend_h(k,2:)),max(tend_h(k,2:)),min(hFluxBot(k,:)),max(hFluxBot(k,:))'
+! 'k, min(tend_h(k,2:)),max(tend_h(k,2:)),min(hFluxTop(k,:)),max(hFluxTop(k,:))'
! do k=1,nVertLevels
! print '(i5,10es10.2)', &
! k,minval(tend_h(k,1:grid % ncellssolve)),maxval(tend_h(k,1:grid % ncellssolve)), &
-! minval(hFluxBot(k,1:grid % ncellssolve)),maxval(hFluxBot(k,1:grid % ncellssolve))
+! minval(hFluxTop(k,1:grid % ncellssolve)),maxval(hFluxTop(k,1:grid % ncellssolve))
! enddo
! print '(a,i5,f10.5)', &
@@ -554,25 +543,24 @@
real (kind=RKIND) :: flux, tracer_edge, r
real (kind=RKIND) :: dist
real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zsurface
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: &
- u,h,hFluxBot, h_edge, zmid, zbot
+ u,h,hFluxTop, h_edge, zMid, zTop
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension(:), allocatable:: hFluxBotCol
- real (kind=RKIND), dimension(:,:), allocatable:: fluxVert, tracerBot
+ real (kind=RKIND), dimension(:), allocatable:: hFluxTopCol
+ real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, tracerTop
real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div
u => s % u % array
h => s % h % array
tracers => s % tracers % array
h_edge => s % h_edge % array
- hFluxBot => tend % wBot % array
- zmid => s % zmid % array
- zbot => s % zbot % array
- zsurface => s % zsurface % array
+ hFluxTop => tend % wTop % array
+ zMid => s % zMid % array
+ zTop => s % zTop % array
tend_tr => tend % tracers % array
@@ -585,7 +573,7 @@
nVertLevels = grid % nVertLevels
!
- ! tracer tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( h \phi u)
+ ! tracer tendency: horizontal advection term -div( h \phi u)
!
tend_tr(:,:,:) = 0.0
if (config_hor_tracer_adv.eq.'centered') then
@@ -641,84 +629,47 @@
!
! tracer tendency: vertical advection term -d/dz( h \phi w)
!
- allocate(tracerBot(num_tracers,0:nVertLevels))
- allocate(hFluxBotCol(0:nVertLevels))
- tracerBot(:,0)=0.0
- tracerBot(:,nVertLevels)=0.0
- hFluxBotCol(0)=0.0
- hFluxBotCol(nVertLevels)=0.0
+ allocate(tracerTop(num_tracers,nVertLevels+1))
+ tracerTop(:,1)=0.0
+ tracerTop(:,nVertLevels+1)=0.0
do iCell=1,grid % nCellsSolve
if (config_vert_tracer_adv.eq.'centered') then
- do k=1,nVertLevels-1
+ do k=2,nVertLevels
do iTracer=1,num_tracers
- tracerBot(iTracer,k) = ( tracers(iTracer,k ,iCell) &
- + tracers(iTracer,k+1,iCell))/2.0
+ 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=1,nVertLevels-1
- if (hFluxBot(k,iCell)>0.0) then
- upwindCell = k+1
- else
+ do k=2,nVertLevels
+ if (hFluxTop(k,iCell)>0.0) then
upwindCell = k
+ else
+ upwindCell = k-1
endif
do iTracer=1,num_tracers
- tracerBot(iTracer,k) = tracers(iTracer,upwindCell,iCell)
+ tracerTop(iTracer,k) = tracers(iTracer,upwindCell,iCell)
end do
end do
endif
- hFluxBotCol(1:nVertLevels)=hFluxBot(1:nVertLevels,iCell)
do k=1,nVertLevels
do iTracer=1,num_tracers
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( hFluxBotCol(k-1)*tracerBot(iTracer,k-1) &
- - hFluxBotCol(k )*tracerBot(iTracer,k ))/h(k,icell)
+ - ( hFluxTop(k ,iCell)*tracerTop(iTracer,k ) &
+ - hFluxTop(k+1,iCell)*tracerTop(iTracer,k+1))/h(k,icell)
end do
end do
enddo ! iCell
- deallocate(tracerBot,hFluxBotcol)
+ deallocate(tracerTop)
!
- ! tracer tendency: vertical mixing d/dz( \kappa_v d\phi/dz))
- !
- allocate(fluxVert(num_tracers,0:nVertLevels))
- fluxVert(:,0) = 0.0
- fluxVert(:,nVertLevels) = 0.0
- do iCell=1,grid % nCellsSolve
- do k=1,nVertLevels-1
- do iTracer=1,num_tracers
- fluxVert(iTracer,k) = config_vert_diffusion &
- * (tracers(iTracer,k,iCell) - tracers(iTracer,k+1,iCell) )&
- / (zMid(k,iCell) -zMid(k+1,iCell))
- enddo
- enddo
-
- k=1
- dist = zSurface(iCell) - zBot(k,iCell)
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- + (fluxVert(iTracer,k-1) - fluxVert(iTracer,k))/dist ! orig
- enddo
-
- do k=2,nVertLevels
- dist = zBot(k-1,iCell) - zBot(k,iCell)
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- + (fluxVert(iTracer,k-1) - fluxVert(iTracer,k))/dist
- enddo
- enddo
-
- enddo ! iCell loop
- deallocate(fluxVert)
-
- !
! tracer tendency: horizontal tracer diffusion
- ! </font>
<font color="black">abla\cdot (\kappa_h </font>
<font color="blue">abla\phi )
+ ! div(\kappa_h </font>
<font color="black">abla\phi )
!
! first compute \kappa_h </font>
<font color="gray">abla\phi at horizontal edges.
allocate(tr_flux(num_tracers,nVertLevels,nEdges))
@@ -760,7 +711,7 @@
end if
end do
- ! add </font>
<font color="black">abla\cdot (\kappa_h </font>
<font color="blue">abla\phi ) to tracer tendency
+ ! add div(\kappa_h </font>
<font color="gray">abla\phi ) to tracer tendency
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
do k = 1,nVertLevels
@@ -772,6 +723,33 @@
enddo
deallocate(tr_flux, tr_div)
+ !
+ ! tracer tendency: vertical mixing d/dz( \kappa_v d\phi/dz))
+ !
+ allocate(fluxVertTop(num_tracers,nVertLevels+1))
+ fluxVertTop(:,1) = 0.0
+ fluxVertTop(:,nVertLevels+1) = 0.0
+ do iCell=1,grid % nCellsSolve
+ do k=2,nVertLevels
+ do iTracer=1,num_tracers
+ fluxVertTop(iTracer,k) = config_vert_diffusion &
+ * (tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell) )&
+ / (zMid(k-1,iCell) -zMid(k,iCell))
+ enddo
+ enddo
+
+ do k=1,nVertLevels
+ dist = zTop(k,iCell) - zTop(k+1,iCell)
+ do iTracer=1,num_tracers
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ + (fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1))/dist
+ enddo
+ enddo
+
+ enddo ! iCell loop
+ deallocate(fluxVertTop)
+
+
! print some diagnostics - for debugging
! print *, 'after vertical mixing',&
! 'iTracer,k, minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))'
@@ -803,19 +781,19 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p, pbot_temp
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p, pTop_temp
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- zSurface, hZLevel, zSurfaceEdge
+ hZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wBot, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wTop, &
circulation, vorticity, ke, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
- zmid, zbot, zmidEdge, zbotEdge, pmid, pbot, MontPot, rho, temperature, salinity, pmidZLevel
+ zMid, zTop, zMidEdge, zTopEdge, p, pTop, MontPot, rho, temperature, salinity, pZLevel
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND), dimension(:,:), allocatable:: div_u
character :: c1*6
@@ -828,7 +806,7 @@
h => s % h % array
u => s % u % array
v => s % v % array
- wBot => s % wBot % array
+ wTop => s % wTop % array
h_edge => s % h_edge % array
circulation => s % circulation % array
vorticity => s % vorticity % array
@@ -841,16 +819,14 @@
gradPVt => s % gradPVt % array
rho => s % rho % array
tracers => s % tracers % array
- zmid => s % zmid % array
- zbot => s % zbot % array
- zmidEdge => s % zmidEdge % array
- zbotEdge => s % zbotEdge % array
- pmid => s % pmid % array
- pmidZLevel => s % pmidZLevel % array
- pbot => s % pbot % array
+ zMid => s % zMid % array
+ zTop => s % zTop % array
+ zMidEdge => s % zMidEdge % array
+ zTopEdge => s % zTopEdge % array
+ p => s % p % array
+ pZLevel => s % pZLevel % array
+ pTop => s % pTop % array
MontPot => s % MontPot % array
- zSurface => s % zSurface % array
- zSurfaceEdge=> s % zSurfaceEdge % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -1107,27 +1083,25 @@
! h_s is ocean topography: elevation above lowest point,
! and is positive. z coordinates are pos upward, and z=0
! is at lowest ocean point.
- zBot(nVertLevels,iCell) = h_s(iCell)
- zMid(nVertLevels,iCell) = zBot(nVertLevels,iCell) + 0.5*h(nVertLevels,iCell)
- do k=nVertLevels-1,1,-1
- zBot(k,iCell) = zBot(k+1,iCell) + h(k+1,iCell)
- zMid(k,iCell) = zBot(k,iCell) + 0.5*h(k,iCell)
+ zTop(nVertLevels+1,iCell) = h_s(iCell)
+ do k=nVertLevels,1,-1
+ zMid(k,iCell) = zTop(k+1,iCell) + 0.5*h(k,iCell)
+ zTop(k,iCell) = zTop(k+1,iCell) + h(k,iCell)
end do
- ! rather than using zBot(0,iCell), I am adding another 2D variable.
- zSurface(iCell) = zBot(1,iCell) + h(1,iCell)
! assume atmospheric pressure at the surface is zero for now.
- pmid(1,iCell) = 0.5*rho(1,iCell)*gravity* h(1,iCell)
- pbot(1,iCell) = rho(1,iCell)*gravity* h(1,iCell)
- MontPot(1,iCell) = gravity * zSurface(iCell)
- do k=2,nVertLevels
+ pTop(1,iCell) = 0.0
+ do k=1,nVertLevels
delta_p = rho(k,iCell)*gravity* h(k,iCell)
- pmid(k,iCell) = pbot(k-1,iCell) + 0.5*delta_p
- pbot(k,iCell) = pbot(k-1,iCell) + delta_p
+ p(k ,iCell) = pTop(k,iCell) + 0.5*delta_p
+ pTop(k+1,iCell) = pTop(k,iCell) + delta_p
+ end do
+ MontPot(1,iCell) = gravity * zTop(1,iCell)
+ do k=2,nVertLevels
! from delta M = p delta / rho
MontPot(k,iCell) = MontPot(k-1,iCell) &
- + pbot(k-1,iCell)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
+ + pTop(k,iCell)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
end do
end do
@@ -1136,31 +1110,32 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if(cell1<=nCells .and. cell2<=nCells) then
- do k=1,nVertLevels
- zBotEdge(k,iEdge) = (zBot(k,cell1)+zBot(k,cell2))/2.0
- zMidEdge(k,iEdge) = (zMid(k,cell1)+zMid(k,cell2))/2.0
- enddo
- zSurfaceEdge(iEdge) = (zSurface(cell1) + zSurface(cell2))/2.0
+ do k=1,nVertLevels
+ zTopEdge(k,iEdge) = (zTop(k,cell1)+zTop(k,cell2))/2.0
+ zMidEdge(k,iEdge) = (zMid(k,cell1)+zMid(k,cell2))/2.0
+ enddo
+ zTopEdge(nVertLevels+1,iEdge) = ( zTop(nVertLevels+1,cell1) &
+ + zTop(nVertLevels+1,cell2))/2.0
endif
enddo
! For z-level model.
! Compute pressure at middle of each level.
- ! pmidZLevel and pmid should only differ at k=1, where pmid is
- ! pressure at middle of layer including SSH, and pmidZLevel is
+ ! pZLevel and p should only differ at k=1, where p is
+ ! pressure at middle of layer including SSH, and pZLevel is
! pressure at a depth of hZLevel(1)/2.
!
do iCell=1,nCells
! compute pressure for z-level coordinates
! assume atmospheric pressure at the surface is zero for now.
- pmidZLevel(1,iCell) = rho(1,iCell)*gravity &
+ pZLevel(1,iCell) = rho(1,iCell)*gravity &
* (h(1,iCell)-0.5*hZLevel(1))
- pbot_temp = rho(1,iCell)*gravity*h(1,iCell)
+ pTop_temp = rho(1,iCell)*gravity*h(1,iCell)
do k=2,nVertLevels
delta_p = rho(k,iCell)*gravity*hZLevel(k)
- pmidZLevel(k,iCell) = pmidZLevel(k-1,iCell) + 0.5*delta_p
- pbot_temp = pbot_temp + delta_p
+ pZLevel(k,iCell) = pZLevel(k-1,iCell) + 0.5*delta_p
+ pTop_temp = pTop_temp + delta_p
end do
end do
@@ -1168,7 +1143,7 @@
! compute vertical velocity
if (config_vert_grid_type.eq.'isopycnal') then
! set vertical velocity to zero in isopycnal case
- wBot=0.0
+ wTop=0.0
elseif (config_vert_grid_type.eq.'zlevel') then
@@ -1200,12 +1175,11 @@
div_u(k,iCell) = div_u(k,iCell) / areaCell(iCell)
end do
- ! change nvertlevels to kmt later
+ ! Vertical velocity at bottom is zero.
! this next line can be set permanently somewhere upon startup.
- wBot(nVertLevels,iCell) = 0.0
-
- do k=nVertLevels,2,-1
- wBot(k-1,iCell) = wBot(k,iCell) - div_u(k,iCell)*h(k,iCell)
+ wTop(nVertLevels+1,iCell) = 0.0
+ do k=nVertLevels,1,-1
+ wTop(k,iCell) = wTop(k+1,iCell) - div_u(k,iCell)*h(k,iCell)
end do
end do
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/mpas_interface.F 2010-05-14 17:06:19 UTC (rev 274)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/mpas_interface.F 2010-05-14 18:09:25 UTC (rev 275)
@@ -27,7 +27,7 @@
! xsad 10-02-09:
! commenting this out until we incorporate the necessary lapack routines into mpas
- ! call init_reconstruct(mesh)
+ call init_reconstruct(mesh)
! xsad 10-02-09 end
! mrp 100316 In order for this to work, we need to pass domain % dminfo as an
</font>
</pre>