Merge branch 'mark-petersen/ocn/Redi-min-layers-diag-terms' (PR #4785)

Turn off Redi diagonal terms for shallow columns The MPAS-Ocean Redi diagonal terms are not guaranteed to produce bounded tracer fields, and in practice produce growing temperatures in a few ocean columns with fewer than 5 vertical cells in particular simulations. Redi is meant for isopycnal mixing in the deep ocean, so not applying Redi diagonal terms in very shallow regions is an acceptable solution. The problematic behavior is described in MPAS-Dev/compass#308. Here we add the flag config_Redi_min_layers_diag_terms. Redi diagonal terms (2 and 3) are turned off from layer 1 through config_Redi_min_layers_diag_terms-1, and on from config_Redi_min_layers_diag_terms through nVertLevels. If config_Redi_min_layers_diag_terms is set to 1, Redi will be applied as before this PR, and simulations should compare bit-for-bit or machine precision to previous simulations. We are adding the default config_Redi_min_layers_diag_terms=6 (no Redi applied in columns with 5 or fewer vertical cells). I found in sensitivity tests that values of 5 and 6 produce no non-monotonic warming in the problematic case of WC14 ocean-only spin-up with no surface forcing and rough bathymetry (the most challenging case, see MPAS-Dev/compass#308) but a value of 4 produces increasing temperatures. Fixes #4784 [NML] [non-BFB]
E3SM-Project · Mar 18, 2022 · 6f7d0ba · 6f7d0ba
2 parents 11f02f2 + 669c0de
commit 6f7d0ba
Show file tree

Hide file tree

Showing 6 changed files with 71 additions and 32 deletions.
diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist
@@ -571,6 +571,7 @@ add_default($nl, 'config_Redi_use_surface_taper');
 add_default($nl, 'config_Redi_N2_based_taper_enable');
 add_default($nl, 'config_Redi_N2_based_taper_min');
 add_default($nl, 'config_Redi_N2_based_taper_limit_term1');
+add_default($nl, 'config_Redi_min_layers_diag_terms');
 
 ############################################
 # Namelist group: GM_eddy_parameterization #

diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section
@@ -114,6 +114,7 @@ add_default($nl, 'config_Redi_use_surface_taper');
 add_default($nl, 'config_Redi_N2_based_taper_enable');
 add_default($nl, 'config_Redi_N2_based_taper_min');
 add_default($nl, 'config_Redi_N2_based_taper_limit_term1');
+add_default($nl, 'config_Redi_min_layers_diag_terms');
 
 ############################################
 # Namelist group: GM_eddy_parameterization #

diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
@@ -152,6 +152,7 @@
 <config_Redi_N2_based_taper_enable>.true.</config_Redi_N2_based_taper_enable>
 <config_Redi_N2_based_taper_min>0.1</config_Redi_N2_based_taper_min>
 <config_Redi_N2_based_taper_limit_term1>.true.</config_Redi_N2_based_taper_limit_term1>
+<config_Redi_min_layers_diag_terms>6</config_Redi_min_layers_diag_terms>
 
 <!-- GM_eddy_parameterization -->
 <config_use_GM>.true.</config_use_GM>

diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
@@ -493,6 +493,14 @@ Valid values: .true. or .false.
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_Redi_min_layers_diag_terms" type="integer"
+	category="Redi_isopycnal_mixing" group="Redi_isopycnal_mixing">
+Redi diagonal terms (2 and 3) are turned off from layer 1 through config_Redi_min_layers_diag_terms-1, and on from config_Redi_min_layers_diag_terms to nVertLevels. The Redi diagonal terms are not guaranteed to produce bounded tracer fields, and in practice produce growing temperatures in a few columns with fewer than 5 vertical cells. Redi is meant for isopycnal mixing in the deep ocean, so not applying Redi diagonal terms in very shallow regions is an acceptable solution.
+
+Valid values: any integer between 0 (all layers on) and nVertLevels (all layers off)
+Default: Defined in namelist_defaults.xml
+</entry>
+
 
 <!-- GM_eddy_parameterization -->
 

diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
@@ -337,9 +337,13 @@
 					possible_values="greater than zero and less than 1"
 		/>
 		<nml_option name="config_Redi_N2_based_taper_limit_term1" type="logical" default_value=".true." units="unitless"
-					description="If true, the N2 lim iting is applied to the horizontal diffusion term"
+					description="If true, the N2 limiting is applied to the horizontal diffusion term"
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_Redi_min_layers_diag_terms" type="integer" default_value="6" units="NA"
+					description="Redi diagonal terms (2 and 3) are turned off from layer 1 through config_Redi_min_layers_diag_terms-1, and on from config_Redi_min_layers_diag_terms to nVertLevels. The Redi diagonal terms are not guaranteed to produce bounded tracer fields, and in practice produce growing temperatures in a few columns with fewer than 5 vertical cells. Redi is meant for isopycnal mixing in the deep ocean, so not applying Redi diagonal terms in very shallow regions is an acceptable solution."
+					possible_values="any integer between 0 (all layers on) and nVertLevels (all layers off)"
+		/>
 	</nml_record>
 	<nml_record name="GM_eddy_parameterization" mode="forward">
 		<nml_option name="config_use_GM" type="logical" default_value=".false." units="NA"

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
@@ -130,6 +130,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg
 
       integer :: iCell, iEdge, cell1, cell2, iCellSelf
       integer :: i, k, iTr, nTracers, nCells, nEdges, nCellsP1
+      logical :: use_Redi_diag_terms
       integer, pointer :: nVertLevels
       integer, dimension(:), pointer :: nCellsArray, nEdgesArray
 
@@ -200,8 +201,14 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg
       !$omp do schedule(runtime) &
       !$omp private(fluxRediZTop, invAreaCell, i, iEdge, cell1, cell2, iCellSelf, r_tmp, coef, k, &
       !$omp         kappaRediEdgeVal, iTr, tracer_turb_flux, flux, flux_term2, flux_term3, &
-      !$omp         dTracerDx)
+      !$omp         dTracerDx, use_Redi_diag_terms)
       do iCell = 1, nCells
+         if (maxLevelCell(iCell) .ge. config_Redi_min_layers_diag_terms) then
+            use_Redi_diag_terms = .true.
+         else
+            use_Redi_diag_terms = .false.
+         endif
+
          redi_term1(:, :, iCell) = 0.0_RKIND
          redi_term2(:, :, iCell) = 0.0_RKIND
          redi_term3(:, :, iCell) = 0.0_RKIND
@@ -240,6 +247,9 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg
                redi_term1(iTr, k, iCell) = redi_term1(iTr, k, iCell) - (1.0_RKIND + term1TaperFactor* &
                                            (kappaRediEdgeVal - 1.0_RKIND))*edgeSignOnCell(i, iCell)*flux*invAreaCell
 
+               if (.not.use_Redi_diag_terms) cycle
+
+               ! Term 2: div( h S dphi/dz)
                flux_term2 = coef*kappaRediEdgeVal*RediKappa(iEdge)*layerThicknessEdge(k, iEdge)* &
                             0.25_RKIND* &
                             (slopeTriadDown(k, 1, iEdge)*(tracers(iTr, k, cell1) - tracers(iTr, k + 1, cell1)) &
@@ -278,6 +288,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg
                   redi_term1(iTr, k, iCell) = redi_term1(iTr, k, iCell) - (1.0_RKIND + term1TaperFactor* &
                                               (kappaRediEdgeVal - 1.0_RKIND))*edgeSignOnCell(i, iCell)*flux*invAreaCell
 
+                  if (.not.use_Redi_diag_terms) cycle
+
                   flux_term2 = coef*RediKappa(iEdge)*kappaRediEdgeVal*layerThicknessEdge(k, iEdge)* &
                                0.25_RKIND* &
                                (slopeTriadUp(k, 1, iEdge)*(tracers(iTr, k - 1, cell1) - tracers(iTr, k, cell1)) &
@@ -311,6 +323,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg
                redi_term1(iTr, k, iCell) = redi_term1(iTr, k, iCell) - (1.0_RKIND + term1TaperFactor* &
                                            (kappaRediEdgeVal - 1.0_RKIND))*edgeSignOnCell(i, iCell)*flux*invAreaCell
 
+               if (.not.use_Redi_diag_terms) cycle
+
                ! For bottom layer, only use triads pointing up:
                flux_term2 = coef*kappaRediEdgeVal*RediKappa(iEdge)*layerThicknessEdge(k, iEdge)* &
                             0.25_RKIND* &
@@ -340,28 +354,31 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg
             end do
          end do ! nEdgesOnCell(iCell)
 
-         ! impose no-flux boundary conditions at top and bottom of column
-         fluxRediZTop(:, 1:minLevelCell(iCell)) = 0.0_RKIND
-         fluxRediZTop(:, maxLevelCell(iCell) + 1) = 0.0_RKIND
-         redi_term3_topOfCell(:, 1:minLevelCell(iCell), iCell) = 0.0_RKIND
-         do k = minLevelCell(iCell) + 1, maxLevelCell(iCell)
-            redi_term3_topOfCell(:, k, iCell) = fluxRediZTop(:, k) * invAreaCell
-         end do
-         redi_term3_topOfCell(:, maxLevelCell(iCell) + 1, iCell) = 0.0_RKIND
+         if ( use_Redi_diag_terms) then
+            ! impose no-flux boundary conditions at top and bottom of column
+            fluxRediZTop(:, 1:minLevelCell(iCell)) = 0.0_RKIND
+            fluxRediZTop(:, maxLevelCell(iCell) + 1) = 0.0_RKIND
+            redi_term3_topOfCell(:, 1:minLevelCell(iCell), iCell) = 0.0_RKIND
+            do k = minLevelCell(iCell) + 1, maxLevelCell(iCell)
+               redi_term3_topOfCell(:, k, iCell) = fluxRediZTop(:, k) * invAreaCell
+            end do
+            redi_term3_topOfCell(:, maxLevelCell(iCell) + 1, iCell) = 0.0_RKIND
 
-         do k = minLevelCell(iCell), maxLevelCell(iCell)
-            do iTr = 1, nTracers
-               ! Add tendency for Term 3: d/dz ( h S grad phi) = ( S grad phi) fluxes
-               ! 2.0 in next line is because a dot product on a C-grid
-               ! requires a factor of 1/2 to average to the cell center.
-               flux_term3 = rediKappaCell(iCell)*2.0_RKIND* &
-                            (RediKappaSfcTaper(k, iCell)*RediKappaScaling(k, iCell)*fluxRediZTop(iTr, k)  &
-                            * invAreaCell - RediKappaSfcTaper(k + 1, iCell)*RediKappaScaling(k + 1, iCell) &
-                            *fluxRediZTop(iTr, k + 1) * invAreaCell)
-
-               redi_term3(iTr, k, iCell) = redi_term3(iTr, k, iCell) + flux_term3
+            do k = minLevelCell(iCell), maxLevelCell(iCell)
+               do iTr = 1, nTracers
+                  ! Add tendency for Term 3: d/dz ( h S grad phi) = ( S grad phi) fluxes
+                  ! 2.0 in next line is because a dot product on a C-grid
+                  ! requires a factor of 1/2 to average to the cell center.
+                  flux_term3 = rediKappaCell(iCell)*2.0_RKIND* &
+                               (RediKappaSfcTaper(k, iCell)*RediKappaScaling(k, iCell)*fluxRediZTop(iTr, k)  &
+                               * invAreaCell - RediKappaSfcTaper(k + 1, iCell)*RediKappaScaling(k + 1, iCell) &
+                               *fluxRediZTop(iTr, k + 1) * invAreaCell)
+
+                  redi_term3(iTr, k, iCell) = redi_term3(iTr, k, iCell) + flux_term3
+               end do
             end do
-         end do
+         end if
+
       end do ! iCell
       !$omp end do
       !$omp end parallel
@@ -416,22 +433,29 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThicknessEdg
          do k = minLevelCell(iCell), maxLevelCell(iCell)
             do iTr = 1, ntracers
                tend(iTr, k, iCell) = tend(iTr, k, iCell) + redi_term1(iTr, k, iCell)
-               tend(iTr, k, iCell) = tend(iTr, k, iCell) + rediKappaCell(iCell)*2.0_RKIND* &
-                                     (RediKappaSfcTaper(k, iCell)*RediKappaScaling(k, iCell)* &
-                                      redi_term3_topOfCell(iTr, k, iCell) - RediKappaSfcTaper(k + 1, iCell) &
-                                      *RediKappaScaling(k + 1, iCell)*redi_term3_topOfCell(iTr, k + 1, iCell))
             end do
          end do
 
-         do i = 1, nEdgesOnCell(iCell)
-            if (cellsOnCell(i, iCell) .eq. nCellsP1) cycle
-            iEdge = edgesOnCell(i, iCell)
-            do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
+         if (maxLevelCell(iCell) .ge. config_Redi_min_layers_diag_terms) then
+            do k = minLevelCell(iCell), maxLevelCell(iCell)
                do iTr = 1, ntracers
-                  tend(iTr, k, iCell) = tend(iTr, k, iCell) + edgeSignOnCell(i, iCell)*invAreaCell*redi_term2_edge(iTr, k, iEdge)
+                  tend(iTr, k, iCell) = tend(iTr, k, iCell) + rediKappaCell(iCell)*2.0_RKIND* &
+                                        (RediKappaSfcTaper(k, iCell)*RediKappaScaling(k, iCell)* &
+                                         redi_term3_topOfCell(iTr, k, iCell) - RediKappaSfcTaper(k + 1, iCell) &
+                                         *RediKappaScaling(k + 1, iCell)*redi_term3_topOfCell(iTr, k + 1, iCell))
                end do
             end do
-         end do
+
+            do i = 1, nEdgesOnCell(iCell)
+               if (cellsOnCell(i, iCell) .eq. nCellsP1) cycle
+               iEdge = edgesOnCell(i, iCell)
+               do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
+                  do iTr = 1, ntracers
+                     tend(iTr, k, iCell) = tend(iTr, k, iCell) + edgeSignOnCell(i, iCell)*invAreaCell*redi_term2_edge(iTr, k, iEdge)
+                  end do
+               end do
+            end do
+         end if
       end do ! iCell
       !$omp end do
       !$omp end parallel