<p><b>duda</b> 2011-06-29 14:16:00 -0600 (Wed, 29 Jun 2011)</p><p>BRANCH COMMIT<br>
<br>
Update sw core w.r.t. trunk.<br>
<br>
<br>
M namelist.input.sw<br>
M src/core_sw/module_time_integration.F<br>
M src/core_sw/Registry<br>
M src/core_sw/module_mpas_core.F<br>
M src/core_sw/module_global_diagnostics.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/namelist.input.sw
===================================================================
--- branches/atmos_physics/namelist.input.sw 2011-06-29 20:02:04 UTC (rev 908)
+++ branches/atmos_physics/namelist.input.sw 2011-06-29 20:16:00 UTC (rev 909)
@@ -5,11 +5,17 @@
config_ntimesteps = 7500
config_output_interval = 500
config_stats_interval = 0
- config_visc = 0.0
+ config_h_ScaleWithMesh = .false.
+ config_h_mom_eddy_visc2 = 0.0
+ config_h_mom_eddy_visc4 = 0.0
+ config_h_tracer_eddy_diff2 = 0.0
+ config_h_tracer_eddy_diff4 = 0.0
config_thickness_adv_order = 2
config_tracer_adv_order = 2
config_positive_definite = .false.
config_monotonic = .false.
+ config_wind_stress = .false.
+ config_bottom_drag = .false.
/
&io
Modified: branches/atmos_physics/src/core_sw/Registry
===================================================================
--- branches/atmos_physics/src/core_sw/Registry 2011-06-29 20:02:04 UTC (rev 908)
+++ branches/atmos_physics/src/core_sw/Registry 2011-06-29 20:16:00 UTC (rev 909)
@@ -7,11 +7,17 @@
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 logical sw_model config_h_ScaleWithMesh false
+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_h_tracer_eddy_diff2 0.0
+namelist real sw_model config_h_tracer_eddy_diff4 0.0
namelist integer sw_model config_thickness_adv_order 2
namelist integer sw_model config_tracer_adv_order 2
namelist logical sw_model config_positive_definite false
namelist logical sw_model config_monotonic false
+namelist logical sw_model config_wind_stress false
+namelist logical sw_model config_bottom_drag false
namelist real sw_model config_apvm_upwinding 0.5
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
@@ -63,6 +69,10 @@
var persistent real zVertex ( nVertices ) 0 iro zVertex mesh - -
var persistent integer indexToVertexID ( nVertices ) 0 iro indexToVertexID mesh - -
+var persistent real meshDensity ( nCells ) 0 iro meshDensity mesh - -
+var persistent real meshScalingDel2 ( nEdges ) 0 ro meshScalingDel2 mesh - -
+var persistent real meshScalingDel4 ( nEdges ) 0 ro meshScalingDel4 mesh - -
+
var persistent integer cellsOnEdge ( TWO nEdges ) 0 iro cellsOnEdge mesh - -
var persistent integer nEdgesOnCell ( nCells ) 0 iro nEdgesOnCell mesh - -
var persistent integer nEdgesOnEdge ( nEdges ) 0 iro nEdgesOnEdge mesh - -
@@ -109,6 +119,7 @@
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 iro boundaryVertex mesh - -
var persistent integer boundaryCell ( nVertLevels nCells ) 0 iro boundaryCell mesh - -
+var persistent real u_src ( nVertLevels nEdges ) 0 iro u_src mesh - -
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 iro u state - -
@@ -142,3 +153,4 @@
var persistent real gradPVt ( nVertLevels nEdges Time ) 2 - gradPVt state - -
var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
var persistent real h_vertex ( nVertLevels nVertices Time ) 2 - h_vertex state - -
+
Modified: branches/atmos_physics/src/core_sw/module_global_diagnostics.F
===================================================================
--- branches/atmos_physics/src/core_sw/module_global_diagnostics.F 2011-06-29 20:02:04 UTC (rev 908)
+++ branches/atmos_physics/src/core_sw/module_global_diagnostics.F 2011-06-29 20:16:00 UTC (rev 909)
@@ -274,7 +274,7 @@
else
open(fileID, file='GlobalIntegrals.txt',POSITION='append')
endif
- write(fileID,'(1i, 100es24.16)') timeIndex, timeIndex*dt, globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, &
+ write(fileID,'(1i0, 100es24.16)') timeIndex, timeIndex*dt, globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, &
globalEnergy, globalCoriolisEnergyTendency, globalKineticEnergyTendency+globalPotentialEnergyTendency, &
globalKineticEnergy, globalPotentialEnergy
close(fileID)
Modified: branches/atmos_physics/src/core_sw/module_mpas_core.F
===================================================================
--- branches/atmos_physics/src/core_sw/module_mpas_core.F 2011-06-29 20:02:04 UTC (rev 908)
+++ branches/atmos_physics/src/core_sw/module_mpas_core.F 2011-06-29 20:16:00 UTC (rev 909)
@@ -55,7 +55,8 @@
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
-
+ call compute_mesh_scaling(mesh)
+
call rbfInterp_initialize(mesh)
call init_reconstruct(mesh)
call reconstruct(block % state % time_levs(1) % state, block % diag, mesh)
@@ -214,4 +215,36 @@
end subroutine mpas_core_finalize
+
+ subroutine compute_mesh_scaling(mesh)
+
+ use grid_types
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: mesh
+
+ integer :: iEdge, cell1, cell2
+ real (kind=RKIND), dimension(:), pointer :: meshDensity, meshScalingDel2, meshScalingDel4
+
+ meshDensity => mesh % meshDensity % array
+ meshScalingDel2 => mesh % meshScalingDel2 % array
+ meshScalingDel4 => mesh % meshScalingDel4 % array
+
+ !
+ ! Compute the scaling factors to be used in the del2 and del4 dissipation
+ !
+ meshScalingDel2(:) = 1.0
+ meshScalingDel4(:) = 1.0
+ if (config_h_ScaleWithMesh) then
+ do iEdge=1,mesh%nEdges
+ cell1 = mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = mesh % cellsOnEdge % array(2,iEdge)
+ meshScalingDel2(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/12.0)
+ meshScalingDel4(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/6.0)
+ end do
+ end if
+
+ end subroutine compute_mesh_scaling
+
end module mpas_core
Modified: branches/atmos_physics/src/core_sw/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_sw/module_time_integration.F 2011-06-29 20:02:04 UTC (rev 908)
+++ branches/atmos_physics/src/core_sw/module_time_integration.F 2011-06-29 20:16:00 UTC (rev 909)
@@ -120,6 +120,16 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, 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 % state % 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 % state % 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
@@ -248,16 +258,23 @@
real (kind=RKIND) :: flux, vorticity_abs, workpv, q, upstream_bias
integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+ real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ meshScalingDel2, meshScalingDel4
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, vorticity, ke, pv_edge, divergence, h_vertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: u_diffusion, visc
+ real (kind=RKIND) :: r, u_diffusion
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
- visc = config_visc
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+ real (kind=RKIND) :: ke_edge
+
h => s % h % array
u => s % u % array
v => s % v % array
@@ -295,26 +312,24 @@
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ u_src => grid % u_src % array
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
+
+
!
! Compute height tendency for each cell
!
tend_h(:,:) = 0.0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend_h(k,cell1) = tend_h(k,cell1) - flux
- end do
- end if
- if (cell2 <= nCells) then
- do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend_h(k,cell2) = tend_h(k,cell2) + flux
- end do
- end if
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell1) = tend_h(k,cell1) - flux
+ tend_h(k,cell2) = tend_h(k,cell2) + flux
+ end do
end do
do iCell=1,grid % nCellsSolve
do k=1,nVertLevels
@@ -334,15 +349,6 @@
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)
@@ -352,12 +358,13 @@
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
#ifdef NCAR_FORMULATION
@@ -386,6 +393,141 @@
end do
#endif
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="blue">abla vorticity
+ ! only valid for visc == constant
+ if (config_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
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = meshScalingDel2(iEdge) * config_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 (config_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
+
+ delsq_u(k,iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+
+ 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)
+
+ u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
+
+ end do
+ end do
+
+ deallocate(delsq_divergence)
+ deallocate(delsq_u)
+ deallocate(delsq_circulation)
+ deallocate(delsq_vorticity)
+
+ end if
+
+ ! Compute u (velocity) tendency from wind stress (u_src)
+ if(config_wind_stress) then
+ do iEdge=1,grid % nEdges
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
+ end do
+ endif
+
+ if (config_bottom_drag) then
+ do iEdge=1,grid % nEdges
+ ! bottom drag is the same as POP:
+ ! -c |u| u where c is unitless and 1.0e-3.
+ ! see POP Reference guide, section 3.4.4.
+ ke_edge = 0.5 * ( ke(1,cellsOnEdge(1,iEdge)) &
+ + ke(1,cellsOnEdge(2,iEdge)))
+
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ - 1.0e-3*u(1,iEdge) &
+ *sqrt(2.0*ke_edge)/h_edge(1,iEdge)
+ end do
+ endif
+
end subroutine compute_tend
@@ -405,7 +547,12 @@
type (mesh_type), intent(in) :: grid
integer :: iCell, iEdge, k, iTracer, cell1, cell2, i
- real (kind=RKIND) :: flux, tracer_edge
+ real (kind=RKIND) :: flux, tracer_edge, r
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
+ integer, dimension(:,:), pointer :: boundaryEdge
+ real (kind=RKIND), dimension(:,:), allocatable :: boundaryMask
+ real (kind=RKIND), dimension(:,:,:), allocatable:: delsq_tracer
+
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
real (kind=RKIND), dimension(:,:,:), pointer :: tracers, tracer_tend
@@ -422,6 +569,7 @@
tracers => s % tracers % array
cellsOnEdge => grid % cellsOnEdge % array
boundaryCell=> grid % boundaryCell % array
+ boundaryEdge=> grid % boundaryEdge % array
areaCell => grid % areaCell % array
tracer_tend => tend % tracers % array
@@ -429,6 +577,7 @@
if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
if (config_tracer_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+
tracer_tend(:,:,:) = 0.0
if (config_tracer_adv_order == 2) then
@@ -556,6 +705,108 @@
endif ! if (config_tracer_adv_order == 2 )
+ !
+ ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="blue">abla \phi)
+ !
+ if ( config_h_tracer_eddy_diff2 > 0.0 ) then
+
+ !
+ ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
+ !
+ allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
+ boundaryMask = 1.0
+ where(boundaryEdge.eq.1) boundaryMask=0.0
+
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ invAreaCell1 = 1.0/areaCell(cell1)
+ invAreaCell2 = 1.0/areaCell(cell2)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1, grid % nTracers
+ ! \kappa_2 </font>
<font color="blue">abla \phi on edge
+ tracer_turb_flux = config_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)
+ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) + flux * invAreaCell1
+ tracer_tend(iTracer,k,cell2) = tracer_tend(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 ( config_h_tracer_eddy_diff4 > 0.0 ) then
+
+ !
+ ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
+ !
+ allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
+ boundaryMask = 1.0
+ where(boundaryEdge.eq.1) boundaryMask=0.0
+
+ allocate(delsq_tracer(grid % nTracers, grid % nVertLevels, grid % nCells+1))
+
+ delsq_tracer(:,:,:) = 0.
+
+ ! first del2: div(h </font>
<font color="blue">abla \phi) at cell center
+ do iEdge=1,grid % nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1, grid % nTracers
+ 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, grid % nCells
+ r = 1.0 / grid % areaCell % array(iCell)
+ do k=1,grid % nVertLevels
+ do iTracer=1,grid % nTracers
+ delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
+ end do
+ end do
+ end do
+
+ ! second del2: div(h </font>
<font color="gray">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 / grid % areaCell % array(cell1)
+ invAreaCell2 = 1.0 / grid % areaCell % array(cell2)
+
+ do k=1,grid % nVertLevels
+ do iTracer=1,grid % nTracers
+ tracer_turb_flux = config_h_tracer_eddy_diff4 * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) / dcEdge(iEdge)
+ flux = dvEdge(iEdge) * tracer_turb_flux
+ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux * invAreaCell1 * boundaryMask(k,iEdge)
+ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux * invAreaCell2 * boundaryMask(k,iEdge)
+ end do
+ enddo
+
+ end do
+
+ deallocate(delsq_tracer)
+ deallocate(boundaryMask)
+
+ end if
+
end subroutine compute_scalar_tend
@@ -865,10 +1116,7 @@
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
- VTX_LOOP: do iVertex = 1,nVertices
- do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) > nCells) cycle VTX_LOOP
- end do
+ do iVertex = 1,nVertices
do k=1,nVertLevels
h_vertex(k,iVertex) = 0.0
do i=1,grid % vertexDegree
@@ -878,7 +1126,7 @@
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex(k,iVertex)
end do
- end do VTX_LOOP
+ end do
!
</font>
</pre>