Skip to content

Commit

Permalink
Merge PR #825 'cbegeman/ocean/new-min-level-cell-feature' into ocean/…
Browse files Browse the repository at this point in the history 
…develop

Add the option for inactive top cells #825

Inactive top cells are introduced to add flexibility to the vertical
coordinate below ice shelves.

This PR introduces phase 1 in the development documented in #811, which
touches only the code needed to run a basic ice shelf test case.

The variable minLevelCell and related minLevelEdge* and minLevelVertex*
variables are introduced by analogy with maxLevel* variables. By
default, there are no inactive top cells and minLevel* = 1. The main
code changes are to k-loop limits and designating surface quantities as
those at minLevel*.
  • Loading branch information
@mark-petersen
mark-petersen committed Mar 22, 2021
2 parents 683ea77 + 8329b9c commit 29a819a
Show file tree
Hide file tree
Showing 29 changed files with 890 additions and 408 deletions.
22 changes: 22 additions & 0 deletions src/core_ocean/Registry.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1412,6 +1412,7 @@
<var name="refBottomDepth"/>
<var name="bottomDepth"/>
<var name="bottomDepthObserved"/>
<var name="minLevelCell"/>
<var name="maxLevelCell"/>
<var name="vertCoordMovementWeights"/>
<var name="edgeMask"/>
Expand Down Expand Up @@ -1506,6 +1507,7 @@
<var name="fVertex"/>
<var name="fCell"/>
<var name="bottomDepth"/>
<var name="minLevelCell"/>
<var name="maxLevelCell"/>
<var name="refBottomDepth" packages="forwardMode;analysisMode"/>
<var name="restingThickness" packages="forwardMode;analysisMode"/>
Expand Down Expand Up @@ -1646,6 +1648,7 @@
<var name="layerThickness"/>
<var name="ssh"/>
<var name="maxLevelEdgeTop"/>
<var name="minLevelEdgeBot"/>
<var name="vertCoordMovementWeights"/>
<var name="edgeMask"/>
<var name="vertexMask"/>
Expand Down Expand Up @@ -2163,21 +2166,40 @@
<var name="coeffs_reconstruct" type="real" dimensions="R3 maxEdges nCells" units="unitless"
description="Coefficients to reconstruct velocity vectors at cells centers."
/>
<var name="minLevelCell" type="integer" dimensions="nCells" units="unitless" default_value="1"
description="Index to the first active ocean cell in each column."
/>
<var name="maxLevelCell" type="integer" dimensions="nCells" units="unitless"
description="Index to the last active ocean cell in each column."
/>
<var name="minLevelEdgeTop" type="integer" dimensions="nEdges" units="unitless"
description="Index to the first edge in a column with at least one active ocean cell on either side of it."
packages="forwardMode;analysisMode"
/>
<var name="maxLevelEdgeTop" type="integer" dimensions="nEdges" units="unitless"
description="Index to the last edge in a column with active ocean cells on both sides of it."
packages="forwardMode;analysisMode"
/>
<var name="minLevelEdgeBot" type="integer" dimensions="nEdges" units="unitless"
description="Index to the first edge in a column with active ocean cells on both sides of it."
packages="forwardMode;analysisMode"
/>
<var name="maxLevelEdgeBot" type="integer" dimensions="nEdges" units="unitless"
description="Index to the last edge in a column with at least one active ocean cell on either side of it."
packages="forwardMode;analysisMode"
/>
<var name="minLevelVertexTop" type="integer" dimensions="nVertices" units="unitless"
description="Index to the first vertex in a column with at least one active ocean cell around it."
packages="forwardMode;analysisMode"
/>
<var name="maxLevelVertexTop" type="integer" dimensions="nVertices" units="unitless"
description="Index to the last vertex in a column with all active cells around it."
packages="forwardMode;analysisMode"
/>
<var name="minLevelVertexBot" type="integer" dimensions="nVertices" units="unitless"
description="Index to the first vertex in a column with all active cells around it."
packages="forwardMode;analysisMode"
/>
<var name="maxLevelVertexBot" type="integer" dimensions="nVertices" units="unitless"
description="Index to the last vertex in a column with at least one active ocean cell around it."
packages="forwardMode;analysisMode"
Expand Down
2 changes: 1 addition & 1 deletion src/core_ocean/mode_analysis/mpas_ocn_analysis_mode.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@ function ocn_analysis_mode_init(domain, startTimeStamp) result(ierr)!{{{

call ocn_init_routines_vert_coord(domain)

call ocn_init_routines_compute_max_level(domain)
call ocn_init_routines_compute_min_max_level(domain)

timeStep = mpas_get_clock_timestep(domain % clock, ierr=err_tmp)
call mpas_get_timeInterval(timeStep, dt=dt)
Expand Down
2 changes: 1 addition & 1 deletion src/core_ocean/mode_forward/mpas_ocn_forward_mode.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -253,7 +253,7 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{

call ocn_init_routines_vert_coord(domain)

call ocn_init_routines_compute_max_level(domain)
call ocn_init_routines_compute_min_max_level(domain)

call mpas_pool_get_config(domain % configs, 'config_time_integrator', config_time_integrator)

Expand Down
72 changes: 37 additions & 35 deletions src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -387,7 +387,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp do schedule(runtime) private(k)
do iCell = 1, nCellsAll
sshNew(iCell) = sshCur(iCell)
do k = 1, maxLevelCell(iCell)
do k = minLevelCell(iCell), maxLevelCell(iCell)
layerThicknessNew(k,iCell) = layerThicknessCur(k,iCell)
end do
end do
Expand All @@ -408,7 +408,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp parallel
!$omp do schedule(runtime) private(k)
do iCell = 1, nCellsAll
do k = 1, maxLevelCell(iCell)
do k = minLevelCell(iCell), maxLevelCell(iCell)
tracersGroupNew(:,k,iCell) = &
tracersGroupCur(:,k,iCell)
end do
Expand Down Expand Up @@ -506,7 +506,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp parallel
!$omp do schedule(runtime) private(k)
do iCell = 1, nCellsAll
do k = 1, maxLevelCell(iCell)
do k = minLevelCell(iCell), maxLevelCell(iCell)
! this is h^{hf}_{n+1}
highFreqThicknessNew(k,iCell) = &
highFreqThicknessCur(k,iCell) + dt* &
Expand Down Expand Up @@ -573,7 +573,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
cell2 = cellsOnEdge(2,iEdge)
! could put this after with uTemp(maxleveledgetop+1:nvertlevels)=0
uTemp = 0.0_RKIND
do k = 1, maxLevelEdgeTop(iEdge)
do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
! normalBaroclinicVelocityNew =
! normalBaroclinicVelocityOld +
Expand All @@ -593,11 +593,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
! on land boundaries maxLevelEdgeTop=0, but we want to
! initialize thicknessSum with a nonzero value to avoid
! a NaN.
normalThicknessFluxSum = layerThickEdge(1,iEdge)* &
uTemp(1)
thicknessSum = layerThickEdge(1,iEdge)
normalThicknessFluxSum = layerThickEdge(minLevelEdgeBot(iEdge),iEdge)* &
uTemp(minLevelEdgeBot(iEdge))
thicknessSum = layerThickEdge(minLevelEdgeBot(iEdge),iEdge)
do k = 2, maxLevelEdgeTop(iEdge)
do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge)
normalThicknessFluxSum = normalThicknessFluxSum + &
layerThickEdge(k,iEdge)*uTemp(k)
thicknessSum = thicknessSum + &
Expand All @@ -606,7 +606,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
barotropicForcing(iEdge) = splitFact* &
normalThicknessFluxSum/thicknessSum/dt
do k = 1, maxLevelEdgeTop(iEdge)
do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
! These two steps are together here:
!{\bf u}'_{k,n+1} =
! {\bf u}'_{k,n} - \Delta t {\overline {\bf G}}
Expand Down Expand Up @@ -1395,11 +1395,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
! on land boundaries maxLevelEdgeTop=0, but I want to
! initialize thicknessSum with a nonzero value to avoid
! a NaN.
normalThicknessFluxSum = layerThickEdge(1,iEdge)* &
uTemp(1)
thicknessSum = layerThickEdge(1,iEdge)
normalThicknessFluxSum = layerThickEdge(minLevelEdgeBot(iEdge),iEdge)* &
uTemp(minLevelEdgeBot(iEdge))
thicknessSum = layerThickEdge(minLevelEdgeBot(iEdge),iEdge)

do k = 2, maxLevelEdgeTop(iEdge)
do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge)
normalThicknessFluxSum = normalThicknessFluxSum + &
layerThickEdge(k,iEdge)*&
uTemp(k)
Expand Down Expand Up @@ -1542,7 +1542,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp parallel
!$omp do schedule(runtime) private(i, k, temp_h, temp)
do iCell = 1, nCellsAll
do k = 1, maxLevelCell(iCell)
do k = minLevelCell(iCell), maxLevelCell(iCell)

! this is h_{n+1}
temp_h = layerThicknessCur(k,iCell) + dt* &
Expand Down Expand Up @@ -1572,7 +1572,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp parallel
!$omp do schedule(runtime) private(k, temp)
do iCell = 1, nCellsAll
do k = 1, maxLevelCell(iCell)
do k = minLevelCell(iCell), maxLevelCell(iCell)

! h^{hf}_{n+1} was computed in Stage 1

Expand Down Expand Up @@ -1628,7 +1628,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp parallel
!$omp do schedule(runtime) private(k)
do iCell = 1, nCellsAll
do k = 1, maxLevelCell(iCell)
do k = minLevelCell(iCell), maxLevelCell(iCell)
! this is h_{n+1}
layerThicknessNew(k,iCell) = &
layerThicknessCur(k,iCell) + &
Expand All @@ -1643,7 +1643,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp parallel
!$omp do schedule(runtime) private(k)
do iEdge = 1, nEdgesAll
do k= 1, maxLevelEdgeTop(iEdge)
do k= minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
activeTracerHorizontalAdvectionEdgeFlux(:,k,iEdge) = &
activeTracerHorizontalAdvectionEdgeFlux(:,k,iEdge) / &
layerThickEdge(k,iEdge)
Expand All @@ -1653,7 +1653,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

!$omp do schedule(runtime) private(k)
do iCell = 1, nCellsAll
do k= 1, maxLevelCell(iCell)
do k= minLevelCell(iCell), maxLevelCell(iCell)
activeTracerHorizontalAdvectionTendency(:,k,iCell) = &
activeTracerHorizontalAdvectionTendency(:,k,iCell) / &
layerThicknessNew(k,iCell)
Expand Down Expand Up @@ -1708,7 +1708,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp parallel
!$omp do schedule(runtime) private(k)
do iCell = 1, nCellsAll
do k = 1, maxLevelCell(iCell)
do k = minLevelCell(iCell), maxLevelCell(iCell)
tracersGroupNew(:,k,iCell) = &
(tracersGroupCur(:,k,iCell) * &
layerThicknessCur(k,iCell) + dt* &
Expand All @@ -1725,7 +1725,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp parallel
!$omp do schedule(runtime) private(k)
do iCell = 1, nCellsAll
do k = 1, maxLevelCell(iCell)
do k = minLevelCell(iCell), maxLevelCell(iCell)
tracersGroupNew(indexSalinity,k,iCell) = &
max(0.001_RKIND, &
tracersGroupNew(indexSalinity,k,iCell))
Expand All @@ -1745,7 +1745,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
do iCell = 1, nCellsAll

! Reset tracer1 to 2 in top n layers
do k = 1, config_reset_debugTracers_top_nLayers
do k = minLevelCell(iCell), minLevelCell(iCell)+config_reset_debugTracers_top_nLayers-1
tracersGroupNew(1,k,iCell) = 2.0_RKIND
end do

Expand All @@ -1757,11 +1757,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
.or.lat>- 2.5.and.lat< 2.5 &
.or.lat> 35.0.and.lat< 40.0 &
.or.lat> 55.0.and.lat< 60.0 ) then
do k = 1, config_reset_debugTracers_top_nLayers
do k = minLevelCell(iCell), minLevelCell(iCell)+config_reset_debugTracers_top_nLayers-1
tracersGroupNew(2,k,iCell) = 2.0_RKIND
end do
else
do k = 1, config_reset_debugTracers_top_nLayers
do k = minLevelCell(iCell), minLevelCell(iCell)+config_reset_debugTracers_top_nLayers-1
tracersGroupNew(2,k,iCell) = 1.0_RKIND
end do
end if
Expand All @@ -1775,11 +1775,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
.or.lat> 10.0.and.lat< 15.0 &
.or.lat> 30.0.and.lat< 35.0 &
.or.lat> 50.0.and.lat< 55.0 ) then
do k = 1, config_reset_debugTracers_top_nLayers
do k = minLevelCell(iCell), minLevelCell(iCell)+config_reset_debugTracers_top_nLayers-1
tracersGroupNew(3,k,iCell) = 2.0_RKIND
end do
else
do k = 1, config_reset_debugTracers_top_nLayers
do k = minLevelCell(iCell), minLevelCell(iCell)+config_reset_debugTracers_top_nLayers-1
tracersGroupNew(3,k,iCell) = 1.0_RKIND
end do
end if
Expand All @@ -1795,7 +1795,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp parallel
!$omp do schedule(runtime) private(k)
do iCell = 1, nCellsAll
do k = 1, maxLevelCell(iCell)
do k = minLevelCell(iCell), maxLevelCell(iCell)

! h^{hf}_{n+1} was computed in Stage 1

Expand Down Expand Up @@ -1828,7 +1828,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp parallel
!$omp do schedule(runtime) private(k)
do iEdge = 1, nEdgesAll
do k = 1, maxLevelEdgeTop(iEdge)
do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
normalVelocityNew(k,iEdge) = &
normalBarotropicVelocityNew(iEdge) + &
2*normalBaroclinicVelocityNew(k,iEdge) - &
Expand Down Expand Up @@ -1956,9 +1956,9 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp do schedule(runtime)
do iCell = 1, nCellsAll
surfaceVelocity(indexSurfaceVelocityZonal,iCell) = &
velocityZonal(1,iCell)
velocityZonal(minLevelCell(iCell),iCell)
surfaceVelocity(indexSurfaceVelocityMeridional,iCell) = &
velocityMeridional(1,iCell)
velocityMeridional(minLevelCell(iCell),iCell)
SSHGradient(indexSSHGradientZonal,iCell) = gradSSHZonal(iCell)
SSHGradient(indexSSHGradientMeridional,iCell) = &
Expand Down Expand Up @@ -2152,6 +2152,8 @@ subroutine ocn_time_integration_split_init(domain)!{{{
call mpas_pool_get_array(meshPool, 'refBottomDepth', refBottomDepth)
call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell)
call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', minLevelEdgeBot)
call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop)
!*** Compute barotropic velocity at first timestep
Expand All @@ -2168,7 +2170,7 @@ subroutine ocn_time_integration_split_init(domain)!{{{
if (config_filter_btr_mode) then
do iCell = 1, nCells
layerThickness(1,iCell) = refBottomDepth(1)
layerThickness(minLevelCell(iCell),iCell) = refBottomDepth(minLevelCell(iCell))
enddo
endif
Expand All @@ -2186,13 +2188,13 @@ subroutine ocn_time_integration_split_init(domain)!{{{
! initialize thicknessSum with a nonzero value to avoid
! a NaN.
layerThicknessEdge1 = 0.5_RKIND* &
(layerThickness(1,cell1) + &
layerThickness(1,cell2) )
(layerThickness(minLevelCell(cell1),cell1) + &
layerThickness(minLevelCell(cell2),cell2) )
normalThicknessFluxSum = layerThicknessEdge1* &
normalVelocity(1,iEdge)
normalVelocity(minLevelEdgeBot(iEdge),iEdge)
layerThicknessSum = layerThicknessEdge1
do k=2, kmax
do k=minLevelEdgeBot(iEdge)+1, kmax
layerThicknessEdge1 = 0.5_RKIND* &
(layerThickness(k,cell1) + &
layerThickness(k,cell2))
Expand All @@ -2209,7 +2211,7 @@ subroutine ocn_time_integration_split_init(domain)!{{{
! normalBaroclinicVelocity = normalVelocity -
! normalBarotropicVelocity
do k = 1, kmax
do k = minLevelEdgeBot(iEdge), kmax
normalBaroclinicVelocity(k,iEdge) = &
normalVelocity(k,iEdge) - &
normalBarotropicVelocity(iEdge)
Expand Down
Loading

0 comments on commit 29a819a

Please sign in to comment.