<p><b>mpetersen@lanl.gov</b> 2013-02-06 09:38:30 -0700 (Wed, 06 Feb 2013)</p><p>Branch commit, diagnistics_revision. Added computation of vertical velocity. This does not match bit-for-bit with previous revision due to change of order of operations in computing wTop, but it matches to the last two digits. The change is the flux computation in wTop, which I consolidated to a single line for clarity.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/diagnostics_revision/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/diagnostics_revision/src/core_ocean/Registry 2013-02-06 13:53:38 UTC (rev 2443)
+++ branches/ocean_projects/diagnostics_revision/src/core_ocean/Registry 2013-02-06 16:38:30 UTC (rev 2444)
@@ -333,6 +333,7 @@
var persistent real MontPot ( nVertLevels nCells Time ) 2 - MontPot state - -
var persistent real pressure ( nVertLevels nCells Time ) 2 - pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 - wTop state - -
+var persistent real vertVelocityTop ( nVertLevelsP1 nCells Time ) 2 o vertVelocityTop state - -
var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 - rhoDisplaced state - -
var persistent real viscosity ( nVertLevels nEdges Time ) 2 o viscosity state - -
Modified: branches/ocean_projects/diagnostics_revision/src/core_ocean/mpas_ocn_diagnostics.F
===================================================================
--- branches/ocean_projects/diagnostics_revision/src/core_ocean/mpas_ocn_diagnostics.F 2013-02-06 13:53:38 UTC (rev 2443)
+++ branches/ocean_projects/diagnostics_revision/src/core_ocean/mpas_ocn_diagnostics.F 2013-02-06 16:38:30 UTC (rev 2444)
@@ -94,7 +94,7 @@
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2, coef_3rd_order, r_tmp, &
invAreaCell1, invAreaCell2, invAreaTri1, invAreaTri2, invLength, h_vertex
- real (kind=RKIND), dimension(:), allocatable:: pTop
+ real (kind=RKIND), dimension(:), allocatable:: pTop, div_hu
real (kind=RKIND), dimension(:), pointer :: &
bottomDepth, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, ssh
@@ -102,9 +102,8 @@
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
Vor_edge, Vor_vertex, Vor_cell, gradVor_n, gradVor_t, divergence, &
- rho, temperature, salinity, kev, kevc, uBolusGM, uTransport
+ rho, temperature, salinity, kev, kevc, uBolusGM, uTransport, vertVelocityTop
real (kind=RKIND), dimension(:,:,:), pointer :: tracers, deriv_two
- real (kind=RKIND), dimension(:,:), allocatable:: div_u
character :: c1*6
h => s % h % array
@@ -131,6 +130,7 @@
zMid => s % zMid % array
ssh => s % ssh % array
tracers => s % tracers % array
+ vertVelocityTop => s % vertVelocityTop % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -198,6 +198,7 @@
circulation(:,:) = 0.0
vorticity(:,:) = 0.0
divergence(:,:) = 0.0
+ vertVelocityTop(:,:)=0.0
ke(:,:) = 0.0
v(:,:) = 0.0
do iVertex = 1, nVertices
@@ -213,7 +214,9 @@
end do
end do
+ allocate(div_hu(nVertLevels))
do iCell = 1, nCells
+ div_hu(:) = 0.0_RKIND
invAreaCell1 = 1.0 / areaCell(iCell)
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(i, iCell)
@@ -221,10 +224,16 @@
r_tmp = dvEdge(iEdge) * u(k, iEdge) * invAreaCell1
divergence(k, iCell) = divergence(k, iCell) - edgeSignOnCell(i, iCell) * r_tmp
+ div_hu(k) = div_hu(k) - h_edge(k, iEdge) * edgeSignOnCell(i, iCell) * r_tmp
ke(k, iCell) = ke(k, iCell) + 0.25 * r_tmp * dcEdge(iEdge) * u(k,iEdge)
end do
end do
+ ! Vertical velocity at bottom (maxLevelCell(iCell)+1) is zero, initialized above.
+ do k=maxLevelCell(iCell),1,-1
+ vertVelocityTop(k,iCell) = vertVelocityTop(k+1,iCell) - div_hu(k)
+ end do
end do
+ deallocate(div_hu)
do iEdge=1,nEdges
! Compute v (tangential) velocities
@@ -463,12 +472,13 @@
!
! routine ocn_wtop
!
-!> \brief Computes vertical velocity
+!> \brief Computes vertical transport
!> \author Mark Petersen
!> \date 23 September 2011
!> \version SVN:$Id$
!> \details
-!> This routine computes the vertical velocity in the top layer for the ocean
+!> This routine computes the vertical transport through the top of each
+!> cell.
!
!-----------------------------------------------------------------------
subroutine ocn_wtop(grid,h,h_edge,u,wTop, err)!{{{
@@ -489,7 +499,7 @@
h_edge !< Input: h interpolated to an edge
real (kind=RKIND), dimension(:,:), intent(in) :: &
- u !< Input: velocity
+ u !< Input: transport
!-----------------------------------------------------------------
!
@@ -498,7 +508,7 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:,:), intent(out) :: &
- wTop !< Output: vertical transport at top edge
+ wTop !< Output: vertical transport at top of cell
integer, intent(out) :: err !< Output: error flag
@@ -544,7 +554,7 @@
if (config_vert_coord_movement.eq.'isopycnal') then
- ! set vertical velocity to zero in isopycnal case
+ ! set vertical transport to zero in isopycnal case
wTop=0.0_RKIND
return
end if
@@ -566,8 +576,7 @@
iEdge = edgesOnCell(i, iCell)
do k = 1, maxLevelEdgeBot(iEdge)
- flux = u(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
- flux = edgeSignOnCell(i, iCell) * flux * invAreaCell
+ flux = h_edge(k, iEdge) * u(k, iEdge) * dvEdge(iEdge) * edgeSignOnCell(i, iCell) * invAreaCell
div_hu(k) = div_hu(k) - flux
div_hu_btr = div_hu_btr - flux
end do
@@ -582,8 +591,7 @@
h_tend_col = h_tend_col / hSum
end if
- ! Vertical velocity through layer interface at top and
- ! bottom is zero.
+ ! Vertical transport through layer interface at top and bottom is zero.
wTop(1,iCell) = 0.0_RKIND
wTop(maxLevelCell(iCell)+1,iCell) = 0.0_RKIND
do k=maxLevelCell(iCell),2,-1
</font>
</pre>