<p><b>mpetersen@lanl.gov</b> 2012-02-14 15:01:52 -0700 (Tue, 14 Feb 2012)</p><p>ztilde branch: Add thickness-weighted divergence, div_hu and div_hu_btr, as new variables, add subroutine for their computation.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/ztilde/namelist.input.ocean
===================================================================
--- branches/ocean_projects/ztilde/namelist.input.ocean 2012-02-14 21:41:33 UTC (rev 1508)
+++ branches/ocean_projects/ztilde/namelist.input.ocean 2012-02-14 22:01:52 UTC (rev 1509)
@@ -22,6 +22,9 @@
config_vert_grid_type = 'isopycnal'
config_pressure_type = 'pressure'
config_rho0 = 1014.65
+ config_div_low_freq_time = 86400.0
+ config_h_restore_time = 86400.0
+ config_h_high_freq_diff = 1.0e-5
/
&split_explicit_ts
config_n_ts_iter = 2
Modified: branches/ocean_projects/ztilde/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/ztilde/src/core_ocean/Registry 2012-02-14 21:41:33 UTC (rev 1508)
+++ branches/ocean_projects/ztilde/src/core_ocean/Registry 2012-02-14 22:01:52 UTC (rev 1509)
@@ -174,6 +174,7 @@
var persistent integer maxLevelVertexTop ( nVertices ) 0 - maxLevelVertexTop mesh - -
var persistent integer maxLevelVertexBot ( nVertices ) 0 - maxLevelVertexBot mesh - -
var persistent real referenceBottomDepth ( nVertLevels ) 0 iro referenceBottomDepth mesh - -
+var persistent real referenceBottomDepthTopOfCell ( nVertLevelsP1 ) 0 - referenceBottomDepthTopOfCell mesh - -
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
% Boundary conditions: read from input, saved in restart and written to output
@@ -255,6 +256,8 @@
var persistent real circulation ( nVertLevels nVertices Time ) 2 - circulation state - -
var persistent real gradPVt ( nVertLevels nEdges Time ) 2 - gradPVt state - -
var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
+var persistent real div_hu ( nVertLevels nCells Time ) 2 - div_hu state - -
+var persistent real div_hu_btr ( nCells Time ) 2 - div_hu_btr state - -
% Globally reduced diagnostic variables: only written to output
var persistent real areaCellGlobal ( Time ) 2 o areaCellGlobal state - -
Modified: branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_mpas_core.F 2012-02-14 21:41:33 UTC (rev 1508)
+++ branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_mpas_core.F 2012-02-14 22:01:52 UTC (rev 1509)
@@ -233,7 +233,8 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(1) % state, mesh)
call mpas_timer_stop("diagnostic solve", initDiagSolveTimer)
- call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(1) % state, mesh, block % tend)
+ call ocn_div_hu(block % state % time_levs(1) % state,block % state % time_levs(1) % state, block % mesh)
+ call ocn_wtop (block % state % time_levs(1) % state,block % state % time_levs(1) % state, mesh, block % tend)
call ocn_compute_mesh_scaling(mesh)
@@ -515,7 +516,8 @@
integer :: iTracer, cell, cell1, cell2
real (kind=RKIND) :: uhSum, hSum, hEdge1
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
+ real (kind=RKIND), dimension(:), pointer :: &
+ referenceBottomDepth, referenceBottomDepthTopOfCell
real (kind=RKIND), dimension(:,:), pointer :: h
integer :: nVertLevels
@@ -525,6 +527,7 @@
h => block % state % time_levs(1) % state % h % array
referenceBottomDepth => block % mesh % referenceBottomDepth % array
+ referenceBottomDepthTopOfCell => block % mesh % referenceBottomDepthTopOfCell % array
nVertLevels = block % mesh % nVertLevels
! mrp 120208 right now hZLevel is in the grid.nc file.
@@ -537,6 +540,12 @@
referenceBottomDepth(k) = referenceBottomDepth(k-1) + block % mesh % hZLevel % array(k)
end do
+ ! TopOfCell needed where zero depth for the very top may be referenced.
+ referenceBottomDepthTopOfCell(1) = 0.0
+ do k = 1,nVertLevels
+ referenceBottomDepthTopOfCell(k+1) = referenceBottomDepth(k)
+ end do
+
! Compute barotropic velocity at first timestep
! This is only done upon start-up.
if (trim(config_time_integration) == 'unsplit_explicit') then
Modified: branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_tendency.F 2012-02-14 21:41:33 UTC (rev 1508)
+++ branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_tendency.F 2012-02-14 22:01:52 UTC (rev 1509)
@@ -66,6 +66,7 @@
ocn_tend_u, &
ocn_tend_scalar, &
ocn_diagnostic_solve, &
+ ocn_div_hu, &
ocn_wtop, &
ocn_fuperp, &
ocn_tendency_init
@@ -988,18 +989,18 @@
!***********************************************************************
!
-! routine ocn_wtop
+! routine ocn_div_hu
!
-!> \brief Computes vertical velocity
-!> \author Doug Jacobsen
-!> \date 23 September 2011
+!> \brief Computes thickness weighted divergence
+!> \author Mark Petesren
+!> \date 14 February 2012
!> \version SVN:$Id$
!> \details
-!> This routine computes the vertical velocity in the top layer for the ocean
+!> This routine computes the thickness weighted divergence
!
!-----------------------------------------------------------------------
- subroutine ocn_wtop(s1,s2, grid, tend)!{{{
+ subroutine ocn_div_hu(s1,s2, grid)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields used in the tendency computations
!
@@ -1013,47 +1014,30 @@
type (state_type), intent(inout) :: s1
type (state_type), intent(inout) :: s2
type (mesh_type), intent(in) :: grid
- type (tend_type), intent(in) :: tend
- ! mrp 110512 could clean this out, remove pointers?
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum
+ integer :: iEdge, iCell, k, cell1, cell2, nCells, nEdges, nVertLevels
+ integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeBot
+ integer, dimension(:,:), pointer :: cellsOnEdge
- integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
+ real (kind=RKIND) :: flux
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, div_hu_btr
+ real (kind=RKIND), dimension(:,:), pointer :: u, h_edge, div_hu
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
- real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
- real (kind=RKIND), dimension(:,:), pointer :: u,h,wTop, h_edge, tend_hHighFreq
- real (kind=RKIND), dimension(:,:), allocatable:: div_hu
- real (kind=RKIND), dimension(:), allocatable:: div_hu_btr, h_tend_col, h_weights
-
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
- verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
- boundaryEdge, boundaryCell
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
- maxLevelVertexBot, maxLevelVertexTop
-
- h => s1 % h % array
h_edge => s1 % h_edge % array
u => s2 % u % array
- wTop => s2 % wTop % array
+ div_hu => s2 % div_hu % array
+ div_hu_btr => s2 % div_hu_btr % array
- tend_hHighFreq => tend % hHighFreq % array
-
areaCell => grid % areaCell % array
cellsOnEdge => grid % cellsOnEdge % array
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
dvEdge => grid % dvEdge % array
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertLevels = grid % nVertLevels
-
- allocate(div_hu(nVertLevels,nCells+1), div_hu_btr(nCells+1), &
- h_tend_col(nVertLevels), h_weights(nVertLevels))
-
!
! Compute div(h^{edge} u) for each cell
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
@@ -1077,6 +1061,60 @@
end do
end do
+ end subroutine ocn_div_hu!}}}
+
+!***********************************************************************
+!
+! routine ocn_wtop
+!
+!> \brief Computes vertical velocity
+!> \author Doug Jacobsen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the vertical velocity in the top layer for the ocean
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_wtop(s1,s2, grid, tend)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields used in the tendency computations
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: s - computed diagnostics
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(inout) :: s1
+ type (state_type), intent(inout) :: s2
+ type (mesh_type), intent(in) :: grid
+ type (tend_type), intent(in) :: tend
+
+ integer :: nCells, nVertLevels, iCell, k
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND) :: hSum
+ real (kind=RKIND), dimension(:), pointer :: div_hu_btr
+ real (kind=RKIND), dimension(:), allocatable:: h_tend_col, h_weights
+ real (kind=RKIND), dimension(:,:), pointer :: u,h,wTop, tend_hHighFreq, div_hu
+
+ nCells = grid % nCells
+ nVertLevels = grid % nVertLevels
+
+ h => s1 % h % array
+ u => s2 % u % array
+ wTop => s2 % wTop % array
+ div_hu => s2 % div_hu % array
+ div_hu_btr => s2 % div_hu_btr % array
+
+ tend_hHighFreq => tend % hHighFreq % array
+
+ maxLevelCell => grid % maxLevelCell % array
+
+ allocate(h_tend_col(nVertLevels), h_weights(nVertLevels))
+
!
! vertical velocity through layer interface
!
@@ -1200,7 +1238,7 @@
endif
- deallocate(div_hu, div_hu_btr, h_tend_col, h_weights)
+ deallocate(h_tend_col, h_weights)
end subroutine ocn_wtop!}}}
Modified: branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-02-14 21:41:33 UTC (rev 1508)
+++ branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-02-14 22:01:52 UTC (rev 1509)
@@ -177,7 +177,7 @@
!call ocn_tend_divergenceLowFreq(block % tend, provis, block % diagnostics, block % mesh)
endif
- ! mrp 111206 put ocn_wtop call at top for ALE
+ call ocn_div_hu(provis, provis, block % mesh)
call ocn_wtop(provis, provis, block % mesh, block % tend)
if (.not.config_implicit_vertical_mix) then
Modified: branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-14 21:41:33 UTC (rev 1508)
+++ branches/ocean_projects/ztilde/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-14 22:01:52 UTC (rev 1509)
@@ -351,6 +351,7 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN Barotropic subcycle loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
do j=1,config_n_btr_subcycles*config_btr_subcycle_loop_factor
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -407,6 +408,7 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Barotropic subcycle: SSH PREDICTOR STEP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
block => domain % blocklist
do while (associated(block))
@@ -432,7 +434,7 @@
sshEdge = 0.5 * (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
- hSum = sshEdge + block % mesh % referenceBottomDepth % array (block % mesh % maxLevelEdgeTop % array(iEdge))
+ hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
flux = ((1.0-config_btr_gam1_uWt1) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ config_btr_gam1_uWt1 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
@@ -546,7 +548,7 @@
+ config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
sshEdge = 0.5 * (sshCell1 + sshCell2)
- hSum = sshEdge + block % mesh % referenceBottomDepth % array (block % mesh % maxLevelEdgeTop % array(iEdge))
+ hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
flux = ((1.0-config_btr_gam3_uWt2) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ config_btr_gam3_uWt2 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
@@ -718,6 +720,7 @@
block => domain % blocklist
do while (associated(block))
+ call ocn_div_hu(block % state % time_levs(1) % state,block % state % time_levs(2) % state, block % mesh)
call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(2) % state, block % mesh,block % tend )
call ocn_tend_h (block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
</font>
</pre>