From 831eda7175ce72f293ac05b625189f083e32dedd Mon Sep 17 00:00:00 2001 From: dengwirda Date: Thu, 27 Apr 2023 18:28:12 -0600 Subject: [PATCH 1/9] Add a quasi-monotone limiter for isopycnal diffusion - Add a new config_Redi_use_quasi_monotone_limiter option for the isopycnal diffusion operator to better control non-monotone departure of tracer values. - Limiter expands on existing strategy of disabling cross-fluxes, but uses a more comprehensive test based on min/max values over stencil of neighbouring cells in the adjacent layers. --- components/mpas-ocean/src/Registry.xml | 4 ++ .../src/shared/mpas_ocn_tracer_hmix_redi.F | 55 +++++++++++++++++-- 2 files changed, 55 insertions(+), 4 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 2091c51429e2..e5e2af1a0034 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -326,6 +326,10 @@ description="If true, the N2 limiting is applied to the horizontal diffusion term" possible_values=".true. or .false." /> + maximumBnd(iTr, k, iCell))) then do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) @@ -463,6 +508,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea deallocate (redi_term3) deallocate (redi_term2_edge) deallocate (redi_term3_topOfCell) + deallocate (minimumBnd) + deallocate (maximumBnd) call mpas_timer_stop("tracer redi") From cb595d200286624a5e34dc10647158f99bef1a06 Mon Sep 17 00:00:00 2001 From: dengwirda Date: Fri, 28 Apr 2023 03:15:15 -0600 Subject: [PATCH 2/9] Add a quasi-monotone limiter for isopycnal diffusion - Scale cross-fluxes toward zero proportional to the degree of non-monotonicity. Enables softer limiter compared to default on/off switch. Control ramp with quasi_monotone_safety_factor. --- components/mpas-ocean/src/Registry.xml | 6 ++- .../src/shared/mpas_ocn_tracer_hmix_redi.F | 37 ++++++++++++------- 2 files changed, 29 insertions(+), 14 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index e5e2af1a0034..d958461b77ad 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -327,9 +327,13 @@ possible_values=".true. or .false." /> + maximumBnd(iTr, k, iCell))) then + (use_quasi_monotone .and. tempTracer < minimumBnd(iTr, k, iCell)) .or. & + (use_quasi_monotone .and. tempTracer > maximumBnd(iTr, k, iCell))) then + + ! scale cross-fluxes down relative to size of bnds violation + over = max(tempTracer - maximumBnd(iTr, k, iCell), & + minimumBnd(iTr, k, iCell) - tempTracer) + + scal = min(1.0_RKIND, max(0.0_RKIND, & + limiter_safety_factor * ( & + 1.0_RKIND - over / abs(tempTracer - tracers(iTr, k, iCell))))) do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) - redi_term2_edge(iTr, k, iEdge) = 0.0_RKIND + redi_term2_edge(iTr, k, iEdge) = & + redi_term2_edge(iTr, k, iEdge) * scal end do - redi_term3_topOfCell(iTr, k, iCell) = 0.0_RKIND - redi_term3_topOfCell(iTr, k + 1, iCell) = 0.0_RKIND + redi_term3_topOfCell(iTr, k, iCell) = & + redi_term3_topOfCell(iTr, k, iCell) * scal + redi_term3_topOfCell(iTr, k + 1, iCell) = & + redi_term3_topOfCell(iTr, k + 1, iCell) * scal if (isActiveTracer) then rediLimiterCount(k, iCell) = rediLimiterCount(k, iCell) + 1 - endif - endif + end if + end if end do end do end do ! iCell From 02a36bd5b589838141518049271b33f5b38ac70e Mon Sep 17 00:00:00 2001 From: dengwirda Date: Sun, 30 Apr 2023 18:42:56 -0600 Subject: [PATCH 3/9] Add a quasi-monotone limiter for isopycnal diffusion - Ensure fluxes are only scaled once by accumulating scaling on cells and then interpolating scaling to edges/levels. - Change default safety_factor to 1.0, based on QU30km spin-ups showing fully monotone behaviour. - Add local variables to OMP directives. --- components/mpas-ocean/src/Registry.xml | 2 +- .../src/shared/mpas_ocn_tracer_hmix_redi.F | 69 ++++++++++++++----- 2 files changed, 51 insertions(+), 20 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index d958461b77ad..a21d1de90a15 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -330,7 +330,7 @@ description="If true, fluxes are reduced to prevent tracers from violating monotonicity. Cross-term fluxes are scaled toward zero to prevent tracers from under/overshooting the min/max values in adjacent cells and layers" possible_values=".true. or .false." /> - 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 index d0396cdeef95..8e71698ad88f 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F @@ -57,6 +57,7 @@ module ocn_tracer_hmix_Redi real(kind=RKIND), dimension(:), allocatable :: rediKappaCell real(kind=RKIND), dimension(:, :, :), allocatable :: redi_term2_edge, redi_term3_topOfCell real(kind=RKIND), dimension(:, :, :), allocatable :: redi_term1, redi_term2, redi_term3 + real(kind=RKIND), dimension(:, :, :), allocatable :: redi_flux_scale real(kind=RKIND), dimension(:, :, :), allocatable :: minimumBnd, maximumBnd !*********************************************************************** @@ -140,7 +141,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea real(kind=RKIND) :: tempTracer, invAreaCell real(kind=RKIND) :: flux, flux_term2, flux_term3, dTracerDx, coef real(kind=RKIND) :: r_tmp, tracer_turb_flux, kappaRediEdgeVal - real(kind=RKIND) :: over, scal, limiter_safety_factor + real(kind=RKIND) :: over, diff, scal, scalUp, scalDn + real(kind=RKIND) :: limiter_safety_factor real(kind=RKIND), dimension(:), pointer :: areaCell, dvEdge, dcEdge real(kind=RKIND), dimension(:), allocatable :: minimumVal real(kind=RKIND), dimension(:, :), allocatable :: fluxRediZTop @@ -180,6 +182,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea allocate (redi_term3_topOfCell(nTracers, nVertLevels + 1, nCells)) allocate (minimumBnd(nTracers, nVertLevels, nCells)) allocate (maximumBnd(nTracers, nVertLevels, nCells)) + allocate (redi_flux_scale(nTracers, nVertLevels, nCells)) ! RediKappa changes every time step if either: ! config_Redi_horizontal_taper == 'RossbyRadius' @@ -213,7 +216,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea !$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, use_Redi_diag_terms) + !$omp dTracerDx, use_Redi_diag_terms, kUp, kDn) do iCell = 1, nCells if (maxLevelCell(iCell) .ge. config_Redi_min_layers_diag_terms) then use_Redi_diag_terms = .true. @@ -221,9 +224,10 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea use_Redi_diag_terms = .false. endif + minimumBnd(:, :, iCell) = 0.0_RKIND + maximumBnd(:, :, iCell) = 0.0_RKIND if (use_quasi_monotone) then - minimumBnd(:, :, iCell) = 0.0_RKIND - maximumBnd(:, :, iCell) = 0.0_RKIND + ! init. min/max bnds in column do k = minLevelCell(iCell), maxLevelCell(iCell) kUp = max(minLevelCell(iCell), k - 1) kDn = min(maxLevelCell(iCell), k + 1) @@ -259,6 +263,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea jCell = cellsOnCell(i, iCell) if (use_quasi_monotone) then + ! expand min/max bnds on edge neighbours do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) kUp = max(minLevelEdgeBot(iEdge), k - 1) kDn = min(maxLevelEdgeTop(iEdge), k + 1) @@ -437,8 +442,10 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea endif !$omp parallel - !$omp do schedule(runtime) private(k, iTr, tempTracer, iEdge) + !$omp do schedule(runtime) & + !$omp private(k, iTr, tempTracer, iEdge, over, scal, diff) do iCell = 1, nCells + redi_flux_scale(:, :, iCell) = 1.0_RKIND do k = minLevelCell(iCell), maxLevelCell(iCell) do iTr = 1, ntracers tempTracer = tracers(iTr, k, iCell) + dt*(redi_term1(iTr, k, iCell) + redi_term2(iTr, k, iCell) + & @@ -452,20 +459,18 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea over = max(tempTracer - maximumBnd(iTr, k, iCell), & minimumBnd(iTr, k, iCell) - tempTracer) + diff = abs(tempTracer - tracers(iTr, k, iCell)) + scal = min(1.0_RKIND, max(0.0_RKIND, & limiter_safety_factor * ( & - 1.0_RKIND - over / abs(tempTracer - tracers(iTr, k, iCell))))) + 1.0_RKIND - over / (diff + 1.0E-16_RKIND)))) - do i = 1, nEdgesOnCell(iCell) - iEdge = edgesOnCell(i, iCell) - redi_term2_edge(iTr, k, iEdge) = & - redi_term2_edge(iTr, k, iEdge) * scal - end do + ! scale to zero if not using monotone limiter + if (.not. use_quasi_monotone) scal = 0.0_RKIND + + ! cell flux scalings, to be interp. onto edge/levs next pass + redi_flux_scale(iTr, k, iCell) = scal - redi_term3_topOfCell(iTr, k, iCell) = & - redi_term3_topOfCell(iTr, k, iCell) * scal - redi_term3_topOfCell(iTr, k + 1, iCell) = & - redi_term3_topOfCell(iTr, k + 1, iCell) * scal if (isActiveTracer) then rediLimiterCount(k, iCell) = rediLimiterCount(k, iCell) + 1 end if @@ -478,7 +483,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea !now go back and reapply all tendencies !$omp parallel - !$omp do schedule(runtime) private(invAreaCell, k, iTr, i, iEdge) + !$omp do schedule(runtime) & + !$omp private(invAreaCell, k, kUp, kDn, iTr, i, iEdge, jCell, scal, scalUp, scalDn) do iCell = 1, nCells invAreaCell = 1.0_RKIND/areaCell(iCell) do k = minLevelCell(iCell), maxLevelCell(iCell) @@ -489,20 +495,44 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea if (maxLevelCell(iCell) .ge. config_Redi_min_layers_diag_terms) then do k = minLevelCell(iCell), maxLevelCell(iCell) + kUp = max(minLevelCell(iCell), k - 1) + kDn = min(maxLevelCell(iCell), k + 1) do iTr = 1, ntracers + ! harmonic-mean; not strictly monotone, but smooth + scalUp = 2.0_RKIND * & + redi_flux_scale(iTr, kUp, iCell) * & + redi_flux_scale(iTr, k, iCell) / ( & + redi_flux_scale(iTr, kUp, iCell) + & + redi_flux_scale(iTr, k, iCell) + 1.0E-16_RKIND) + + scalDn = 2.0_RKIND * & + redi_flux_scale(iTr, kDn, iCell) * & + redi_flux_scale(iTr, k, iCell) / ( & + redi_flux_scale(iTr, kDn, iCell) + & + redi_flux_scale(iTr, k, iCell) + 1.0E-16_RKIND) + tend(iTr, k, iCell) = tend(iTr, k, iCell) + rediKappaCell(iCell)*2.0_RKIND* & (RediKappaScaling(k, iCell)* & - redi_term3_topOfCell(iTr, k, iCell) - & - RediKappaScaling(k + 1, iCell)*redi_term3_topOfCell(iTr, k + 1, iCell)) + scalUp*redi_term3_topOfCell(iTr, k, iCell) - & + RediKappaScaling(k + 1, iCell)* & + scalDn*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) + jCell = cellsOnCell(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) + scal = 2.0_RKIND * & + redi_flux_scale(iTr, k, iCell) * & + redi_flux_scale(iTr, k, jCell) / ( & + redi_flux_scale(iTr, k, iCell) + & + redi_flux_scale(iTr, k, jCell) + 1.0E-16_RKIND) + + tend(iTr, k, iCell) = tend(iTr, k, iCell) + & + scal*edgeSignOnCell(i, iCell)*invAreaCell*redi_term2_edge(iTr, k, iEdge) end do end do end do @@ -521,6 +551,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea deallocate (redi_term3_topOfCell) deallocate (minimumBnd) deallocate (maximumBnd) + deallocate (redi_flux_scale) call mpas_timer_stop("tracer redi") From dfe41e3e42015e07394061f176c7f72b0157606d Mon Sep 17 00:00:00 2001 From: dengwirda Date: Fri, 5 May 2023 05:28:29 -0600 Subject: [PATCH 4/9] Fix errors in isopycnal diffusion flux terms - Update use of triad slope in Redi term-3. - Fix edge-to-cell remapping weights for Redi term-3. - Switch to module-level config. variables. --- .../src/shared/mpas_ocn_tracer_hmix_redi.F | 90 +++++++++++-------- 1 file changed, 52 insertions(+), 38 deletions(-) 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 index 8e71698ad88f..60340a28624b 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F @@ -54,6 +54,8 @@ module ocn_tracer_hmix_Redi !-------------------------------------------------------------------- real(kind=RKIND) :: term1TaperFactor + real(kind=RKIND) :: limiter_safety_factor + logical :: use_quasi_monotone real(kind=RKIND), dimension(:), allocatable :: rediKappaCell real(kind=RKIND), dimension(:, :, :), allocatable :: redi_term2_edge, redi_term3_topOfCell real(kind=RKIND), dimension(:, :, :), allocatable :: redi_term1, redi_term2, redi_term3 @@ -131,18 +133,18 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea integer :: iCell, jCell, iEdge, cell1, cell2, iCellSelf integer :: i, k, kUp, kDn, iTr, nTracers, nCells, nEdges, nCellsP1 - logical :: use_Redi_diag_terms, use_quasi_monotone + logical :: use_Redi_diag_terms integer, pointer :: nVertLevels integer, dimension(:), pointer :: nCellsArray, nEdgesArray - integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop, nEdgesOnCell, minLevelCell, maxLevelCell + integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop, & + nEdgesOnCell, minLevelCell, maxLevelCell integer, dimension(:, :), pointer :: cellsOnEdge, edgesOnCell, edgeSignOnCell, cellsOnCell real(kind=RKIND) :: tempTracer, invAreaCell real(kind=RKIND) :: flux, flux_term2, flux_term3, dTracerDx, coef real(kind=RKIND) :: r_tmp, tracer_turb_flux, kappaRediEdgeVal real(kind=RKIND) :: over, diff, scal, scalUp, scalDn - real(kind=RKIND) :: limiter_safety_factor real(kind=RKIND), dimension(:), pointer :: areaCell, dvEdge, dcEdge real(kind=RKIND), dimension(:), allocatable :: minimumVal real(kind=RKIND), dimension(:, :), allocatable :: fluxRediZTop @@ -207,11 +209,10 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea nCells = nCellsArray(2) nCellsP1 = nCellsArray(size(nCellsArray)) + 1 - use_quasi_monotone = config_Redi_use_quasi_monotone_limiter - limiter_safety_factor = config_Redi_quasi_monotone_safety_factor - - ! Term 1: this is the "standard" horizontal del2 term, but with RediKappa coefficient. - ! \kappa_2 \nabla \phi on edge + ! Term 1: the "standard" horizontal del^2 term, but with RediKappa coefficient. + ! Term 2: a cross-flux: kappa div( h S d/dz psi ) + ! Term 3: a cross-flux: kappa d/dz div( S grad psi ) + ! Term 4: the vertical flux (done in vmix) d/dz kappa S^2 d/dz psi !$omp parallel !$omp do schedule(runtime) & !$omp private(fluxRediZTop, invAreaCell, i, iEdge, cell1, cell2, iCellSelf, r_tmp, coef, k, & @@ -252,6 +253,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea ! Check if neighboring cell exists if (cellsOnCell(i, iCell) .eq. nCellsP1) cycle iEdge = edgesOnCell(i, iCell) + jCell = cellsOnCell(i, iCell) cell1 = cellsOnEdge(1, iEdge) cell2 = cellsOnEdge(2, iEdge) if (cell1 == iCell) then @@ -260,8 +262,6 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea iCellSelf = 2 endif - jCell = cellsOnCell(i, iCell) - if (use_quasi_monotone) then ! expand min/max bnds on edge neighbours do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) @@ -289,6 +289,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2)) do iTr = 1, nTracers + ! term 1: div( h grad psi ) tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1) flux = layerThickEdgeMean(k, iEdge)*tracer_turb_flux*r_tmp*RediKappa(iEdge) redi_term1(iTr, k, iCell) = redi_term1(iTr, k, iCell) - (1.0_RKIND + term1TaperFactor* & @@ -296,7 +297,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea if (.not.use_Redi_diag_terms) cycle - ! Term 2: div( h S dphi/dz) + ! term 2: div( h S d/dz psi ) flux_term2 = coef*kappaRediEdgeVal*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* & 0.25_RKIND* & (slopeTriadDown(k, 1, iEdge)*(tracers(iTr, k, cell1) - tracers(iTr, k + 1, cell1)) & @@ -316,6 +317,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea redi_term2(iTr, k, iCell) = redi_term2(iTr, k, iCell) - edgeSignOnCell(i, iCell)*flux_term2*invAreaCell redi_term2_edge(iTr, k, iEdge) = -flux_term2 + + ! term 3: zero at surface end do do k = minLevelEdgeBot(iEdge) + 1, maxLevelEdgeTop(iEdge) - 1 @@ -324,17 +327,15 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2)) do iTr = 1, nTracers - ! \kappa_2 \nabla \phi on edge + ! term 1: div( h grad psi ) tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1) - - ! div(h \kappa_2 \nabla \phi) at cell center flux = layerThickEdgeMean(k, iEdge)*tracer_turb_flux*r_tmp*RediKappa(iEdge) - 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 d/dz psi ) flux_term2 = coef*RediKappa(iEdge)*kappaRediEdgeVal*layerThickEdgeMean(k, iEdge)* & 0.25_RKIND* & (slopeTriadUp(k, 1, iEdge)*(tracers(iTr, k - 1, cell1) - tracers(iTr, k, cell1)) & @@ -349,10 +350,11 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea invAreaCell redi_term2_edge(iTr, k, iEdge) = -flux_term2 - dTracerDx = (tracers(iTr, k, cell2) - tracers(iTr, k, cell1))*dvEdge(iEdge)*0.25_RKIND - fluxRediZTop(iTr, k) = fluxRediZTop(iTr, k) + & - 0.5_RKIND*(slopeTriadUp(k, iCellSelf, iEdge) + slopeTriadDown(k - 1, iCellSelf, iEdge)) & - *dTracerDx + ! term 3 : d/dz div( S grad psi ) + ! includes weights for edge-to-cell remap + fluxRediZTop(iTr, k) = fluxRediZTop(iTr, k) + 0.125_RKIND*dvEdge(iEdge)*( & + slopeTriadDown(k - 1, iCellSelf, iEdge)*(tracers(iTr, k - 1, cell2) - tracers(iTr, k - 1, cell1)) + & + slopeTriadUp(k, iCellSelf, iEdge)*(tracers(iTr, k, cell2) - tracers(iTr, k, cell1))) end do end do @@ -361,6 +363,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2)) do iTr = 1, nTracers + ! term 1: div( h grad psi ) tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1) flux = layerThickEdgeMean(k, iEdge)*tracer_turb_flux*r_tmp*RediKappa(iEdge) redi_term1(iTr, k, iCell) = redi_term1(iTr, k, iCell) - (1.0_RKIND + term1TaperFactor* & @@ -368,6 +371,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea if (.not.use_Redi_diag_terms) cycle + ! term 2: div( h S d/dz psi ) ! For bottom layer, only use triads pointing up: flux_term2 = coef*kappaRediEdgeVal*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* & 0.25_RKIND* & @@ -390,32 +394,32 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea *invAreaCell redi_term2_edge(iTr, k, iEdge) = -flux_term2 - dTracerDx = (tracers(iTr, k, cell2) - tracers(iTr, k, cell1))*dvEdge(iEdge)*0.25_RKIND - fluxRediZTop(iTr, k) = fluxRediZTop(iTr, k) + & - 0.5_RKIND*(slopeTriadUp(k, iCellSelf, iEdge) + slopeTriadDown(k - 1, iCellSelf, iEdge)) & - *dTracerDx + ! term 3: d/dz div( S grad psi ) + ! includes weights for edge-to-cell remap + fluxRediZTop(iTr, k) = fluxRediZTop(iTr, k) + 0.125_RKIND*dvEdge(iEdge)*( & + slopeTriadDown(k - 1, iCellSelf, iEdge)*(tracers(iTr, k - 1, cell2) - tracers(iTr, k - 1, cell1)) + & + slopeTriadUp(k, iCellSelf, iEdge)*(tracers(iTr, k, cell2) - tracers(iTr, k, cell1))) end do end do ! nEdgesOnCell(iCell) 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 + fluxRediZTop(:, maxLevelCell(iCell) + 1) = 0.0_RKIND redi_term3_topOfCell(:, maxLevelCell(iCell) + 1, iCell) = 0.0_RKIND + ! take d/dz and remap div( S grad psi ) from edge to cells 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* & - (RediKappaScaling(k, iCell)*fluxRediZTop(iTr, k) & - * invAreaCell - RediKappaScaling(k + 1, iCell) & - *fluxRediZTop(iTr, k + 1) * invAreaCell) + flux_term3 = rediKappaCell(iCell) * ( & + RediKappaScaling(k, iCell)* & + redi_term3_topOfCell(iTr, k, iCell) - & + RediKappaScaling(k + 1, iCell)* & + redi_term3_topOfCell(iTr, k + 1, iCell)) redi_term3(iTr, k, iCell) = redi_term3(iTr, k, iCell) + flux_term3 end do @@ -448,8 +452,9 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea redi_flux_scale(:, :, iCell) = 1.0_RKIND do k = minLevelCell(iCell), maxLevelCell(iCell) do iTr = 1, ntracers - tempTracer = tracers(iTr, k, iCell) + dt*(redi_term1(iTr, k, iCell) + redi_term2(iTr, k, iCell) + & - redi_term3(iTr, k, iCell))/layerThickness(k, iCell) + tempTracer = tracers(iTr, k, iCell) + dt*( & + redi_term1(iTr, k, iCell) + redi_term2(iTr, k, iCell) + & + redi_term3(iTr, k, iCell))/layerThickness(k, iCell) if (tempTracer < minimumVal(iTr) .or. & (use_quasi_monotone .and. tempTracer < minimumBnd(iTr, k, iCell)) .or. & @@ -494,6 +499,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea end do if (maxLevelCell(iCell) .ge. config_Redi_min_layers_diag_terms) then + ! term 3: d/dz div( S grad psi ) do k = minLevelCell(iCell), maxLevelCell(iCell) kUp = max(minLevelCell(iCell), k - 1) kDn = min(maxLevelCell(iCell), k + 1) @@ -511,14 +517,16 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea redi_flux_scale(iTr, kDn, iCell) + & redi_flux_scale(iTr, k, iCell) + 1.0E-16_RKIND) - tend(iTr, k, iCell) = tend(iTr, k, iCell) + rediKappaCell(iCell)*2.0_RKIND* & - (RediKappaScaling(k, iCell)* & - scalUp*redi_term3_topOfCell(iTr, k, iCell) - & - RediKappaScaling(k + 1, iCell)* & - scalDn*redi_term3_topOfCell(iTr, k + 1, iCell)) + tend(iTr, k, iCell) = tend(iTr, k, iCell) + & + rediKappaCell(iCell) * ( & + RediKappaScaling(k, iCell)* & + scalUp*redi_term3_topOfCell(iTr, k, iCell) - & + RediKappaScaling(k + 1, iCell)* & + scalDn*redi_term3_topOfCell(iTr, k + 1, iCell)) end do end do + ! term 2: div( h S d/dz psi ) do i = 1, nEdgesOnCell(iCell) if (cellsOnCell(i, iCell) .eq. nCellsP1) cycle iEdge = edgesOnCell(i, iCell) @@ -634,6 +642,12 @@ subroutine ocn_tracer_hmix_Redi_init(domain, err)!{{{ call mpas_dmpar_finalize(domain%dminfo) end if + ! quasi bounds-preserving limiter + use_quasi_monotone = config_Redi_use_quasi_monotone_limiter + + ! safety factor for limiter: 0 <= fac <= 1 + limiter_safety_factor = config_Redi_quasi_monotone_safety_factor + ! Initialize horizontal taper if (config_Redi_horizontal_taper == 'none' .or. & config_Redi_horizontal_taper == 'RossbyRadius') then From b8c3c99d25237fb377b44efd5fd35c5ec46cd0cd Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Wed, 20 Sep 2023 10:24:25 -0500 Subject: [PATCH 5/9] initialize redi_flux_scale --- .../src/shared/mpas_ocn_tracer_hmix_redi.F | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) 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 index 60340a28624b..ba4a50cf6438 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F @@ -192,6 +192,18 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea ! in that case, RediKappa is updated every time step in mpas_ocn_gm.F ! For static RediKappa, it is set on init and remains unchanged. + !$omp parallel + !$omp do schedule(runtime) private(k,iTr) + do iCell=1,nCells + do k=1,nVertLevels + do iTr = 1,nTracers + redi_flux_scale(iTr,k,iCell) = 1.0_RKIND + end do + end do + end do + !$omp end do + !$omp end parallel + nCells = nCellsArray(2) !$omp parallel !$omp do schedule(runtime) private(i, iEdge) @@ -449,7 +461,6 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea !$omp do schedule(runtime) & !$omp private(k, iTr, tempTracer, iEdge, over, scal, diff) do iCell = 1, nCells - redi_flux_scale(:, :, iCell) = 1.0_RKIND do k = minLevelCell(iCell), maxLevelCell(iCell) do iTr = 1, ntracers tempTracer = tracers(iTr, k, iCell) + dt*( & From 45bdbab584f36f85e611e34610e4822831edd3b1 Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Thu, 21 Sep 2023 17:55:53 -0500 Subject: [PATCH 6/9] Make bld files consistent with Registry changes --- components/mpas-ocean/bld/build-namelist | 2 ++ components/mpas-ocean/bld/build-namelist-section | 2 ++ .../namelist_files/namelist_defaults_mpaso.xml | 2 ++ .../namelist_files/namelist_definition_mpaso.xml | 16 ++++++++++++++++ 4 files changed, 22 insertions(+) diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index c7030bd5fc5a..2dfa6a9d1e92 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -560,6 +560,8 @@ 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_min_layers_diag_terms'); add_default($nl, 'config_Redi_horizontal_taper'); add_default($nl, 'config_Redi_horizontal_ramp_min'); diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index a47c8cc74656..98cfe66feea2 100644 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -99,6 +99,8 @@ 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_min_layers_diag_terms'); add_default($nl, 'config_Redi_horizontal_taper'); add_default($nl, 'config_Redi_horizontal_ramp_min'); diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index 08964af1cff9..0ca032f74e2b 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -144,6 +144,8 @@ .true. 0.1 .true. +.false. +1.0 6 15 'ramp' diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml index 352303431ec7..f3da8bfe4c2d 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -447,6 +447,22 @@ Valid values: .true. or .false. Default: Defined in namelist_defaults.xml + +If true, fluxes are reduced to prevent tracers from violating monotonicity. Cross-term fluxes are scaled toward zero to prevent tracers from under/overshooting the min/max values in adjacent cells and layers + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +A safety factor applied to flux scaling when monotonicity is violated. Smaller values scale fluxes toward zero more aggressively. + +Valid values: A value between 0 and 1 +Default: Defined in namelist_defaults.xml + + 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. From aca6daa3787cb23fb1640ab14edafd0c43f25603 Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Tue, 26 Sep 2023 09:09:29 -0500 Subject: [PATCH 7/9] Set E3SM default to config_Redi_use_quasi_monotone_limiter=true --- .../mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index 0ca032f74e2b..3b64fce70785 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -144,7 +144,7 @@ .true. 0.1 .true. -.false. +.true. 1.0 6 15 From f22e3f3e58a8a1ba62f27fff44809b0e56356a62 Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Tue, 26 Sep 2023 12:45:34 -0500 Subject: [PATCH 8/9] Set default as config_Redi_quasi_monotone_safety_factor=0.9 --- .../mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml | 2 +- components/mpas-ocean/src/Registry.xml | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index 3b64fce70785..d86a29620d5e 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -145,7 +145,7 @@ 0.1 .true. .true. -1.0 +0.9 6 15 'ramp' diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index a21d1de90a15..fdfa0cb2dea8 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -326,11 +326,11 @@ description="If true, the N2 limiting is applied to the horizontal diffusion term" possible_values=".true. or .false." /> - - From ec67ee79227883e68b2b8b469f43e2312e708d81 Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Wed, 25 Oct 2023 13:57:10 -0500 Subject: [PATCH 9/9] Add missing variable to omp directive --- components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 index ba4a50cf6438..227b27592ac4 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F @@ -229,7 +229,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea !$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, use_Redi_diag_terms, kUp, kDn) + !$omp dTracerDx, use_Redi_diag_terms, jCell, kUp, kDn) do iCell = 1, nCells if (maxLevelCell(iCell) .ge. config_Redi_min_layers_diag_terms) then use_Redi_diag_terms = .true.