<p><b>mpetersen@lanl.gov</b> 2010-09-01 10:18:45 -0600 (Wed, 01 Sep 2010)</p><p>Merge dissipation branch to trunk. This adds del2 and del4<br>
dissipation operators for momentum and tracers. Dissipation branch<br>
was created by Todd Ringler, reviewed by Mark Petersen. <br>
<br>
Review process: Read code, added clarifying comments. Tested in the<br>
x3 global domain, 25 z-levels with flat bottom. Looked at diffusion<br>
of a tracer blob with del2 and del4 operators separately. Qualitative<br>
images look good, with no sign of boundary problems. Global<br>
volume-weighted tracer sum is conserved to 15 digits in all tests.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/namelist.input.ocean
===================================================================
--- trunk/mpas/namelist.input.ocean 2010-09-01 16:04:24 UTC (rev 488)
+++ trunk/mpas/namelist.input.ocean 2010-09-01 16:18:45 UTC (rev 489)
@@ -1,11 +1,10 @@
&sw_model
config_test_case = 0
config_time_integration = 'RK4'
- config_dt = 60.0
- config_ntimesteps = 1440000
- config_output_interval = 14400
- config_stats_interval = 1440
- config_visc = 1.0e5
+ config_dt = 90.0
+ config_ntimesteps = 1920000
+ config_output_interval = 19200
+ config_stats_interval = 1920
/
&io
@@ -17,29 +16,32 @@
&restart
config_restart_interval = 115200
config_do_restart = .false.
- config_restart_time = 1036800.0
+ config_restart_time = 31104000
/
&grid
- config_vert_grid_type = 'zlevel'
- config_rho0 = 1028
+ config_vert_grid_type = 'isopycnal'
+ config_rho0 = 1015
/
&hmix
- config_hor_diffusion = 1.0e4
+ config_h_mom_eddy_visc2 = 0.0
+ config_h_mom_eddy_visc4 = 5.0e8
+ config_h_tracer_eddy_diff2 = 10.0
+ config_h_tracer_eddy_diff4 = 0.0
/
&vmix
- config_vert_visc_type = 'tanh'
- config_vert_diff_type = 'tanh'
+ config_vert_visc_type = 'const'
+ config_vert_diff_type = 'const'
config_vmixTanhViscMax = 2.5e-1
config_vmixTanhViscMin = 1.0e-4
config_vmixTanhDiffMax = 2.5e-2
config_vmixTanhDiffMin = 1.0e-5
config_vmixTanhZMid = -100
config_vmixTanhZWidth = 100
- config_vert_viscosity = 2.5e-4
- config_vert_diffusion = 2.5e-5
+ config_vert_viscosity = 1.0e-4
+ config_vert_diffusion = 1.0e-4
/
&advection
- config_hor_tracer_adv = 'upwind'
- config_vert_tracer_adv = 'upwind'
+ config_hor_tracer_adv = 'centered'
+ config_vert_tracer_adv = 'centered'
/
Modified: trunk/mpas/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2010-09-01 16:04:24 UTC (rev 488)
+++ trunk/mpas/src/core_ocean/Registry 2010-09-01 16:18:45 UTC (rev 489)
@@ -7,7 +7,6 @@
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 character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
@@ -16,7 +15,10 @@
namelist real restart config_restart_time 172800.0
namelist character grid config_vert_grid_type isopycnal
namelist real grid config_rho0 1028
-namelist real hmix config_hor_diffusion 2000.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_diff2 0.0
+namelist real hmix config_h_tracer_eddy_diff4 0.0
namelist character vmix config_vert_visc_type const
namelist character vmix config_vert_diff_type const
namelist real vmix config_vert_viscosity 2.5e-4
@@ -30,7 +32,6 @@
namelist character advection config_hor_tracer_adv 'centered'
namelist character advection config_vert_tracer_adv 'centered'
-
#
# dim type name_in_file name_in_code
#
@@ -128,6 +129,7 @@
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 - -
+var real h_vertex ( nVertLevels nVertices Time ) o h_vertex - -
var real ke ( nVertLevels nCells Time ) o ke - -
var real ke_edge ( nVertLevels nEdges Time ) o ke_edge - -
var real pv_vertex ( nVertLevels nVertices Time ) o pv_vertex - -
Modified: trunk/mpas/src/core_ocean/module_test_cases.F
===================================================================
--- trunk/mpas/src/core_ocean/module_test_cases.F 2010-09-01 16:04:24 UTC (rev 488)
+++ trunk/mpas/src/core_ocean/module_test_cases.F 2010-09-01 16:18:45 UTC (rev 489)
@@ -28,10 +28,10 @@
! mrp 100507: for diagnostic output
integer :: iTracer
real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
- hZLevel, zMidZLevel, zTopZLevel
+ hZLevel, zMidZLevel, zTopZLevel, latCell,LonCell
real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- real (kind=RKIND) :: delta_rho
+ real (kind=RKIND) :: delta_rho, pi, latCenter, lonCenter, dist
integer :: nCells, nEdges, nVertices, nVertLevels
! mrp 100507 end: for diagnostic output
@@ -106,6 +106,8 @@
u_src => block_ptr % mesh % u_src % array
xCell => block_ptr % mesh % xCell % array
yCell => block_ptr % mesh % yCell % array
+ latCell => block_ptr % mesh % latCell % array
+ lonCell => block_ptr % mesh % lonCell % array
hZLevel => block_ptr % mesh % hZLevel % array
zMidZLevel => block_ptr % mesh % zMidZLevel % array
@@ -116,6 +118,13 @@
nVertices = block_ptr % mesh % nVertices
nVertLevels = block_ptr % mesh % nVertLevels
+ pi=3.1415
+ ! Tracer blob in Central Pacific, away from boundaries:
+ !latCenter=pi/16; lonCenter=9./8.*pi
+
+ ! Tracer blob in Central Pacific, near boundaries:
+ latCenter=pi*2./16; lonCenter=13./16.*pi
+
if (config_vert_grid_type.eq.'zlevel') then
! These should eventually be in an input file. For now
! I just read them in from h(:,1).
@@ -138,19 +147,27 @@
! Set tracers, if not done in grid.nc file
!tracers = 0.0
do iCell = 1,nCells
+ dist = sqrt( (latCell(iCell)-latCenter)**2 + (lonCell(iCell)-lonCenter)**2)
do iLevel = 1,nVertLevels
! for 20 layer test
! tracers(index_temperature,iLevel,iCell) = 5.0 ! temperature
! tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
! for x3, 25 layer test
- tracers(index_temperature,iLevel,iCell) = 10.0 ! temperature
- tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
+ !tracers(index_temperature,iLevel,iCell) = 10.0 ! temperature
+ !tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
- tracers(index_tracer1,iLevel,iCell) = 1.0
- tracers(index_tracer2,iLevel,iCell) = &
- (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
+ ! tracers(index_tracer1,iLevel,iCell) = 1.0
+ ! tracers(index_tracer2,iLevel,iCell) = &
+ ! (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
+ ! Tracer blob
+ !if (dist.lt.pi/16) then
+ ! tracers(index_tracer1,iLevel,iCell) = 1.0
+ !!else
+ ! tracers(index_tracer1,iLevel,iCell) = 0.0
+ !endif
+
rho(iLevel,iCell) = 1000.0*( 1.0 &
- 2.5e-4*tracers(index_temperature,iLevel,iCell) &
+ 7.6e-4*tracers(index_salinity,iLevel,iCell))
Modified: trunk/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_ocean/module_time_integration.F 2010-09-01 16:04:24 UTC (rev 488)
+++ trunk/mpas/src/core_ocean/module_time_integration.F 2010-09-01 16:18:45 UTC (rev 489)
@@ -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)
!
@@ -370,12 +389,16 @@
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 +447,138 @@
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="blue">abla \xi)
+ ! velocity tendency: del2 dissipation, </font>
<font color="black">u_2 </font>
<font color="blue">abla^2 u
+ ! computed as </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity )
+ ! strictly only valid for h_mom_eddy_visc2 == constant
!
- ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
- ! only valid for visc == constant
- do iEdge=1,grid % nEdgesSolve
+ 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
+
+ ! Here -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ ! is - </font>
<font color="blue">abla vorticity pointing from vertex 2 to vertex 1, or equivalently
+ ! + k \times </font>
<font color="blue">abla vorticity pointing from cell1 to cell2.
+
+ 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
+
+ !
+ ! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="blue">abla^4 u
+ ! computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity
+ ! applied recursively.
+ ! strictly only valid for h_mom_eddy_visc4 == constant
+ !
+ 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
+
+ ! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="blue">abla vorticity
+ 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
+
+ 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
+
+ ! vorticity using </font>
<font color="blue">abla^2 u
+ delsq_circulation(:,:) = 0.0
+ do iEdge=1,nEdges
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
+ - dcEdge(iEdge) * delsq_u(k,iEdge)
+ delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
+ + 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
+
+ ! Divergence using </font>
<font color="blue">abla^2 u
+ 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
+
+ ! Compute - \kappa </font>
<font color="blue">abla^4 u
+ ! as </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="gray">abla^2 u) )
+ 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 = ( 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 +587,6 @@
end do
tend_u(k,iEdge) = tend_u(k,iEdge) &
+ q &
- + u_diffusion &
- ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
end do
@@ -471,13 +608,6 @@
- 1.0e-3*u(nVertLevels,iEdge) &
*sqrt(2.0*ke_edge(nVertLevels,iEdge))
- ! mrp 100603 The following method is more straight forward,
- ! that the above method of computing ke_edge, but I have
- ! not verified that v is working correctly yet.
- !tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- ! - 1.0e-3*u(nVertLevels,iEdge) &
- ! *sqrt(u(nVertLevels,iEdge)**2 + v(nVertLevels,iEdge)**2)
-
! old bottom drag, just linear friction
! du/dt = u/tau where tau=100 days.
!tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
@@ -546,6 +676,7 @@
integer :: iCell, iEdge, k, iTracer, cell1, cell2, upwindCell,&
nEdges, nCells, nVertLevels
+ real (kind=RKIND) :: h_tracer_eddy_diff2, h_tracer_eddy_diff4, invAreaCell1, invAreaCell2, tracer_turb_flux
real (kind=RKIND) :: flux, tracer_edge, r
real (kind=RKIND) :: dist
real (kind=RKIND), dimension(:), pointer :: &
@@ -554,14 +685,15 @@
u,h,wTop, h_edge, zMid, zTop
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
+ integer, dimension(:,:), pointer :: boundaryEdge
type (dm_info) :: dminfo
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
integer, dimension(:,:), pointer :: cellsOnEdge
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:: fluxVertTop, tracerTop, boundaryMask
+ real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer
real (kind=RKIND), dimension(:), allocatable:: vertDiffTop
@@ -580,11 +712,16 @@
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
zTopZLevel => grid % zTopZLevel % array
+ boundaryEdge => grid % boundaryEdge % array
nEdges = grid % nEdges
nCells = grid % nCells
nVertLevels = grid % nVertLevels
+
+ h_tracer_eddy_diff2 = config_h_tracer_eddy_diff2
+ h_tracer_eddy_diff4 = config_h_tracer_eddy_diff4
+
!
! tracer tendency: horizontal advection term -div( h \phi u)
!
@@ -681,61 +818,118 @@
deallocate(tracerTop)
!
- ! tracer tendency: horizontal tracer diffusion
- ! div(h \kappa_h </font>
<font color="blue">abla\phi )
+ ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="red">abla \phi)
!
- ! 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_diff2 > 0.0 ) then
+
+ !
+ ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
+ !
+ allocate(boundaryMask(nVertLevels, nEdges+1))
+ boundaryMask = 1.0
+ where(boundaryEdge.eq.1) boundaryMask=0.0
+
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ invAreaCell1 = 1.0/areaCell(cell1)
+ invAreaCell2 = 1.0/areaCell(cell2)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1,num_tracers
+ ! \kappa_2 </font>
<font color="blue">abla \phi on edge
+ tracer_turb_flux = h_tracer_eddy_diff2 &
+ *( tracers(iTracer,k,cell2) &
+ - tracers(iTracer,k,cell1))/dcEdge(iEdge)
+
+ ! div(h \kappa_2 </font>
<font color="blue">abla \phi) at cell center
+ flux = dvEdge (iEdge) * h_edge(k,iEdge) &
+ * tracer_turb_flux * boundaryMask(k, iEdge)
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) + flux * invAreaCell1
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) - flux * invAreaCell2
+ end do
+ end do
+
+ end do
+
+ deallocate(boundaryMask)
+
+ end if
+
+ !
+ ! tracer tendency: del4 horizontal tracer diffusion, &
+ ! div(h \kappa_4 </font>
<font color="black">abla [div(h </font>
<font color="blue">abla \phi)])
+ !
+ if ( h_tracer_eddy_diff4 > 0.0 ) then
+
+ !
+ ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
+ !
+ allocate(boundaryMask(nVertLevels, nEdges+1))
+ boundaryMask = 1.0
+ where(boundaryEdge.eq.1) boundaryMask=0.0
+
+ allocate(delsq_tracer(num_tracers,nVertLevels, nCells+1))
+
+ delsq_tracer(:,:,:) = 0.
+
+ ! first del2: div(h </font>
<font color="red">abla \phi) at cell center
+ 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) * boundaryMask(k,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) * boundaryMask(k,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
+
+ ! second del2: div(h </font>
<font color="red">abla [delsq_tracer]) at cell center
+ do iEdge=1,grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1,num_tracers
+ tracer_turb_flux = h_tracer_eddy_diff4 &
+ *( 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 * invAreaCell1 * boundaryMask(k,iEdge)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) &
+ + flux * invAreaCell2 * boundaryMask(k,iEdge)
+
+ 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>