From 328f9a4064dc3168e3972a27e7e8cea2ee91c94d Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Tue, 8 Sep 2020 13:38:48 -0700 Subject: [PATCH 1/8] Implement ocn_mesh in tracer_advection_std --- .../shared/mpas_ocn_tracer_advection.F | 6 +- .../shared/mpas_ocn_tracer_advection_std.F | 57 ++++++------------- 2 files changed, 18 insertions(+), 45 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_advection.F b/src/core_ocean/shared/mpas_ocn_tracer_advection.F index d50bac8ba6..0dbc5d1a84 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_advection.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_advection.F @@ -149,10 +149,8 @@ subroutine ocn_tracer_advection_tend(tracers, normalThicknessFlux, w, & diagnosticsPool, computeBudgets) else - call ocn_tracer_advection_std_tend(tracers, advCoefs, advCoefs3rd, & - nAdvCellsForEdge, advCellsForEdge, normalThicknessFlux, w, layerThickness, & - layerThickness, dt, meshPool, tend, maxLevelCell, maxLevelEdgeTop, & - highOrderAdvectionMask, edgeSignOnCell) + call ocn_tracer_advection_std_tend(tracers, normalThicknessFlux, w, layerThickness, & + layerThickness, dt, tend) endif call mpas_timer_stop("tracer adv") diff --git a/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F b/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F index a34a88e1f4..73c67e0bc4 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F @@ -26,6 +26,7 @@ module ocn_tracer_advection_std use mpas_pool_routines use mpas_io_units use mpas_threading + use ocn_mesh use mpas_tracer_advection_helpers @@ -63,37 +64,23 @@ module ocn_tracer_advection_std !> Both horizontal and vertical. ! !----------------------------------------------------------------------- - subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdvCellsForEdge, advCellsForEdge, &!{{{ - normalThicknessFlux, w, layerThickness, verticalCellSize, dt, meshPool, & - tend, maxLevelCell, maxLevelEdgeTop, & - highOrderAdvectionMask, edgeSignOnCell) + subroutine ocn_tracer_advection_std_tend(tracers, normalThicknessFlux, w, layerThickness, & + verticalCellSize, dt, tend)!{{{ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input: current tracer values - real (kind=RKIND), dimension(:,:), intent(in) :: adv_coefs !< Input: Advection coefficients for 2nd order advection - real (kind=RKIND), dimension(:,:), intent(in) :: adv_coefs_3rd !< Input: Advection coeffs for blending in 3rd/4th order - integer, dimension(:), intent(in) :: nAdvCellsForEdge !< Input: Number of advection cells for each edge - integer, dimension(:,:), intent(in) :: advCellsForEdge !< Input: List of advection cells for each edge real (kind=RKIND), dimension(:,:), intent(in) :: normalThicknessFlux !< Input: Thichness weighted velocitiy real (kind=RKIND), dimension(:,:), intent(in) :: w !< Input: Vertical velocity real (kind=RKIND), dimension(:,:), intent(in) :: layerThickness !< Input: Thickness real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of a cell real (kind=RKIND), intent(in) :: dt !< Input: Timestep - type (mpas_pool_type), intent(in) :: meshPool !< Input: Mesh information real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer tendency - integer, dimension(:), pointer :: maxLevelCell !< Input: Index to max level at cell center - integer, dimension(:), pointer :: maxLevelEdgeTop !< Input: Index to max level at edge with non-land cells on both sides - integer, dimension(:,:), pointer :: highOrderAdvectionMask !< Input: Mask for high order advection - integer, dimension(:, :), pointer :: edgeSignOnCell !< Input: Sign for flux from edge on each cell. integer :: i, iCell, iEdge, k, iTracer, cell1, cell2 integer :: nVertLevels, num_tracers - integer, pointer :: nCells, nEdges, nCellsSolve, maxEdges - integer, dimension(:), pointer :: nEdgesOnCell - integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, edgesOnCell real (kind=RKIND) :: tracer_weight, invAreaCell1 real (kind=RKIND) :: verticalWeightK, verticalWeightKm1 - real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, verticalDivergenceFactor + real (kind=RKIND), dimension(:), pointer :: verticalDivergenceFactor real (kind=RKIND), dimension(:,:), allocatable :: tracer_cur real (kind=RKIND), dimension(:,:), allocatable :: high_order_horiz_flux @@ -102,34 +89,22 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv real (kind=RKIND), parameter :: eps = 1.e-10_RKIND ! Get dimensions - call mpas_pool_get_dimension(meshPool, 'nCells', nCells) - call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) - call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) - call mpas_pool_get_dimension(meshPool, 'maxEdges', maxEdges) nVertLevels = size(tracers,dim=2) num_tracers = size(tracers,dim=1) - ! Initialize pointers - call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) - call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) - call mpas_pool_get_array(meshPool, 'areaCell', areaCell) - call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) - allocate(verticalDivergenceFactor(nVertLevels)) verticalDivergenceFactor = 1.0_RKIND - allocate(high_order_horiz_flux(nVertLevels, nEdges), & - tracer_cur(nVertLevels, nCells), & - high_order_vert_flux(nVertLevels, nCells)) + allocate(high_order_horiz_flux(nVertLevels, nEdgesAll), & + tracer_cur(nVertLevels, nCellsAll), & + high_order_vert_flux(nVertLevels, nCellsAll)) ! Loop over tracers. One tracer is advected at a time. It is copied into a temporary array in order to improve locality do iTracer = 1, num_tracers ! Initialize variables for use in this iTracer iteration !$omp parallel !$omp do schedule(runtime) - do iCell = 1, nCells + do iCell = 1, nCellsAll tracer_cur(:, iCell) = tracers(iTracer, :, iCell) high_order_vert_flux(:, iCell) = 0.0_RKIND @@ -137,14 +112,14 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv !$omp end do !$omp do schedule(runtime) - do iEdge = 1, nEdges + do iEdge = 1, nEdgesAll high_order_horiz_flux(:, iEdge) = 0.0_RKIND end do !$omp end do ! Compute the high order vertical flux. Also determine bounds on tracer_cur. !$omp do schedule(runtime) private(k, verticalWeightK, verticalWeightKm1) - do iCell = 1, nCells + do iCell = 1, nCellsAll k = max(1, min(maxLevelCell(iCell), 2)) verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell)) verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell)) @@ -175,13 +150,13 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv ! Compute the high order horizontal flux !$omp do schedule(runtime) private(cell1, cell2, k, tracer_weight, i, iCell) - do iEdge = 1, nEdges + do iEdge = 1, nEdgesAll cell1 = cellsOnEdge(1, iEdge) cell2 = cellsOnEdge(2, iEdge) ! Compute 2nd order fluxes where needed. do k = 1, maxLevelEdgeTop(iEdge) - tracer_weight = iand(highOrderAdvectionMask(k, iEdge)+1, 1) * (dvEdge(iEdge) * 0.5_RKIND) & + tracer_weight = iand(int(highOrderAdvectionMask(k, iEdge)+1), 1) * (dvEdge(iEdge) * 0.5_RKIND) & * normalThicknessFlux(k, iEdge) high_order_horiz_flux(k, iEdge) = high_order_horiz_flux(k, iedge) + tracer_weight & @@ -192,8 +167,8 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv do i = 1, nAdvCellsForEdge(iEdge) iCell = advCellsForEdge(i,iEdge) do k = 1, maxLevelCell(iCell) - tracer_weight = highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef3rdOrder & - * sign(1.0_RKIND,normalThicknessFlux(k,iEdge))*adv_coefs_3rd(i,iEdge)) + tracer_weight = highOrderAdvectionMask(k, iEdge) * (advCoefs(i,iEdge) + coef3rdOrder & + * sign(1.0_RKIND,normalThicknessFlux(k,iEdge))*advCoefs3rd(i,iEdge)) tracer_weight = normalThicknessFlux(k,iEdge)*tracer_weight high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) + tracer_weight * tracer_cur(k,iCell) @@ -204,7 +179,7 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv ! Accumulate the scaled high order horizontal tendencies !$omp do schedule(runtime) private(invAreaCell1, i, iEdge, k) - do iCell = 1, nCells + do iCell = 1, nCellsAll invAreaCell1 = 1.0_RKIND / areaCell(iCell) do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) @@ -218,7 +193,7 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv ! Accumulate the scaled high order vertical tendencies. !$omp do schedule(runtime) private(k) - do iCell = 1, nCellsSolve + do iCell = 1, nCellsOwned do k = 1,maxLevelCell(iCell) tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + verticalDivergenceFactor(k) * (high_order_vert_flux(k+1, iCell) & - high_order_vert_flux(k, iCell)) From 92902ef17ac42be82fe66cbb03aa41a9d7de237e Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Wed, 9 Sep 2020 07:20:17 -0700 Subject: [PATCH 2/8] Use ocn_mesh in tracer adv mono --- .../shared/mpas_ocn_tracer_advection.F | 5 -- .../shared/mpas_ocn_tracer_advection_mono.F | 87 +++++-------------- 2 files changed, 21 insertions(+), 71 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_advection.F b/src/core_ocean/shared/mpas_ocn_tracer_advection.F index 0dbc5d1a84..75d8543d05 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_advection.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_advection.F @@ -141,11 +141,6 @@ subroutine ocn_tracer_advection_tend(tracers, normalThicknessFlux, w, & if (monotonicOn) then call ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & normalThicknessFlux, w, dt, & - advCoefs, advCoefs3rd, & - nAdvCellsForEdge, advCellsForEdge,& - maxLevelCell, maxLevelEdgeTop, & - highOrderAdvectionMask, & - edgeSignOnCell, meshPool, & diagnosticsPool, computeBudgets) else diff --git a/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F b/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F index c9faabbe40..3fd13a57f1 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F @@ -33,6 +33,7 @@ module ocn_tracer_advection_mono use mpas_tracer_advection_helpers use ocn_constants + use ocn_mesh implicit none private @@ -90,13 +91,8 @@ module ocn_tracer_advection_mono !------------------------------------------------------------------------------- subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & - normalThicknessFlux, w, dt, & - advCoefs, advCoefs3rd, & - nAdvCellsForEdge, advCellsForEdge,& - maxLevelCell, maxLevelEdgeTop, & - highOrderAdvectionMask, & - edgeSignOnCell, meshPool, & - diagnosticsPool, computeBudgets)!{{{ + normalThicknessFlux, w, dt, diagnosticsPool, & + computeBudgets)!{{{ !*** Input/Output parameters @@ -115,24 +111,6 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & w !< [in] Vertical velocity real (kind=RKIND), intent(in) :: & dt !< [in] Timestep - real (kind=RKIND), dimension(:,:), intent(in) :: & - advCoefs !< [in] Coefficients for 2nd order advection - real (kind=RKIND), dimension(:,:), intent(in) :: & - advCoefs3rd !< [in] Coeffs for missing 3rd/4th order advection - integer, dimension(:), intent(in) :: & - nAdvCellsForEdge !< [in] Number of advection cells for each edge - integer, dimension(:,:), intent(in) :: & - advCellsForEdge !< [in] List of advection cells for each edge - integer, dimension(:), intent(in) :: & - maxLevelCell !< [in] Index to max level at cell center - integer, dimension(:), intent(in) :: & - maxLevelEdgeTop !< [in] Indx max edge lvl with both non-land cells - integer, dimension(:,:), intent(in) :: & - highOrderAdvectionMask !< [in] Mask for high order advection - integer, dimension(:,:), intent(in) :: & - edgeSignOnCell !< [in] Sign for flux from edge on each cell - type (mpas_pool_type), intent(in) :: & - meshPool !< [in] Mesh information type (mpas_pool_type), intent(in) :: & diagnosticsPool !< [in] pool for traceradvection budget term logical, intent(in) :: & @@ -149,19 +127,6 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & iTracer, &! tracer index numTracers ! total number of tracers - ! pointers for mesh variables retrieved from mesh pool - integer, dimension(:), pointer, contiguous :: & - nCellsArray, &! number of cells for each halo depth - nEdgesArray, &! number of edges for each halo depth - nEdgesOnCell ! number of edges on each cell - integer, dimension(:,:), pointer, contiguous :: & - cellsOnEdge, &! indices of cells on edges - cellsOnCell, &! indices of cells on cells - edgesOnCell ! indices of edges on cells - real (kind=RKIND), dimension(:), pointer, contiguous :: & - dvEdge, &! edge metric - areaCell ! cell area - real (kind=RKIND) :: & signedFactor, &! temp factor including flux sign tracerNew, &! updated tracer @@ -198,16 +163,6 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & #ifdef _ADV_TIMERS call mpas_timer_start('startup') #endif - ! Get mesh data - call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) - call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) - call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) - call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) - call mpas_pool_get_array(meshPool, 'areaCell', areaCell) - call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) - if (computeBudgets) then call mpas_pool_get_array(diagnosticsPool, & 'activeTracerVerticalAdvectionTopFlux', & @@ -226,8 +181,8 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & ! Get dimensions nVertLevels = size(tracers,dim=2) numTracers = size(tracers,dim=1) - nCells = nCellsArray(size(nCellsArray)) - nEdges = nEdgesArray(size(nEdgesArray)) + nCells = nCellsAll + nEdges = nEdgesAll ! allocate temporary arrays allocate(tracerCur (nVertLevels ,nCells+1), & @@ -286,7 +241,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & call mpas_timer_start('cell init') #endif ! reset nCells to include all cells - nCells = nCellsArray( size(nCellsArray) ) + nCells = nCellsAll ! Extract current tracer and change index order to improve locality !$omp parallel @@ -306,7 +261,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & #endif ! set nCells to first halo level - nCells = nCellsArray( 2 ) + nCells = nCellsHalo( 1 ) ! Determine bounds on tracer (tracerMin and tracerMax) from ! surrounding cells for later limiting. @@ -334,7 +289,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & !$omp end parallel ! Need all the edges around the 1 halo cells and owned cells - nEdges = nEdgesArray( 3 ) + nEdges = nEdgesHalo( 2 ) ! Compute the high order horizontal flux @@ -375,7 +330,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & ! Remove low order flux from the high order flux ! Store left over high order flux in highOrderFlx array do k = 1, maxLevelEdgeTop(iEdge) - tracerWeight = iand(highOrderAdvectionMask(k, iEdge)+1, 1) & + tracerWeight = iand(int(highOrderAdvectionMask(k, iEdge)+1), 1) & * (dvEdge(iEdge) * 0.5_RKIND) & * normalThicknessFlux(k, iEdge) @@ -399,7 +354,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & call mpas_timer_start('scale factor build') #endif ! Initialize flux arrays for all cells - nCells = nCellsArray(size(nCellsArray)) + nCells = nCellsAll !$omp parallel !$omp do schedule(runtime) @@ -414,7 +369,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & !$omp end parallel ! Need one halo of cells around owned cells - nCells = nCellsArray( 2 ) + nCells = nCellsHalo( 1 ) !$omp parallel !$omp do schedule(runtime) & @@ -476,7 +431,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & call mpas_timer_start('rescale horiz fluxes') #endif ! Need all of the edges around owned cells - nEdges = nEdgesArray( 2 ) + nEdges = nEdgesHalo( 1 ) ! rescale the high order horizontal fluxes !$omp parallel !$omp do schedule(runtime) private(cell1, cell2, k) @@ -497,7 +452,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & call mpas_timer_start('flux accumulate') #endif - nCells = nCellsArray( 1 ) + nCells = nCellsOwned ! Accumulate the scaled high order vertical tendencies, and the upwind tendencies !$omp parallel !$omp do schedule(runtime) private(invAreaCell1, signedFactor, i, iEdge, k) @@ -537,7 +492,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & !$omp parallel !$omp do schedule(runtime) private(k) - do iEdge = 1,nEdgesArray( 2 ) + do iEdge = 1,nEdgesHalo( 1 ) do k = 1,maxLevelEdgeTop(iEdge) ! Save u*h*T flux on edge for analysis. This variable will be ! divided by h at the end of the time step. @@ -548,7 +503,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & !$omp end do !$omp do schedule(runtime) private(k) - do iCell = 1, nCellsArray( 1 ) + do iCell = 1, nCellsOwned do k = 1, maxLevelCell(iCell) activeTracerHorizontalAdvectionTendency(iTracer,k,iCell) = & workTend(k,iCell) @@ -560,7 +515,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & end if ! computeBudgets if (monotonicityCheck) then - nCells = nCellsArray( 1 ) + nCells = nCellsOwned ! Check tracer values against local min,max to detect ! non-monotone values and write warning if found !$omp parallel @@ -599,7 +554,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & #ifdef _ADV_TIMERS call mpas_timer_start('cell init') #endif - nCells = nCellsArray( size(nCellsArray) ) + nCells = nCellsAll ! Initialize variables for use in this iTracer iteration !$omp parallel !$omp do schedule(runtime) private(k) @@ -622,7 +577,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & #endif ! Need all owned and 1 halo cells - nCells = nCellsArray( 2 ) + nCells = nCellsHalo( 1 ) ! Determine bounds on tracerCur from neighbor values for limiting !$omp parallel @@ -778,7 +733,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & call mpas_timer_start('scale factor build') #endif ! Need one halo of cells around owned cells - nCells = nCellsArray( 2 ) + nCells = nCellsHalo( 1 ) !$omp parallel !$omp do schedule(runtime) & !$omp private(k, tracerMinNew, tracerMaxNew, tracerUpwindNew, scaleFactor) @@ -816,7 +771,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & call mpas_timer_start('flux accumulate') #endif - nCells = nCellsArray( 1 ) + nCells = nCellsOwned ! Accumulate the scaled high order vertical tendencies ! and the upwind tendencies !$omp parallel @@ -870,7 +825,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & end if ! computeBudgets if (monotonicityCheck) then - nCells = nCellsArray( 1 ) + nCells = nCellsOwned ! Check for monotonicity of new tracer value !$omp parallel !$omp do schedule(runtime) private(k, tracerNew) From 270905e95be8705bc2d6ff0dbbce48beb9122c0f Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Wed, 9 Sep 2020 09:33:27 -0700 Subject: [PATCH 3/8] Tracer adv mono cpu optimization --- src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F b/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F index 3fd13a57f1..10d72a366c 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F @@ -152,6 +152,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & activeTracerVerticalAdvectionTendency, & activeTracerHorizontalAdvectionEdgeFlux, & activeTracerVerticalAdvectionTopFlux + integer :: kmax1 real (kind=RKIND), parameter :: & eps = 1.e-10_RKIND ! small number to avoid numerical difficulties @@ -269,13 +270,14 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & !$omp parallel !$omp do schedule(runtime) private(k, i, cell2, kmax) do iCell = 1, nCells - do k=1, maxLevelCell(iCell) + kmax1 = maxLevelCell(iCell) + do k=1, kmax1 tracerMin(k,iCell) = tracerCur(k,iCell) tracerMax(k,iCell) = tracerCur(k,iCell) end do do i = 1, nEdgesOnCell(iCell) cell2 = cellsOnCell(i,iCell) - kmax = min(maxLevelCell(iCell), & + kmax = min(kmax1, & maxLevelCell(cell2)) do k=1, kmax tracerMax(k,iCell) = max(tracerMax(k,iCell), & From 559c37271d1ff64c238ffca68aaff77bd6bfe203 Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Wed, 9 Sep 2020 11:53:11 -0700 Subject: [PATCH 4/8] New cpu optimization to tracer adv mono --- .../shared/mpas_ocn_tracer_advection_mono.F | 21 +++++++++++-------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F b/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F index 10d72a366c..d0180f7111 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F @@ -154,6 +154,8 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & activeTracerVerticalAdvectionTopFlux integer :: kmax1 + real (kind=RKIND) :: dvEdgeTmp + real (kind=RKIND), parameter :: & eps = 1.e-10_RKIND ! small number to avoid numerical difficulties @@ -331,10 +333,10 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & ! Also compute low order upwind horizontal flux (monotonic and diffused) ! Remove low order flux from the high order flux ! Store left over high order flux in highOrderFlx array + dvEdgeTmp = dvEdge(iEdge) * 0.5_RKIND do k = 1, maxLevelEdgeTop(iEdge) tracerWeight = iand(int(highOrderAdvectionMask(k, iEdge)+1), 1) & - * (dvEdge(iEdge) * 0.5_RKIND) & - * normalThicknessFlux(k, iEdge) + * dvEdgeTmp * normalThicknessFlux(k, iEdge) lowOrderFlx(k,iEdge) = dvEdge(iEdge) * & (max(0.0_RKIND,normalThicknessFlux(k,iEdge))*tracerCur(k,cell1) & @@ -493,14 +495,15 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & if (computeBudgets) then !$omp parallel - !$omp do schedule(runtime) private(k) + !$omp do schedule(runtime) private(k, dvEdgeTmp) do iEdge = 1,nEdgesHalo( 1 ) - do k = 1,maxLevelEdgeTop(iEdge) - ! Save u*h*T flux on edge for analysis. This variable will be - ! divided by h at the end of the time step. - activeTracerHorizontalAdvectionEdgeFlux(iTracer,k,iEdge) = & - (lowOrderFlx(k,iEdge) + highOrderFlx(k,iEdge))/dvEdge(iEdge) - enddo + dvEdgeTmp = 1.0_RKIND / dvEdge(iEdge) + do k = 1,maxLevelEdgeTop(iEdge) + ! Save u*h*T flux on edge for analysis. This variable will be + ! divided by h at the end of the time step. + activeTracerHorizontalAdvectionEdgeFlux(iTracer,k,iEdge) = & + (lowOrderFlx(k,iEdge) + highOrderFlx(k,iEdge))*dvEdgeTmp + enddo enddo !$omp end do From 2fd8a1c831d23358090157bb7d2b4b94ae0446ee Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Thu, 10 Sep 2020 08:00:29 -0700 Subject: [PATCH 5/8] Add missing dependencies in Makefile --- src/core_ocean/mode_init/Makefile | 2 ++ src/core_ocean/shared/Makefile | 8 +++++--- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/core_ocean/mode_init/Makefile b/src/core_ocean/mode_init/Makefile index 4f0e87502e..946f4bcff4 100644 --- a/src/core_ocean/mode_init/Makefile +++ b/src/core_ocean/mode_init/Makefile @@ -83,6 +83,8 @@ mpas_ocn_init_hurricane.o: $(UTILS) mpas_ocn_init_tidal_boundary.o: $(UTILS) +mpas_ocn_init_dam_break.o: $(UTILS) + #mpas_ocn_init_TEMPLATE.o: $(UTILS) clean: diff --git a/src/core_ocean/shared/Makefile b/src/core_ocean/shared/Makefile index 3db7fa8c08..cec8807a3f 100644 --- a/src/core_ocean/shared/Makefile +++ b/src/core_ocean/shared/Makefile @@ -82,7 +82,7 @@ mpas_ocn_mesh.o: mpas_ocn_thick_ale.o: mpas_ocn_constants.o mpas_ocn_config.o -mpas_ocn_time_average_coupled.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_tracer_ecosys.o +mpas_ocn_time_average_coupled.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_thick_hadv.o: mpas_ocn_constants.o mpas_ocn_config.o @@ -120,9 +120,9 @@ mpas_ocn_tracer_hmix_del4.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_tracer_advection.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_tracer_advection_mono.o mpas_ocn_tracer_advection_std.o -mpas_ocn_tracer_advection_mono.o: mpas_ocn_constants.o mpas_ocn_config.o +mpas_ocn_tracer_advection_mono.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o -mpas_ocn_tracer_advection_std.o: mpas_ocn_constants.o mpas_ocn_config.o +mpas_ocn_tracer_advection_std.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_tracer_hmix_redi.o: mpas_ocn_constants.o mpas_ocn_config.o @@ -198,6 +198,8 @@ mpas_ocn_wetting_drying.o: mpas_ocn_diagnostics.o mpas_ocn_gm.o mpas_ocn_tidal_potential_forcing.o: mpas_ocn_constants.o mpas_ocn_config.o +mpas_ocn_vel_tidal_potential.o: mpas_ocn_constants.o + mpas_ocn_vel_pressure_grad.o: mpas_ocn_constants.o clean: From c467fb3d3180ff09368df9858fedc8e5b126361b Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Thu, 10 Sep 2020 10:55:27 -0600 Subject: [PATCH 6/8] Add missing OpenMP private vars Add the newly-added variables in tracer_advection_mono to the OpenMP private clauses. --- src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F b/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F index d0180f7111..de56295c75 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_advection_mono.F @@ -270,7 +270,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & ! surrounding cells for later limiting. !$omp parallel - !$omp do schedule(runtime) private(k, i, cell2, kmax) + !$omp do schedule(runtime) private(k, i, cell2, kmax, kmax1) do iCell = 1, nCells kmax1 = maxLevelCell(iCell) do k=1, kmax1 @@ -300,7 +300,7 @@ subroutine ocn_tracer_advection_mono_tend(tend, tracers, layerThickness, & !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, k, wgtTmp, sgnTmp, flxTmp, i, iCell, coef1, coef3, & - !$omp tracerWeight) + !$omp tracerWeight, dvEdgeTmp) do iEdge = 1, nEdges cell1 = cellsOnEdge(1, iEdge) cell2 = cellsOnEdge(2, iEdge) From c7418470a65214090a4cdfd43b83554920640c58 Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Thu, 10 Sep 2020 14:09:44 -0600 Subject: [PATCH 7/8] CPU optimization for tracer adv std --- src/core_ocean/shared/mpas_ocn_tracer_advection_std.F | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F b/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F index 73c67e0bc4..67132ecbe2 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_advection_std.F @@ -88,6 +88,8 @@ subroutine ocn_tracer_advection_std_tend(tracers, normalThicknessFlux, w, layerT real (kind=RKIND), parameter :: eps = 1.e-10_RKIND + real (kind=RKIND) :: coef1 + ! Get dimensions nVertLevels = size(tracers,dim=2) num_tracers = size(tracers,dim=1) @@ -154,9 +156,11 @@ subroutine ocn_tracer_advection_std_tend(tracers, normalThicknessFlux, w, layerT cell1 = cellsOnEdge(1, iEdge) cell2 = cellsOnEdge(2, iEdge) + coef1 = dvEdge(iEdge) * 0.5_RKIND + ! Compute 2nd order fluxes where needed. do k = 1, maxLevelEdgeTop(iEdge) - tracer_weight = iand(int(highOrderAdvectionMask(k, iEdge)+1), 1) * (dvEdge(iEdge) * 0.5_RKIND) & + tracer_weight = iand(int(highOrderAdvectionMask(k, iEdge)+1), 1) * coef1 & * normalThicknessFlux(k, iEdge) high_order_horiz_flux(k, iEdge) = high_order_horiz_flux(k, iedge) + tracer_weight & From a6bd07523985d0bb95b7e1842263f72038512e88 Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Thu, 10 Sep 2020 15:13:52 -0600 Subject: [PATCH 8/8] Add ocn_mesh to tracer advection driver --- src/core_ocean/shared/mpas_ocn_tendency.F | 2 +- .../shared/mpas_ocn_tracer_advection.F | 33 ++----------------- 2 files changed, 4 insertions(+), 31 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 7b6977ec54..1a9aae5d40 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -814,7 +814,7 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, diagnosticsPool, me ! Tendency for tracer budget is stored within the tracer adv ! routine call ocn_tracer_advection_tend(tracerGroup, normalThicknessFlux, vertAleTransportTop, layerThickness, & - dt, meshPool, scratchPool, diagnosticsPool, tracerGroupTend, & + dt, diagnosticsPool, tracerGroupTend, & trim(groupItr % memberName)) ! diff --git a/src/core_ocean/shared/mpas_ocn_tracer_advection.F b/src/core_ocean/shared/mpas_ocn_tracer_advection.F index 75d8543d05..d55d797edf 100644 --- a/src/core_ocean/shared/mpas_ocn_tracer_advection.F +++ b/src/core_ocean/shared/mpas_ocn_tracer_advection.F @@ -34,6 +34,7 @@ module ocn_tracer_advection use ocn_constants use ocn_config + use ocn_mesh implicit none private @@ -67,9 +68,8 @@ module ocn_tracer_advection !------------------------------------------------------------------------------- subroutine ocn_tracer_advection_tend(tracers, normalThicknessFlux, w, & - layerThickness, dt, meshPool, & - scratchPool, diagnosticsPool, tend, & - tracerGroupName)!{{{ + layerThickness, dt, diagnosticsPool, & + tend, tracerGroupName)!{{{ !*** Input/Output parameters @@ -88,10 +88,6 @@ subroutine ocn_tracer_advection_tend(tracers, normalThicknessFlux, w, & layerThickness !< [in] Thickness field real (kind=RKIND), intent(in) :: & dt !< [in] Time step - type (mpas_pool_type), intent(in) :: & - meshPool !< [in] Mesh information - type (mpas_pool_type), intent(in) :: & - scratchPool !< [in] Scratch fields type (mpas_pool_type), intent(in) :: & diagnosticsPool !< [in] Diagnostic fields character (len=*), intent(in) :: & @@ -102,18 +98,6 @@ subroutine ocn_tracer_advection_tend(tracers, normalThicknessFlux, w, & logical :: & computeBudgets ! flag to compute active tracer budget - real (kind=RKIND), dimension(:,:), pointer, contiguous :: & - advCoefs, advCoefs3rd ! advection coefficients - - integer, dimension(:), pointer, contiguous :: & - maxLevelCell, &! index of max level at each cell - maxLevelEdgeTop, &! max level at edge with both cells active - nAdvCellsForEdge ! number of advective cells for each edge - integer, dimension(:,:), pointer, contiguous :: &! - highOrderAdvectionMask, &! mask for higher order contributions - edgeSignOnCell, &! sign at cell edge for fluxes - advCellsForEdge ! index of advective cells for each edge - ! end of preamble !---------------- ! begin code @@ -123,17 +107,6 @@ subroutine ocn_tracer_advection_tend(tracers, normalThicknessFlux, w, & call mpas_timer_start("tracer adv") - ! extract pool variables - call mpas_pool_get_array(meshPool, 'advCoefs', advCoefs) - call mpas_pool_get_array(meshPool, 'advCoefs3rd', advCoefs3rd) - call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'highOrderAdvectionMask', & - highOrderAdvectionMask) - call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell) - call mpas_pool_get_array(meshPool, 'nAdvCellsForEdge', nAdvCellsForEdge) - call mpas_pool_get_array(meshPool, 'advCellsForEdge', advCellsForEdge) - ! determine whether active tracer budgets should be computed computeBudgets = (budgetDiagsOn .and. tracerGroupName == 'activeTracers')