<p><b>mpetersen@lanl.gov</b> 2012-02-08 11:43:33 -0700 (Wed, 08 Feb 2012)</p><p>Rewrote all occurances of hZLevel, zMidZLevel, zTopZLevel using the new referenceBottomDepth variable. Right now the code expects hZLevel in the grid.nc or restart.nc file, and them computes referenceBottomDepth. This can be changed when basin.F is revised to create the referenceBottomDepth variable. Tested for bit-for-bit matching.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/Registry 2012-02-08 16:56:13 UTC (rev 1482)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/Registry 2012-02-08 18:43:33 UTC (rev 1483)
@@ -168,9 +168,8 @@
var persistent integer maxLevelEdgeBot ( nEdges ) 0 - maxLevelEdgeBot mesh - -
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 hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
-var persistent real zMidZLevel ( nVertLevels ) 0 - zMidZLevel mesh - -
-var persistent real zTopZLevel ( nVertLevelsP1 ) 0 - zTopZLevel mesh - -
% Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_equation_of_state_jm.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2012-02-08 16:56:13 UTC (rev 1482)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2012-02-08 18:43:33 UTC (rev 1483)
@@ -95,7 +95,7 @@
real (kind=RKIND), dimension(:), pointer :: &
- zMidZLevel, pRefEOS
+ referenceBottomDepth, pRefEOS
real (kind=RKIND), dimension(:,:), intent(inout) :: &
rho
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
@@ -196,7 +196,7 @@
nCells = grid % nCells
maxLevelCell => grid % maxLevelCell % array
nVertLevels = grid % nVertLevels
- zMidZLevel => grid % zMidZLevel % array
+ referenceBottomDepth => grid % referenceBottomDepth % array
! Jackett and McDougall
@@ -205,20 +205,24 @@
smin = 0.0 ! valid salinity, in psu
smax = 42.0
- ! This could be put in a startup routine.
- ! Note I am using zMidZlevel, so pressure on top level does
- ! not include SSH contribution. I am not sure if that matters.
-
! This function computes pressure in bars from depth in meters
! using a mean density derived from depth-dependent global
! average temperatures and salinities from Levitus 1994, and
! integrating using hydrostatic balance.
allocate(pRefEOS(nVertLevels),p(nVertLevels),p2(nVertLevels))
- do k = 1,nVertLevels
- depth = -zMidZLevel(k)
- pRefEOS(k) = 0.059808*(exp(-0.025*depth) - 1.0) &
- + 0.100766*depth + 2.28405e-7*depth**2
+
+ ! This could be put in the init routine.
+ ! Note I am using referenceBottomDepth, so pressure on top level does
+ ! not include SSH contribution. I am not sure if that matters, but
+ ! POP does it the same way.
+ depth = 0.5*referenceBottomDepth(1)
+ pRefEOS(1) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ + 0.100766*depth + 2.28405e-7*depth**2
+ do k = 2,nVertLevels
+ depth = 0.5*(referenceBottomDepth(k)+referenceBottomDepth(k-1))
+ pRefEOS(k) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ + 0.100766*depth + 2.28405e-7*depth**2
enddo
! If k_displaced=0, in-situ density is returned (no displacement)
Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_mpas_core.F 2012-02-08 16:56:13 UTC (rev 1482)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_mpas_core.F 2012-02-08 18:43:33 UTC (rev 1483)
@@ -118,6 +118,13 @@
call mpas_dmpar_abort(dminfo)
endif
+ if (config_filter_btr_mode.and. &
+ config_vert_grid_type.ne.'zlevel')then
+ print *, 'filter_btr_mode has only been tested with'// &
+ ' config_vert_grid_type=zlevel.'
+ call mpas_dmpar_abort(dminfo)
+ endif
+
!
! Initialize core
!
@@ -486,59 +493,46 @@
end subroutine mpas_timestep!}}}
-subroutine ocn_init_z_level(domain)!{{{
-! Initialize maxLevel and bouncary grid variables.
+ subroutine ocn_init_z_level(domain)!{{{
+ ! Initialize maxLevel and bouncary grid variables.
- use mpas_grid_types
- use mpas_configure
+ use mpas_grid_types
+ use mpas_configure
- implicit none
+ implicit none
- type (domain_type), intent(inout) :: domain
+ type (domain_type), intent(inout) :: domain
- integer :: i, iCell, iEdge, iVertex, k
- type (block_type), pointer :: block
+ integer :: i, iCell, iEdge, iVertex, k
+ type (block_type), pointer :: block
- integer :: iTracer, cell, cell1, cell2
- real (kind=RKIND) :: uhSum, hSum, hEdge1
- real (kind=RKIND), dimension(:), pointer :: &
- hZLevel, zTopZLevel, zMidZLevel
- real (kind=RKIND), dimension(:,:), pointer :: h
- integer :: nVertLevels
+ integer :: iTracer, cell, cell1, cell2
+ real (kind=RKIND) :: uhSum, hSum, hEdge1
+ real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
+ real (kind=RKIND), dimension(:,:), pointer :: h
+ integer :: nVertLevels
- ! Initialize z-level grid variables from h, read in from input file.
- block => domain % blocklist
- do while (associated(block))
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
- h => block % state % time_levs(1) % state % h % array
- hZLevel => block % mesh % hZLevel % array
- zTopZLevel => block % mesh % zTopZLevel % array
- zMidZLevel => block % mesh % zMidZLevel % array
- nVertLevels = block % mesh % nVertLevels
+ h => block % state % time_levs(1) % state % h % array
+ referenceBottomDepth => block % mesh % referenceBottomDepth % array
+ nVertLevels = block % mesh % nVertLevels
- ! hZLevel should be in the grid.nc or restart.nc file.
+ ! mrp 120208 right now hZLevel is in the grid.nc file.
+ ! We would like to transition to using referenceBottomDepth
+ ! as the defining variable instead, and will transition soon.
+ ! When the transition is done, hZLevel can be removed from
+ ! registry and the following four lines deleted.
+ referenceBottomDepth(1) = block % mesh % hZLevel % array(1)
+ do k = 2,nVertLevels
+ referenceBottomDepth(k) = referenceBottomDepth(k-1) + block % mesh % hZLevel % array(k)
+ end do
- ! zTopZLevel is used for the z-coordinate of the bottom of
- ! the water column. This is used to compute zMid from bottom
- ! up at each timestep.
-
- ! zMidZLevel is used in the EOS computation to compute pressure.
- zTopZLevel(1) = 0.0
- do k = 1,nVertLevels
- zMidZLevel(k) = zTopZLevel(k)-0.5*hZLevel(k)
- zTopZLevel(k+1) = zTopZLevel(k)- hZLevel(k)
- end do
-
- ! Compute ssh from h at level 1.
- do iCell=1,block % mesh % nCells
- block % state % time_levs(1) % state % ssh % array(iCell) &
- = block % state % time_levs(1) % state % h % array(1,iCell) &
- - block % mesh % hZLevel % array(1)
- enddo
-
! Compute barotropic velocity at first timestep
! This is only done upon start-up.
- if (trim(config_time_integration) == 'unsplit_explicit') then
+ if (trim(config_time_integration) == 'unsplit_explicit') then
block % state % time_levs(1) % state % uBtr % array(:) = 0.0
block % state % time_levs(1) % state % uBcl % array(:,:) &
@@ -549,9 +543,7 @@
if (config_filter_btr_mode) then
do iCell=1,block % mesh % nCells
block % state % time_levs(1) % state % h % array(1,iCell) &
- = block % mesh % hZLevel % array(1)
-
- block % state % time_levs(1) % state % ssh % array(iCell) = 0.0
+ = block % mesh % referenceBottomDepth % array(1)
enddo
endif
@@ -560,8 +552,8 @@
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
! uBtr = sum(u)/sum(h) on each column
- ! ocn_diagnostic_solve has not yet been called, so compute hEdge
- ! just for this edge.
+ ! ocn_diagnostic_solve has not yet been called, so compute hEdge
+ ! just for this edge.
! hSum is initialized outside the loop because on land boundaries
! maxLevelEdgeTop=0, but I want to initialize hSum with a
@@ -611,11 +603,10 @@
endif
block => block % next
+ end do
- end do
+ end subroutine ocn_init_z_level!}}}
-end subroutine ocn_init_z_level!}}}
-
subroutine ocn_compute_max_level(domain)!{{{
! Initialize maxLevel and bouncary grid variables.
Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_tendency.F 2012-02-08 16:56:13 UTC (rev 1482)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_tendency.F 2012-02-08 18:43:33 UTC (rev 1483)
@@ -437,7 +437,7 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, &
- zTopZLevel, ssh
+ referenceBottomDepth, ssh
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
@@ -485,7 +485,7 @@
areaTriangle => grid % areaTriangle % array
h_s => grid % h_s % array
fVertex => grid % fVertex % array
- zTopZLevel => grid % zTopZLevel % array
+ referenceBottomDepth => grid % referenceBottomDepth % array
deriv_two => grid % deriv_two % array
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
@@ -828,8 +828,10 @@
! Compute zMid, the z-coordinate of the middle of the layer.
! This is used for the rho g grad z momentum term.
+ ! Note the negative sign, since referenceBottomDepth is positive
+ ! and z-coordinates are negative below the surface.
k = maxLevelCell(iCell)
- zMid(k:nVertLevels,iCell) = zTopZLevel(k+1) + 0.5*h(k,iCell)
+ zMid(k:nVertLevels,iCell) = -referenceBottomDepth(k) + 0.5*h(k,iCell)
do k=maxLevelCell(iCell)-1, 1, -1
zMid(k,iCell) = zMid(k+1,iCell) &
@@ -847,9 +849,11 @@
do iCell=1,nCells
! Start at the bottom where we know the depth, and go up.
! The bottom depth for this cell is
- ! zTopZLevel(maxLevelCell(iCell)+1)
+ ! referenceBottomDepth(maxLevelCell(iCell)).
+ ! Note the negative sign, since referenceBottomDepth is positive
+ ! and z-coordinates are negative below the surface.
- ssh(iCell) = zTopZLevel(maxLevelCell(iCell)+1) &
+ ssh(iCell) = -referenceBottomDepth(maxLevelCell(iCell)) &
+ sum(h(1:maxLevelCell(iCell),iCell))
end do
Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-08 16:56:13 UTC (rev 1482)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-08 18:43:33 UTC (rev 1483)
@@ -430,19 +430,13 @@
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
-
-!TDR: seems to me that we do not want to use hZlevel to compute total depth at edges.
-!TDR: seems like before we start this process we should compute the sum of h_edge(at time n) and subtract ssh(at time n), call this depthAtRestEdge
-!TDR: then here we sum depthAtRestEdge with sshEdge to compute flux. depthAtRest{Edge,Cell} could get moved in the the "grid.nc" file at some point
-!TDR: this way we remove any reference to a field (like hZlevel) that is associated with a specific vertical grid
-
sshEdge = 0.5 * (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
- hSum = sum(block % mesh % hZLevel % array (1:block % mesh % maxLevelEdgeTop % array(iEdge)))
+ hSum = sshEdge + block % mesh % referenceBottomDepth % array (block % mesh % maxLevelEdgeTop % array(iEdge))
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)) &
- * (hSum + sshEdge)
+ * hSum
block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
@@ -552,11 +546,11 @@
+ config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
sshEdge = 0.5 * (sshCell1 + sshCell2)
- hSum = sum(block % mesh % hZLevel % array (1:block % mesh % maxLevelEdgeTop % array(iEdge)))
+ hSum = sshEdge + block % mesh % referenceBottomDepth % array (block % mesh % maxLevelEdgeTop % array(iEdge))
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)) &
- * (sshEdge + hSum)
+ * hSum
block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-02-08 16:56:13 UTC (rev 1482)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-02-08 18:43:33 UTC (rev 1483)
@@ -248,8 +248,6 @@
else
! for Ri<0 and explicit vertical mix,
! use maximum diffusion allowed by CFL criterion
- ! mrp 110324 efficiency note: for z-level, could use fixed
- ! grid array hMeanTopZLevel and compute maxdiff on startup.
vertViscTopOfEdge(k,iEdge) = &
((h_edge(k-1,iEdge)+h_edge(k,iEdge))/2.0)**2/config_dt/4.0
end if
@@ -357,8 +355,6 @@
else
! for Ri<0 and explicit vertical mix,
! use maximum diffusion allowed by CFL criterion
- ! mrp 110324 efficiency note: for z-level, could use fixed
- ! grid array hMeanTopZLevel and compute maxdiff on startup.
vertDiffTopOfCell(k,iCell) = &
((h(k-1,iCell)+h(k,iCell))/2.0)**2/config_dt/4.0
end if
Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F 2012-02-08 16:56:13 UTC (rev 1482)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F 2012-02-08 18:43:33 UTC (rev 1483)
@@ -177,21 +177,22 @@
integer :: k, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: zTopZLevel
+ real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
err = 0
if(.not.tanhViscOn) return
nVertLevels = grid % nVertLevels
- zTopZLevel => grid % zTopZLevel % array
+ referenceBottomDepth => grid % referenceBottomDepth % array
- ! zTopZLevel is used here for simplicity. Using zMid and h, which
+ ! referenceBottomDepth is used here for simplicity. Using zMid and h, which
! vary in time, would give the exact location of the top, but it
! would only change the diffusion value very slightly.
- do k=1,nVertLevels+1
- vertViscTopOfEdge(k,:) = -(config_max_visc_tanh-config_min_visc_tanh)/2.0 &
- *tanh(-(zTopZLevel(k)-config_ZMid_tanh) &
+ vertViscTopOfEdge = 0.0
+ do k=2,nVertLevels
+ vertViscTopOfEdge(k,:) = -(config_max_visc_tanh-config_min_visc_tanh)/2.0 &
+ *tanh((referenceBottomDepth(k-1)+config_ZMid_tanh) &
/config_zWidth_tanh) &
+ (config_max_visc_tanh+config_min_visc_tanh)/2
end do
@@ -249,21 +250,22 @@
integer :: k, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: zTopZLevel
+ real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
err = 0
if(.not.tanhDiffOn) return
nVertLevels = grid % nVertLevels
- zTopZLevel => grid % zTopZLevel % array
+ referenceBottomDepth => grid % referenceBottomDepth % array
- ! zTopZLevel is used here for simplicity. Using zMid and h, which
+ ! referenceBottomDepth is used here for simplicity. Using zMid and h, which
! vary in time, would give the exact location of the top, but it
! would only change the diffusion value very slightly.
- do k=1,nVertLevels+1
+ vertDiffTopOfCell = 0.0
+ do k=2,nVertLevels
vertDiffTopOfCell(k,:) = -(config_max_diff_tanh-config_min_diff_tanh)/2.0 &
- *tanh(-(zTopZLevel(k)-config_ZMid_tanh) &
+ *tanh((referenceBottomDepth(k-1)+config_ZMid_tanh) &
/config_zWidth_tanh) &
+ (config_max_diff_tanh+config_min_diff_tanh)/2
end do
</font>
</pre>