<p><b>qchen3@fsu.edu</b> 2013-04-02 14:14:47 -0600 (Tue, 02 Apr 2013)</p><p>BRANCH COMMIT<br>
<br>
The horizontal components of the Redi diffusion is implemented. Code compiles, but not tested.<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-04-02 19:43:18 UTC (rev 2709)
+++ branches/ocean_projects/gm_implement/src/core_ocean/Registry 2013-04-02 20:14:47 UTC (rev 2710)
@@ -395,10 +395,15 @@
var scratch real grad_rho ( nVertLevels nEdges ) 0 - grad_rho scratch - -
var scratch real grad_rho_interface ( nVertLevelsP1 nEdges ) 0 - grad_rho_interface scratch - -
+var scratch real grad_hTracer ( nVertLevels nEdges ) 0 - grad_hTracer scratch - -
+var scratch real grad_hTracer_interface ( nVertLevelsP1 nEdges ) 0 - grad_hTracer_interface scratch - -
+var scratch real grad_hTracer_sloped ( nVertLevelsP1 nCells ) 0 - grad_hTracer_sloped scratch - -
var scratch real grad_rho_constZ ( nVertLevelsP1 nEdges ) 0 - grad_rho_constZ scratch - -
var scratch real rhoz ( nVertLevelsP1 nCells ) 0 - rhoz scratch - -
var scratch real rhoz_edge ( nVertLevelsP1 nEdges ) 0 - rhoz_edge scratch - -
var scratch real rhoz_center ( nVertLevels nCells ) 0 - rhoz_center scratch - -
+var scratch real dhTracerdZ ( nVertLevelsP1 nCells ) 0 - dhTracerdZ scratch - -
+var scratch real dhTracerdZ_edge ( nVertLevelsP1 nEdges ) 0 - dhTracerdZ_edge scratch - -
var scratch real grad_zMid ( nVertLevels nEdges ) 0 - grad_zMid scratch - -
var scratch real grad_zMid_interface ( nVertLevelsP1 nEdges ) 0 - grad_zMid_interface scratch - -
var scratch real gmStreamFunc ( nVertLevelsP1 nEdges ) 0 - gmStreamFunc scratch - -
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-04-02 19:43:18 UTC (rev 2709)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tendency.F 2013-04-02 20:14:47 UTC (rev 2710)
@@ -253,17 +253,18 @@
!> This routine computes tracer tendencies for the ocean
!
!-----------------------------------------------------------------------
- subroutine ocn_tend_tracer(tend, s, d, grid, dt)!{{{
+ subroutine ocn_tend_tracer(tend, s, d, grid, scratch, dt)!{{{
implicit none
type (tend_type), intent(inout) :: tend !< Input/Output: Tendency structure
type (state_type), intent(in) :: s !< Input: State information
type (diagnostics_type), intent(in) :: d !< Input: Diagnostic information
type (mesh_type), intent(in) :: grid !< Input: Grid information
+ type (scratch_type), intent(inout) :: scratch !< Input: Grid information
real (kind=RKIND), intent(in) :: dt !< Input: Time step
real (kind=RKIND), dimension(:,:), pointer :: &
- uTransport, h,wTop, h_edge, vertDiffTopOfCell, tend_h, uh
+ uTransport, h,wTop, h_edge, vertDiffTopOfCell, tend_h, uh, zMid, slopeRelative
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
@@ -277,6 +278,8 @@
tracers => s % tracers % array
h_edge => s % h_edge % array
vertDiffTopOfCell => d % vertDiffTopOfCell % array
+ zMid => s % zMid % array
+ slopeRelative => s % slopeRelative % array
tend_tr => tend % tracers % array
tend_h => tend % h % array
@@ -309,7 +312,7 @@
! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="gray">abla \phi)
!
call mpas_timer_start("hmix", .false., tracerHmixTimer)
- call ocn_tracer_hmix_tend(grid, h_edge, tracers, tend_tr, err)
+ call ocn_tracer_hmix_tend(grid, scratch, h, h_edge, zMid, tracers, slopeRelative, tend_tr, err)
call mpas_timer_stop("hmix", tracerHmixTimer)
!
Modified: branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-04-02 19:43:18 UTC (rev 2709)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-04-02 20:14:47 UTC (rev 2710)
@@ -163,7 +163,7 @@
call ocn_filter_btr_mode_tend_u(block % tend, block % provis, block % mesh)
endif
- call ocn_tend_tracer(block % tend, block % provis, block % diagnostics, block % mesh, dt)
+ call ocn_tend_tracer(block % tend, block % provis, block % diagnostics, block % mesh, block % scratch, dt)
block => block % next
end do
call mpas_timer_stop("RK4-tendency computations")
Modified: branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_time_integration_split.F 2013-04-02 19:43:18 UTC (rev 2709)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_time_integration_split.F 2013-04-02 20:14:47 UTC (rev 2710)
@@ -793,7 +793,7 @@
block => domain % blocklist
do while (associated(block))
- call ocn_tend_tracer(block % tend, block % state % time_levs(2) % state, block % diagnostics, block % mesh, dt)
+ call ocn_tend_tracer(block % tend, block % state % time_levs(2) % state, block % diagnostics, block % mesh, block % scratch, dt)
block => block % next
end do
Modified: branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tracer_hmix.F
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tracer_hmix.F 2013-04-02 19:43:18 UTC (rev 2709)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tracer_hmix.F 2013-04-02 20:14:47 UTC (rev 2710)
@@ -73,7 +73,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_hmix_tend(grid, h_edge, tracers, tend, err)!{{{
+ subroutine ocn_tracer_hmix_tend(grid, scratch, h, h_edge, zMid, tracers, slopeRelative, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -84,9 +84,21 @@
real (kind=RKIND), dimension(:,:), intent(in) :: &
h_edge !< Input: thickness at edge
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness at edge
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ zMid !< Input: thickness at edge
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ slopeRelative !< Input: thickness at edge
+
type (mesh_type), intent(in) :: &
grid !< Input: grid information
+ type (scratch_type), intent(inout) :: &
+ scratch !< Input: grid information
+
real (kind=RKIND), dimension(:,:,:), intent(in) :: &
tracers !< Input: tracer quantities
@@ -126,7 +138,7 @@
if(.not.tracerHmixOn) return
call mpas_timer_start("del2", .false., del2Timer)
- call ocn_tracer_hmix_del2_tend(grid, h_edge, tracers, tend, err1)
+ call ocn_tracer_hmix_del2_tend(grid, scratch, h, h_edge, zMid, tracers, slopeRelative, tend, err1)
call mpas_timer_stop("del2", del2Timer)
call mpas_timer_start("del4", .false., del4Timer)
call ocn_tracer_hmix_del4_tend(grid, h_edge, tracers, tend, err2)
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-04-02 19:43:18 UTC (rev 2709)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2013-04-02 20:14:47 UTC (rev 2710)
@@ -67,7 +67,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_hmix_del2_tend(grid, h_edge, tracers, slopeRelative, tend, err)!{{{
+ subroutine ocn_tracer_hmix_del2_tend(grid, scratch, h, h_edge, zMid, tracers, slopeRelative, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -76,14 +76,23 @@
!-----------------------------------------------------------------
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
+ type (scratch_type), intent(inout) :: &
+ scratch !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness at cell center
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h_edge !< Input: thickness at edge
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ zMid !< Input: height of cell center and layer center
+
real (kind=RKIND), dimension(:,:,:), intent(in) :: &
tracers !< Input: tracer quantities
@@ -115,15 +124,24 @@
integer, dimension(:,:), allocatable :: boundaryMask
- integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnCell
+ integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnCell, maxLevelCell
integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask, edgesOnCell, edgeSignOnCell
real (kind=RKIND) :: invAreaCell1, invAreaCell2
- real (kind=RKIND) :: tracer_turb_flux, flux, r_tmp
+ real (kind=RKIND) :: tracer_turb_flux, flux, r_tmp, h1, h2
real (kind=RKIND), dimension(:), pointer :: areaCell, dvEdge, dcEdge
real (kind=RKIND), dimension(:), pointer :: meshScalingDel2
+ real (kind=RKIND), dimension(:,:), pointer :: grad_hTracer, grad_hTracer_interface, grad_hTracer_sloped, dhTracerdZ, &
+ dhTracerdZ_edge
+
+ call mpas_allocate_scratch_field(scratch % grad_hTracer, .True.)
+ call mpas_allocate_scratch_field(scratch % grad_hTracer_interface, .True.)
+ call mpas_allocate_scratch_field(scratch % grad_hTracer_sloped, .True.)
+ call mpas_allocate_scratch_field(scratch % dhTracerdZ, .True.)
+ call mpas_allocate_scratch_field(scratch % dhTracerdZ_edge, .True.)
+
!-----------------------------------------------------------------
!
! call relevant routines for computing tendencies
@@ -142,6 +160,7 @@
num_tracers = size(tracers, dim=1)
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelCell => grid % maxLevelCell % array
cellsOnEdge => grid % cellsOnEdge % array
edgeMask => grid % edgeMask % array
areaCell => grid % areaCell % array
@@ -153,6 +172,12 @@
edgesOnCell => grid % edgesOnCell % array
edgeSignOnCell => grid % edgeSignOnCell % array
+ grad_hTracer => scratch % grad_hTracer % array
+ grad_hTracer_interface => scratch % grad_hTracer_interface % array
+ grad_hTracer_sloped => scratch % grad_hTracer_sloped % array
+ dhTracerdZ => scratch % dhTracerdZ % array
+ dhTracerdZ_edge => scratch % dhTracerdZ_edge % array
+
!
! compute a boundary mask to enforce insulating boundary conditions in the horizontal
!
@@ -183,32 +208,65 @@
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
+ !
+ ! NEED TO COUPLE TRACER WITH THICKNESS!!!!
+ !
+
+
+ !
+ ! 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 iTracer = 1, num_tracers
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))
+ dhTracerdZ(k,iCell) = (h(k-1,iCell)*tracers(iTracer,k-1,iCell) - h(k,iCell)*tracers(iTracer,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
+ ! Approximation of dhTracerdZ 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
+ ! Essentially, this enforces the boundary condition (d tracer)/dz = 0 at the top and bottom.
+ dhTracerdZ(1,iCell) = 0.0
+ dhTracerdZ(maxLevelCell(iCell)+1,iCell) = 0.0
end do
- ! Interpolate dTracersdZ to cell edge and layer interface
+ ! Compute tracer gradient (grad_hTracer) along the constant coordinate surface.
+ ! The computed variables lives at cell edge and layer center
do iEdge = 1, nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+ grad_hTracer(k,iEdge) = (h(k,cell2)*tracers(iTracer,k,cell2) - h(k,cell1)*tracers(iTracer,k,cell1)) / dcEdge(iEdge)
+ end do
+ end do
+
+ ! Interpolate dhTracerdZ 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))
+ dhTracerdZ_edge(k,iEdge) = 0.5 * (dhTracerdZ(k,cell1) + dhTracerdZ(k,cell2))
end do
end do
+ ! Interpolate grad_hTracer to layer interface
+ do iEdge = 1, nEdges
+ do k = 2, maxLevelEdgeTop(iEdge)
+ h1 = h_edge(k-1,iEdge)
+ h2 = h_edge(k,iEdge)
+
+ ! Using second-order interpolation below
+ grad_hTracer_interface(k,iEdge) = (h2 * grad_hTracer(k-1,iEdge) + h1 * grad_hTracer(k,iEdge)) / (h1 + h2)
+ end do
+
+ ! Approximation of values 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 concentration.
+ grad_hTracer_interface(1,iEdge) = grad_hTracer(1,iEdge)
+ grad_hTracer_interface(maxLevelEdgeTop(iEdge)+1,iEdge) = grad_hTracer(maxLevelEdgeTop(iEdge),iEdge)
+ end do
+
! Compute </font>
<font color="gray">abla\cdot(slopeRelative d\phi/dz)
do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
@@ -216,18 +274,41 @@
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
+ flux = 0.5*dvEdge(iEdge)*(slopeRelative(k,iEdge)*dhTracerdZ_edge(k,iEdge) + slopeRelative(k+1,iEdge)*dhTracerdZ_edge(k+1,iEdge))
+ tend(iTracer,k,cell1) = tend(iTracer,k,cell1) + config_Redi_kappa * flux * invAreaCell1
+ tend(iTracer,k,cell2) = tend(iTracer,k,cell2) - config_Redi_kappa * flux * invAreaCell2
end do
-
end do
+ ! Compute d(slopeRelative\cdot</font>
<font color="blue">abla\phi)/dz
+ do iEdge = 1, nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ invAreaCell1 = 1./areaCell(cell1)
+ invAreaCell2 = 1./areaCell(cell2)
+
+ do k = 1, maxLevelEdgeTop(iEdge) + 1
+ r_tmp = 0.5 * slopeRelative(k,iEdge) * grad_hTracer_interface(k,iEdge) * dcEdge(iEdge) * dvEdge(iEdge)
+ grad_hTracer_sloped(k,cell1) = grad_hTracer_sloped(k,cell1) + r_tmp * invAreaCell1
+ grad_hTracer_sloped(k,cell2) = grad_hTracer_sloped(k,cell2) + r_tmp * invAreaCell2
+ end do
+ end do
+
+ do iCell = 1, nCells
+ do k = 1, maxLevelCell(iCell)
+ tend(iTracer,k,iCell) = tend(iTracer,k,iCell) + config_Redi_kappa * (grad_hTracer_sloped(k,iCell) - grad_hTracer_sloped(k+1,iCell)) / h(k,iCell)
+ end do
+ end do
+
+ end do
+
+ call mpas_deallocate_scratch_field(scratch % grad_hTracer)
+ call mpas_deallocate_scratch_field(scratch % grad_hTracer_interface)
+ call mpas_deallocate_scratch_field(scratch % grad_hTracer_sloped)
+ call mpas_deallocate_scratch_field(scratch % dhTracerdZ)
+ call mpas_deallocate_scratch_field(scratch % dhTracerdZ_edge)
+
!--------------------------------------------------------------------
end subroutine ocn_tracer_hmix_del2_tend!}}}
</font>
</pre>