<p><b>qchen3@fsu.edu</b> 2013-03-15 10:28:51 -0600 (Fri, 15 Mar 2013)</p><p>BRANCH COMMIT<br>
<br>
In the gm_implement branch, the computations of slopeRelative and k33 are now<br>
completed. The program compiles, but no execution of the code has been attempted.<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-15 15:49:27 UTC (rev 2621)
+++ branches/ocean_projects/gm_implement/src/core_ocean/Registry 2013-03-15 16:28:51 UTC (rev 2622)
@@ -56,6 +56,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 logical Rayleigh_damping config_Rayleigh_friction .false.
namelist real Rayleigh_damping config_Rayleigh_damping_coeff 0.0
@@ -301,6 +302,8 @@
var persistent real uBolusGMZ ( nVertLevels nEdges Time ) 2 - uBolusGMZ state - -
var persistent real uBolusGMZonal ( nVertLevels nEdges Time ) 2 o uBolusGMZonal state - -
var persistent real uBolusGMMeridional ( nVertLevels nEdges Time ) 2 o uBolusGMMeridional state - -
+var persistent real slopeRelative ( nVertLevelsP1 nEdges Time ) 2 o slopeRelative state - -
+var persistent real k33 ( nVertLevelsP1 nCells Time ) 2 o k33 state - -
var persistent real hEddyFlux ( nVertLevels nEdges Time ) 2 - hEddyFlux state - -
var persistent real h_kappa ( nVertLevels nEdges Time ) 2 - h_kappa state - -
var persistent real h_kappa_q ( nVertLevels nEdges Time ) 2 - h_kappa_q state - -
@@ -381,9 +384,18 @@
var scratch real scratch_edge_var1 ( nVertLevels nEdges ) 0 - scratch_edge_var1 scratch - -
var scratch real scratch_edge_var2 ( nVertLevels nEdges ) 0 - scratch_edge_var2 scratch - -
var scratch real scratch_edge_var3 ( nVertLevels nEdges ) 0 - scratch_edge_var3 scratch - -
-var scratch real scratch_edge_var4 ( nVertLevels nEdges ) 0 - scratch_edge_var3 scratch - -
-var scratch real scratch_edge_var5 ( nVertLevels nEdges ) 0 - scratch_edge_var3 scratch - -
+var scratch real scratch_edge_var4 ( nVertLevels nEdges ) 0 - scratch_edge_var4 scratch - -
+var scratch real scratch_edge_var5 ( nVertLevels nEdges ) 0 - scratch_edge_var5 scratch - -
var scratch real scratch_vertex_var1 ( nVertLevels nVertices ) 0 - scratch_vertex_var1 scratch - -
var scratch real scratch_vertex_var2 ( nVertLevels nVertices ) 0 - scratch_vertex_var2 scratch - -
var scratch real scratch_vertex_var3 ( nVertLevels nVertices ) 0 - scratch_vertex_var3 scratch - -
+
+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 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 grad_zMid ( nVertLevels nEdges ) 0 - grad_zMid scratch - -
+var scratch real grad_zMid_interface ( nVertLevelsP1 nEdges ) 0 - grad_zMid_interface scratch - -
+
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-15 15:49:27 UTC (rev 2621)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_gm.F 2013-03-15 16:28:51 UTC (rev 2622)
@@ -37,36 +37,61 @@
type(scratch_type), intent(inout) :: scratch
real(kind=RKIND), dimension(:,:), pointer :: rho, zMid, uBolusGM, hEddyFlux, h_edge, &
- grad_rho, grad_rho_interface, rhoz_edge, rhoz, grad_zMid, grad_zMid_interface
+ grad_rho, grad_rho_interface, rhoz_edge, rhoz, grad_zMid, grad_zMid_interface, &
+ slopeRelative, k33, rhoz_center
+ real(kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge
+ integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell
+ integer, dimension(:,:), pointer :: cellsOnEdge
+ integer :: k, iEdge, nEdges, cell1, cell2, iCell, nCells
+ real(kind=RKIND) :: h1, h2, areaEdge
- integer, dimension(:), pointer :: maxLevelEdgeTop
- integer :: k, iEdge, nEdges
- real(kind=RKIND) :: h1, h2
+ call mpas_allocate_scratch_field(scratch % grad_rho, .True.)
+ call mpas_allocate_scratch_field(scratch % grad_rho_interface, .True.)
+ call mpas_allocate_scratch_field(scratch % rhoz, .True.)
+ call mpas_allocate_scratch_field(scratch % rhoz_edge, .True.)
+ call mpas_allocate_scratch_field(scratch % rhoz_center, .True.)
+ call mpas_allocate_scratch_field(scratch % grad_zMid, .True.)
+ call mpas_allocate_scratch_field(scratch % grad_zMid_interface, .True.)
- call mpas_allocate_scratch_field(sratch % scratch_edge_var1, .True.)
- call mpas_allocate_scratch_field(sratch % scratch_edge_var2, .True.)
- call mpas_allocate_scratch_field(sratch % scratch_edge_var3, .True.)
- call mpas_allocate_scratch_field(sratch % scratch_edge_var4, .True.)
- call mpas_allocate_scratch_field(sratch % scratch_edge_var5, .True.)
- call mpas_allocate_scratch_field(sratch % scratch_cell_var1, .True.)
-
rho => s % rho % array
zMid => s % zMid % array
uBolusGM => s % uBolusGM % array
+ slopeRelative => s % slopeRelative % array
+ k33 => s % k33 % array
h_edge => s % h_edge % array
hEddyFlux => s % hEddyFlux % array
zMid => s % zMid % array
- grad_rho => scratch % scratch_edge_var1 % array
- grad_rho_interface => scratch % scratch_edge_var2 % array
- rhoz_edge => scratch % scratch_edge_var3 % array
- grad_zMid => scratch % scratch_edge_var4 % array
- grad_zMid_interface => scratch % scratch_edge_var5 % array
- rhoz => scratch % scratch_cell_var1 % array
+ grad_rho => scratch % grad_rho % array
+ grad_rho_interface => scratch % grad_rho_interface % array
+ rhoz => scratch % rhoz % array
+ rhoz_edge => scratch % rhoz_edge % array
+ rhoz_center => scratch % rhoz_center % array
+ grad_zMid => scratch % grad_zMid % array
+ grad_zMid_interface => scratch % grad_zMid_interface % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelCell => grid % maxLevelCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ areaCell => grid % areaCell % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
nEdges = grid % nEdges
+ nCells = grid % nCells
+ ! Assign a huge value to the scratch variables which may manifest itself when
+ ! there is a bug.
+ grad_rho(:,:) = huge(0D0)
+ grad_rho_interface(:,:) = huge(0D0)
+ rhoz(:,:) = huge(0D0)
+ rhoz_edge(:,:) = huge(0D0)
+ rhoz_center(:,:) = huge(0D0)
+ grad_zMid(:,:) = huge(0D0)
+ grad_zMid_interface(:,:) = huge(0D0)
+
+ slopeRelative(:,:) = huge(0D0)
+ k33(:,:) = huge(0D0)
+
! Compute density gradient (grad_rho) and gradient of zMid (grad_zMid)
! along the constant coordinate surface.
! The computed variables lives at cell edge and layer center
@@ -82,9 +107,15 @@
! Compute vertical derivative of density (rhoz) at cell center and layer interface
do iCell = 1, nCells
- do k = 2, nVertLevelCell(iCell)
+ do k = 2, maxLevelCell(iCell)
rhoz(k,iCell) = (rho(k-1,iCell) - rho(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 density.
+ ! Essentially, this enforces the boundary condition (d rho)/dz = 0 at the top and bottom.
+ rhoz(1,iCell) = 0.0
+ rhoz(maxLevelCell(iCell)+1,iCell) = 0.0
end do
! Interpolate grad_rho and grad_zMid to layer interface
@@ -98,18 +129,57 @@
grad_zMid_interface(k,iEdge) = (h2 * grad_zMid(k-1,iEdge) + h1 * grad_zMid(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 density.
+ grad_rho_interface(1,iEdge) = grad_rho(1,iEdge)
+ grad_rho_interface(maxLevelEdgeTop(iEdge)+1,iEdge) = grad_rho(maxLevelEdgeTop(iEdge),iEdge)
+ grad_zMid_interface(1,iEdge) = grad_zMid(1,iEdge)
+ grad_zMid_interface(maxLevelEdgeTop(iEdge)+1,iEdge) = grad_zMid(maxLevelEdgeTop(iEdge),iEdge)
+
end do
! Interpolate rhoz to cell edge and layer interface
do iEdge = 1, nEdges
- do k = 2, maxLevelEdgeTop(iEdge)
- cell1 = cellsOnEdge(iEdge)
- cell2 = cellsOnEdge(iEdge)
+ do k = 1, maxLevelEdgeTop(iEdge)+1
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
rhoz_edge(k,iEdge) = 0.5 * (rhoz(k,cell1) + rhoz(k,cell2))
end do
end do
+ ! Interpolate rhoz to cell center and layer center
+ do iCell = 1, nCells
+ do k = 1, maxLevelCell(iCell)
+ rhoz_center(k,iCell) = 0.5 * (rhoz(k,iCell) + rhoz(k+1,iCell))
+ end do
+ end do
+ ! Compute slopeRelative at cell edge and layer interface
+ ! QC Note: We may need to overwrite the invalid or unrealistically large values in the variable.
+ slopeRelative(:,:) = -(1.0 - config_diapycnal_diff) * grad_rho_interface(:,:) / rhoz_edge(:,:)
+ slopeRelative(:,:) = slopeRelative(:,:) - config_diapycnal_diff * grad_zMid_interface(:,:)
+
+ ! Compute k33 at cell center and layer interface
+ k33(:,:) = 0.0
+ do iEdge = 1, nEdges
+ do k = 1, maxLevelEdgeTop(iEdge)+1
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ areaEdge = dcEdge(iEdge) * dvEdge(iEdge)
+ k33(k,cell1) = k33(k,cell1) + 0.25 * areaEdge * grad_rho_interface(k,iEdge)**2
+ k33(k,cell2) = k33(k,cell2) + 0.25 * areaEdge * grad_rho_interface(k,iEdge)**2
+ end do
+ end do
+
+ ! QC Note: We may need to overwrite invalid or unrealistically large values of k33.
+ do iCell = 1,nCells
+ do k = 1, maxLevelCell(iCell)+1
+ k33(k,iCell) = k33(k,iCell) / areaCell(iCell) / rhoz(k,iCell)**2 + config_diapycnal_diff
+ end do
+ end do
+
+
call ocn_gm_compute_hEddyFlux(s, grid)
if (config_vert_coord_movement .EQ. 'isopycnal') then
@@ -127,8 +197,15 @@
end if
- call mpas_deallocate_scratch_field(scratch % grad_rho, .True.)
+ call mpas_deallocate_scratch_field(scratch % grad_rho)
+ call mpas_deallocate_scratch_field(scratch % grad_rho_interface)
+ call mpas_deallocate_scratch_field(scratch % rhoz)
+ call mpas_deallocate_scratch_field(scratch % rhoz_edge)
+ call mpas_deallocate_scratch_field(scratch % rhoz_center)
+ call mpas_deallocate_scratch_field(scratch % grad_zMid)
+ call mpas_deallocate_scratch_field(scratch % grad_zMid_interface)
+
end subroutine ocn_gm_compute_uBolus!}}}
subroutine ocn_gm_compute_slopeRelative(s, grid, scratch)
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-03-15 15:49:27 UTC (rev 2621)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-03-15 16:28:51 UTC (rev 2622)
@@ -209,7 +209,7 @@
! Compute uBolusGM and then tracer transport velocity
if (config_h_kappa .GE. epsilon(0D0)) then
- call ocn_gm_compute_uBolus(block % provis, block % mesh)
+ call ocn_gm_compute_uBolus(block % provis, block % mesh, block % scratch)
else
block % provis % uBolusGM % array(:,:) = 0.0
end if
@@ -310,7 +310,7 @@
! Compute uBolusGM and the tracer transport velocity
if (config_h_kappa .GE. epsilon(0D0)) then
- call ocn_gm_compute_uBolus(block % state % time_levs(2) % state, block % mesh)
+ call ocn_gm_compute_uBolus(block % state % time_levs(2) % state, block % mesh, block % scratch)
else
block % state % time_levs(2) % state % uBolusGM % array(:,:) = 0.0
end if
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-03-15 15:49:27 UTC (rev 2621)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_time_integration_split.F 2013-03-15 16:28:51 UTC (rev 2622)
@@ -880,7 +880,7 @@
! Compute uBolusGM
! QC note: Not sure if uBolus needs to be computed here. If not, then the next few lines can be deleted.
!if (config_h_kappa .GE. epsilon(0D0)) then
- ! call ocn_gm_compute_uBolus(block % state % time_levs(2) % state, block % mesh)
+ ! call ocn_gm_compute_uBolus(block % state % time_levs(2) % state, block % mesh, block % scratch)
!else
! block % state % time_levs(2) % state % uBolusGM % array(:,:) = 0.0
!end if
@@ -989,7 +989,7 @@
! Compute uBolusGM
if (config_h_kappa .GE. epsilon(0D0)) then
- call ocn_gm_compute_uBolus(block % state % time_levs(2) % state, block % mesh)
+ call ocn_gm_compute_uBolus(block % state % time_levs(2) % state, block % mesh, block % scratch)
else
block % state % time_levs(2) % state % uBolusGM % array(:,:) = 0.0
end if
</font>
</pre>