<p><b>ringler@lanl.gov</b> 2009-10-27 13:49:10 -0600 (Tue, 27 Oct 2009)</p><p><br>
added option for del2 mixing on velocity<br>
u_diffusion is a RHS source term in the velocity tendency equation<br>
coefficient visc is in Registry and namelist, default is zero<br>
added diagnostic field (as output) called divergence, this is computed in compute_solve_diagnostics<br>
<br>
NOTE: The del2 source term is insided "LANL_FORMULATION"<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/swmodel/Registry/Registry
===================================================================
--- trunk/swmodel/Registry/Registry 2009-10-27 16:21:59 UTC (rev 63)
+++ trunk/swmodel/Registry/Registry 2009-10-27 19:49:10 UTC (rev 64)
@@ -6,6 +6,7 @@
namelist real sw_model config_dt 172.8
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
+namelist real sw_model config_visc 0.0
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
@@ -79,6 +80,7 @@
# Diagnostic fields: only written to output
var real v ( nVertLevels nEdges Time ) o v
+var real divergence ( nVertLevels nCells Time ) o divergence
var real vorticity ( nVertLevels nVertices Time ) o vorticity
var real pv_edge ( nVertLevels nEdges Time ) o pv_edge
var real h_edge ( nVertLevels nEdges Time ) o h_edge
Modified: trunk/swmodel/namelist.input
===================================================================
--- trunk/swmodel/namelist.input 2009-10-27 16:21:59 UTC (rev 63)
+++ trunk/swmodel/namelist.input 2009-10-27 19:49:10 UTC (rev 64)
@@ -4,6 +4,7 @@
config_dt = 172.8
config_ntimesteps = 7500
config_output_interval = 500
+ config_visc = 0.0
/
&restart
Modified: trunk/swmodel/src/module_time_integration.F
===================================================================
--- trunk/swmodel/src/module_time_integration.F 2009-10-27 16:21:59 UTC (rev 63)
+++ trunk/swmodel/src/module_time_integration.F 2009-10-27 19:49:10 UTC (rev 64)
@@ -242,17 +242,21 @@
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge
+ circulation, vorticity, ke, pv_edge, divergence
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ real (kind=RKIND) :: u_diffusion, visc
+ visc = config_visc
+
h => s % h % array
u => s % u % array
v => s % v % array
h_edge => s % h_edge % array
circulation => s % circulation % array
vorticity => s % vorticity % array
+ divergence => s % divergence % array
ke => s % ke % array
pv_edge => s % pv_edge % array
@@ -318,19 +322,31 @@
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
do k=1,nVertLevels
+
+ !
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="gray">abla vorticity
+ ! only valid for visc == constant
+ !
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = visc * u_diffusion
+
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(j,iEdge)
workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
end do
- tend_u(k,iEdge) = q - &
- (ke(k,cell2) - ke(k,cell1) + &
- gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
- ) / &
- dcEdge(iEdge)
+ tend_u(k,iEdge) = &
+ q &
+ + u_diffusion &
+ - ( ke(k,cell2) - ke(k,cell1) + &
+ gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
+ ) / dcEdge(iEdge)
end do
end do
#endif
@@ -444,9 +460,10 @@
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt
+ circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ real (kind=RKIND) :: r
h => s % h % array
@@ -458,6 +475,7 @@
tend_u => s % u % array
circulation => s % circulation % array
vorticity => s % vorticity % array
+ divergence => s % divergence % array
ke => s % ke % array
pv_edge => s % pv_edge % array
pv_vertex => s % pv_vertex % array
@@ -524,6 +542,31 @@
end do
!
+ ! Compute the divergence at each cell center
+ !
+ divergence(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0) then
+ do k=1,nVertLevels
+ divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+ enddo
+ endif
+ if(cell2 > 0) then
+ do k=1,nVertLevels
+ divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+ enddo
+ end if
+ end do
+ do iCell = 1,nCells
+ r = 1.0 / areaCell(iCell)
+ do k = 1,nVertLevels
+ divergence(k,iCell) = divergence(k,iCell) * r
+ enddo
+ enddo
+
+ !
! Compute kinetic energy in each cell
!
ke(:,:) = 0.0
</font>
</pre>