<p><b>duda</b> 2010-01-06 15:17:15 -0700 (Wed, 06 Jan 2010)</p><p>Add del4 mixing for velocity and theta, with eddy viscosities<br>
specified via namelist parameters config_h_mom_eddy_visc4 and <br>
config_h_theta_eddy_visc4.<br>
<br>
M src/module_time_integration.F<br>
M Registry/Registry<br>
M namelist.input<br>
</p><hr noshade><pre><font color="gray">Modified: branches/hyd_model/Registry/Registry
===================================================================
--- branches/hyd_model/Registry/Registry 2010-01-06 03:33:00 UTC (rev 106)
+++ branches/hyd_model/Registry/Registry 2010-01-06 22:17:15 UTC (rev 107)
@@ -7,8 +7,10 @@
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
namelist real sw_model config_h_mom_eddy_visc2 0.0
+namelist real sw_model config_h_mom_eddy_visc4 0.0
namelist real sw_model config_v_mom_eddy_visc2 0.0
namelist real sw_model config_h_theta_eddy_visc2 0.0
+namelist real sw_model config_h_theta_eddy_visc4 0.0
namelist real sw_model config_v_theta_eddy_visc2 0.0
namelist integer sw_model config_number_of_sub_steps 4
namelist integer sw_model config_theta_adv_order 2
Modified: branches/hyd_model/namelist.input
===================================================================
--- branches/hyd_model/namelist.input 2010-01-06 03:33:00 UTC (rev 106)
+++ branches/hyd_model/namelist.input 2010-01-06 22:17:15 UTC (rev 107)
@@ -6,8 +6,10 @@
config_output_interval = 24
config_number_of_sub_steps = 4
config_h_mom_eddy_visc2 = 0.0
+ config_h_mom_eddy_visc4 = 0.0
config_v_mom_eddy_visc2 = 0.0
config_h_theta_eddy_visc2 = 0.0
+ config_h_theta_eddy_visc4 = 0.0
config_v_theta_eddy_visc2 = 0.0
config_theta_adv_order = 2
config_scalar_adv_order = 2
Modified: branches/hyd_model/src/module_time_integration.F
===================================================================
--- branches/hyd_model/src/module_time_integration.F 2010-01-06 03:33:00 UTC (rev 106)
+++ branches/hyd_model/src/module_time_integration.F 2010-01-06 22:17:15 UTC (rev 107)
@@ -118,21 +118,29 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % cqu % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % h % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % h % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pressure % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pressure % array(:,:), &
block % mesh % nVertLevels+1, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % geopotential % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % geopotential % array(:,:), &
block % mesh % nVertLevels+1, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % alpha % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % alpha % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pv_edge % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % pv_edge % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ if (config_h_mom_eddy_visc4 > 0.0) then
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % divergence % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % vorticity % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ end if
block => block % next
end do
@@ -497,8 +505,8 @@
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, upstream_bias
integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND) :: h_mom_eddy_visc2, v_mom_eddy_visc2
- real (kind=RKIND) :: h_theta_eddy_visc2, v_theta_eddy_visc2
+ real (kind=RKIND) :: h_mom_eddy_visc2, v_mom_eddy_visc2, h_mom_eddy_visc4
+ real (kind=RKIND) :: h_theta_eddy_visc2, v_theta_eddy_visc2, h_theta_eddy_visc4
real (kind=RKIND) :: u_diffusion
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, p_s
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
@@ -509,12 +517,14 @@
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND), dimension( grid % nVertLevels + 1 ) :: wdtn, wdun
- real (kind=RKIND) :: theta_edge, theta_turb_flux, z1, z2, z3, z4, zm, z0, zp
+ real (kind=RKIND) :: theta_edge, theta_turb_flux, z1, z2, z3, z4, zm, z0, zp, r
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:), pointer :: rdnu, rdnw, fnm, fnp
-! real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: delsq_theta
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_theta, delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
h => s % h % array
u => s % u % array
@@ -553,10 +563,13 @@
nCells = grid % nCells
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
+ nVertices = grid % nVertices
h_mom_eddy_visc2 = config_h_mom_eddy_visc2
+ h_mom_eddy_visc4 = config_h_mom_eddy_visc4
v_mom_eddy_visc2 = config_v_mom_eddy_visc2
h_theta_eddy_visc2 = config_h_theta_eddy_visc2
+ h_theta_eddy_visc4 = config_h_theta_eddy_visc4
v_theta_eddy_visc2 = config_v_theta_eddy_visc2
@@ -620,7 +633,6 @@
! horizontal mixing for u
!
if ( h_mom_eddy_visc2 > 0.0 ) then
-
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -640,7 +652,104 @@
tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
end do
end do
+ end if
+ if ( h_mom_eddy_visc4 > 0.0 ) then
+
+ allocate(delsq_divergence(nVertLevels, nCells))
+ allocate(delsq_u(nVertLevels, nEdges))
+ allocate(delsq_circulation(nVertLevels, nVertices))
+ allocate(delsq_vorticity(nVertLevels, nVertices))
+
+ delsq_u(:,:) = 0.0
+
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ if (cell1 > 0 .and. cell2 > 0) then
+ do k=1,nVertLevels
+
+ !
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="blue">abla vorticity
+ ! only valid for h_mom_eddy_visc4 == constant
+ !
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+
+ delsq_u(k,iEdge) = delsq_u(k,iEdge) + u_diffusion
+ end do
+ end if
+ end do
+
+ delsq_circulation(:,:) = 0.0
+ do iEdge=1,nEdges
+ if (verticesOnEdge(1,iEdge) > 0) then
+ do k=1,nVertLevels
+ delsq_circulation(k,verticesOnEdge(1,iEdge)) = delsq_circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * delsq_u(k,iEdge)
+ end do
+ end if
+ if (verticesOnEdge(2,iEdge) > 0) then
+ do k=1,nVertLevels
+ delsq_circulation(k,verticesOnEdge(2,iEdge)) = delsq_circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * delsq_u(k,iEdge)
+ end do
+ end if
+ end do
+ do iVertex=1,nVertices
+ r = 1.0 / areaTriangle(iVertex)
+ do k=1,nVertLevels
+ delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
+ end do
+ end do
+
+ delsq_divergence(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0) then
+ do k=1,nVertLevels
+ delsq_divergence(k,cell1) = delsq_divergence(k,cell1) + delsq_u(k,iEdge)*dvEdge(iEdge)
+ end do
+ end if
+ if(cell2 > 0) then
+ do k=1,nVertLevels
+ delsq_divergence(k,cell2) = delsq_divergence(k,cell2) - delsq_u(k,iEdge)*dvEdge(iEdge)
+ end do
+ end if
+ end do
+ do iCell = 1,nCells
+ r = 1.0 / areaCell(iCell)
+ do k = 1,nVertLevels
+ delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
+ end do
+ end do
+
+ 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 h_mom_eddy_visc4 == constant
+ !
+ u_diffusion = ( delsq_divergence(k,cell2) - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( delsq_vorticity(k,vertex2) - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+
+ tend_u(k,iEdge) = tend_u(k,iEdge) - h_mom_eddy_visc4 * u_diffusion
+ end do
+ end do
+
+ deallocate(delsq_divergence)
+ deallocate(delsq_u)
+ deallocate(delsq_circulation)
+ deallocate(delsq_vorticity)
+
end if
@@ -699,7 +808,6 @@
tend_theta(:,:) = 0.
-! delsq_theta(:,:) = 0.
!
! horizontal mixing for theta - we could combine this with advection directly (i.e. as a turbulent flux),
@@ -717,14 +825,57 @@
flux = dvEdge (iEdge) * h_edge(k,iEdge) * theta_turb_flux
tend_theta(k,cell1) = tend_theta(k,cell1) + flux
tend_theta(k,cell2) = tend_theta(k,cell2) - flux
+ end do
-! delsq_theta(k,cell1) = delsq_theta(k,iCell) + dvEdge(iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
-! delsq_theta(k,cell2) = delsq_theta(k,iCell) - dvEdge(iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ end if
+ end do
+
+ end if
+
+ if ( h_theta_eddy_visc4 > 0.0 ) then
+
+ allocate(delsq_theta(nVertLevels, nCells))
+
+ delsq_theta(:,:) = 0.
+
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+
+ do k=1,grid % nVertLevels
+ delsq_theta(k,cell1) = delsq_theta(k,cell1) + dvEdge(iEdge)*h_edge(k,iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
+ delsq_theta(k,cell2) = delsq_theta(k,cell2) - dvEdge(iEdge)*h_edge(k,iEdge)*(theta(k,cell2) - theta(k,cell1))/dcEdge(iEdge)
end do
end if
end do
+ do iCell = 1, nCells
+ r = 1.0 / areaCell(iCell)
+ do k=1,nVertLevels
+ delsq_theta(k,iCell) = delsq_theta(k,iCell) * r
+ end do
+ end do
+
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+
+ do k=1,grid % nVertLevels
+ theta_turb_flux = h_theta_eddy_visc4*prandtl*(delsq_theta(k,cell2) - delsq_theta(k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * theta_turb_flux
+
+ tend_theta(k,cell1) = tend_theta(k,cell1) - flux
+ tend_theta(k,cell2) = tend_theta(k,cell2) + flux
+ end do
+
+ end if
+ end do
+
+ deallocate(delsq_theta)
+
end if
</font>
</pre>