Merge branch 'cbegeman/ocn/bugfix-inactive-top-cells' into next (PR #…

…5246) This PR fixes 2 bugs that affect the inactive top cells feature (i.e., when minLevelCell > 1). One in ocn_diagnostic_solve_ssh and the other affects the sea ice salinity flux in ocn_surface_bulk_forcing. No changes in the solution should occur when there are no inactive top cells. This PR also includes a cleanup of Redi. The original clunky workaround was designed to avoid non-bfb changes in BGC tracers, and this revision introduces those potentially non-bfb changes as @maltrud determined that they are insignificant. Fixes #5245 [BFB]
E3SM-Project · Nov 7, 2022 · a1e57b4 · a1e57b4
2 parents 350203c + 719a98f
commit a1e57b4
Show file tree

Hide file tree

Showing 6 changed files with 14 additions and 14 deletions.
diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_vertical_grids.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_vertical_grids.F
@@ -791,6 +791,7 @@ subroutine ocn_compute_Haney_number(domain, iErr)
               dzVert2 = zMid(k-1,c2)-zMid(k,c2)
               dzEdgeK = zMid(k-1,c2)-zMid(k-1,c1)
             end if
+            if ( dzVert1+dzVert2 == 0.0_RKIND ) cycle
             dzEdgeKp1 = zMid(k,c2)-zMid(k,c1)
 
             rx1 = abs(dzEdgeK+dzEdgeKp1)/(dzVert1+dzVert2)

diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile
@@ -189,7 +189,7 @@ mpas_ocn_transport_tests.o: mpas_ocn_config.o
 
 mpas_ocn_effective_density_in_land_ice.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o
 
-mpas_ocn_forcing_restoring.o: mpas_ocn_constants.o mpas_ocn_config.o
+mpas_ocn_forcing_restoring.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o
 
 mpas_ocn_tracer_surface_restoring.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_framework_forcing.o mpas_ocn_diagnostics_variables.o
 

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F
@@ -2587,7 +2587,7 @@ subroutine ocn_diagnostic_solve_ssh(forcingPool, ssh, seaIcePressure, &
          cell1 = cellsOnEdge(1, iEdge)
          cell2 = cellsOnEdge(2, iEdge)
 
-         gradSSH(iEdge) = edgeMask(1, iEdge) * ( pressureAdjustedSSH(cell2) - pressureAdjustedSSH(cell1) ) / dcEdge(iEdge)
+         gradSSH(iEdge) = edgeMask(minLevelEdgeBot(iEdge), iEdge) * ( pressureAdjustedSSH(cell2) - pressureAdjustedSSH(cell1) ) / dcEdge(iEdge)
       end do
 #ifndef MPAS_OPENACC
       !$omp end do

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F b/components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F
@@ -519,8 +519,8 @@ subroutine ocn_surface_bulk_forcing_active_tracers(meshPool, forcingPool, tracer
 
         ! Negative seaIceSalinityFlux is an extraction of salt from the ocean
         ! So, we negate seaIceSalinityFlux when determining how much salt this flux needs.
-        requiredSalt = - seaIceSalinityFlux(iCell) * sflux_factor * dt / layerThickness(1, iCell)
-        allowedSalt = min( 4.0_RKIND, tracerGroup(index_salinity_flux, 1, iCell) )
+        requiredSalt = - seaIceSalinityFlux(iCell) * sflux_factor * dt / layerThickness(minLevelCell(iCell), iCell)
+        allowedSalt = min( 4.0_RKIND, tracerGroup(index_salinity_flux, minLevelCell(iCell), iCell) )
 
         if ( allowedSalt < requiredSalt ) then
            tracersSurfaceFluxRemoved(index_salinity_flux, iCell) = tracersSurfaceFluxRemoved(index_salinity_flux, iCell)  &

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F
@@ -235,10 +235,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
             r_tmp = dvEdge(iEdge)/dcEdge(iEdge)
             coef = dvEdge(iEdge)
 
-            k = minLevelEdgeBot(iEdge)
-            kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + &
-                                           RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))
             k = 1
+            if (minLevelEdgeBot(iEdge) .ne. 1) k = minLevelEdgeBot(iEdge)
             kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + &
                                            RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))
 
@@ -259,13 +257,13 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
                              /(zMid(k, cell2) - zMid(k + 1, cell2)))
                if (minLevelCell(cell1) < minLevelEdgeBot(iEdge)) then
                   flux_term2 = flux_term2 + coef*RediKappa(iEdge)* &
-                               layerThickEdgeMean(minLevelEdgeBot(iEdge), iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(minLevelEdgeBot(iEdge), 1, iEdge)* &
-                               (tracers(iTr, minLevelEdgeBot(iEdge) - 1, cell1) - tracers(iTr, minLevelEdgeBot(iEdge), cell1))/(zMid(minLevelEdgeBot(iEdge) - 1, cell1) - zMid(minLevelEdgeBot(iEdge), cell1))
+                               layerThickEdgeMean(k, iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(k, 1, iEdge)* &
+                               (tracers(iTr, k - 1, cell1) - tracers(iTr, k, cell1))/(zMid(k - 1, cell1) - zMid(k, cell1))
                endif
                if (minLevelCell(cell2) < minLevelEdgeBot(iEdge)) then
                   flux_term2 = flux_term2 + coef*RediKappa(iEdge)* &
-                               layerThickEdgeMean(minLevelEdgeBot(iEdge), iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(minLevelEdgeBot(iEdge), 2, iEdge)* &
-                               (tracers(iTr, minLevelEdgeBot(iEdge) - 1, cell2) - tracers(iTr, minLevelEdgeBot(iEdge), cell2))/(zMid(minLevelEdgeBot(iEdge) - 1, cell2) - zMid(minLevelEdgeBot(iEdge), cell2))
+                               layerThickEdgeMean(k, iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(k, 2, iEdge)* &
+                               (tracers(iTr, k - 1, cell2) - tracers(iTr, k, cell2))/(zMid(k - 1, cell2) - zMid(k, cell2))
                endif
 
                redi_term2(iTr, k, iCell) = redi_term2(iTr, k, iCell) - edgeSignOnCell(i, iCell)*flux_term2*invAreaCell

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_surface_restoring.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_surface_restoring.F
@@ -30,6 +30,7 @@ module ocn_tracer_surface_restoring
    use ocn_config
    use ocn_framework_forcing
    use ocn_diagnostics_variables
+   use ocn_mesh
 
    implicit none
 
@@ -135,11 +136,11 @@ subroutine ocn_tracer_surface_restoring_compute(groupName, nTracers, nCells, tra
 
       err = 0
 
-      iLevel = 1  ! base surface flux restoring on tracer fields in the top layer
-
       !$omp parallel
-      !$omp do schedule(runtime) private(iTracer)
+      !$omp do schedule(runtime) private(iTracer, iLevel)
       do iCell=1,nCells
+        iLevel = minLevelCell(iCell)  ! base surface flux restoring on tracer fields in the top layer
+
         do iTracer=1,nTracers
 
           ! For monthly salinity restoring, tracersSurfaceRestoringValue contains the zero-mean deltaS