From b031624fa7cd95aa18784c1b500a89f77deac128 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Tue, 5 Sep 2023 23:10:56 -0500 Subject: [PATCH 1/6] Changes redi tapering and slopes This changes tapering on the redi terms to be fully on the slopes themselves and fixes a bug where k33 slopes are tapered by the square and not linearly. A 'safe slope' calculation is also introduced to bound slopes from -1 to 1. --- components/mpas-ocean/src/Registry.xml | 12 ++- .../shared/mpas_ocn_diagnostics_variables.F | 12 ++- .../mpas-ocean/src/shared/mpas_ocn_gm.F | 92 ++++++++++--------- .../mpas-ocean/src/shared/mpas_ocn_tendency.F | 5 +- .../src/shared/mpas_ocn_tracer_hmix.F | 9 +- .../src/shared/mpas_ocn_tracer_hmix_redi.F | 53 +++++------ 6 files changed, 101 insertions(+), 82 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 2091c51429e2..c9d4802932d1 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -3008,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" /> - + + 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 @@ -229,30 +245,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) @@ -271,8 +263,14 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & nCells = nCellsArray(3) !$omp parallel !$omp do schedule(runtime) & - !$omp private(invAreaCell, k33Norm, k, dzTop, dTdzTop, dSdzTop, i, iEdge, cell1, cell2, & - !$omp iCellSelf, dcEdgeInv, areaEdge, drhoDT, drhoDS, dTdx, dSdx, drhoDx, & +! !$omp private(invAreaCell, k33Norm, k, dzTop, dTdzTop, dSdzTop, i, iEdge, cell1, cell2, & +! !$omp iCellSelf, dcEdgeInv, areaEdge, drhoDT, drhoDS, dTdx, dSdx, drhoDx, & +! !$omp slopeTaperUp, slopeTaperDown, sfcTaperUp, sfcTaperDown, sfcTaper) + !$omp private(invAreaCell, dzTop, dTdzTop, dSdzTop, & + !$omp k, i, iEdge, cell1, cell2, & + !$omp iCellSelf, dcEdgeInv, drhoDTTop, drhoDTBot, & + !$omp drhoDSTop, drhoDSBot, dTdxTop, dTdxBot, dSdxTop, dSdxBot, & + !$omp drhoDxTop, drhoDxBot, & !$omp slopeTaperUp, slopeTaperDown, sfcTaperUp, sfcTaperDown, sfcTaper) do iCell = 1, nCells invAreaCell = 1.0_RKIND/areaCell(iCell) @@ -318,17 +316,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 +361,34 @@ 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 + RediKappaScaling(1:minLevelCell(iCell),iCell) = 0.0_RKIND + RediKappaScaling(maxLevelCell(iCell):nVertLevels,iCell) = 0.0_RKIND 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 +906,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 +920,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 diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F index a806c57412ef..2719e1526930 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F +++ b/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 diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix.F index c0e2cdadf16a..7b10b012f88c 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix.F @@ -80,8 +80,8 @@ module ocn_tracer_hmix !----------------------------------------------------------------------- subroutine ocn_tracer_hmix_tend(meshPool, layerThickEdgeMean, layerThickness, zMid, tracers, & - RediKappa, slopeTriadUp, slopeTriadDown, dt, isActiveTracer, & - RediKappaScaling, rediLimiterCount, tend, err)!{{{ + RediKappa, slopeTriadUp, slopeTriadDown, limiterUp, limiterDown, & + dt, isActiveTracer, rediLimiterCount, tend, err)!{{{ !----------------------------------------------------------------- @@ -98,8 +98,6 @@ subroutine ocn_tracer_hmix_tend(meshPool, layerThickEdgeMean, layerThickness, zM layerThickEdgeMean, &!< Input: mean thickness at edge layerThickness, &!< Input: thickness at centers zMid !< Input: Z coordinate at the center of a cell - real (kind=RKIND), dimension(:,:), intent(in), optional :: & - RediKappaScaling !< Input: vertical structure of GM/Redi limiter real (kind=RKIND), dimension(:,:,:), intent(in) :: & slopeTriadUp, slopeTriadDown, & @@ -154,7 +152,8 @@ subroutine ocn_tracer_hmix_tend(meshPool, layerThickEdgeMean, layerThickness, zM call ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMean, zMid, tracers, & RediKappa, dt, isActiveTracer, & slopeTriadUp, slopeTriadDown, & - RediKappaScaling, rediLimiterCount, tend, err1) + limiterUp, limiterDown, rediLimiterCount, & + tend, err1) endif err = ior(err1, err2) call mpas_timer_stop("tracer hmix") 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 83dd0986a239..602fcc2db9ad 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 @@ -75,7 +75,7 @@ module ocn_tracer_hmix_Redi subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMean, zMid, tracers, & RediKappa, dt, isActiveTracer, slopeTriadUp, slopeTriadDown, & - RediKappaScaling, rediLimiterCount, tend, err)!{{{ + limiterUp, limiterDown, rediLimiterCount, tend, err)!{{{ !----------------------------------------------------------------- ! @@ -91,11 +91,10 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea real(kind=RKIND), dimension(:, :), intent(in) :: & layerThickEdgeMean, &!< Input: mean thickness at edge zMid, &!< Input: Z coordinate at the center of a cell - RediKappaScaling, & layerThickness real(kind=RKIND), dimension(:, :, :), intent(in) :: & - slopeTriadUp, slopeTriadDown, & + slopeTriadUp, slopeTriadDown, limiterUp, limiterDown, & tracers !< Input: tracer quantities real(kind=RKIND), intent(in) :: dt @@ -237,10 +236,12 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea 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)) + + kappaRediEdgeVal = 0.25_RKIND*(limiterUp(k,1,iEdge) + limiterUp(k,2,iEdge) + & + limiterDown(k,1,iEdge) + limiterDown(k,2,iEdge)) 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* & @@ -248,8 +249,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea if (.not.use_Redi_diag_terms) cycle - ! Term 2: div( h S dphi/dz) - flux_term2 = coef*kappaRediEdgeVal*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* & + ! term 2: div( h S d/dz psi ) + flux_term2 = coef*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* & 0.25_RKIND* & (slopeTriadDown(k, 1, iEdge)*(tracers(iTr, k, cell1) - tracers(iTr, k + 1, cell1)) & /(zMid(k, cell1) - zMid(k + 1, cell1)) & @@ -257,37 +258,37 @@ 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(k, iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(k, 1, iEdge)* & + layerThickEdgeMean(k, iEdge)*0.25_RKIND*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(k, iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadUp(k, 2, iEdge)* & + layerThickEdgeMean(k, iEdge)*0.25_RKIND*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 redi_term2_edge(iTr, k, iEdge) = -flux_term2 + + ! term 3: zero at surface end do do k = minLevelEdgeBot(iEdge) + 1, maxLevelEdgeTop(iEdge) - 1 - kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + & - RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2)) + kappaRediEdgeVal = 0.25_RKIND*(limiterUp(k,1,iEdge) + limiterUp(k,2,iEdge) + & + limiterDown(k,1,iEdge) + limiterDown(k,2,iEdge)) 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 - flux_term2 = coef*RediKappa(iEdge)*kappaRediEdgeVal*layerThickEdgeMean(k, iEdge)* & + ! term 2: div( h S d/dz psi ) + flux_term2 = coef*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* & 0.25_RKIND* & (slopeTriadUp(k, 1, iEdge)*(tracers(iTr, k - 1, cell1) - tracers(iTr, k, cell1)) & /(zMid(k - 1, cell1) - zMid(k, cell1)) & @@ -309,10 +310,11 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea end do k = maxLevelEdgeTop(iEdge) - kappaRediEdgeVal = 0.25_RKIND*(RediKappaScaling(k, cell1) + RediKappaScaling(k, cell2) + & - RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2)) + kappaRediEdgeVal = 0.25_RKIND*(limiterUp(k,1,iEdge) + limiterUp(k,2,iEdge) + & + limiterDown(k,1,iEdge) + limiterDown(k,2,iEdge)) 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* & @@ -321,7 +323,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea if (.not.use_Redi_diag_terms) cycle ! For bottom layer, only use triads pointing up: - flux_term2 = coef*kappaRediEdgeVal*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* & + flux_term2 = coef*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* & 0.25_RKIND* & (slopeTriadUp(k, 1, iEdge)*(tracers(iTr, k - 1, cell1) - tracers(iTr, k, cell1)) & /(zMid(k - 1, cell1) - zMid(k, cell1)) & @@ -330,12 +332,12 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea if (maxLevelCell(cell1) > maxLevelEdgeTop(iEdge)) then flux_term2 = flux_term2 + coef*RediKappa(iEdge)* & - layerThickEdgeMean(k, iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadDown(k, 1, iEdge)* & + layerThickEdgeMean(k, iEdge)*0.25_RKIND*slopeTriadDown(k, 1, iEdge)* & (tracers(iTr, k, cell1) - tracers(iTr, k + 1, cell1))/(zMid(k, cell1) - zMid(k + 1, cell1)) endif if (maxLevelCell(cell2) > maxLevelEdgeTop(iEdge)) then flux_term2 = flux_term2 + coef*RediKappa(iEdge)* & - layerThickEdgeMean(k, iEdge)*0.25_RKIND*kappaRediEdgeVal*slopeTriadDown(k, 2, iEdge)* & + layerThickEdgeMean(k, iEdge)*0.25_RKIND*slopeTriadDown(k, 2, iEdge)* & (tracers(iTr, k, cell2) - tracers(iTr, k + 1, cell2))/(zMid(k, cell2) - zMid(k + 1, cell2)) endif redi_term2(iTr, k, iCell) = redi_term2(iTr, k, iCell) - edgeSignOnCell(i, iCell)*flux_term2 & @@ -365,8 +367,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea ! 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) & + * invAreaCell - & *fluxRediZTop(iTr, k + 1) * invAreaCell) redi_term3(iTr, k, iCell) = redi_term3(iTr, k, iCell) + flux_term3 @@ -435,9 +437,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea do k = minLevelCell(iCell), maxLevelCell(iCell) do iTr = 1, ntracers 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)) + (redi_term3_topOfCell(iTr, k, iCell) - & + redi_term3_topOfCell(iTr, k + 1, iCell)) end do end do From 135dea2b4e281c25579cf0640b261b06dbcc21ec Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Tue, 5 Sep 2023 23:55:48 -0500 Subject: [PATCH 2/6] fixes build errors --- components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix.F | 2 +- components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix.F index 7b10b012f88c..a40cc58933cd 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix.F @@ -100,7 +100,7 @@ subroutine ocn_tracer_hmix_tend(meshPool, layerThickEdgeMean, layerThickness, zM zMid !< Input: Z coordinate at the center of a cell real (kind=RKIND), dimension(:,:,:), intent(in) :: & - slopeTriadUp, slopeTriadDown, & + slopeTriadUp, slopeTriadDown, limiterUp, limiterDown, & tracers !< Input: tracer quantities real (kind=RKIND), intent(in) :: dt 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 602fcc2db9ad..bd93fbe9ee6c 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 @@ -369,7 +369,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea flux_term3 = rediKappaCell(iCell)*2.0_RKIND* & (fluxRediZTop(iTr, k) & * invAreaCell - & - *fluxRediZTop(iTr, k + 1) * invAreaCell) + fluxRediZTop(iTr, k + 1) * invAreaCell) redi_term3(iTr, k, iCell) = redi_term3(iTr, k, iCell) + flux_term3 end do From fbb9ce21b72399b68ecfe6a6bbd140587de263e5 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Wed, 6 Sep 2023 13:25:01 -0500 Subject: [PATCH 3/6] changes redi namelist options the redi n2 based limiter options are removed from the namelist and the term1 limiter name is changed to reference tapering not just the n2 based tapering for clarity --- components/mpas-ocean/src/Registry.xml | 10 +--------- .../mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F | 2 +- 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index c9d4802932d1..2616766e514d 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/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." /> - - - 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 bd93fbe9ee6c..7a2b8f58b673 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 @@ -506,7 +506,7 @@ subroutine ocn_tracer_hmix_Redi_init(domain, err)!{{{ if (.not. config_use_Redi) return - if (config_Redi_N2_based_taper_limit_term1) then + if (config_Redi_limit_term1) then term1TaperFactor = 1.0_RKIND else term1TaperFactor = 0.0_RKIND From a0fdeca71121d5685c2ff58c9fa2384f2bf489c1 Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Thu, 14 Sep 2023 11:00:37 -0500 Subject: [PATCH 4/6] Fixed OpenMP --- components/mpas-ocean/src/shared/mpas_ocn_gm.F | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_gm.F b/components/mpas-ocean/src/shared/mpas_ocn_gm.F index c27b87f892f7..fe8d9353bbe3 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_gm.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_gm.F @@ -263,14 +263,8 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & nCells = nCellsArray(3) !$omp parallel !$omp do schedule(runtime) & -! !$omp private(invAreaCell, k33Norm, k, dzTop, dTdzTop, dSdzTop, i, iEdge, cell1, cell2, & -! !$omp iCellSelf, dcEdgeInv, areaEdge, drhoDT, drhoDS, dTdx, dSdx, drhoDx, & -! !$omp slopeTaperUp, slopeTaperDown, sfcTaperUp, sfcTaperDown, sfcTaper) - !$omp private(invAreaCell, dzTop, dTdzTop, dSdzTop, & - !$omp k, i, iEdge, cell1, cell2, & - !$omp iCellSelf, dcEdgeInv, drhoDTTop, drhoDTBot, & - !$omp drhoDSTop, drhoDSBot, dTdxTop, dTdxBot, dSdxTop, dSdxBot, & - !$omp drhoDxTop, drhoDxBot, & + !$omp private(invAreaCell, k33Norm, k, dzTop, dTdzTop, dSdzTop, i, iEdge, cell1, cell2, & + !$omp iCellSelf, dcEdgeInv, areaEdge, drhoDT, drhoDS, dTdx, dSdx, drhoDx, & !$omp slopeTaperUp, slopeTaperDown, sfcTaperUp, sfcTaperDown, sfcTaper) do iCell = 1, nCells invAreaCell = 1.0_RKIND/areaCell(iCell) From fac4c9cb983b44af4dd4f103b83b95571b1671c0 Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Fri, 22 Sep 2023 13:24:06 -0500 Subject: [PATCH 5/6] Update bld files to match Registry changes --- components/mpas-ocean/bld/build-namelist | 4 +--- .../mpas-ocean/bld/build-namelist-section | 4 +--- .../namelist_files/namelist_defaults_mpaso.xml | 4 +--- .../namelist_definition_mpaso.xml | 18 +----------------- 4 files changed, 4 insertions(+), 26 deletions(-) diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index c7030bd5fc5a..9cc8f4dd1414 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -557,9 +557,7 @@ 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_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'); diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index a47c8cc74656..bf9d47935e66 100644 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -96,9 +96,7 @@ 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_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'); 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..3fbbfa0c096d 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -141,9 +141,7 @@ 0.01 .true. .true. -.true. -0.1 -.true. +.true. 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..868f9c63ddb5 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/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 - -If true apply the N2 limiting of Danabasoglu and Marshall 2007 - -Valid values: .true. or .false. -Default: Defined in namelist_defaults.xml - - - -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 - - - If true, the N2 limiting is applied to the horizontal diffusion term From 534209c13c1b923769f0e6c283ca055553955b57 Mon Sep 17 00:00:00 2001 From: Mark Petersen Date: Mon, 25 Sep 2023 15:02:26 -0500 Subject: [PATCH 6/6] Remove RediKappaScaling creation and initialization --- .../mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F | 6 ------ components/mpas-ocean/src/shared/mpas_ocn_gm.F | 3 --- 2 files changed, 9 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F index d618647930d1..6a64f86c2a56 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -121,7 +121,6 @@ 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 @@ -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', & @@ -798,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, & @@ -1042,7 +1038,6 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc RediKappa, & !$acc gmBolusKappa, & !$acc gmKappaScaling, & - !$acc RediKappaScaling, & !$acc rediLimiterCount, & !$acc RediKappaSfcTaper, & !$acc RediHorizontalTaper, & @@ -1240,7 +1235,6 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ nullify(RediKappa, & gmBolusKappa, & gmKappaScaling, & - RediKappaScaling, & rediLimiterCount, & RediKappaSfcTaper, & RediHorizontalTaper, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_gm.F b/components/mpas-ocean/src/shared/mpas_ocn_gm.F index fe8d9353bbe3..6f4442103fa7 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_gm.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_gm.F @@ -230,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 @@ -379,8 +378,6 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, & end do ! nEdgesOnCell(iCell) ! ! Normalize k33 - RediKappaScaling(1:minLevelCell(iCell),iCell) = 0.0_RKIND - RediKappaScaling(maxLevelCell(iCell):nVertLevels,iCell) = 0.0_RKIND do k = minLevelCell(iCell)+1, maxLevelCell(iCell) k33(k, iCell) = k33(k, iCell)/k33Norm(k) end do