<p><b>duda</b> 2011-04-06 18:18:11 -0600 (Wed, 06 Apr 2011)</p><p>BRANCH COMMIT<br>
<br>
Add correction to the computation of ke in non-hydrostatic <br>
model to avoid Hollingsworth instability.<br>
<br>
Note: Once we have subTriangleAreasOnVertex in our regular<br>
mesh files, the corrections can be made more exact.<br>
<br>
<br>
M module_time_integration.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F        2011-04-06 00:07:21 UTC (rev 785)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_time_integration.F        2011-04-07 00:18:11 UTC (rev 786)
@@ -3028,7 +3028,12 @@
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ !WCS-instability
+ logical, parameter :: hollingsworth=.true.
+ real (kind=RKIND), allocatable, dimension(:,:) :: ke_vertex
+ real (kind=RKIND) :: ke_fact
+
h => s % rho % array
u => s % u % array
v => diag % v % array
@@ -3147,7 +3152,51 @@
end do
end do
+ !WCS-instability
+ ! Compute ke at cell vertices - AG's new KE construction, part 1
+ ! *** approximation here because we don't have inner triangle areas
!
+ if (hollingsworth) then
+ allocate (ke_vertex(nVertLevels,nVertices))
+ do iVertex=1,nVertices
+ do k=1,nVertLevels
+! ke_vertex(k,iVertex) = ( subTriangleAreasOnVertex(1,iVertex)*u(k,EdgesOnVertex(1,iVertex))**2.0 &
+! + subTriangleAreasOnVertex(2,iVertex)*u(k,EdgesOnVertex(2,iVertex))**2.0 &
+! + subTriangleAreasOnVertex(3,iVertex)*u(k,EdgesOnVertex(3,iVertex))**2.0 &
+! ) / AreaTriangle(iVertex)
+
+ ke_vertex(k,iVertex) = ( dcEdge(EdgesOnVertex(1,iVertex))*dvEdge(EdgesOnVertex(1,iVertex))*u(k,EdgesOnVertex(1,iVertex))**2.0 &
+ +dcEdge(EdgesOnVertex(2,iVertex))*dvEdge(EdgesOnVertex(2,iVertex))*u(k,EdgesOnVertex(2,iVertex))**2.0 &
+ +dcEdge(EdgesOnVertex(3,iVertex))*dvEdge(EdgesOnVertex(3,iVertex))*u(k,EdgesOnVertex(3,iVertex))**2.0 &
+ ) * 0.25 / AreaTriangle(iVertex)
+
+ end do
+ end do
+
+ ! adjust ke at cell vertices - AG's new KE construction, part 2
+ !
+
+ ke_fact = 1.0 - .375
+
+ do iCell=1,nCells
+ do k=1,nVertLevels
+ ke(k,iCell) = ke_fact*ke(k,iCell)
+ end do
+ end do
+
+ do iVertex = 1, nVertices
+ do i=1,grid % vertexDegree
+ iCell = cellsOnVertex(i,iVertex)
+ do k = 1,nVertLevels
+ ke(k,iCell) = ke(k,iCell) + (1.-ke_fact)*kiteAreasOnVertex(i, iVertex) * ke_vertex(k, iVertex) / areaCell(iCell)
+ end do
+ end do
+ end do
+ deallocate (ke_vertex)
+ end if
+ !END of WCS-instability
+
+ !
! Compute v (tangential) velocities
!
v(:,:) = 0.0
</font>
</pre>