Zero out normalVelocity in RK4

E3SM-Project · Aug 29, 2022 · 3757046 · 3757046
1 parent 3c9f1c5
commit 3757046
Showing 1 changed file with 16 additions and 5 deletions.
diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F
@@ -1135,7 +1135,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight,
 
       logical, pointer :: config_prescribe_velocity, config_prescribe_thickness, config_use_GM
 
-      integer, pointer :: nCells, nEdges
+      integer, pointer :: nCells, nEdges, nVertLevels
       integer :: iCell, iEdge, k
 
       type (mpas_pool_type), pointer :: statePool, tendPool, meshPool, scratchPool
@@ -1173,6 +1173,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight,
 
       call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells)
       call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges)
+      call mpas_pool_get_dimension(block % dimensions, 'nVertLevels', nVertLevels)
 
       call mpas_pool_get_subpool(block % structs, 'state', statePool)
       call mpas_pool_get_subpool(block % structs, 'tend', tendPool)
@@ -1213,6 +1214,9 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight,
                                            * normalVelocityTend(k, iEdge)
             normalVelocityProvis(k, iEdge) = normalVelocityProvis(k, iEdge) * (1.0_RKIND - wettingVelocity(k, iEdge))
          end do
+         do k = maxLevelEdgeTop(iEdge)+1, nVertLevels
+            normalVelocityProvis(k, iEdge) = 0.0_RKIND
+         end do
       end do
       !$omp end do
 
@@ -1384,7 +1388,7 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{
       real (kind=RKIND), intent(in) :: rkWeight
       integer, intent(out) :: err
 
-      integer, pointer :: nCells, nEdges
+      integer, pointer :: nCells, nEdges, nVertLevels
       integer :: iCell, iEdge, k
 
       type (mpas_pool_type), pointer :: statePool, tendPool, meshPool
@@ -1397,7 +1401,7 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{
 
       real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroupNew, tracersGroupTend
 
-      integer, dimension(:), pointer :: minLevelCell, maxLevelCell
+      integer, dimension(:), pointer :: minLevelCell, maxLevelCell, maxLevelEdgeTop
 
       logical, pointer :: config_use_tracerGroup
       type (mpas_pool_iterator_type) :: groupItr
@@ -1408,6 +1412,7 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{
 
       call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells)
       call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges)
+      call mpas_pool_get_dimension(block % dimensions, 'nVertLevels', nVertLevels)
 
       call mpas_pool_get_subpool(block % structs, 'state', statePool)
       call mpas_pool_get_subpool(block % structs, 'tend', tendPool)
@@ -1434,6 +1439,7 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{
 
       call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell)
       call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
+      call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop)
 
       !$omp parallel
       !$omp do schedule(runtime) private(k)
@@ -1446,8 +1452,13 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{
 
       !$omp do schedule(runtime)
       do iEdge = 1, nEdges
-         normalVelocityNew(:, iEdge) = normalVelocityNew(:, iEdge) + rkWeight * normalVelocityTend(:, iEdge)
-         normalVelocityNew(:, iEdge) = normalVelocityNew(:, iEdge) * (1.0_RKIND - wettingVelocity(:, iEdge))
+         do k = 1, maxLevelEdgeTop(iEdge)
+           normalVelocityNew(k, iEdge) = normalVelocityNew(k, iEdge) + rkWeight * normalVelocityTend(k, iEdge)
+           normalVelocityNew(k, iEdge) = normalVelocityNew(k, iEdge) * (1.0_RKIND - wettingVelocity(k, iEdge))
+         end do
+         do k = maxLevelEdgeTop(iEdge)+1, nVertLevels
+            normalVelocityNew(k, iEdge) = 0.0_RKIND
+         end do
       end do
       !$omp end do
       !$omp end parallel