<p><b>qchen3@fsu.edu</b> 2013-03-29 11:24:47 -0600 (Fri, 29 Mar 2013)</p><p>BRANCH COMMIT<br>
<br>
For gm_implement. Not compiled.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/gm_implement/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/Registry        2013-03-29 16:32:28 UTC (rev 2681)
+++ branches/ocean_projects/gm_implement/src/core_ocean/Registry        2013-03-29 17:24:47 UTC (rev 2682)
@@ -57,6 +57,7 @@
namelist real standard_GM config_h_kappa 0.0
namelist real standard_GM config_h_kappa_q 0.0
namelist real standard_GM config_diapycnal_diff 0.0000001
+namelist real standard_GM config_Redi_kappa 0.0
namelist real standard_GM config_gravWaveSpeed_trunc 0.3
namelist logical Rayleigh_damping config_Rayleigh_friction .false.
Modified: branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_gm.F
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_gm.F        2013-03-29 16:32:28 UTC (rev 2681)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_gm.F        2013-03-29 17:24:47 UTC (rev 2682)
@@ -247,7 +247,7 @@
end do
-
+ ! Deallocate scratch variables
call mpas_deallocate_scratch_field(scratch % grad_rho)
call mpas_deallocate_scratch_field(scratch % grad_rho_interface)
call mpas_deallocate_scratch_field(scratch % rhoz)
Modified: branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tendency.F        2013-03-29 16:32:28 UTC (rev 2681)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tendency.F        2013-03-29 17:24:47 UTC (rev 2682)
@@ -289,9 +289,7 @@
if(config_disable_tr_all_tend) return
allocate(uh(grid % nVertLevels, grid % nEdges+1))
- !
- ! QC Comment (3/15/12): need to make sure that uTransport is the right
- ! transport velocity for the tracer.
+
do iEdge = 1, grid % nEdges
do k = 1, grid % nVertLevels
uh(k, iEdge) = uTransport(k, iEdge) * h_edge(k, iEdge)
Modified: branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tracer_hmix_del2.F
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tracer_hmix_del2.F        2013-03-29 16:32:28 UTC (rev 2681)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tracer_hmix_del2.F        2013-03-29 17:24:47 UTC (rev 2682)
@@ -67,7 +67,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_hmix_del2_tend(grid, h_edge, tracers, tend, err)!{{{
+ subroutine ocn_tracer_hmix_del2_tend(grid, h_edge, tracers, slopeRelative, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -78,6 +78,9 @@
real (kind=RKIND), dimension(:,:), intent(in) :: &
h_edge !< Input: thickness at edge
+ real (kind=RKIND), dimension(:,:), slopeRelative(in) :: &
+ slopeRelative !< Input: the relative slope at cell edge and layer interface
+
type (mesh_type), intent(in) :: &
grid !< Input: grid information
@@ -154,6 +157,8 @@
! compute a boundary mask to enforce insulating boundary conditions in the horizontal
!
do iCell = 1, nCells
+
+ ! Compute the horizontal diffusion
invAreaCell1 = 1.0 / areaCell(iCell)
do i = 1, nEdgesOncell(iCell)
iEdge = edgesOnCell(i, iCell)
@@ -175,8 +180,54 @@
end do
end do
+
end do
+ !
+ ! Compute the extra terms arising due to mismatch between the constant coordiate surfaces and the
+ ! isopycnal surfaces.
+ !
+ ! Compute vertical derivative of tracers at cell center and layer interface
+ do iCell = 1, nCells
+ do k = 2, maxLevelCell(iCell)
+ dTracersdZ(:, k,iCell) = (tracers(:,k-1,iCell) - tracers(:,k,iCell)) / (zMid(k-1,iCell) - zMid(k,iCell))
+ end do
+
+ ! Approximation of rhoz on the top and bottom interfaces through the idea of having
+ ! ghost cells above the top and below the bottom layers of the same depths and tracer density.
+ ! Essentially, this enforces the boundary condition (d tracers)/dz = 0 at the top and bottom.
+ dTracersdZ(:,1,iCell) = 0.0
+ dTracersdZ(:,maxLevelCell(iCell)+1,iCell) = 0.0
+ end do
+
+ ! Interpolate dTracersdZ to cell edge and layer interface
+ do iEdge = 1, nEdges
+ do k = 1, maxLevelEdgeTop(iEdge)+1
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ dTracersdZ_edge(:,k,iEdge) = 0.5 * (dTracersdZ(:,k,cell1) + dTracersdZ(:,k,cell2))
+ end do
+ end do
+
+ ! Compute </font>
<font color="gray">abla\cdot(slopeRelative d\phi/dz)
+ do iEdge = 1, nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ invAreaCell1 = 1./areaCell(cell1)
+ invAreaCell2 = 1./areaCell(cell2)
+
+ r_temp = 0.5 * config_Redi_kappa * dvEdge(iEdge)
+
+ do k = 1, maxLevelEdgeTop(iEdge)
+ do iTracer=1,num_tracers
+ flux = r_temp*(slopeRelative(k,iEdge)*dTracersdZ_edge(iTracer,k,iEdge) + slopeRelative(k+1,iEdge)*dTracersdZ_edge(iTracer,k+1,iEdge))
+ tend(iTracer,k,cell1) = tend(iTracer,k,cell1) + flux * invAreaCell1
+ tend(iTracer,k,cell2) = tend(iTracer,k,cell2) - flux * invAreaCell2
+ end do
+ end do
+
+ end do
+
!--------------------------------------------------------------------
end subroutine ocn_tracer_hmix_del2_tend!}}}
@@ -213,7 +264,7 @@
if ( config_tracer_del2 > 0.0 ) then
del2On = .true.
- eddyDiff2 = config_tracer_del2
+ eddyDiff2 = config_tracer_del2 + config_Redi_kappa
endif
if(.not.config_use_tracer_del2) del2on = .false.
</font>
</pre>