Merge branch 'vanroekel/ocean/change-redi-slopes-and-limiter' into ne…

…xt (PR #5947)

Improve slope calculation and fix slope limiting for redi mixing

In current MPAS code, slopes are altered directly, e.g. when S > Scrit
it is set to zero. Griffies et al 1998 (Appendix C) equates this to
slope clipping. The correct approach is to taper the redi kappa
coefficient for each individual flux instead. The cross terms ("term2"
and "term 3") end up being identical as they are d/dz(kappa S grad phi)
so tapering kappa or S are identical. However for the vertical it is
d/dz(kappa |S|^2 d phi / dz) so tapering S is different from tapering
kappa. Also term 1 (laplacian diffusion) was not being tapered with
slopes creating a cross isopycnal flux in steeply sloping regions.
These issues are fixed here.

In addition a new slope calculation from @dengwirda is introduced to
prevent the slopes from getting very large when grad phi or d/dz(phi)
get small. This function limits the slope to be between -1 and 1.

[NML]
[CC]

components/mpas-ocean/bld/build-namelist

@@ -557,11 +557,9 @@ add_default($nl, 'config_Redi_constant_kappa');
 add_default($nl, 'config_Redi_maximum_slope');
 add_default($nl, 'config_Redi_use_slope_taper');
 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_use_quasi_monotone_limiter');
 add_default($nl, 'config_Redi_quasi_monotone_safety_factor');
+add_default($nl, 'config_Redi_limit_term1');
 add_default($nl, 'config_Redi_min_layers_diag_terms');
 add_default($nl, 'config_Redi_horizontal_taper');
 add_default($nl, 'config_Redi_horizontal_ramp_min');

components/mpas-ocean/bld/build-namelist-section

@@ -96,11 +96,9 @@ add_default($nl, 'config_Redi_constant_kappa');
 add_default($nl, 'config_Redi_maximum_slope');
 add_default($nl, 'config_Redi_use_slope_taper');
 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_use_quasi_monotone_limiter');
 add_default($nl, 'config_Redi_quasi_monotone_safety_factor');
+add_default($nl, 'config_Redi_limit_term1');
 add_default($nl, 'config_Redi_min_layers_diag_terms');
 add_default($nl, 'config_Redi_horizontal_taper');
 add_default($nl, 'config_Redi_horizontal_ramp_min');

components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml

@@ -146,11 +146,9 @@
 <config_Redi_maximum_slope>0.01</config_Redi_maximum_slope>
 <config_Redi_use_slope_taper>.true.</config_Redi_use_slope_taper>
 <config_Redi_use_surface_taper>.true.</config_Redi_use_surface_taper>
-<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_use_quasi_monotone_limiter>.true.</config_Redi_use_quasi_monotone_limiter>
 <config_Redi_quasi_monotone_safety_factor>0.9</config_Redi_quasi_monotone_safety_factor>
+<config_Redi_limit_term1>.true.</config_Redi_limit_term1>
 <config_Redi_min_layers_diag_terms>6</config_Redi_min_layers_diag_terms>
 <config_Redi_min_layers_diag_terms ocn_grid="ARRM10to60E2r1">15</config_Redi_min_layers_diag_terms>
 <config_Redi_horizontal_taper>'ramp'</config_Redi_horizontal_taper>

components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml

@@ -423,23 +423,7 @@ Valid values: .true. or .false.
 Default: Defined in namelist_defaults.xml
 </entry>
 
-<entry id="config_Redi_N2_based_taper_enable" type="logical"
-	category="Redi_isopycnal_mixing" group="Redi_isopycnal_mixing">
-If true apply the N2 limiting of Danabasoglu and Marshall 2007
-
-Valid values: .true. or .false.
-Default: Defined in namelist_defaults.xml
-</entry>
-
-<entry id="config_Redi_N2_based_taper_min" type="real"
-	category="Redi_isopycnal_mixing" group="Redi_isopycnal_mixing">
-minimum taper value for the N2 limiting of Danabasoglu and Marshall 2007
-
-Valid values: greater than zero and less than 1
-Default: Defined in namelist_defaults.xml
-</entry>
-
-<entry id="config_Redi_N2_based_taper_limit_term1" type="logical"
+<entry id="config_Redi_limit_term1" type="logical"
 	category="Redi_isopycnal_mixing" group="Redi_isopycnal_mixing">
 If true, the N2 limiting is applied to the horizontal diffusion term
 

components/mpas-ocean/src/Registry.xml

@@ -314,15 +314,7 @@
 					description="If true, Redi slope is tapered near sfc based on Large et al 1997"
 					possible_values=".true. or .false."
 		/>
-		<nml_option name="config_Redi_N2_based_taper_enable" type="logical" default_value=".true."
-					description="If true apply the N2 limiting of Danabasoglu and Marshall 2007"
-					possible_values=".true. or .false."
-		/>
-		<nml_option name="config_Redi_N2_based_taper_min" type="real" default_value="0.1"
-					description="minimum taper value for the N2 limiting of Danabasoglu and Marshall 2007"
-					possible_values="greater than zero and less than 1"
-		/>
-		<nml_option name="config_Redi_N2_based_taper_limit_term1" type="logical" default_value=".true."
+		<nml_option name="config_Redi_limit_term1" type="logical" default_value=".true."
 					description="If true, the N2 limiting is applied to the horizontal diffusion term"
 					possible_values=".true. or .false."
 		/>
@@ -3016,10 +3008,6 @@
 			description="spatially and depth varying GM kappa.  The scaling is based on the Brunt Vaisala Frequency relative to a maximum value below the mixed layer, follows from Danabasoglu and Marshall 2007.  If config_GM_closure is not set to N2_dependent the scaling value is set to 1 everywhere."
 			packages="gm"
 		/>
-		<var name="RediKappaScaling" type="real" dimensions="nVertLevelsP1 nCells Time" units="1"
-			description="Scaling coefficient for GM kappa. Varies from 0 to 1.  The scaling is based on the Brunt Vaisala Frequency relative to a maximum value below the mixed layer, follows from Danabasoglu and Marshall 2007.  If config_Redi_N2_based_taper_enable = .false. the scaling is set to 1 everywhere."
-			packages="gm"
-		/>
 		<var name="RediKappaSfcTaper" type="real" dimensions="nVertLevelsP1 nCells Time" units="1"
 			description="Scaling coefficient for Redi kappa. Varies from 0 to 1."
 			packages="gm"
@@ -3159,6 +3147,14 @@
 			 description="Magnitude of slope of isopycnal surface, using triad through this cell and edge, angled up. Uses expansion of equation of state."
 			 packages="forwardMode;analysisMode"
 		/>
+		<var name="limiterUp" type="real" dimensions="nVertLevels TWO nEdges Time" units="non-dimensional" default_value="0.0"
+			 description="Magnitude of slope of isopycnal surface, using triad through this cell and edge, angled up. Uses expansion of equation of state."
+			 packages="forwardMode;analysisMode"
+		/>
+		<var name="limiterDown" type="real" dimensions="nVertLevels TWO nEdges Time" units="non-dimensional" default_value="0.0"
+			 description="Magnitude of slope of isopycnal surface, using triad through this cell and edge, angled up. Uses expansion of equation of state."
+			 packages="forwardMode;analysisMode"
+		/>
 		<var name="k33" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2 s^-1"
 			 description="The (3,3) entry of the Redi diffusion tensor. Added to the model vertical diffusion. Defined at the top of cell k"
 			 packages="forwardMode;analysisMode"

components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F

@@ -121,13 +121,12 @@ module ocn_diagnostics_variables
    real (kind=RKIND), dimension(:,:), pointer :: velocityZonal
 
    real(kind=RKIND), dimension(:,:), pointer :: gmKappaScaling
-   real(kind=RKIND), dimension(:,:), pointer :: RediKappaScaling
    real(kind=RKIND), dimension(:,:), pointer :: RediKappaSfcTaper
    real(kind=RKIND), dimension(:,:), pointer :: k33
    real(kind=RKIND), dimension(:,:), pointer :: gradDensityEddy
    real(kind=RKIND), dimension(:), pointer   :: gmBolusKappa
    real(kind=RKIND), dimension(:), pointer   :: cGMphaseSpeed
-   real(kind=RKIND), dimension(:,:,:), pointer :: slopeTriadUp, slopeTriadDown
+   real(kind=RKIND), dimension(:,:,:), pointer :: limiterUp, limiterDown, slopeTriadUp, slopeTriadDown
 
    integer, dimension(:), pointer :: indMLD
    real(kind=RKIND), dimension(:), pointer :: dThreshMLD
@@ -265,8 +264,6 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
       if (config_use_GM.or.config_use_Redi) then
          call mpas_pool_get_array(diagnosticsPool, 'RediKappa', &
                    RediKappa)
-         call mpas_pool_get_array(diagnosticsPool, 'RediKappaScaling', &
-                   RediKappaScaling)
          call mpas_pool_get_array(diagnosticsPool, 'RediKappaSfcTaper', &
                    RediKappaSfcTaper)
          call mpas_pool_get_array(diagnosticsPool, 'rediLimiterCount', &
@@ -420,6 +417,10 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
                    slopeTriadUp)
          call mpas_pool_get_array(diagnosticsPool, 'slopeTriadDown', &
                    slopeTriadDown)
+         call mpas_pool_get_array(diagnosticsPool, 'limiterUp', &
+                   limiterUp)
+         call mpas_pool_get_array(diagnosticsPool, 'limiterDown', &
+                   limiterDown)
 
          ! GM-related fields
          call mpas_pool_get_array(diagnosticsPool, 'cGMphaseSpeed', &
@@ -719,6 +720,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
          !$acc                   SSHGradient,                     &
          !$acc                   circulation,                     &
          !$acc                   slopeTriadUp,                    &
+         !$acc                   limiterUp,                       &
          !$acc                   gradSSHZonal,                    &
          !$acc                   cGMphaseSpeed,                   &
          !$acc                   velocityZonal,                   &
@@ -729,6 +731,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
          !$acc                   layerThickEdgeFlux,              &
          !$acc                   layerThickEdgeMean,              &
          !$acc                   slopeTriadDown,                  &
+         !$acc                   limiterDown,                     &
          !$acc                   normRelVortEdge,                 &
          !$acc                   surfaceVelocity,                 &
          !$acc                   vertVelocityTop,                 &
@@ -792,7 +795,6 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
          !$acc                   RediKappa,         &
          !$acc                   gmBolusKappa,      &
          !$acc                   gmKappaScaling,    &
-         !$acc                   RediKappaScaling,  &
          !$acc                   rediLimiterCount,  &
          !$acc                   RediKappaSfcTaper, &
          !$acc                   RediHorizontalTaper,  &
@@ -961,6 +963,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
          !$acc                   SSHGradient,                     &
          !$acc                   circulation,                     &
          !$acc                   slopeTriadUp,                    &
+         !$acc                   limiterUp,                       &
+         !$acc                   limiterDown,                     &
          !$acc                   gradSSHZonal,                    &
          !$acc                   cGMphaseSpeed,                   &
          !$acc                   velocityZonal,                   &
@@ -1034,7 +1038,6 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
          !$acc                   RediKappa,         &
          !$acc                   gmBolusKappa,      &
          !$acc                   gmKappaScaling,    &
-         !$acc                   RediKappaScaling,  &
          !$acc                   rediLimiterCount,  &
          !$acc                   RediKappaSfcTaper, &
          !$acc                   RediHorizontalTaper,  &
@@ -1161,6 +1164,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
                  SSHGradient,                     &
                  circulation,                     &
                  slopeTriadUp,                    &
+                 limiterUp,                       &
+                 limiterDown,                     &
                  gradSSHZonal,                    &
                  cGMphaseSpeed,                   &
                  velocityZonal,                   &
@@ -1230,7 +1235,6 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
          nullify(RediKappa,         &
                  gmBolusKappa,      &
                  gmKappaScaling,    &
-                 RediKappaScaling,  &
                  rediLimiterCount,  &
                  RediKappaSfcTaper, &
                  RediHorizontalTaper, &

components/mpas-ocean/src/shared/mpas_ocn_gm.F

@@ -62,6 +62,22 @@ module ocn_gm
 
 contains
 
+   function safe_slope(num, den) result(slope)
+
+      real(kind=RKIND), intent(in) :: num, den
+      real(kind=RKIND) :: slope
+
+      ! compute isopycnal slopes safely wrt. finite precision
+      ! small slope assumption, so don't want |s| > 1
+
+      if (abs(den) > abs(num)) then
+         slope = num / den
+      else
+         slope = sign(1.0_RKIND, den * num)
+      end if
+
+   end function
+
 !***********************************************************************
 !
 !  routine ocn_GM_compute_Bolus_velocity
@@ -214,7 +230,6 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
 
       !$omp do schedule(runtime)
       do iCell = 1, nCells
-         RediKappaScaling(:, iCell) = 1.0_RKIND
          RediKappaSfcTaper(:, iCell) = 1.0_RKIND
       end do
       !$omp end do
@@ -229,30 +244,6 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
          allocate(dSdzTop(nVertLevels + 1))
          allocate(k33Norm(nVertLevels + 1))
 
-         if ( config_Redi_N2_based_taper_enable ) then
-
-            !$omp parallel
-            !$omp do schedule(runtime) private(maxLocation, k, BruntVaisalaFreqTopEdge, maxN)
-            do iCell = 1, nCells
-               k = min(maxLevelCell(iCell) - 1, max(minLevelCell(iCell), indMLD(iCell)))
-               maxN = max(BruntVaisalaFreqTop(k, iCell), 0.0_RKIND)
-               do while (BruntVaisalaFreqTop(k + 1, iCell) > maxN .and. k < maxLevelCell(iCell) - 1)
-                 k = k + 1
-                 maxN = max(maxN,max(BruntVaisalaFreqTop(k, iCell), 0.0_RKIND))
-               enddo
-
-               maxLocation = k
-               do k = maxLocation, maxLevelCell(iCell)
-                  BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTop(k, iCell), 0.0_RKIND)
-                  RediKappaScaling(k, iCell) = min(max(config_Redi_N2_based_taper_min,                  &
-                                                       BruntVaisalaFreqTopEdge/(maxN + 1.0E-10_RKIND)), &
-                                                       1.0_RKIND)
-               end do
-            end do
-            !$omp end do
-            !$omp end parallel
-         end if
-
          if (config_Redi_use_surface_taper) then
             !$omp parallel
             !$omp do schedule (runtime) private(zMLD, k)
@@ -318,17 +309,17 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
                          *dcEdgeInv
                   drhoDx = drhoDT*dTdx + drhoDS*dSdx
 
+
+                  sfcTaperUp = drhoDT*dTdzTop(k) + drhoDS*dSdzTop(k)
+                  sfcTaperDown = drhoDT*dTdzTop(k+1) + drhoDS*dSdzTop(k+1)
+
                   ! Always compute *Up on the top cell and *Down on the bottom
                   ! cell, even though they are never used. This avoids an if
                   ! statement or separate computation for top and bottom.
                   slopeTriadUp(k, iCellSelf, iEdge) = &
-                     -drhoDx/ &
-                     (drhoDT*dTdzTop(k) &
-                      + drhoDS*dSdzTop(k) + 1E-15_RKIND)
+                     safe_slope( -drhoDx, sfcTaperUp )
                   slopeTriadDown(k, iCellSelf, iEdge) = &
-                     -drhoDx/ &
-                     (drhoDT*dTdzTop(k + 1) &
-                      + drhoDS*dSdzTop(k + 1) + 1E-15_RKIND)
+                      safe_slope( -drhoDx, sfcTaperDown )
 
                   ! set taper of slope ('F' function from Danabasoglu and McWilliams 95)
                   if (abs(slopeTriadDown(k, iCellSelf, iEdge)) > 0.6_RKIND*config_redi_maximum_slope) then
@@ -363,30 +354,32 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
                   sfcTaperUp = 1.0_RKIND + sfcTaperFactor*(sfcTaper - 1.0_RKIND)
                   sfcTaperDown = 1.0_RKIND + sfcTaperFactor*(sfcTaper - 1.0_RKIND)
 
-                  slopeTriadUp(k, iCellSelf, iEdge) = &
-                     slopeTaperUp*sfcTaperUp*slopeTriadUp(k, iCellSelf, iEdge)
-                  slopeTriadDown(k, iCellSelf, iEdge) = &
-                     slopeTaperDown*sfcTaperDown*slopeTriadDown(k, iCellSelf, iEdge)
-
-                  ! Griffies 1998 eqn 34
+!                  ! Griffies 1998 eqn 34
                   if (k > minLevelCell(iCell)) then
                      k33(k, iCell) = k33(k, iCell) +  &
-                                     areaEdge*dzTop(k)*slopeTriadUp(k, iCellSelf, iEdge)**2
+                                     slopeTaperUp*sfcTaperUp*areaEdge*dzTop(k)*slopeTriadUp(k, iCellSelf, iEdge)**2
                      k33Norm(k) = k33Norm(k) + areaEdge*dzTop(k)
                   end if
 
-                  !Taper Redi by tapering the slopes
+!                  !Taper Redi by tapering the slopes
                   k33(k + 1, iCell) = k33(k + 1, iCell) +  &
-                                      areaEdge*dzTop(k + 1)*slopeTriadDown(k, iCellSelf, iEdge)**2
+                                      slopeTaperDown*sfcTaperDown*areaEdge*dzTop(k + 1)*slopeTriadDown(k, iCellSelf, iEdge)**2
                   k33Norm(k + 1) = k33Norm(k + 1) + areaEdge*dzTop(k + 1)
 
+                  limiterUp(k,iCellSelf,iEdge) = slopeTaperUp*sfcTaperUp
+                  limiterDown(k,iCellSelf,iEdge) = slopeTaperDown*sfcTaperDown
+
+                  slopeTriadUp(k, iCellSelf, iEdge) = &
+                     slopeTaperUp*sfcTaperUp*slopeTriadUp(k, iCellSelf, iEdge)
+                  slopeTriadDown(k, iCellSelf, iEdge) = &
+                     slopeTaperDown*sfcTaperDown*slopeTriadDown(k, iCellSelf, iEdge)
 
                end do ! maxLevelEdgeTop(iEdge)
             end do ! nEdgesOnCell(iCell)
 
-            ! Normalize k33
+!            ! Normalize k33
             do k = minLevelCell(iCell)+1, maxLevelCell(iCell)
-               k33(k, iCell) = k33(k, iCell)/k33Norm(k)*RediKappaScaling(k, iCell)
+               k33(k, iCell) = k33(k, iCell)/k33Norm(k)
             end do
             k33(1:minLevelCell(iCell), iCell) = 0.0_RKIND
             k33(maxLevelCell(iCell) + 1, iCell) = 0.0_RKIND
@@ -904,7 +897,7 @@ subroutine ocn_GM_init(domain, err)!{{{
       real(kind=RKIND), dimension(:), pointer :: dcEdge
 
       integer :: iEdge
-      integer, pointer :: nEdges
+      integer, pointer :: nCells, nEdges
       integer, pointer :: nVertLevels
       real(kind=RKIND), pointer :: sphere_radius
       real(kind=RKIND), dimension(:), pointer   :: latEdge, fEdge
@@ -918,10 +911,12 @@ subroutine ocn_GM_init(domain, err)!{{{
       do while (associated(block))
          call mpas_pool_get_subpool(block%structs, 'mesh', meshPool)
          call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
+         call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
          call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges)
          call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
          call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)
          call mpas_pool_get_array(meshPool, 'fEdge', fEdge)
+
          if (config_Redi_use_slope_taper) then
             slopeTaperFactor = 1.0_RKIND
          else

components/mpas-ocean/src/shared/mpas_ocn_tendency.F

@@ -1141,8 +1141,9 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, &
          call ocn_tracer_hmix_tend(meshPool, layerThickEdgeMean, &
                                    layerThickness, zMid, tracerGroup, &
                                    RediKappa, slopeTriadUp, &
-                                   slopeTriadDown, dt, isActiveTracer, &
-                                   RediKappaScaling, rediLimiterCount, &
+                                   slopeTriadDown, limiterUp, &
+                                   limiterDown, dt, isActiveTracer, &
+                                   rediLimiterCount, &
                                    tracerGroupTend, err)
 
          if (computeBudgets) then