<p><b>mpetersen@lanl.gov</b> 2010-11-02 17:03:02 -0600 (Tue, 02 Nov 2010)</p><p>BRANCH COMMIT: Merge MontPot for isopycnal mode and pZLevel for z-level mode into a single p variable.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-11-02 22:11:00 UTC (rev 590)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-11-02 23:03:02 UTC (rev 591)
@@ -168,11 +168,7 @@
var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
var persistent real p ( nVertLevels nCells Time ) 2 o p state - -
-var persistent real pZLevel ( nVertLevels nCells Time ) 2 - pZLevel state - -
-var persistent real MontPot ( nVertLevels nCells Time ) 2 o MontPot state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 o wTop state - -
-var persistent real pgrad ( nVertLevels nCells Time ) 2 o pgrad state - -
-var persistent real pgrad_edge ( nVertLevels nEdges Time ) 2 o pgrad_edge state - -
# Other diagnostic variables: neither read nor written to any files
var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
@@ -180,7 +176,6 @@
var persistent real gradPVt ( nVertLevels nEdges Time ) 2 - gradPVt state - -
var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
-# xsad 10-02-05:
# Globally reduced diagnostic variables: only written to output
var persistent real areaCellGlobal ( Time ) 2 o areaCellGlobal state - -
var persistent real areaEdgeGlobal ( Time ) 2 o areaEdgeGlobal state - -
@@ -189,5 +184,4 @@
var persistent real volumeCellGlobal ( Time ) 2 o volumeCellGlobal state - -
var persistent real volumeEdgeGlobal ( Time ) 2 o volumeEdgeGlobal state - -
var persistent real CFLNumberGlobal ( Time ) 2 o CFLNumberGlobal state - -
-# xsad 10-02-05 end
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/module_global_diagnostics.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/module_global_diagnostics.F 2010-11-02 22:11:00 UTC (rev 590)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_global_diagnostics.F 2010-11-02 23:03:02 UTC (rev 591)
@@ -37,7 +37,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, p, MontPot, wTop, rho, tracerTemp
+ pv_cell, gradPVn, gradPVt, p, wTop, rho, tracerTemp
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
@@ -81,7 +81,6 @@
gradPVn => state % gradPVn % array
gradPVt => state % gradPVt % array
p => state % p % array
- MontPot => state % MontPot % array
variableIndex = 0
! h
@@ -159,12 +158,6 @@
p(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
- ! MontPot
- variableIndex = variableIndex + 1
- call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
- MontPot(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
- verticalSumMaxes(variableIndex))
-
! wTop vertical velocity
variableIndex = variableIndex + 1
call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels+1, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
@@ -292,10 +285,6 @@
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
- ! MontPot
- variableIndex = variableIndex + 1
- averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
-
! wTop vertical velocity
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F 2010-11-02 22:11:00 UTC (rev 590)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F 2010-11-02 23:03:02 UTC (rev 591)
@@ -275,9 +275,9 @@
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, wTop, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, p, wTop, &
tend_h, tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
- divergence, MontPot, pZLevel
+ divergence
type (dm_info) :: dminfo
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
@@ -307,8 +307,7 @@
ke => s % ke % array
ke_edge => s % ke_edge % array
pv_edge => s % pv_edge % array
- MontPot => s % MontPot % array
- pZLevel => s % pZLevel % array
+ p => s % p % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -435,7 +434,7 @@
cell2 = cellsOnEdge(2,iEdge)
do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- - (MontPot(k,cell2) - MontPot(k,cell1))/dcEdge(iEdge)
+ - (p(k,cell2) - p(k,cell1))/dcEdge(iEdge)
end do
enddo
elseif (config_vert_grid_type.eq.'zlevel') then
@@ -444,8 +443,8 @@
cell2 = cellsOnEdge(2,iEdge)
do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- - rho0Inv*( pZLevel(k,cell2) &
- - pZLevel(k,cell1) )/dcEdge(iEdge)
+ - rho0Inv*( p(k,cell2) &
+ - p(k,cell1) )/dcEdge(iEdge)
end do
enddo
endif
@@ -1113,7 +1112,7 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p, rho0Inv
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, pTop, rho0Inv
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
@@ -1122,13 +1121,12 @@
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
hZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wTop, &
- circulation, vorticity, ke, ke_edge, pgrad, pgrad_edge, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, p, wTop, &
+ circulation, vorticity, ke, ke_edge, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
- MontPot, rho, temperature, salinity, pZLevel
+ rho, temperature, salinity
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND), dimension(:,:), allocatable:: div_u
- real (kind=RKIND), dimension(:), allocatable:: pTop
character :: c1*6
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
@@ -1159,10 +1157,7 @@
gradPVt => s % gradPVt % array
rho => s % rho % array
tracers => s % tracers % array
- pZLevel => s % pZLevel % array
- MontPot => s % MontPot % array
- pgrad => s % pgrad % array
- pgrad_edge => s % pgrad_edge % array
+ p => s % p % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -1505,50 +1500,46 @@
! This section must be after computing rho
!
if (config_vert_grid_type.eq.'isopycnal') then
- allocate(pTop(nVertLevels+1))
! For Isopycnal model.
! Compute pressure at top of each layer, and then
! Montgomery Potential.
do iCell=1,nCells
- ! assume atmospheric pressure at the surface is zero for now.
- pTop(1) = 0.0
- do k=1,nVertLevels
- delta_p = rho(k,iCell)*gravity* h(k,iCell)
- pTop(k+1) = pTop(k) + delta_p
- end do
+ ! assume atmospheric pressure at the surface is zero for now.
+ pTop = 0.0
+ ! For isopycnal mode, p is the Montgomery Potential.
+ ! At top layer it is g*SSH, where SSH may be off by a
+ ! constant (ie, h_s can be relative to top or bottom)
+ p(1,iCell) = gravity &
+ * (h_s(iCell) + sum(h(1:nVertLevels,iCell)))
- ! MontPot at top layer is g*SSH, where SSH may be off by a
- ! constant i the horizontal.
- MontPot(1,iCell) = gravity &
- * (h_s(iCell) + sum(h(1:nVertLevels,iCell)))
- do k=2,nVertLevels
- ! from delta M = p delta / rho
- MontPot(k,iCell) = MontPot(k-1,iCell) &
- + pTop(k)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
- end do
+ do k=2,nVertLevels
+ pTop = pTop + rho(k-1,iCell)*gravity* h(k-1,iCell)
+ ! from delta M = p delta / rho
+ p(k,iCell) = p(k-1,iCell) &
+ + pTop*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
+ end do
+
end do
- deallocate(pTop)
elseif (config_vert_grid_type.eq.'zlevel') then
! For z-level model.
! Compute pressure at middle of each level.
- ! 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.
+ ! At k=1, where p is pressure at a depth of hZLevel(1)/2, not
+ ! pressure at middle of layer including SSH.
do iCell=1,nCells
! compute pressure for z-level coordinates
! assume atmospheric pressure at the surface is zero for now.
- pZLevel(1,iCell) = rho(1,iCell)*gravity &
+ p(1,iCell) = rho(1,iCell)*gravity &
* (h(1,iCell)-0.5*hZLevel(1))
do k=2,maxLevelCell(iCell)
- pZLevel(k,iCell) = pZLevel(k-1,iCell) &
+ p(k,iCell) = p(k-1,iCell) &
+ 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
+ rho(k ,iCell)*hZLevel(k ))
end do
@@ -1598,39 +1589,6 @@
endif
- !
- ! Compute pressure gradient in each cell
- !
- ! This is for output and testing only. We could eventually
- ! remove the variables pgrad and pgrad_edge.
- rho0Inv = 1.0/config_rho0
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
- pgrad_edge(k,iEdge) = - rho0Inv*( pZLevel(k,cell2) &
- - pZLevel(k,cell1) )/dcEdge(iEdge)
- end do
- end do
-
- pgrad(:,:) = 0.0
- do iCell=1,nCells
- do i=1,nEdgesOnCell(iCell)
- iEdge = edgesOnCell(i,iCell)
- do k=1,maxLevelEdgeTop(iEdge)
- ! mrp first attempt, I think it is wrong:
- !pgrad(k,iCell) = pgrad(k,iCell) + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * abs(pgrad_edge(k,iEdge))
- pgrad(k,iCell) = pgrad(k,iCell) + abs(pgrad_edge(k,iEdge))
- end do
- end do
- do k=1,maxLevelCell(iCell)
- ! mrp first attempt, I think it is wrong:
-! pgrad(k,iCell) = pgrad(k,iCell) / areaCell(iCell)
- pgrad(k,iCell) = pgrad(k,iCell) / nEdgesOnCell(iCell)
- end do
- end do
-
-
end subroutine compute_solve_diagnostics
</font>
</pre>