<p><b>ringler@lanl.gov</b> 2010-08-05 21:51:58 -0600 (Thu, 05 Aug 2010)</p><p><br>
del2 and del4 options added for velocity and tracers<br>
<br>
this version compiles, but has not been tested.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/dissipation/src/core_ocean/Registry
===================================================================
--- branches/dissipation/src/core_ocean/Registry 2010-08-05 23:38:23 UTC (rev 465)
+++ branches/dissipation/src/core_ocean/Registry 2010-08-06 03:51:58 UTC (rev 466)
@@ -7,7 +7,10 @@
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
namelist integer sw_model config_stats_interval 100
-namelist real sw_model config_visc 0.0
+namelist real hmix config_h_mom_eddy_visc2 0.0
+namelist real hmix config_h_mom_eddy_visc4 0.0
+namelist real hmix config_h_tracer_eddy_visc2 0.0
+namelist real hmix config_h_tracer_eddy_visc4 0.0
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
Modified: branches/dissipation/src/core_ocean/module_time_integration.F
===================================================================
--- branches/dissipation/src/core_ocean/module_time_integration.F 2010-08-05 23:38:23 UTC (rev 465)
+++ branches/dissipation/src/core_ocean/module_time_integration.F 2010-08-06 03:51:58 UTC (rev 466)
@@ -122,6 +122,16 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % intermediate_step(PROVIS) % 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
@@ -252,8 +262,9 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
integer :: nCells, nEdges, nVertices, nVertLevels
+ real (kind=RKIND) :: h_mom_eddy_visc2, h_mom_eddy_visc4
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
- upstream_bias, wTopEdge, rho0Inv
+ upstream_bias, wTopEdge, rho0Inv, r
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
zMidZLevel, zTopZLevel
@@ -270,6 +281,11 @@
real (kind=RKIND) :: u_diffusion
real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge, vertViscTop
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
@@ -319,6 +335,9 @@
u_src => grid % u_src % array
+ h_mom_eddy_visc2 = config_h_mom_eddy_visc2
+ h_mom_eddy_visc4 = config_h_mom_eddy_visc4
+
!
! height tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( hu)
!
@@ -369,13 +388,18 @@
end do
endif ! coordinate type
+
!
+ ! velocity tendency: start accumulating tendency terms
+ !
+ tend_u(:,:) = 0.0
+
+ !
! velocity tendency: vertical advection term -w du/dz
!
allocate(w_dudzTopEdge(nVertLevels+1))
w_dudzTopEdge(1) = 0.0
w_dudzTopEdge(nVertLevels+1) = 0.0
- tend_u(:,:) = 0.0
if (config_vert_grid_type.eq.'zlevel') then
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
@@ -424,23 +448,126 @@
endif
!
- ! velocity tendency: -q(h u^\perp) - </font>
<font color="red">abla K
- ! +</font>
<font color="black">u_h(</font>
<font color="black">abla \delta + {\bf k}\times </font>
<font color="red">abla \xi)
+ ! velocity tendency: del2 and/or del4 dissipation
!
- ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
- ! only valid for visc == constant
+ if ( h_mom_eddy_visc2 > 0.0 ) then
+ 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="blue">abla vorticity
+ ! only valid for h_mom_eddy_visc2 == constant
+ !
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = h_mom_eddy_visc2 * u_diffusion
+
+ 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+1))
+ allocate(delsq_u(nVertLevels, nEdges+1))
+ allocate(delsq_circulation(nVertLevels, nVertices+1))
+ allocate(delsq_vorticity(nVertLevels, nVertices+1))
+
+ 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)
+
+ 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 do
+
+ delsq_circulation(:,:) = 0.0
+ do iEdge=1,nEdges
+ do k=1,nVertLevels
+ delsq_circulation(k,verticesOnEdge(1,iEdge)) = delsq_circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * delsq_u(k,iEdge)
+ delsq_circulation(k,verticesOnEdge(2,iEdge)) = delsq_circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * delsq_u(k,iEdge)
+ end do
+ 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)
+ do k=1,nVertLevels
+ delsq_divergence(k,cell1) = delsq_divergence(k,cell1) + delsq_u(k,iEdge)*dvEdge(iEdge)
+ delsq_divergence(k,cell2) = delsq_divergence(k,cell2) - delsq_u(k,iEdge)*dvEdge(iEdge)
+ end do
+ 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
+
+ !
+ ! velocity tendency: nonlinear Coriolis term and grad of kinetic energy
+ !
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
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = config_visc * u_diffusion
-
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(j,iEdge)
@@ -449,7 +576,6 @@
end do
tend_u(k,iEdge) = tend_u(k,iEdge) &
+ q &
- + u_diffusion &
- ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
end do
@@ -546,6 +672,7 @@
integer :: iCell, iEdge, k, iTracer, cell1, cell2, upwindCell,&
nEdges, nCells, nVertLevels
+ real (kind=RKIND) :: h_tracer_eddy_visc2, h_tracer_eddy_visc4, areaCell1, areaCell2, tracer_turb_flux
real (kind=RKIND) :: flux, tracer_edge, r
real (kind=RKIND) :: dist
real (kind=RKIND), dimension(:), pointer :: &
@@ -561,7 +688,7 @@
real (kind=RKIND), dimension(:), pointer :: &
zTopZLevel
real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, tracerTop
- real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div
+ real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer
real (kind=RKIND), dimension(:), allocatable:: vertDiffTop
@@ -585,6 +712,9 @@
nCells = grid % nCells
nVertLevels = grid % nVertLevels
+ h_tracer_eddy_visc2 = config_h_tracer_eddy_visc2
+ h_tracer_eddy_visc4 = config_h_tracer_eddy_visc4
+
!
! tracer tendency: horizontal advection term -div( h \phi u)
!
@@ -681,61 +811,81 @@
deallocate(tracerTop)
!
- ! tracer tendency: horizontal tracer diffusion
- ! div(h \kappa_h </font>
<font color="red">abla\phi )
+ ! tracer tendency: horizontal tracer diffusion, del2 and/or del4
!
- ! first compute \kappa_h </font>
<font color="red">abla\phi at horizontal edges.
- allocate(tr_flux(num_tracers,nVertLevels,nEdges))
- tr_flux(:,:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
+
+ if ( h_tracer_eddy_visc2 > 0.0 ) then
+
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ areaCell1 = areaCell(cell1)
+ areaCell2 = areaCell(cell2)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1,num_tracers
+ tracer_turb_flux = h_tracer_eddy_visc2*(tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * h_edge(k,iEdge) * tracer_turb_flux
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) + flux/areaCell1
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) - flux/areaCell2
+ enddo
+ end do
+
+ end do
+
+ end if
+
+ if ( h_tracer_eddy_visc4 > 0.0 ) then
+
+ allocate(delsq_tracer(num_tracers,nVertLevels, nCells+1))
+
+ delsq_tracer(:,:,:) = 0.
+
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1,num_tracers
+ delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) + dvEdge(iEdge)*h_edge(k,iEdge)*(tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1))/dcEdge(iEdge)
+ delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) - dvEdge(iEdge)*h_edge(k,iEdge)*(tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1))/dcEdge(iEdge)
+ end do
+ end do
+
+ end do
+
+ do iCell = 1, nCells
+ r = 1.0 / areaCell(iCell)
+ do k=1,nVertLevels
do iTracer=1,num_tracers
- tr_flux(iTracer,k,iEdge) = h_edge(k,iEdge)*config_hor_diffusion * &
- (Tracers(iTracer,k,cell2) - Tracers(iTracer,k,cell1))/dcEdge(iEdge)
+ delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
+ end do
+ end do
+ end do
+
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ areaCell1 = areaCell(cell1)
+ areaCell2 = areaCell(cell2)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1,num_tracers
+ tracer_turb_flux = h_tracer_eddy_visc4*(delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * tracer_turb_flux
+
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell1
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell2
+
+ end do
enddo
- enddo
- endif
- enddo
- ! Compute the divergence, div(h \kappa_h </font>
<font color="red">abla\phi) at cell centers
- allocate(tr_div(num_tracers,nVertLevels,nCells))
- tr_div(:,:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,nVertLevels
- do iTracer=1,num_tracers
- tr_div(iTracer,k,cell1) = tr_div(iTracer,k,cell1) &
- + tr_flux(iTracer,k,iEdge)*dvEdge(iEdge)
- enddo
- enddo
- endif
- if (cell2 <= nCells) then
- do k=1,nVertLevels
- do iTracer=1,num_tracers
- tr_div(iTracer,k,cell2) = tr_div(iTracer,k,cell2) &
- - tr_flux(iTracer,k,iEdge)*dvEdge(iEdge)
- enddo
- enddo
- end if
- end do
+ end do
- ! add div(h \kappa_h </font>
<font color="red">abla\phi ) to tracer tendency
- do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k = 1,nVertLevels
- do iTracer=1,num_tracers
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- + tr_div(iTracer,k,iCell)*r
- enddo
- enddo
- enddo
- deallocate(tr_flux, tr_div)
+ deallocate(delsq_tracer)
+ end if
+
!
! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
!
</font>
</pre>