<p><b>dwj07@fsu.edu</b> 2012-10-19 15:02:11 -0600 (Fri, 19 Oct 2012)</p><p><br>
        -- BRANCH COMMIT --<br>
<br>
        Cleaner version of ocn_wtop at Mark's request.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tendency.F        2012-10-19 20:39:30 UTC (rev 2230)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tendency.F        2012-10-19 21:02:11 UTC (rev 2231)
@@ -869,8 +869,8 @@
real (kind=RKIND), dimension(:), pointer :: &
dvEdge, areaCell, zstarWeight
real (kind=RKIND), dimension(:,:), pointer :: uTransport,h,wTop, h_edge
- real (kind=RKIND), dimension(:,:), allocatable:: div_hu
- real (kind=RKIND), dimension(:), allocatable:: div_hu_btr, h_tend_col
+ real (kind=RKIND), dimension(:), allocatable:: div_hu, h_tend_col
+ real (kind=RKIND) :: div_hu_btr
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
@@ -898,69 +898,57 @@
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
- allocate(div_hu(nVertLevels,nCells+1), div_hu_btr(nCells+1), &
- h_tend_col(nVertLevels))
+ if (config_vert_grid_type.eq.'isopycnal') then
+ ! set vertical velocity to zero in isopycnal case
+ wTop=0.0_RKIND
+ return
+ end if
+
+ allocate(div_hu(nVertLevels), h_tend_col(nVertLevels))
+
!
! Compute div(h^{edge} u) for each cell
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
!
- div_hu(:,:) = 0.0_RKIND
do iCell=1,nCells
+ div_hu(:) = 0.0_RKIND
+ div_hu_btr = 0.0_RKIND
+ hSum = 0.0_RKIND
invAreaCell = 1.0_RKIND / areaCell(iCell)
+
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(i, iCell)
do k = 1, maxLevelEdgeBot(iEdge)
flux = uTransport(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
- div_hu(k, iCell) = div_hu(k, iCell) - edgeSignOnCell(i, iCell) * flux * invAreaCell
+ flux = edgeSignOnCell(i, iCell) * flux * invAreaCell
+ div_hu(k) = div_hu(k) - flux
+ div_hu_btr = div_hu_btr - flux
end do
end do
- end do
- do iCell=1,nCells
- div_hu_btr(iCell) = 0.0_RKIND
- do k=1,maxLevelCell(iCell)
- div_hu_btr(iCell) = div_hu_btr(iCell) + div_hu(k,iCell)
- end do
- end do
+ do k = 1, maxLevelCell(iCell)
+ h_tend_col(k) = - zstarWeight(k) * h(k, iCell) * div_hu_btr
+ hSum = hSum + zstarWeight(k) * h(k, iCell)
+ end do
- !
- ! vertical velocity through layer interface
- !
- !dwj: 10/25/2011 - Need to explore isopycnal vs zlevel flags
- if (config_vert_grid_type.eq.'isopycnal') then
- ! set vertical velocity to zero in isopycnal case
- wTop=0.0_RKIND
+ if(hSum > 0.0) then
+ h_tend_col = h_tend_col / hSum
+ end if
- else ! zlevel or zstar type vertical grid
-
- do iCell=1,nCells
-
- hSum = 0.0_RKIND
- do k=1,maxLevelCell(iCell)
- h_tend_col(k) = - zstarWeight(k)*h(k,iCell)*div_hu_btr(iCell)
- hSum = hSum + zstarWeight(k)*h(k,iCell)
- end do
- if(hSum > 0.0) then
- h_tend_col = h_tend_col / hSum
- else
- end if
-
- ! Vertical velocity 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
- wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k,iCell) - h_tend_col(k)
- end do
+ ! Vertical velocity 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
+ wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k) - h_tend_col(k)
end do
+ end do
- endif
+ deallocate(div_hu, h_tend_col)
- deallocate(div_hu, div_hu_btr, h_tend_col)
-
end subroutine ocn_wtop!}}}
!***********************************************************************
</font>
</pre>