From a92a5214ecbefa18881ca765797e70ee488cf618 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 27 Sep 2020 11:31:18 +0200 Subject: [PATCH 01/22] Add PPR to the build --- components/mpas-ocean/src/Makefile | 13 ++++++++++--- components/mpas-ocean/src/build_options.mk | 1 + components/mpas-ocean/src/ocean.cmake | 13 +++++++++++-- 3 files changed, 22 insertions(+), 5 deletions(-) diff --git a/components/mpas-ocean/src/Makefile b/components/mpas-ocean/src/Makefile index 6ace38223a57..822ed12d5409 100644 --- a/components/mpas-ocean/src/Makefile +++ b/components/mpas-ocean/src/Makefile @@ -8,7 +8,7 @@ endif FRAMEWORK_DIR = ../../../mpas-framework/src OCEAN_SHARED_INCLUDES = -I$(FRAMEWORK_DIR)/framework -I$(FRAMEWORK_DIR)/external/esmf_time_f90 -I$(FRAMEWORK_DIR)/operators OCEAN_SHARED_INCLUDES += -I$(PWD)/shared -I$(PWD)/analysis_members -I$(PWD)/mode_forward -I$(PWD)/mode_analysis -OCEAN_SHARED_INCLUDES += -I$(PWD)/BGC -I$(PWD)/MARBL/include -I$(PWD)/cvmix/src/shared -I$(PWD)/gotm/build -I$(PWD)/gotm/build/modules -I$(PWD)/SHTNS +OCEAN_SHARED_INCLUDES += -I$(PWD)/BGC -I$(PWD)/MARBL/include -I$(PWD)/cvmix/src/shared -I$(PWD)/gotm/build -I$(PWD)/gotm/build/modules -I$(PWD)/SHTNS -I$(PWD)/ppr/src ifneq "$(EXCLUDE_INIT_MODE)" "true" OCEAN_SHARED_INCLUDES += -I$(PWD)/mode_init endif @@ -22,7 +22,7 @@ ifneq "$(EXCLUDE_INIT_MODE)" "true" OCEAN_SHARED_INCLUDES += -I$(PWD)/mode_init endif -all: shared libcvmix analysis_members libBGC libMARBL libgotm libfftw libshtns +all: shared libcvmix analysis_members libBGC libMARBL libgotm libfftw libshtns libppr (cd mode_forward; $(MAKE) FCINCLUDES="$(FCINCLUDES) $(OCEAN_SHARED_INCLUDES)" all ) (cd mode_analysis; $(MAKE) FCINCLUDES="$(FCINCLUDES) $(OCEAN_SHARED_INCLUDES)" all ) ifneq "$(EXCLUDE_INIT_MODE)" "true" @@ -131,7 +131,14 @@ libfftw: fi \ fi -shared: libcvmix libBGC libMARBL libgotm libshtns libfftw +libppr: + if [ -e ppr/.git ]; then \ + (cd ppr/src; $(CPP) $(CPPFLAGS) $(CPPINCLUDES) ppr_1d.F90 > ppr_1d.F; $(FC) $(FFLAGS) -c ppr_1d.F -o ppr_1d.o) \ + else \ + (pwd ; echo "Missing core_ocean/ppr/.git, did you forget to 'git submodule update --init --recursive' ?"; exit 1) \ + fi + +shared: libcvmix libBGC libMARBL libgotm libshtns libfftw libppr (cd shared; $(MAKE) FCINCLUDES="$(FCINCLUDES) $(OCEAN_SHARED_INCLUDES)") analysis_members: libcvmix shared libgotm diff --git a/components/mpas-ocean/src/build_options.mk b/components/mpas-ocean/src/build_options.mk index 2670b7e163a1..8e854d2db2ff 100644 --- a/components/mpas-ocean/src/build_options.mk +++ b/components/mpas-ocean/src/build_options.mk @@ -10,6 +10,7 @@ FCINCLUDES += -I$(ROOT_DIR)/cvmix/src/shared FCINCLUDES += -I$(ROOT_DIR)/BGC FCINCLUDES += -I$(ROOT_DIR)/MARBL/include FCINCLUDES += -I$(ROOT_DIR)/gotm/build/modules +FCINCLUDES += -I$(ROOT_DIR)/ppr/src override CPPFLAGS += -DCORE_OCEAN report_builds: diff --git a/components/mpas-ocean/src/ocean.cmake b/components/mpas-ocean/src/ocean.cmake index b0dabd4f25e5..8422b44b05b7 100644 --- a/components/mpas-ocean/src/ocean.cmake +++ b/components/mpas-ocean/src/ocean.cmake @@ -173,6 +173,15 @@ set(MARBL_FILES core_ocean/MARBL/src/marbl_timing_mod.F90 core_ocean/MARBL/src/marbl_utils_mod.F90 ) + +# Add PPR +if (NOT EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/core_ocean/ppr/.git) + message(FATAL_ERROR "Missing core_ocean/ppr/.git, did you forget to 'git submodule update --init --recursive' ?") +endif() +set(PPR_FILES + core_ocean/ppr/src/ppr_1d.F90 +) + # Add GOTM if (NOT EXISTS core_ocean/gotm/.git) message(FATAL "Missing core_ocean/gotm/.git, did you forget to 'git submodule update --init --recursive' ?") @@ -220,8 +229,8 @@ set(GOTM_FILES core_ocean/gotm/src/turbulence/variances.F90 ) -list(APPEND RAW_SOURCES ${CVMIX_FILES} ${BGC_FILES} ${MARBL_FILES} ${GOTM_FILES}) -list(APPEND NO_PREPROCESS ${CVMIX_FILES} ${BGC_FILES} ${MARBL_FILES} ${GOTM_FILES}) +list(APPEND RAW_SOURCES ${CVMIX_FILES} ${BGC_FILES} ${MARBL_FILES} ${GOTM_FILES} ${PPR_FILES}) +list(APPEND NO_PREPROCESS ${CVMIX_FILES} ${BGC_FILES} ${MARBL_FILES} ${GOTM_FILES} ${PPR_FILES}) # Add analysis members list(APPEND RAW_SOURCES From 825f2bc5bc263b65ad782be9b9e239e4b2982da9 Mon Sep 17 00:00:00 2001 From: dengwirda Date: Fri, 18 Dec 2020 11:11:33 -0700 Subject: [PATCH 02/22] Add (very rough) vertical remap to RK4 via PPR --- .../mpas_ocn_time_integration_rk4.F | 245 ++++++++++++++++++ 1 file changed, 245 insertions(+) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F index 844b5bc2483b..cb55cfd85c78 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F @@ -43,6 +43,8 @@ module ocn_time_integration_rk4 use ocn_surface_land_ice_fluxes use ocn_transport_tests + use ppr_1d + implicit none private save @@ -631,6 +633,22 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_timer_stop("RK4-implicit vert mix halos") + + + !! Hack in remapping here? + !! just do this for "impermeable-interfaces" for now... + + if (config_vert_coord_movement.eq.'impermeable_interfaces') then + block => domain % blocklist + do while(associated(block)) + call ocn_remap_vert_state(block, err) + + block => block % next + end do + end if + + + call mpas_timer_stop("RK4-implicit vert mix") block => domain % blocklist @@ -1647,6 +1665,233 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{ end subroutine ocn_time_integrator_rk4_cleanup!}}} + + + subroutine ocn_remap_vert_state(block, err) + + type (block_type), intent(in) :: block + integer, intent(out) :: err + + integer, pointer :: nCells, nEdges, nVertLevels + integer :: nLevels, nVars, nDoFs, nTracers + integer :: iCell, jCell, iEdge, iTrac, k + real (kind=RKIND) :: scale + + type (mpas_pool_type), pointer :: meshPool, verticalMeshPool + type (mpas_pool_type), pointer :: statePool + type (mpas_pool_type), pointer :: tracersPool + type (mpas_pool_iterator_type) :: groupItr + + type(rmap_work) :: work + type(rmap_opts) :: opts + type(rcon_ends), dimension(:), allocatable :: bcUpper, bcLower + + real (kind=RKIND), dimension(:, :, :), allocatable :: dataNow, dataNew + real (kind=RKIND), dimension(:, :), allocatable :: heightCellNow, heightEdgeNow + real (kind=RKIND), dimension(:, :), allocatable :: heightCellNew, heightEdgeNew + + integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot + integer, dimension(:, :), pointer :: cellsOnEdge + real (kind=RKIND), dimension(:), pointer :: bottomDepth + real (kind=RKIND), dimension(:, :), pointer :: layerThickness, restingThickness + real (kind=RKIND), dimension(:, :), pointer :: normalVelocity + real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroup + + call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) + call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges) + + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) + call mpas_pool_get_subpool(block % structs, 'state', statePool) + + call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) + + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', maxLevelEdgeBot) + call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + + call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) + + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 2) + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 2) + + + !---------------------------------------------------------------------------- + !- SETUP VERTICAL GRIDS: + !- heightCellNow, heightEdgeNow: depth at layer interfaces (lagrangian grid) + !- heightCellNew, heighyEdgeNew: depth at layer interfaces (new target grid) + + !- this should all go into a separate "regrid" routine + !- (eventually) this would handle regridding to "other" coordinates, rather + !- than just the z* approach here. + + allocate(heightCellNow(nVertLevels + 1, nCells)) + allocate(heightCellNew(nVertLevels + 1, nCells)) + allocate(heightEdgeNow(nVertLevels + 1, nEdges)) + allocate(heightEdgeNew(nVertLevels + 1, nEdges)) + + do iCell = 1, nCells + heightCellNow(:, iCell) = -bottomDepth(iCell) + heightCellNew(:, iCell) = -bottomDepth(iCell) + + !- reconstruct "now" heights + + do k = maxLevelCell(iCell), 1, -1 + heightCellNow(k, iCell) = & + heightCellNow(k + 1, iCell) + layerThickness(k, iCell) + end do + + !- calc. z*-style expansions + + scale = & + (heightCellNow(1, iCell) + bottomDepth(iCell)) / bottomDepth(iCell) + + !- reconstruct "new" heights + + heightCellNew(1, iCell) = heightCellNow(1, iCell) + do k = maxLevelCell(iCell), 2, -1 + heightCellNew(k, iCell) = & + heightCellNew(k + 1, iCell) + restingThickness(k, iCell) * scale + end do + end do + + !- There's currently a silly issue in PPR: it assumes the "coordinates" for + !- the data being remapped increases with index. Here, we have the opposite + !- (z becomes more -ve with k)... + !- I'll make some updates to PPR to handle this, but for now, just flip the + !- signs to make PPR happy. This shouldn't have any actual impact. + + heightCellNow = -1.0_RKIND * heightCellNow + heightCellNew = -1.0_RKIND * heightCellNew + + + do iEdge = 1, nEdges + iCell = cellsOnEdge(1, iEdge) + jCell = cellsOnEdge(2, iEdge) + + heightEdgeNow(:, iEdge) = -bottomDepth(iCell) + heightEdgeNew(:, iEdge) = -bottomDepth(iCell) + + !- this is just a quick hack for the edge heights + !- really should use layerThicknessEdge here, but it's unclear if this has + !- been properly initiallised at this point... + !- layerThicknessEdge could involve upwinding, so need to be careful! + + do k = maxLevelEdgeTop(iEdge) + 1, 1, -1 + heightEdgeNow(k, iEdge) = & + 0.5_RKIND * (heightCellNow(k, iCell) + heightCellNow(k, jCell)) + heightEdgeNew(k, iEdge) = & + 0.5_RKIND * (heightCellNew(k, iCell) + heightCellNew(k, jCell)) + end do + end do + + + !---------------------------------------------------------------------------- + !- ACTUAL REMAPPING FROM HERE + + !- just simple hard-coded options for PPR at this stage... + + nVars = 1 ! how many tracers to be remapped concurrently per column + nDoFs = 1 ! always 1 for MPAS, it'd be >1 for DG methods (with more DoF per layer) + + allocate(dataNow(nDoFs, nVars, nVertLevels)) + allocate(dataNew(nDoFs, nVars, nVertLevels)) + allocate(bcUpper(nVars)) + allocate(bcLower(nVars)) + + dataNow = 0.0_RKIND + dataNew = 0.0_RKIND + + opts%edge_meth = p3e_method ! 3rd-order edge interp. + opts%cell_meth = ppm_method ! PPM method + opts%cell_lims = mono_limit ! monotone slope limits + + bcUpper%bcopt = bcon_loose ! "loose" => extrapolate + bcLower%bcopt = bcon_loose + + call work%init(nVertLevels + 1, nVars, opts) ! this is internal workspace for PPR + + !- Remap all edge-centred variables from heightEdgeNow to heightEdgeNew + + do iEdge = 1, nEdges + nLevels = maxLevelEdgeTop(iEdge) + 1 + + if (nLevels .lt. 2) cycle + + dataNow(1, 1, :) = normalVelocity(:, iEdge) + !also do normalTransportVelocity, etc...? + + call rmap1d(nLevels, nLevels, nVars ,nDoFs, & + heightEdgeNow(:, iEdge), & + heightEdgeNew(:, iEdge), & + dataNow, dataNew, bcUpper, bcLower, work, opts) + + normalVelocity(:, iEdge) = dataNew(1, 1, :) + end do + + !- Remap all cell-centred variables from heightCellNow to heightCellNew + + call mpas_pool_begin_iteration(tracersPool) + do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) + if ( groupItr % memberType == MPAS_POOL_FIELD ) then + call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroup, 2) + if ( associated(tracersGroup) ) then + + !- we can remap all of the tracers in one go via PPR, but for now + !- just a simple loop over each one at a time for clarity... + + nTracers = size(tracersGroup, dim=1) + do iCell = 1, nCells + do iTrac = 1, nTracers + nLevels = maxLevelCell(iCell) + 1 + + if (nLevels .lt. 2) cycle + + dataNow(1, 1, :) = tracersGroup(iTrac, :, iCell) + + call rmap1d(nLevels, nLevels, nVars ,nDoFs, & + heightCellNow(:, iCell), & + heightCellNew(:, iCell), & + dataNow, dataNew, bcUpper, bcLower, work, opts) + + tracersGroup(iTrac, :, iCell) = dataNew(1, 1, :) + end do + end do + + end if + end if + end do + + !- Assign new layerThicknesses + + do iCell = 1, nCells + do k = maxLevelCell(iCell), 1, -1 + layerThickness(k, iCell) = & + heightCellNew(k + 1, iCell) - heightCellNew(k, iCell) + end do + end do + + !- dealloc any local data, etc + + call work%free() + + deallocate(heightCellNow) + deallocate(heightEdgeNow) + deallocate(heightCellNew) + deallocate(heightEdgeNew) + + deallocate(dataNow) + deallocate(dataNew) + deallocate(bcUpper) + deallocate(bcLower) + + end subroutine ocn_remap_vert_state + + + end module ocn_time_integration_rk4 ! vim: foldmethod=marker From c75ab8cf15994dfdd962558d914d2af27e53818b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 10 Feb 2021 13:15:38 -0700 Subject: [PATCH 03/22] Add config options relevant to VLR: - Add config option to turn remap on/off - Add config option for advection, remap order - Add config option for target grid - Only call vert_transport_top if advection_method is flux-form - Construct opts for remapping using config options --- components/mpas-ocean/src/Registry.xml | 14 +- .../mpas_ocn_time_integration_rk4.F | 167 ++++++++++++------ 2 files changed, 119 insertions(+), 62 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 9e9128f39dda..8211abe45e51 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1000,12 +1000,16 @@ - - + domain % blocklist do while(associated(block)) call ocn_remap_vert_state(block, err) @@ -906,17 +915,20 @@ subroutine ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, & call mpas_pool_get_array(provisStatePool, 'normalVelocity', normalVelocityProvis, 1) call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) + ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) - else - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) + if ( configVertAdvMethod == 1) then + if (associated(highFreqThicknessProvis)) then + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err, highFreqThicknessProvis) + else + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err) + endif endif call ocn_tend_vel(domain, tendPool, provisStatePool, forcingPool, 1, dminfo, dt) @@ -975,16 +987,18 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh end if ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) - else - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) + if ( configVertAdvMethod == 1) then + if (associated(highFreqThicknessProvis)) then + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err, highFreqThicknessProvis) + else + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err) + endif endif call ocn_tend_thick(tendPool, forcingPool) @@ -1045,16 +1059,18 @@ subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeig end if ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) - else - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) + if ( configVertAdvMethod == 1) then + if (associated(highFreqThicknessProvis)) then + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err, highFreqThicknessProvis) + else + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err) + endif endif if (config_filter_btr_mode) then @@ -1108,30 +1124,34 @@ subroutine ocn_time_integrator_rk4_compute_tends(domain, block, dt, rkWeight, dm call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) - else - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkWeight, & - vertAleTransportTop, err) + if ( configVertAdvMethod == 1) then + if (associated(highFreqThicknessProvis)) then + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & + sshCur, rkWeight, & + vertAleTransportTop, err, highFreqThicknessProvis) + else + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & + sshCur, rkWeight, & + vertAleTransportTop, err) + endif endif call ocn_tend_vel(domain, tendPool, provisStatePool, forcingPool, 1, dminfo, dt) - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) - else - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkWeight, & - vertAleTransportTop, err) + if ( configVertAdvMethod == 1) then + if (associated(highFreqThicknessProvis)) then + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkWeight, & + vertAleTransportTop, err, highFreqThicknessProvis) + else + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkWeight, & + vertAleTransportTop, err) + endif endif call ocn_tend_thick(tendPool, forcingPool) @@ -1695,8 +1715,14 @@ subroutine ocn_remap_vert_state(block, err) real (kind=RKIND), dimension(:), pointer :: bottomDepth real (kind=RKIND), dimension(:, :), pointer :: layerThickness, restingThickness real (kind=RKIND), dimension(:, :), pointer :: normalVelocity + real (kind=RKIND), dimension(:, :), pointer :: highFreqThickness,lowFreqDivergence real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroup + integer,pointer :: config_vert_remap_order + + err = 0 + + call mpas_pool_get_config(block % configs, 'config_vert_remap_order', config_vert_remap_order) call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges) @@ -1717,7 +1743,25 @@ subroutine ocn_remap_vert_state(block, err) call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 2) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 2) - + if (config_use_freq_filtered_thickness) then + call mpas_pool_get_array(statePool, 'highFreqThickness', highFreqThickness, 2) + call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergence, 2) + end if + + ! Options for remapping + ! TODO move to initialization routine + if ( config_vert_remap_order == 1 ) then + opts%cell_meth = pcm_method ! PCM method + elseif ( config_vert_remap_order == 2 ) then + opts%cell_meth = plm_method ! PLM method + elseif ( config_vert_remap_order == 3 ) then + opts%cell_meth = ppm_method ! PPM method + opts%edge_meth = p3e_method ! 3rd-order edge interp. + elseif ( config_vert_remap_order == 5 ) then + opts%cell_meth = pqm_method ! PPM method + opts%edge_meth = p5e_method ! 3rd-order edge interp. + endif + opts%cell_lims = mono_limit ! monotone slope limits TODO: consider config_monotonic !---------------------------------------------------------------------------- !- SETUP VERTICAL GRIDS: @@ -1805,8 +1849,17 @@ subroutine ocn_remap_vert_state(block, err) dataNow = 0.0_RKIND dataNew = 0.0_RKIND - opts%edge_meth = p3e_method ! 3rd-order edge interp. - opts%cell_meth = ppm_method ! PPM method + if ( config_vert_remap_order == 1 ) then + opts%cell_meth = pcm_method ! PCM method + elseif ( config_vert_remap_order == 2 ) then + opts%cell_meth = plm_method ! PLM method + elseif ( config_vert_remap_order == 3 ) then + opts%cell_meth = ppm_method ! PPM method + opts%edge_meth = p3e_method ! 3rd-order edge interp. + elseif ( config_vert_remap_order == 5 ) then + opts%cell_meth = pqm_method ! PPM method + opts%edge_meth = p5e_method ! 3rd-order edge interp. + endif opts%cell_lims = mono_limit ! monotone slope limits bcUpper%bcopt = bcon_loose ! "loose" => extrapolate From bc649ff0a8c6dd3df277d10edfd936b500d40dc7 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 10 Feb 2021 12:13:55 -0700 Subject: [PATCH 04/22] Cleanup rough remap routine --- .../mpas_ocn_time_integration_rk4.F | 189 +++++++++++------- 1 file changed, 122 insertions(+), 67 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F index f58f2caed509..799b712d22d1 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F @@ -53,9 +53,9 @@ module ocn_time_integration_rk4 ! ! Public parameters ! - + integer :: configVertAdvMethod - + !-------------------------------------------------------------------- !-------------------------------------------------------------------- @@ -202,6 +202,10 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersNew + ! Remap variables + real (kind=RKIND) :: totalThickness + real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessRemap + ! Get config options call mpas_pool_get_config(domain % configs, 'config_mom_del4', config_mom_del4) call mpas_pool_get_config(domain % configs, 'config_filter_btr_mode', config_filter_btr_mode) @@ -223,7 +227,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_pool_get_config(domain % configs, 'config_tidal_forcing_type', config_tidal_forcing_type) call mpas_pool_get_config(domain % configs, 'config_vert_advection_method', config_vert_advection_method) - + if ( config_vert_advection_method == 'flux-form' ) then configVertAdvMethod = 1 elseif ( config_vert_advection_method == 'remap' ) then @@ -650,14 +654,46 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ if ( configVertAdvMethod == 2 ) then block => domain % blocklist do while(associated(block)) - call ocn_remap_vert_state(block, err) + call mpas_pool_get_subpool(block % structs, 'state', statePool) + call mpas_pool_get_subpool(block % structs, 'tend', tendPool) + call mpas_pool_get_array(tendPool, 'layerThickness', layerThicknessTend) + call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessNew, 2) + call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) + call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) + call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityCur, 1) + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityNew, 2) + + allocate(layerThicknessRemap(nVertLevels, nCells)) + layerThicknessRemap = layerThicknessNew + + !- calc. z*-style expansions + !$omp parallel + !$omp do schedule(runtime) private(k,totalThickness) + do iCell = 1, nCells + totalThickness = 0.0_RKIND + do k = maxLevelCell(iCell), 1, -1 + totalThickness = totalThickness + layerThicknessNew(k,iCell) + end do + do k = maxLevelCell(iCell), 1, -1 + layerThicknessRemap(k, iCell) = restingThickness(k, iCell) * & + totalThickness / bottomDepth(iCell) + end do + end do + !$omp end do + !$omp end parallel + + call ocn_remap_vert_state(block, layerThicknessRemap, err) + deallocate(layerThicknessRemap) block => block % next end do end if - call mpas_timer_stop("RK4-implicit vert mix") block => domain % blocklist @@ -1687,15 +1723,15 @@ end subroutine ocn_time_integrator_rk4_cleanup!}}} - subroutine ocn_remap_vert_state(block, err) + subroutine ocn_remap_vert_state(block, layerThicknessNew, err) type (block_type), intent(in) :: block + real (kind=RKIND), dimension(:, :), intent(in) :: layerThicknessNew integer, intent(out) :: err integer, pointer :: nCells, nEdges, nVertLevels integer :: nLevels, nVars, nDoFs, nTracers integer :: iCell, jCell, iEdge, iTrac, k - real (kind=RKIND) :: scale type (mpas_pool_type), pointer :: meshPool, verticalMeshPool type (mpas_pool_type), pointer :: statePool @@ -1780,7 +1816,7 @@ subroutine ocn_remap_vert_state(block, err) do iCell = 1, nCells heightCellNow(:, iCell) = -bottomDepth(iCell) heightCellNew(:, iCell) = -bottomDepth(iCell) - + !- reconstruct "now" heights do k = maxLevelCell(iCell), 1, -1 @@ -1788,19 +1824,16 @@ subroutine ocn_remap_vert_state(block, err) heightCellNow(k + 1, iCell) + layerThickness(k, iCell) end do - !- calc. z*-style expansions - - scale = & - (heightCellNow(1, iCell) + bottomDepth(iCell)) / bottomDepth(iCell) - !- reconstruct "new" heights + !- Assign new layerThicknesses heightCellNew(1, iCell) = heightCellNow(1, iCell) do k = maxLevelCell(iCell), 2, -1 heightCellNew(k, iCell) = & - heightCellNew(k + 1, iCell) + restingThickness(k, iCell) * scale + heightCellNew(k + 1, iCell) + layerThicknessNew(k, iCell) end do end do + layerThickness = layerThicknessNew !- There's currently a silly issue in PPR: it assumes the "coordinates" for !- the data being remapped increases with index. Here, we have the opposite @@ -1824,13 +1857,16 @@ subroutine ocn_remap_vert_state(block, err) !- been properly initiallised at this point... !- layerThicknessEdge could involve upwinding, so need to be careful! - do k = maxLevelEdgeTop(iEdge) + 1, 1, -1 + do k = maxLevelEdgeTop(iEdge), 1, -1 heightEdgeNow(k, iEdge) = & - 0.5_RKIND * (heightCellNow(k, iCell) + heightCellNow(k, jCell)) + 0.5_RKIND * (heightCellNow(k, iCell) + heightCellNow(k, jCell)) heightEdgeNew(k, iEdge) = & - 0.5_RKIND * (heightCellNew(k, iCell) + heightCellNew(k, jCell)) + 0.5_RKIND * (heightCellNew(k, iCell) + heightCellNew(k, jCell)) end do end do + !- Assign new layerThicknesses + write(*,*) 'set layerThicknessNew to layerThicknessRemap' + layerThickness = layerThicknessNew !---------------------------------------------------------------------------- @@ -1838,7 +1874,11 @@ subroutine ocn_remap_vert_state(block, err) !- just simple hard-coded options for PPR at this stage... - nVars = 1 ! how many tracers to be remapped concurrently per column + if (associated(highFreqThickness) .and. associated(lowFreqDivergence)) then + nVars = 3 ! how many tracers to be remapped concurrently per column + else + nVars = 1 ! how many tracers to be remapped concurrently per column + endif nDoFs = 1 ! always 1 for MPAS, it'd be >1 for DG methods (with more DoF per layer) allocate(dataNow(nDoFs, nVars, nVertLevels)) @@ -1849,98 +1889,113 @@ subroutine ocn_remap_vert_state(block, err) dataNow = 0.0_RKIND dataNew = 0.0_RKIND - if ( config_vert_remap_order == 1 ) then - opts%cell_meth = pcm_method ! PCM method - elseif ( config_vert_remap_order == 2 ) then - opts%cell_meth = plm_method ! PLM method - elseif ( config_vert_remap_order == 3 ) then - opts%cell_meth = ppm_method ! PPM method - opts%edge_meth = p3e_method ! 3rd-order edge interp. - elseif ( config_vert_remap_order == 5 ) then - opts%cell_meth = pqm_method ! PPM method - opts%edge_meth = p5e_method ! 3rd-order edge interp. - endif - opts%cell_lims = mono_limit ! monotone slope limits - bcUpper%bcopt = bcon_loose ! "loose" => extrapolate bcLower%bcopt = bcon_loose call work%init(nVertLevels + 1, nVars, opts) ! this is internal workspace for PPR !- Remap all edge-centred variables from heightEdgeNow to heightEdgeNew + write(*,*) 'remap normalVelocity' - do iEdge = 1, nEdges - nLevels = maxLevelEdgeTop(iEdge) + 1 + if (associated(highFreqThickness) .and. associated(lowFreqDivergence)) then + do iEdge = 1, nEdges + nLevels = maxLevelEdgeTop(iEdge) + 1 - if (nLevels .lt. 2) cycle + if (nLevels .lt. 2) cycle - dataNow(1, 1, :) = normalVelocity(:, iEdge) - !also do normalTransportVelocity, etc...? + dataNow(1, 1, :) = normalVelocity(:, iEdge) + dataNow(1, 2, :) = highFreqThickness(:, iEdge) + dataNow(1, 3, :) = lowFreqDivergence(:, iEdge) - call rmap1d(nLevels, nLevels, nVars ,nDoFs, & - heightEdgeNow(:, iEdge), & - heightEdgeNew(:, iEdge), & - dataNow, dataNew, bcUpper, bcLower, work, opts) + call rmap1d(nLevels, nLevels, nVars ,nDoFs, & + heightEdgeNow(:, iEdge), & + heightEdgeNew(:, iEdge), & + dataNow, dataNew, bcUpper, bcLower, work, opts) - normalVelocity(:, iEdge) = dataNew(1, 1, :) - end do + normalVelocity (1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 1, 1:maxLevelEdgeTop(iEdge)) + highFreqThickness(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 2, 1:maxLevelEdgeTop(iEdge)) + lowFreqDivergence(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 3, 1:maxLevelEdgeTop(iEdge)) + + end do + else + do iEdge = 1, nEdges + nLevels = maxLevelEdgeTop(iEdge) + 1 + + if (nLevels .lt. 2) cycle + + dataNow(1, 1, :) = normalVelocity(:, iEdge) + + call rmap1d(nLevels, nLevels, nVars ,nDoFs, & + heightEdgeNow(:, iEdge), & + heightEdgeNew(:, iEdge), & + dataNow, dataNew, bcUpper, bcLower, work, opts) + + normalVelocity(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 1, 1:maxLevelEdgeTop(iEdge)) + + end do + end if + deallocate(dataNow,dataNew) + deallocate(bcUpper,bcLower) + call work%free() !- Remap all cell-centred variables from heightCellNow to heightCellNew + write(*,*) 'remap Tracers' call mpas_pool_begin_iteration(tracersPool) do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) if ( groupItr % memberType == MPAS_POOL_FIELD ) then call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroup, 2) if ( associated(tracersGroup) ) then - + !- we can remap all of the tracers in one go via PPR, but for now !- just a simple loop over each one at a time for clarity... nTracers = size(tracersGroup, dim=1) + + allocate(dataNow(nDoFs, nTracers, nVertLevels)) + allocate(dataNew(nDoFs, nTracers, nVertLevels)) + + call work%init(nVertLevels + 1, nTracers, opts) ! this is internal workspace for PPR + + allocate(bcUpper(nTracers)) + allocate(bcLower(nTracers)) + ! TODO replace hardcoded bcon_loose with config option + bcUpper%bcopt = bcon_loose ! "loose" => extrapolate + bcLower%bcopt = bcon_loose + do iCell = 1, nCells do iTrac = 1, nTracers nLevels = maxLevelCell(iCell) + 1 if (nLevels .lt. 2) cycle - dataNow(1, 1, :) = tracersGroup(iTrac, :, iCell) - - call rmap1d(nLevels, nLevels, nVars ,nDoFs, & - heightCellNow(:, iCell), & - heightCellNew(:, iCell), & - dataNow, dataNew, bcUpper, bcLower, work, opts) + dataNow(1, iTrac, :) = tracersGroup(iTrac, :, iCell) + end do + call rmap1d(nLevels, nLevels, nTracers,nDoFs, & + heightCellNow(:, iCell), & + heightCellNew(:, iCell), & + dataNow, dataNew, bcUpper, bcLower, work, opts) - tracersGroup(iTrac, :, iCell) = dataNew(1, 1, :) + do iTrac = 1, nTracers + tracersGroup(iTrac, 1:maxLevelCell(iCell), iCell) = dataNew(1, iTrac, 1:maxLevelCell(iCell)) end do + end do - + deallocate(dataNow,dataNew) + deallocate(bcUpper,bcLower) + call work%free() + end if end if end do - !- Assign new layerThicknesses - - do iCell = 1, nCells - do k = maxLevelCell(iCell), 1, -1 - layerThickness(k, iCell) = & - heightCellNew(k + 1, iCell) - heightCellNew(k, iCell) - end do - end do - !- dealloc any local data, etc - call work%free() - deallocate(heightCellNow) deallocate(heightEdgeNow) deallocate(heightCellNew) deallocate(heightEdgeNew) - deallocate(dataNow) - deallocate(dataNew) - deallocate(bcUpper) - deallocate(bcLower) - end subroutine ocn_remap_vert_state From c8b81a518752f3f4dd22a696bcdd7af05e9a7fe0 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 7 Jul 2021 09:26:33 -0500 Subject: [PATCH 05/22] Move VLR from RK4 to split-explicit and a new module --- .../src/mode_forward/mpas_ocn_forward_mode.F | 4 + .../mpas_ocn_time_integration_rk4.F | 453 ++---------------- .../mpas_ocn_time_integration_split.F | 50 +- components/mpas-ocean/src/ocean.cmake | 1 + components/mpas-ocean/src/shared/Makefile | 1 + .../src/shared/mpas_ocn_vertical_remap.F | 353 ++++++++++++++ 6 files changed, 455 insertions(+), 407 deletions(-) create mode 100644 components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F index 1375ccca8a76..001a881131aa 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F @@ -83,6 +83,8 @@ module ocn_forward_mode use ocn_vmix use ocn_vmix_gotm, only: ocn_vmix_gotm_finalize + use ocn_vertical_remap + use ocn_forcing use ocn_time_varying_forcing @@ -462,6 +464,8 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{ ierr = ior(ierr, err_tmp) call ocn_tendency_init(err_tmp) ierr = ior(ierr,err_tmp) + call ocn_vertical_remap_init(err_tmp) + ierr = ior(ierr,err_tmp) if (ierr /= 0) then call mpas_log_write( & diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F index 799b712d22d1..844b5bc2483b 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F @@ -43,8 +43,6 @@ module ocn_time_integration_rk4 use ocn_surface_land_ice_fluxes use ocn_transport_tests - use ppr_1d - implicit none private save @@ -53,9 +51,6 @@ module ocn_time_integration_rk4 ! ! Public parameters ! - - integer :: configVertAdvMethod - !-------------------------------------------------------------------- !-------------------------------------------------------------------- @@ -202,10 +197,6 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersNew - ! Remap variables - real (kind=RKIND) :: totalThickness - real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessRemap - ! Get config options call mpas_pool_get_config(domain % configs, 'config_mom_del4', config_mom_del4) call mpas_pool_get_config(domain % configs, 'config_filter_btr_mode', config_filter_btr_mode) @@ -225,16 +216,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_pool_get_config(domain % configs, 'config_zero_drying_velocity', config_zero_drying_velocity) call mpas_pool_get_config(domain % configs, 'config_use_tidal_forcing', config_use_tidal_forcing) call mpas_pool_get_config(domain % configs, 'config_tidal_forcing_type', config_tidal_forcing_type) - call mpas_pool_get_config(domain % configs, 'config_vert_advection_method', config_vert_advection_method) - - if ( config_vert_advection_method == 'flux-form' ) then - configVertAdvMethod = 1 - elseif ( config_vert_advection_method == 'remap' ) then - configVertAdvMethod = 2 - endif - - ! ! Initialize time_levs(2) with state at current time ! Initialize first RK state @@ -649,51 +631,6 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_timer_stop("RK4-implicit vert mix halos") - - - if ( configVertAdvMethod == 2 ) then - block => domain % blocklist - do while(associated(block)) - call mpas_pool_get_subpool(block % structs, 'state', statePool) - call mpas_pool_get_subpool(block % structs, 'tend', tendPool) - call mpas_pool_get_array(tendPool, 'layerThickness', layerThicknessTend) - call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessNew, 2) - call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) - call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) - call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) - call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) - call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) - call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityCur, 1) - call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityNew, 2) - - allocate(layerThicknessRemap(nVertLevels, nCells)) - layerThicknessRemap = layerThicknessNew - - !- calc. z*-style expansions - !$omp parallel - !$omp do schedule(runtime) private(k,totalThickness) - do iCell = 1, nCells - totalThickness = 0.0_RKIND - do k = maxLevelCell(iCell), 1, -1 - totalThickness = totalThickness + layerThicknessNew(k,iCell) - end do - do k = maxLevelCell(iCell), 1, -1 - layerThicknessRemap(k, iCell) = restingThickness(k, iCell) * & - totalThickness / bottomDepth(iCell) - end do - end do - !$omp end do - !$omp end parallel - - call ocn_remap_vert_state(block, layerThicknessRemap, err) - - deallocate(layerThicknessRemap) - block => block % next - end do - end if - - call mpas_timer_stop("RK4-implicit vert mix") block => domain % blocklist @@ -951,20 +888,17 @@ subroutine ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, & call mpas_pool_get_array(provisStatePool, 'normalVelocity', normalVelocityProvis, 1) call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) - ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if ( configVertAdvMethod == 1) then - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) - else - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) - endif + if (associated(highFreqThicknessProvis)) then + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err, highFreqThicknessProvis) + else + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err) endif call ocn_tend_vel(domain, tendPool, provisStatePool, forcingPool, 1, dminfo, dt) @@ -1023,18 +957,16 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh end if ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if ( configVertAdvMethod == 1) then - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) - else - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) - endif + if (associated(highFreqThicknessProvis)) then + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err, highFreqThicknessProvis) + else + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err) endif call ocn_tend_thick(tendPool, forcingPool) @@ -1095,18 +1027,16 @@ subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeig end if ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if ( configVertAdvMethod == 1) then - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) - else - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) - endif + if (associated(highFreqThicknessProvis)) then + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err, highFreqThicknessProvis) + else + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkSubstepWeight, & + vertAleTransportTop, err) endif if (config_filter_btr_mode) then @@ -1160,34 +1090,30 @@ subroutine ocn_time_integrator_rk4_compute_tends(domain, block, dt, rkWeight, dm call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if ( configVertAdvMethod == 1) then - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) - else - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkWeight, & - vertAleTransportTop, err) - endif + if (associated(highFreqThicknessProvis)) then + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & + sshCur, rkWeight, & + vertAleTransportTop, err, highFreqThicknessProvis) + else + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & + sshCur, rkWeight, & + vertAleTransportTop, err) endif call ocn_tend_vel(domain, tendPool, provisStatePool, forcingPool, 1, dminfo, dt) - if ( configVertAdvMethod == 1) then - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) - else - call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkWeight, & - vertAleTransportTop, err) - endif + if (associated(highFreqThicknessProvis)) then + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkWeight, & + vertAleTransportTop, err, highFreqThicknessProvis) + else + call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, & + layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & + sshCur, rkWeight, & + vertAleTransportTop, err) endif call ocn_tend_thick(tendPool, forcingPool) @@ -1721,285 +1647,6 @@ subroutine ocn_time_integrator_rk4_cleanup(block, dt, err)!{{{ end subroutine ocn_time_integrator_rk4_cleanup!}}} - - - subroutine ocn_remap_vert_state(block, layerThicknessNew, err) - - type (block_type), intent(in) :: block - real (kind=RKIND), dimension(:, :), intent(in) :: layerThicknessNew - integer, intent(out) :: err - - integer, pointer :: nCells, nEdges, nVertLevels - integer :: nLevels, nVars, nDoFs, nTracers - integer :: iCell, jCell, iEdge, iTrac, k - - type (mpas_pool_type), pointer :: meshPool, verticalMeshPool - type (mpas_pool_type), pointer :: statePool - type (mpas_pool_type), pointer :: tracersPool - type (mpas_pool_iterator_type) :: groupItr - - type(rmap_work) :: work - type(rmap_opts) :: opts - type(rcon_ends), dimension(:), allocatable :: bcUpper, bcLower - - real (kind=RKIND), dimension(:, :, :), allocatable :: dataNow, dataNew - real (kind=RKIND), dimension(:, :), allocatable :: heightCellNow, heightEdgeNow - real (kind=RKIND), dimension(:, :), allocatable :: heightCellNew, heightEdgeNew - - integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot - integer, dimension(:, :), pointer :: cellsOnEdge - real (kind=RKIND), dimension(:), pointer :: bottomDepth - real (kind=RKIND), dimension(:, :), pointer :: layerThickness, restingThickness - real (kind=RKIND), dimension(:, :), pointer :: normalVelocity - real (kind=RKIND), dimension(:, :), pointer :: highFreqThickness,lowFreqDivergence - real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroup - - integer,pointer :: config_vert_remap_order - - err = 0 - - call mpas_pool_get_config(block % configs, 'config_vert_remap_order', config_vert_remap_order) - call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) - call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges) - - call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) - call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) - call mpas_pool_get_subpool(block % structs, 'state', statePool) - - call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) - - call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) - call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', maxLevelEdgeBot) - call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - - call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) - - call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 2) - call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 2) - if (config_use_freq_filtered_thickness) then - call mpas_pool_get_array(statePool, 'highFreqThickness', highFreqThickness, 2) - call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergence, 2) - end if - - ! Options for remapping - ! TODO move to initialization routine - if ( config_vert_remap_order == 1 ) then - opts%cell_meth = pcm_method ! PCM method - elseif ( config_vert_remap_order == 2 ) then - opts%cell_meth = plm_method ! PLM method - elseif ( config_vert_remap_order == 3 ) then - opts%cell_meth = ppm_method ! PPM method - opts%edge_meth = p3e_method ! 3rd-order edge interp. - elseif ( config_vert_remap_order == 5 ) then - opts%cell_meth = pqm_method ! PPM method - opts%edge_meth = p5e_method ! 3rd-order edge interp. - endif - opts%cell_lims = mono_limit ! monotone slope limits TODO: consider config_monotonic - - !---------------------------------------------------------------------------- - !- SETUP VERTICAL GRIDS: - !- heightCellNow, heightEdgeNow: depth at layer interfaces (lagrangian grid) - !- heightCellNew, heighyEdgeNew: depth at layer interfaces (new target grid) - - !- this should all go into a separate "regrid" routine - !- (eventually) this would handle regridding to "other" coordinates, rather - !- than just the z* approach here. - - allocate(heightCellNow(nVertLevels + 1, nCells)) - allocate(heightCellNew(nVertLevels + 1, nCells)) - allocate(heightEdgeNow(nVertLevels + 1, nEdges)) - allocate(heightEdgeNew(nVertLevels + 1, nEdges)) - - do iCell = 1, nCells - heightCellNow(:, iCell) = -bottomDepth(iCell) - heightCellNew(:, iCell) = -bottomDepth(iCell) - - !- reconstruct "now" heights - - do k = maxLevelCell(iCell), 1, -1 - heightCellNow(k, iCell) = & - heightCellNow(k + 1, iCell) + layerThickness(k, iCell) - end do - - !- reconstruct "new" heights - !- Assign new layerThicknesses - - heightCellNew(1, iCell) = heightCellNow(1, iCell) - do k = maxLevelCell(iCell), 2, -1 - heightCellNew(k, iCell) = & - heightCellNew(k + 1, iCell) + layerThicknessNew(k, iCell) - end do - end do - layerThickness = layerThicknessNew - - !- There's currently a silly issue in PPR: it assumes the "coordinates" for - !- the data being remapped increases with index. Here, we have the opposite - !- (z becomes more -ve with k)... - !- I'll make some updates to PPR to handle this, but for now, just flip the - !- signs to make PPR happy. This shouldn't have any actual impact. - - heightCellNow = -1.0_RKIND * heightCellNow - heightCellNew = -1.0_RKIND * heightCellNew - - - do iEdge = 1, nEdges - iCell = cellsOnEdge(1, iEdge) - jCell = cellsOnEdge(2, iEdge) - - heightEdgeNow(:, iEdge) = -bottomDepth(iCell) - heightEdgeNew(:, iEdge) = -bottomDepth(iCell) - - !- this is just a quick hack for the edge heights - !- really should use layerThicknessEdge here, but it's unclear if this has - !- been properly initiallised at this point... - !- layerThicknessEdge could involve upwinding, so need to be careful! - - do k = maxLevelEdgeTop(iEdge), 1, -1 - heightEdgeNow(k, iEdge) = & - 0.5_RKIND * (heightCellNow(k, iCell) + heightCellNow(k, jCell)) - heightEdgeNew(k, iEdge) = & - 0.5_RKIND * (heightCellNew(k, iCell) + heightCellNew(k, jCell)) - end do - end do - !- Assign new layerThicknesses - write(*,*) 'set layerThicknessNew to layerThicknessRemap' - layerThickness = layerThicknessNew - - - !---------------------------------------------------------------------------- - !- ACTUAL REMAPPING FROM HERE - - !- just simple hard-coded options for PPR at this stage... - - if (associated(highFreqThickness) .and. associated(lowFreqDivergence)) then - nVars = 3 ! how many tracers to be remapped concurrently per column - else - nVars = 1 ! how many tracers to be remapped concurrently per column - endif - nDoFs = 1 ! always 1 for MPAS, it'd be >1 for DG methods (with more DoF per layer) - - allocate(dataNow(nDoFs, nVars, nVertLevels)) - allocate(dataNew(nDoFs, nVars, nVertLevels)) - allocate(bcUpper(nVars)) - allocate(bcLower(nVars)) - - dataNow = 0.0_RKIND - dataNew = 0.0_RKIND - - bcUpper%bcopt = bcon_loose ! "loose" => extrapolate - bcLower%bcopt = bcon_loose - - call work%init(nVertLevels + 1, nVars, opts) ! this is internal workspace for PPR - - !- Remap all edge-centred variables from heightEdgeNow to heightEdgeNew - write(*,*) 'remap normalVelocity' - - if (associated(highFreqThickness) .and. associated(lowFreqDivergence)) then - do iEdge = 1, nEdges - nLevels = maxLevelEdgeTop(iEdge) + 1 - - if (nLevels .lt. 2) cycle - - dataNow(1, 1, :) = normalVelocity(:, iEdge) - dataNow(1, 2, :) = highFreqThickness(:, iEdge) - dataNow(1, 3, :) = lowFreqDivergence(:, iEdge) - - call rmap1d(nLevels, nLevels, nVars ,nDoFs, & - heightEdgeNow(:, iEdge), & - heightEdgeNew(:, iEdge), & - dataNow, dataNew, bcUpper, bcLower, work, opts) - - normalVelocity (1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 1, 1:maxLevelEdgeTop(iEdge)) - highFreqThickness(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 2, 1:maxLevelEdgeTop(iEdge)) - lowFreqDivergence(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 3, 1:maxLevelEdgeTop(iEdge)) - - end do - else - do iEdge = 1, nEdges - nLevels = maxLevelEdgeTop(iEdge) + 1 - - if (nLevels .lt. 2) cycle - - dataNow(1, 1, :) = normalVelocity(:, iEdge) - - call rmap1d(nLevels, nLevels, nVars ,nDoFs, & - heightEdgeNow(:, iEdge), & - heightEdgeNew(:, iEdge), & - dataNow, dataNew, bcUpper, bcLower, work, opts) - - normalVelocity(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 1, 1:maxLevelEdgeTop(iEdge)) - - end do - end if - deallocate(dataNow,dataNew) - deallocate(bcUpper,bcLower) - call work%free() - - !- Remap all cell-centred variables from heightCellNow to heightCellNew - write(*,*) 'remap Tracers' - - call mpas_pool_begin_iteration(tracersPool) - do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) - if ( groupItr % memberType == MPAS_POOL_FIELD ) then - call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroup, 2) - if ( associated(tracersGroup) ) then - - !- we can remap all of the tracers in one go via PPR, but for now - !- just a simple loop over each one at a time for clarity... - - nTracers = size(tracersGroup, dim=1) - - allocate(dataNow(nDoFs, nTracers, nVertLevels)) - allocate(dataNew(nDoFs, nTracers, nVertLevels)) - - call work%init(nVertLevels + 1, nTracers, opts) ! this is internal workspace for PPR - - allocate(bcUpper(nTracers)) - allocate(bcLower(nTracers)) - ! TODO replace hardcoded bcon_loose with config option - bcUpper%bcopt = bcon_loose ! "loose" => extrapolate - bcLower%bcopt = bcon_loose - - do iCell = 1, nCells - do iTrac = 1, nTracers - nLevels = maxLevelCell(iCell) + 1 - - if (nLevels .lt. 2) cycle - - dataNow(1, iTrac, :) = tracersGroup(iTrac, :, iCell) - end do - call rmap1d(nLevels, nLevels, nTracers,nDoFs, & - heightCellNow(:, iCell), & - heightCellNew(:, iCell), & - dataNow, dataNew, bcUpper, bcLower, work, opts) - - do iTrac = 1, nTracers - tracersGroup(iTrac, 1:maxLevelCell(iCell), iCell) = dataNew(1, iTrac, 1:maxLevelCell(iCell)) - end do - - end do - deallocate(dataNow,dataNew) - deallocate(bcUpper,bcLower) - call work%free() - - end if - end if - end do - - !- dealloc any local data, etc - - deallocate(heightCellNow) - deallocate(heightEdgeNow) - deallocate(heightCellNew) - deallocate(heightEdgeNew) - - end subroutine ocn_remap_vert_state - - - end module ocn_time_integration_rk4 ! vim: foldmethod=marker diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index 09e290ad5f69..af7e7b06bc5e 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -40,6 +40,7 @@ module ocn_time_integration_split use ocn_equation_of_state use ocn_vmix + use ocn_vertical_remap use ocn_time_average_coupled use ocn_effective_density_in_land_ice @@ -83,10 +84,11 @@ module ocn_time_integration_split ! outer timestep iteration integer :: & - neededHalos, &! number of halo levels needed - haloDecrement, &! number of halos to invalidate each cycle - numTSIterations, &! number of outer timestep iterations - nBtrSubcycles ! number of barotropic subcycles + configVertAdvMethod, &! type of vertical advection method + neededHalos, &! number of halo levels needed + haloDecrement, &! number of halos to invalidate each cycle + numTSIterations, &! number of outer timestep iterations + nBtrSubcycles ! number of barotropic subcycles integer :: ssh_sal_on @@ -263,6 +265,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersNew + ! Remap variables + real (kind=RKIND) :: totalThickness + real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessRemap + real (kind=RKIND), dimension(:,:), pointer :: restingThickness + ! These are public variables used from the diagnostics module ! indexSurfaceVelocityZonal, indexSurfaceVelocityMeridional ! indexSSHGradientZonal, indexSSHGradientMeridional @@ -2431,6 +2438,35 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call mpas_timer_stop("se implicit vert mix") + if ( configVertAdvMethod == 2 ) then + + call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) + + allocate(layerThicknessRemap(nVertLevels, nCellsAll)) + layerThicknessRemap = layerThicknessNew + + ! calc. z*-style expansions + !$omp parallel + !$omp do schedule(runtime) private(k,totalThickness) + do iCell = 1, nCellsAll + totalThickness = 0.0_RKIND + do k = maxLevelCell(iCell), 1, -1 + totalThickness = totalThickness + layerThicknessNew(k,iCell) + end do + do k = maxLevelCell(iCell), 1, -1 + layerThicknessRemap(k, iCell) = restingThickness(k, iCell) * & + totalThickness / bottomDepth(iCell) + end do + + end do + !$omp end do + !$omp end parallel + call ocn_remap_vert_state(block, layerThicknessRemap, err) + + deallocate(layerThicknessRemap) + + end if + call mpas_timer_start('se fini') #ifdef MPAS_OPENACC @@ -2721,6 +2757,12 @@ subroutine ocn_time_integration_split_init(domain)!{{{ useVelocityCorrection = 0.0_RKIND endif + if ( config_vert_advection_method == 'flux-form' ) then + configVertAdvMethod = 1 + elseif ( config_vert_advection_method == 'remap' ) then + configVertAdvMethod = 2 + endif + !*** Compute some halo needs to reduce the number !*** of halo updates during subcycling neededHalos = 1 + config_n_btr_cor_iter diff --git a/components/mpas-ocean/src/ocean.cmake b/components/mpas-ocean/src/ocean.cmake index 8422b44b05b7..422116ef79ab 100644 --- a/components/mpas-ocean/src/ocean.cmake +++ b/components/mpas-ocean/src/ocean.cmake @@ -54,6 +54,7 @@ list(APPEND RAW_SOURCES core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F core_ocean/shared/mpas_ocn_vel_pressure_grad.F core_ocean/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F + core_ocean/shared/mpas_ocn_vertical_remap.F core_ocean/shared/mpas_ocn_vmix.F core_ocean/shared/mpas_ocn_vmix_coefs_redi.F core_ocean/shared/mpas_ocn_vmix_cvmix.F diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index 197069ea6e6d..db0a51226922 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -23,6 +23,7 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_vel_forcing_surface_stress.o \ mpas_ocn_vel_forcing_explicit_bottom_drag.o \ mpas_ocn_vel_pressure_grad.o \ + mpas_ocn_vertical_remap.o \ mpas_ocn_vmix.o \ mpas_ocn_vmix_coefs_redi.o \ mpas_ocn_vmix_cvmix.o \ diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F new file mode 100644 index 000000000000..8f23b91a2980 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -0,0 +1,353 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vertical_remap +! +!> \brief MPAS ocean vertical Lagrangian remapping +!> \author Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis +!> \date July 2021 +!> \details +!> This module contains the vertical remapping routine. +! +!----------------------------------------------------------------------- + +module ocn_vertical_remap + + use mpas_derived_types + use mpas_pool_routines + use mpas_constants + use mpas_dmpar + use mpas_threading + use mpas_vector_reconstruction + use mpas_spline_interpolation + use mpas_timer + + use ocn_constants + use ocn_config + + use ppr_1d + + ! TODO add use ocn_mesh + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + + type(rmap_opts) :: opts + + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_remap_vert_state + public :: ocn_vertical_remap_init + + contains + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_remap_vert_state +! +!> \brief MPAS ocean vertical Lagrangian remapping +!> \author Mark Petersen, Doug Jacobsen, Todd Ringler +!> \date September 2011 +!> \details +!> This routine remaps state variables from current layerThickness to +!> a new layerThickness +! +!----------------------------------------------------------------------- + + subroutine ocn_remap_vert_state(block, layerThicknessNew, err) + + type (block_type), intent(in) :: block + real (kind=RKIND), dimension(:, :), intent(in) :: layerThicknessNew + integer, intent(out) :: err + + integer, pointer :: nCells, nEdges, nVertLevels + integer :: nLevels, nVars, nDoFs, nTracers + integer :: iCell, jCell, iEdge, iTrac, k + + type (mpas_pool_type), pointer :: meshPool, verticalMeshPool + type (mpas_pool_type), pointer :: statePool + type (mpas_pool_type), pointer :: tracersPool + type (mpas_pool_iterator_type) :: groupItr + + type(rmap_work) :: work + type(rcon_ends), dimension(:), allocatable :: bcUpper, bcLower + + real (kind=RKIND), dimension(:, :, :), allocatable :: dataNow, dataNew + real (kind=RKIND), dimension(:, :), allocatable :: heightCellNow, heightEdgeNow + real (kind=RKIND), dimension(:, :), allocatable :: heightCellNew, heightEdgeNew + + integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot + integer, dimension(:, :), pointer :: cellsOnEdge + real (kind=RKIND), dimension(:), pointer :: bottomDepth + real (kind=RKIND), dimension(:, :), pointer :: layerThickness, restingThickness + real (kind=RKIND), dimension(:, :), pointer :: normalVelocity + real (kind=RKIND), dimension(:, :), pointer :: highFreqThickness,lowFreqDivergence + real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroup + + err = 0 + + call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) + call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges) + + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) + call mpas_pool_get_subpool(block % structs, 'state', statePool) + + call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) + + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', maxLevelEdgeBot) + call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + + call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) + + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 2) + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 2) + if (config_use_freq_filtered_thickness) then + call mpas_pool_get_array(statePool, 'highFreqThickness', highFreqThickness, 2) + call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergence, 2) + end if + + !---------------------------------------------------------------------------- + !- SETUP VERTICAL GRIDS: + !- heightCellNow, heightEdgeNow: depth at layer interfaces (lagrangian grid) + !- heightCellNew, heighyEdgeNew: depth at layer interfaces (new target grid) + + !- this should all go into a separate "regrid" routine + !- (eventually) this would handle regridding to "other" coordinates, rather + !- than just the z* approach here. + + allocate(heightCellNow(nVertLevels + 1, nCells)) + allocate(heightCellNew(nVertLevels + 1, nCells)) + allocate(heightEdgeNow(nVertLevels + 1, nEdges)) + allocate(heightEdgeNew(nVertLevels + 1, nEdges)) + + do iCell = 1, nCells + heightCellNow(:, iCell) = -bottomDepth(iCell) + heightCellNew(:, iCell) = -bottomDepth(iCell) + + !- reconstruct "now" heights + + do k = maxLevelCell(iCell), 1, -1 + heightCellNow(k, iCell) = & + heightCellNow(k + 1, iCell) + layerThickness(k, iCell) + end do + + !- reconstruct "new" heights + !- Assign new layerThicknesses + + heightCellNew(1, iCell) = heightCellNow(1, iCell) + do k = maxLevelCell(iCell), 2, -1 + heightCellNew(k, iCell) = & + heightCellNew(k + 1, iCell) + layerThicknessNew(k, iCell) + end do + end do + layerThickness = layerThicknessNew + + !- There's currently a silly issue in PPR: it assumes the "coordinates" for + !- the data being remapped increases with index. Here, we have the opposite + !- (z becomes more -ve with k)... + !- I'll make some updates to PPR to handle this, but for now, just flip the + !- signs to make PPR happy. This shouldn't have any actual impact. + + heightCellNow = -1.0_RKIND * heightCellNow + heightCellNew = -1.0_RKIND * heightCellNew + + + do iEdge = 1, nEdges + iCell = cellsOnEdge(1, iEdge) + jCell = cellsOnEdge(2, iEdge) + + heightEdgeNow(:, iEdge) = -bottomDepth(iCell) + heightEdgeNew(:, iEdge) = -bottomDepth(iCell) + + !- this is just a quick hack for the edge heights + !- really should use layerThicknessEdge here, but it's unclear if this has + !- been properly initiallised at this point... + !- layerThicknessEdge could involve upwinding, so need to be careful! + + do k = maxLevelEdgeTop(iEdge), 1, -1 + heightEdgeNow(k, iEdge) = & + 0.5_RKIND * (heightCellNow(k, iCell) + heightCellNow(k, jCell)) + heightEdgeNew(k, iEdge) = & + 0.5_RKIND * (heightCellNew(k, iCell) + heightCellNew(k, jCell)) + end do + end do + + !---------------------------------------------------------------------------- + !- ACTUAL REMAPPING FROM HERE + + !- just simple hard-coded options for PPR at this stage... + + if (associated(highFreqThickness) .and. associated(lowFreqDivergence)) then + nVars = 3 ! how many tracers to be remapped concurrently per column + else + nVars = 1 ! how many tracers to be remapped concurrently per column + endif + nDoFs = 1 ! always 1 for MPAS, it'd be >1 for DG methods (with more DoF per layer) + + allocate(dataNow(nDoFs, nVars, nVertLevels)) + allocate(dataNew(nDoFs, nVars, nVertLevels)) + allocate(bcUpper(nVars)) + allocate(bcLower(nVars)) + + dataNow = 0.0_RKIND + dataNew = 0.0_RKIND + + bcUpper%bcopt = bcon_loose ! "loose" => extrapolate + bcLower%bcopt = bcon_loose + + call work%init(nVertLevels + 1, nVars, opts) ! this is internal workspace for PPR + + !- Remap all edge-centred variables from heightEdgeNow to heightEdgeNew + if (associated(highFreqThickness) .and. associated(lowFreqDivergence)) then + do iEdge = 1, nEdges + nLevels = maxLevelEdgeTop(iEdge) + 1 + + if (nLevels .lt. 2) cycle + + dataNow(1, 1, :) = normalVelocity(:, iEdge) + dataNow(1, 2, :) = highFreqThickness(:, iEdge) + dataNow(1, 3, :) = lowFreqDivergence(:, iEdge) + + call rmap1d(nLevels, nLevels, nVars ,nDoFs, & + heightEdgeNow(:, iEdge), & + heightEdgeNew(:, iEdge), & + dataNow, dataNew, bcUpper, bcLower, work, opts) + + normalVelocity (1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 1, 1:maxLevelEdgeTop(iEdge)) + highFreqThickness(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 2, 1:maxLevelEdgeTop(iEdge)) + lowFreqDivergence(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 3, 1:maxLevelEdgeTop(iEdge)) + + end do + else + do iEdge = 1, nEdges + nLevels = maxLevelEdgeTop(iEdge) + 1 + + if (nLevels .lt. 2) cycle + + dataNow(1, 1, :) = normalVelocity(:, iEdge) + + call rmap1d(nLevels, nLevels, nVars ,nDoFs, & + heightEdgeNow(:, iEdge), & + heightEdgeNew(:, iEdge), & + dataNow, dataNew, bcUpper, bcLower, work, opts) + + normalVelocity(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 1, 1:maxLevelEdgeTop(iEdge)) + + end do + end if + deallocate(dataNow,dataNew) + deallocate(bcUpper,bcLower) + call work%free() + + !- Remap all cell-centred variables from heightCellNow to heightCellNew + + call mpas_pool_begin_iteration(tracersPool) + do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) + if ( groupItr % memberType == MPAS_POOL_FIELD ) then + call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroup, 2) + if ( associated(tracersGroup) ) then + + !- we can remap all of the tracers in one go via PPR, but for now + !- just a simple loop over each one at a time for clarity... + + nTracers = size(tracersGroup, dim=1) + + allocate(dataNow(nDoFs, nTracers, nVertLevels)) + allocate(dataNew(nDoFs, nTracers, nVertLevels)) + + call work%init(nVertLevels + 1, nTracers, opts) ! this is internal workspace for PPR + + allocate(bcUpper(nTracers)) + allocate(bcLower(nTracers)) + ! TODO replace hardcoded bcon_loose with config option + bcUpper%bcopt = bcon_loose ! "loose" => extrapolate + bcLower%bcopt = bcon_loose + + do iCell = 1, nCells + do iTrac = 1, nTracers + nLevels = maxLevelCell(iCell) + 1 + + if (nLevels .lt. 2) cycle + + dataNow(1, iTrac, :) = tracersGroup(iTrac, :, iCell) + end do + call rmap1d(nLevels, nLevels, nTracers,nDoFs, & + heightCellNow(:, iCell), & + heightCellNew(:, iCell), & + dataNow, dataNew, bcUpper, bcLower, work, opts) + + do iTrac = 1, nTracers + tracersGroup(iTrac, 1:maxLevelCell(iCell), iCell) = dataNew(1, iTrac, 1:maxLevelCell(iCell)) + end do + + end do + + deallocate(dataNow,dataNew) + deallocate(bcUpper,bcLower) + call work%free() + + end if + end if + end do + + !- dealloc any local data, etc + + deallocate(heightCellNow) + deallocate(heightEdgeNow) + deallocate(heightCellNew) + deallocate(heightEdgeNew) + + end subroutine ocn_remap_vert_state + + subroutine ocn_vertical_remap_init(err) + + integer, intent(out) :: err !< Output: Error flag + + err = 0 + + if ( config_time_integrator == 'RK4' ) then + CALL mpas_log_write('Vertical remap not supported for RK4', & + MPAS_LOG_CRIT) + endif + + ! Options for remapping + if ( config_vert_remap_order == 1 ) then + opts%cell_meth = pcm_method ! PCM method + elseif ( config_vert_remap_order == 2 ) then + opts%cell_meth = plm_method ! PLM method + elseif ( config_vert_remap_order == 3 ) then + opts%cell_meth = ppm_method ! PPM method + opts%edge_meth = p3e_method ! 3rd-order edge interp. + elseif ( config_vert_remap_order == 5 ) then + opts%cell_meth = pqm_method ! PPM method + opts%edge_meth = p5e_method ! 3rd-order edge interp. + endif + opts%cell_lims = mono_limit ! monotone slope limits TODO: consider config_monotonic + + end subroutine ocn_vertical_remap_init + +end module ocn_vertical_remap +! vim: foldmethod=marker From 7f95f874154d43bae84b4278944ae9ac5ce12efe Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 9 Jul 2021 13:21:13 -0500 Subject: [PATCH 06/22] Treatment of verticalVelocityTop for remap cases: - vertAleTransportTop = 0 so no vertical velocity during momentum solve - Add computation of vertVelocity from lagrangian layerThickness to diagnostics routines --- components/mpas-ocean/src/Registry.xml | 6 ++ .../src/driver/mpas_ocn_core_interface.F | 11 +++ .../mpas_ocn_time_integration_split.F | 12 ++- .../src/shared/mpas_ocn_diagnostics.F | 94 ++++++++++++++++++- 4 files changed, 119 insertions(+), 4 deletions(-) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 8211abe45e51..1a368f8a09e2 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1482,6 +1482,7 @@ + @@ -2181,6 +2182,11 @@ description="GOTM: turbulent length scale defined at the cell center (horizontally) and top (vertically)" packages="gotmPKG" /> + + \brief Computes diagnostic variables for velocity +!> \author Carolyn Begeman +!> \date October 2021 +!> \details +!> This routine computes the Eulerian vertical velocity valid for +!> Vertical Lagrangian Remapping +! +!----------------------------------------------------------------------- + subroutine ocn_diagnostic_solve_vertVel_remap(layerThickness, & + layerThicknessLag, dt, vertVelocityTop)!{{{ + + use ocn_mesh + + implicit none + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness ! layerThickness after remapping at timeLevelIn + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThicknessLag ! layerThickness before remapping at timeLevelIn + real (kind=RKIND), intent(in) :: dt + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + vertVelocityTop + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iCell, k + real (kind=RKIND), dimension(:), allocatable :: & + heightCell, heightCellLag + + allocate(heightCell(nVertLevels+1)) + allocate(heightCellLag(nVertLevels+1)) + + !$omp parallel + !$omp do schedule(runtime) private(k, heightCell, heightCellLag) + do iCell = 1, nCellsAll + + heightCell(:) = -1.0_RKIND * bottomDepth(iCell) + heightCellLag(:) = -1.0_RKIND * bottomDepth(iCell) + + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 + ! reconstruct heights + heightCellLag(k) = heightCellLag(k + 1) + layerThicknessLag(k, iCell) + heightCell(k) = heightCell(k + 1) + layerThickness(k, iCell) + end do + + vertVelocityTop(:, iCell) = 0.0_RKIND + + do k = minLevelCell(iCell), maxLevelCell(iCell) + ! compute vertical velocity at the top of each Cell + vertVelocityTop(k,iCell) = (heightCellLag(k) - heightCell(k))/dt + end do + end do + !$omp end do + !$omp end parallel + + deallocate(heightCell) + deallocate(heightCellLag) + + end subroutine ocn_diagnostic_solve_vertVel_remap!}}} + !*********************************************************************** ! ! routine ocn_diagnostic_solve_surface_pressure @@ -2209,7 +2298,8 @@ subroutine ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, oldLayerT err = 0 - if (config_vert_coord_movement.eq.'impermeable_interfaces') then + if (config_vert_coord_movement.eq.'impermeable_interfaces' .or. & + config_vert_advection_method.eq.'remap') then vertAleTransportTop=0.0_RKIND !$acc update device(vertAleTransportTop) return From beae6b8e241f2f4d281fd64506116b9d01f92443 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 9 Jul 2021 14:13:38 -0500 Subject: [PATCH 07/22] Add output file for debugging VLR --- components/mpas-ocean/src/Registry.xml | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 1a368f8a09e2..c1833cf2df68 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1692,6 +1692,25 @@ + + + + + + + + + + + Date: Tue, 13 Jul 2021 15:12:16 -0500 Subject: [PATCH 08/22] Fixup vert remap coord, incl minLevelCell --- .../mpas_ocn_time_integration_split.F | 4 +- .../src/shared/mpas_ocn_vertical_remap.F | 162 ++++++++++-------- 2 files changed, 91 insertions(+), 75 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index f620e4515553..d44a0e253afb 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -2459,10 +2459,10 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ !$omp do schedule(runtime) private(k,totalThickness) do iCell = 1, nCellsAll totalThickness = 0.0_RKIND - do k = maxLevelCell(iCell), 1, -1 + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 totalThickness = totalThickness + layerThicknessNew(k,iCell) end do - do k = maxLevelCell(iCell), 1, -1 + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 layerThicknessRemap(k, iCell) = restingThickness(k, iCell) * & totalThickness / bottomDepth(iCell) end do diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F index 8f23b91a2980..18f32f0b2bf2 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -30,11 +30,10 @@ module ocn_vertical_remap use ocn_constants use ocn_config + use ocn_mesh use ppr_1d - ! TODO add use ocn_mesh - implicit none private save @@ -45,6 +44,7 @@ module ocn_vertical_remap ! type(rmap_opts) :: opts + integer :: bc_upper, bc_lower !-------------------------------------------------------------------- @@ -78,11 +78,11 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) real (kind=RKIND), dimension(:, :), intent(in) :: layerThicknessNew integer, intent(out) :: err - integer, pointer :: nCells, nEdges, nVertLevels - integer :: nLevels, nVars, nDoFs, nTracers - integer :: iCell, jCell, iEdge, iTrac, k + integer :: nCells, nEdges + integer :: nLayers, nLevels, nVars, nDoFs, nTracers + integer :: iCell, jCell, iEdge, iTrac, k, kTop, kBot - type (mpas_pool_type), pointer :: meshPool, verticalMeshPool + type (mpas_pool_type), pointer :: verticalMeshPool type (mpas_pool_type), pointer :: statePool type (mpas_pool_type), pointer :: tracersPool type (mpas_pool_iterator_type) :: groupItr @@ -94,9 +94,6 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) real (kind=RKIND), dimension(:, :), allocatable :: heightCellNow, heightEdgeNow real (kind=RKIND), dimension(:, :), allocatable :: heightCellNew, heightEdgeNew - integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot - integer, dimension(:, :), pointer :: cellsOnEdge - real (kind=RKIND), dimension(:), pointer :: bottomDepth real (kind=RKIND), dimension(:, :), pointer :: layerThickness, restingThickness real (kind=RKIND), dimension(:, :), pointer :: normalVelocity real (kind=RKIND), dimension(:, :), pointer :: highFreqThickness,lowFreqDivergence @@ -104,22 +101,15 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) err = 0 - call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) - call mpas_pool_get_dimension(block % dimensions, 'nEdges', nEdges) - - call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + ! Remapping currently only supports one block for the whole domain rather + ! than sub-blocks + call mpas_log_write('Entering vertical remapping') + call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) call mpas_pool_get_subpool(block % structs, 'state', statePool) call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) - call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) - call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', maxLevelEdgeBot) - call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 2) @@ -129,6 +119,9 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergence, 2) end if + nCells = nCellsAll + nEdges = nEdgesAll + !---------------------------------------------------------------------------- !- SETUP VERTICAL GRIDS: !- heightCellNow, heightEdgeNow: depth at layer interfaces (lagrangian grid) @@ -143,51 +136,48 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) allocate(heightEdgeNow(nVertLevels + 1, nEdges)) allocate(heightEdgeNew(nVertLevels + 1, nEdges)) + call mpas_log_write(' -- cells') + !- There's currently a silly issue in PPR: it assumes the "coordinates" for + !- the data being remapped increases with index. Here, we have the opposite + !- (z becomes more -ve with k)... + !- I'll make some updates to PPR to handle this, but for now, just flip the + !- signs to make PPR happy. This shouldn't have any actual impact. do iCell = 1, nCells - heightCellNow(:, iCell) = -bottomDepth(iCell) - heightCellNew(:, iCell) = -bottomDepth(iCell) + heightCellNow(nVertLevels+1:maxLevelCell(iCell)+1, iCell) = bottomDepth(iCell) + heightCellNew(nVertLevels+1:maxLevelCell(iCell)+1, iCell) = bottomDepth(iCell) !- reconstruct "now" heights - do k = maxLevelCell(iCell), 1, -1 + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 heightCellNow(k, iCell) = & - heightCellNow(k + 1, iCell) + layerThickness(k, iCell) + heightCellNow(k + 1, iCell) - layerThickness(k, iCell) end do !- reconstruct "new" heights !- Assign new layerThicknesses - heightCellNew(1, iCell) = heightCellNow(1, iCell) - do k = maxLevelCell(iCell), 2, -1 + heightCellNew(minLevelCell(iCell), iCell) = heightCellNow(minLevelCell(iCell), iCell) + do k = maxLevelCell(iCell), minLevelCell(iCell)+1, -1 heightCellNew(k, iCell) = & - heightCellNew(k + 1, iCell) + layerThicknessNew(k, iCell) + heightCellNew(k + 1, iCell) - layerThicknessNew(k, iCell) end do end do layerThickness = layerThicknessNew - !- There's currently a silly issue in PPR: it assumes the "coordinates" for - !- the data being remapped increases with index. Here, we have the opposite - !- (z becomes more -ve with k)... - !- I'll make some updates to PPR to handle this, but for now, just flip the - !- signs to make PPR happy. This shouldn't have any actual impact. - - heightCellNow = -1.0_RKIND * heightCellNow - heightCellNew = -1.0_RKIND * heightCellNew - - + call mpas_log_write(' -- edges') do iEdge = 1, nEdges iCell = cellsOnEdge(1, iEdge) jCell = cellsOnEdge(2, iEdge) - heightEdgeNow(:, iEdge) = -bottomDepth(iCell) - heightEdgeNew(:, iEdge) = -bottomDepth(iCell) + heightEdgeNow(:, iEdge) = bottomDepth(iCell) + heightEdgeNew(:, iEdge) = bottomDepth(iCell) !- this is just a quick hack for the edge heights !- really should use layerThicknessEdge here, but it's unclear if this has !- been properly initiallised at this point... !- layerThicknessEdge could involve upwinding, so need to be careful! - do k = maxLevelEdgeTop(iEdge), 1, -1 + do k = maxLevelEdgeTop(iEdge), minLevelEdgeBot(iEdge), -1 heightEdgeNow(k, iEdge) = & 0.5_RKIND * (heightCellNow(k, iCell) + heightCellNow(k, jCell)) heightEdgeNew(k, iEdge) = & @@ -207,54 +197,65 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) endif nDoFs = 1 ! always 1 for MPAS, it'd be >1 for DG methods (with more DoF per layer) - allocate(dataNow(nDoFs, nVars, nVertLevels)) - allocate(dataNew(nDoFs, nVars, nVertLevels)) + allocate(dataNow(nDoFs, nVars, nVertLevels+1)) + allocate(dataNew(nDoFs, nVars, nVertLevels+1)) allocate(bcUpper(nVars)) allocate(bcLower(nVars)) dataNow = 0.0_RKIND dataNew = 0.0_RKIND - bcUpper%bcopt = bcon_loose ! "loose" => extrapolate - bcLower%bcopt = bcon_loose + bcUpper%bcopt = bc_upper ! "loose" => extrapolate + bcLower%bcopt = bc_lower call work%init(nVertLevels + 1, nVars, opts) ! this is internal workspace for PPR !- Remap all edge-centred variables from heightEdgeNow to heightEdgeNew if (associated(highFreqThickness) .and. associated(lowFreqDivergence)) then do iEdge = 1, nEdges - nLevels = maxLevelEdgeTop(iEdge) + 1 + + kTop = minLevelEdgeBot(iEdge) + kBot = maxLevelEdgeTop(iEdge) + + nLayers = kBot - kTop + 1 + nLevels = nLayers + 1 if (nLevels .lt. 2) cycle - dataNow(1, 1, :) = normalVelocity(:, iEdge) - dataNow(1, 2, :) = highFreqThickness(:, iEdge) - dataNow(1, 3, :) = lowFreqDivergence(:, iEdge) + dataNow(1, 1, 1:nLayers) = normalVelocity (kTop:kBot, iEdge) + dataNow(1, 2, 1:nLayers) = highFreqThickness(kTop:kBot, iEdge) + dataNow(1, 3, 1:nLayers) = lowFreqDivergence(kTop:kBot, iEdge) call rmap1d(nLevels, nLevels, nVars ,nDoFs, & - heightEdgeNow(:, iEdge), & - heightEdgeNew(:, iEdge), & - dataNow, dataNew, bcUpper, bcLower, work, opts) - - normalVelocity (1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 1, 1:maxLevelEdgeTop(iEdge)) - highFreqThickness(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 2, 1:maxLevelEdgeTop(iEdge)) - lowFreqDivergence(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 3, 1:maxLevelEdgeTop(iEdge)) - + heightEdgeNow(kTop:kBot+1, iEdge), & + heightEdgeNew(kTop:kBot+1, iEdge), & + dataNow, dataNew, & + bcUpper, bcLower, work, opts) + + normalVelocity (kTop:kBot, iEdge) = dataNew(1, 1, 1:nLayers) + highFreqThickness(kTop:kBot, iEdge) = dataNew(1, 2, 1:nLayers) + lowFreqDivergence(kTop:kBot, iEdge) = dataNew(1, 3, 1:nLayers) end do else do iEdge = 1, nEdges - nLevels = maxLevelEdgeTop(iEdge) + 1 + + kTop = minLevelEdgeBot(iEdge) + kBot = maxLevelEdgeTop(iEdge) + + nLayers = kBot - kTop + 1 + nLevels = nLayers + 1 if (nLevels .lt. 2) cycle - dataNow(1, 1, :) = normalVelocity(:, iEdge) + dataNow(1, 1, 1:nLayers) = normalVelocity(kTop:kBot, iEdge) call rmap1d(nLevels, nLevels, nVars ,nDoFs, & - heightEdgeNow(:, iEdge), & - heightEdgeNew(:, iEdge), & - dataNow, dataNew, bcUpper, bcLower, work, opts) + heightEdgeNow(kTop:kBot+1, iEdge), & + heightEdgeNew(kTop:kBot+1, iEdge), & + dataNow, dataNew, & + bcUpper, bcLower, work, opts) - normalVelocity(1:maxLevelEdgeTop(iEdge), iEdge) = dataNew(1, 1, 1:maxLevelEdgeTop(iEdge)) + normalVelocity (kTop:kBot, iEdge) = dataNew(1, 1, 1:nLayers) end do end if @@ -282,27 +283,34 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) allocate(bcUpper(nTracers)) allocate(bcLower(nTracers)) - ! TODO replace hardcoded bcon_loose with config option - bcUpper%bcopt = bcon_loose ! "loose" => extrapolate - bcLower%bcopt = bcon_loose + bcUpper%bcopt = bc_upper ! "loose" => extrapolate + bcLower%bcopt = bc_lower do iCell = 1, nCells + + kTop = minLevelCell(iCell) + kBot = maxLevelCell(iCell) + + nLayers = kBot - kTop + 1 + nLevels = nLayers + 1 + + if (nLevels .lt. 2) cycle + do iTrac = 1, nTracers - nLevels = maxLevelCell(iCell) + 1 - if (nLevels .lt. 2) cycle + dataNow(1, iTrac, 1:nLayers) = tracersGroup(iTrac, kTop:kBot, iCell) - dataNow(1, iTrac, :) = tracersGroup(iTrac, :, iCell) end do call rmap1d(nLevels, nLevels, nTracers,nDoFs, & - heightCellNow(:, iCell), & - heightCellNew(:, iCell), & - dataNow, dataNew, bcUpper, bcLower, work, opts) + heightCellNow(kTop:kBot+1, iCell), & + heightCellNew(kTop:kBot+1, iCell), & + dataNow, dataNew, & + bcUpper, bcLower, work, opts) do iTrac = 1, nTracers - tracersGroup(iTrac, 1:maxLevelCell(iCell), iCell) = dataNew(1, iTrac, 1:maxLevelCell(iCell)) + tracersGroup(iTrac, kTop:kBot, iCell) = dataNew(1, iTrac, 1:nLayers) end do - + end do deallocate(dataNow,dataNew) @@ -346,6 +354,14 @@ subroutine ocn_vertical_remap_init(err) opts%edge_meth = p5e_method ! 3rd-order edge interp. endif opts%cell_lims = mono_limit ! monotone slope limits TODO: consider config_monotonic + ! TODO replace hardcoded bcon_loose with config option + !select case (config_vert_remap_bc_upper) + !case ('extrapolate') + bc_upper = bcon_loose + !case ('dirichlet') + !bc_upper = bcon_value + bc_lower = bcon_loose + !bc_lower = bcon_value end subroutine ocn_vertical_remap_init From eb1cfd1f8d27d25db93073b405bd9cfd6d98f5c8 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 23 Jul 2021 12:38:59 -0600 Subject: [PATCH 09/22] Add omp calls and timers to VLR --- .../mpas_ocn_time_integration_split.F | 8 +- .../src/shared/mpas_ocn_vertical_remap.F | 80 ++++++++++++++----- 2 files changed, 65 insertions(+), 23 deletions(-) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index d44a0e253afb..bb9de5fe95fc 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -2442,7 +2442,9 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call mpas_timer_stop("se implicit vert mix") if ( configVertAdvMethod == 2 ) then - + + call mpas_timer_start("vertical remap") + call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) call mpas_pool_get_array(statePool, 'layerThicknessLag', layerThicknessLagNew, 2) @@ -2472,7 +2474,9 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call ocn_remap_vert_state(block, layerThicknessRemap, err) deallocate(layerThicknessRemap) - + + call mpas_timer_stop("vertical remap") + end if call mpas_timer_start('se fini') diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F index 18f32f0b2bf2..7330e16f16bf 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -101,9 +101,12 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) err = 0 + itimestepLastRemap = itimestepLastRemap + 1 + if (itimestepLastRemap < config_vert_remap_interval) return + itimestepLastRemap = 0 + ! Remapping currently only supports one block for the whole domain rather ! than sub-blocks - call mpas_log_write('Entering vertical remapping') call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) call mpas_pool_get_subpool(block % structs, 'state', statePool) @@ -136,15 +139,16 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) allocate(heightEdgeNow(nVertLevels + 1, nEdges)) allocate(heightEdgeNew(nVertLevels + 1, nEdges)) - call mpas_log_write(' -- cells') !- There's currently a silly issue in PPR: it assumes the "coordinates" for !- the data being remapped increases with index. Here, we have the opposite !- (z becomes more -ve with k)... !- I'll make some updates to PPR to handle this, but for now, just flip the !- signs to make PPR happy. This shouldn't have any actual impact. + !$omp parallel + !$omp do schedule(runtime) private(k) do iCell = 1, nCells - heightCellNow(nVertLevels+1:maxLevelCell(iCell)+1, iCell) = bottomDepth(iCell) - heightCellNew(nVertLevels+1:maxLevelCell(iCell)+1, iCell) = bottomDepth(iCell) + heightCellNow(:, iCell) = bottomDepth(iCell) + heightCellNew(:, iCell) = bottomDepth(iCell) !- reconstruct "now" heights @@ -162,9 +166,12 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) heightCellNew(k + 1, iCell) - layerThicknessNew(k, iCell) end do end do + !$omp end do + !$omp end parallel layerThickness = layerThicknessNew - call mpas_log_write(' -- edges') + !$omp parallel + !$omp do schedule(runtime) private(k, iCell, jCell) do iEdge = 1, nEdges iCell = cellsOnEdge(1, iEdge) jCell = cellsOnEdge(2, iEdge) @@ -184,34 +191,39 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) 0.5_RKIND * (heightCellNew(k, iCell) + heightCellNew(k, jCell)) end do end do + !$omp end do + !$omp end parallel - !---------------------------------------------------------------------------- - !- ACTUAL REMAPPING FROM HERE + !---------------------------------------------------------------------------- + ! ACTUAL REMAPPING FROM HERE - !- just simple hard-coded options for PPR at this stage... + nDoFs = 1 ! always 1 for MPAS, it'd be >1 for DG methods (with more DoF per layer) + ! how many variables to be remapped concurrently per column if (associated(highFreqThickness) .and. associated(lowFreqDivergence)) then - nVars = 3 ! how many tracers to be remapped concurrently per column + nVars = 3 else - nVars = 1 ! how many tracers to be remapped concurrently per column + nVars = 1 endif - nDoFs = 1 ! always 1 for MPAS, it'd be >1 for DG methods (with more DoF per layer) - + allocate(dataNow(nDoFs, nVars, nVertLevels+1)) allocate(dataNew(nDoFs, nVars, nVertLevels+1)) allocate(bcUpper(nVars)) allocate(bcLower(nVars)) + ! set boundary conditions + bcUpper%bcopt = bc_upper + bcLower%bcopt = bc_lower + dataNow = 0.0_RKIND dataNew = 0.0_RKIND - bcUpper%bcopt = bc_upper ! "loose" => extrapolate - bcLower%bcopt = bc_lower - call work%init(nVertLevels + 1, nVars, opts) ! this is internal workspace for PPR - !- Remap all edge-centred variables from heightEdgeNow to heightEdgeNew + ! Remap all edge-centred variables from heightEdgeNow to heightEdgeNew if (associated(highFreqThickness) .and. associated(lowFreqDivergence)) then + !$omp parallel + !$omp do schedule(runtime) private(iEdge, kTop, kBot, nLayers, nLevels, dataNow, dataNew) do iEdge = 1, nEdges kTop = minLevelEdgeBot(iEdge) @@ -226,17 +238,23 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) dataNow(1, 2, 1:nLayers) = highFreqThickness(kTop:kBot, iEdge) dataNow(1, 3, 1:nLayers) = lowFreqDivergence(kTop:kBot, iEdge) + call mpas_timer_start("ppr") call rmap1d(nLevels, nLevels, nVars ,nDoFs, & heightEdgeNow(kTop:kBot+1, iEdge), & heightEdgeNew(kTop:kBot+1, iEdge), & dataNow, dataNew, & bcUpper, bcLower, work, opts) + call mpas_timer_stop("ppr") normalVelocity (kTop:kBot, iEdge) = dataNew(1, 1, 1:nLayers) highFreqThickness(kTop:kBot, iEdge) = dataNew(1, 2, 1:nLayers) lowFreqDivergence(kTop:kBot, iEdge) = dataNew(1, 3, 1:nLayers) end do + !$omp end do + !$omp end parallel else + !$omp parallel + !$omp do schedule(runtime) private(iEdge, kTop, kBot, nLayers, nLevels, dataNow, dataNew) do iEdge = 1, nEdges kTop = minLevelEdgeBot(iEdge) @@ -249,21 +267,25 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) dataNow(1, 1, 1:nLayers) = normalVelocity(kTop:kBot, iEdge) + call mpas_timer_start("ppr") call rmap1d(nLevels, nLevels, nVars ,nDoFs, & heightEdgeNow(kTop:kBot+1, iEdge), & heightEdgeNew(kTop:kBot+1, iEdge), & dataNow, dataNew, & bcUpper, bcLower, work, opts) + call mpas_timer_stop("ppr") normalVelocity (kTop:kBot, iEdge) = dataNew(1, 1, 1:nLayers) end do + !$omp end do + !$omp end parallel end if deallocate(dataNow,dataNew) deallocate(bcUpper,bcLower) call work%free() - !- Remap all cell-centred variables from heightCellNow to heightCellNew + ! Remap all cell-centred variables from heightCellNow to heightCellNew call mpas_pool_begin_iteration(tracersPool) do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) @@ -271,9 +293,6 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroup, 2) if ( associated(tracersGroup) ) then - !- we can remap all of the tracers in one go via PPR, but for now - !- just a simple loop over each one at a time for clarity... - nTracers = size(tracersGroup, dim=1) allocate(dataNow(nDoFs, nTracers, nVertLevels)) @@ -283,9 +302,11 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) allocate(bcUpper(nTracers)) allocate(bcLower(nTracers)) - bcUpper%bcopt = bc_upper ! "loose" => extrapolate + bcUpper%bcopt = bc_upper bcLower%bcopt = bc_lower + !$omp parallel + !$omp do schedule(runtime) private(iCell, iTrac, kTop, kBot, nLayers, nLevels, dataNow, dataNew) do iCell = 1, nCells kTop = minLevelCell(iCell) @@ -301,17 +322,21 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) dataNow(1, iTrac, 1:nLayers) = tracersGroup(iTrac, kTop:kBot, iCell) end do + call mpas_timer_start("ppr") call rmap1d(nLevels, nLevels, nTracers,nDoFs, & heightCellNow(kTop:kBot+1, iCell), & heightCellNew(kTop:kBot+1, iCell), & dataNow, dataNew, & bcUpper, bcLower, work, opts) + call mpas_timer_stop("ppr") do iTrac = 1, nTracers tracersGroup(iTrac, kTop:kBot, iCell) = dataNew(1, iTrac, 1:nLayers) end do end do + !$omp end do + !$omp end parallel deallocate(dataNow,dataNew) deallocate(bcUpper,bcLower) @@ -330,6 +355,19 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) end subroutine ocn_remap_vert_state +!*********************************************************************** +! +! routine ocn_vertical_remap_init +! +!> \brief Initializes ocean vertical remapping +!> \author Carolyn Begeman +!> \date July 2021 +!> \details +!> This routine initializes parameters required for vertical Lagrangian +!> remapping +! +!----------------------------------------------------------------------- + subroutine ocn_vertical_remap_init(err) integer, intent(out) :: err !< Output: Error flag From 33698dbbe46352ed89c66184efb19586db0729fe Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 23 Jul 2021 12:39:23 -0600 Subject: [PATCH 10/22] Use layerThicknessEdge from diags for remapping --- components/mpas-ocean/src/shared/Makefile | 2 + .../src/shared/mpas_ocn_diagnostics.F | 1 + .../src/shared/mpas_ocn_vertical_remap.F | 98 +++++++++++-------- 3 files changed, 58 insertions(+), 43 deletions(-) diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index db0a51226922..0a4bed12fe3c 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -218,6 +218,8 @@ mpas_ocn_vel_pressure_grad.o: mpas_ocn_constants.o mpas_ocn_vel_tidal_potential.o: mpas_ocn_diagnostics_variables.o +mpas_ocn_vertical_remap.o: mpas_ocn_diagnostics.o mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_mesh.o + clean: $(RM) *.o *.i *.mod *.f90 diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index feb602cd3f2b..97d5bd35029c 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -55,6 +55,7 @@ module ocn_diagnostics !-------------------------------------------------------------------- public :: ocn_diagnostic_solve, & + ocn_diagnostic_solve_layerThicknessEdge, & ocn_relativeVorticity_circulation, & ocn_vert_transport_velocity_top, & ocn_fuperp, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F index 7330e16f16bf..1dd29557d734 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -31,6 +31,7 @@ module ocn_vertical_remap use ocn_constants use ocn_config use ocn_mesh + use ocn_diagnostics use ppr_1d @@ -80,7 +81,7 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) integer :: nCells, nEdges integer :: nLayers, nLevels, nVars, nDoFs, nTracers - integer :: iCell, jCell, iEdge, iTrac, k, kTop, kBot + integer :: cell1, cell2, iCell, jCell, iEdge, iTrac, k, kTop, kBot type (mpas_pool_type), pointer :: verticalMeshPool type (mpas_pool_type), pointer :: statePool @@ -91,8 +92,11 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) type(rcon_ends), dimension(:), allocatable :: bcUpper, bcLower real (kind=RKIND), dimension(:, :, :), allocatable :: dataNow, dataNew - real (kind=RKIND), dimension(:, :), allocatable :: heightCellNow, heightEdgeNow - real (kind=RKIND), dimension(:, :), allocatable :: heightCellNew, heightEdgeNew + real (kind=RKIND), dimension(:, :), allocatable :: & + layerThickEdge, & ! layerThickness at edges (lagrangian grid) + layerThickEdgeNew, & ! layerThickness at edges (new target grid) + heightCellNow, heightEdgeNow, & ! depth at layer interfaces (lagrangian grid) + heightCellNew, heightEdgeNew ! depth at layer interfaces (new target grid) real (kind=RKIND), dimension(:, :), pointer :: layerThickness, restingThickness real (kind=RKIND), dimension(:, :), pointer :: normalVelocity @@ -107,7 +111,6 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) ! Remapping currently only supports one block for the whole domain rather ! than sub-blocks - call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) call mpas_pool_get_subpool(block % structs, 'state', statePool) @@ -125,41 +128,35 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) nCells = nCellsAll nEdges = nEdgesAll - !---------------------------------------------------------------------------- - !- SETUP VERTICAL GRIDS: - !- heightCellNow, heightEdgeNow: depth at layer interfaces (lagrangian grid) - !- heightCellNew, heighyEdgeNew: depth at layer interfaces (new target grid) - - !- this should all go into a separate "regrid" routine - !- (eventually) this would handle regridding to "other" coordinates, rather - !- than just the z* approach here. + ! SETUP VERTICAL GRIDS: allocate(heightCellNow(nVertLevels + 1, nCells)) allocate(heightCellNew(nVertLevels + 1, nCells)) allocate(heightEdgeNow(nVertLevels + 1, nEdges)) allocate(heightEdgeNew(nVertLevels + 1, nEdges)) + allocate(layerThickEdge(nVertLevels, nEdges)) + allocate(layerThickEdgeNew(nVertLevels, nEdges)) - !- There's currently a silly issue in PPR: it assumes the "coordinates" for - !- the data being remapped increases with index. Here, we have the opposite - !- (z becomes more -ve with k)... - !- I'll make some updates to PPR to handle this, but for now, just flip the - !- signs to make PPR happy. This shouldn't have any actual impact. + + ! There's currently a silly issue in PPR: it assumes the "coordinates" for + ! the data being remapped increases with index. Here, we have the opposite + ! (z becomes more -ve with k)... + ! I'll make some updates to PPR to handle this, but for now, just flip the + ! signs to make PPR happy. This shouldn't have any actual impact. + !$omp parallel !$omp do schedule(runtime) private(k) do iCell = 1, nCells heightCellNow(:, iCell) = bottomDepth(iCell) heightCellNew(:, iCell) = bottomDepth(iCell) - !- reconstruct "now" heights - + ! reconstruct "now" heights do k = maxLevelCell(iCell), minLevelCell(iCell), -1 heightCellNow(k, iCell) = & heightCellNow(k + 1, iCell) - layerThickness(k, iCell) end do - !- reconstruct "new" heights - !- Assign new layerThicknesses - + ! reconstruct "new" heights heightCellNew(minLevelCell(iCell), iCell) = heightCellNow(minLevelCell(iCell), iCell) do k = maxLevelCell(iCell), minLevelCell(iCell)+1, -1 heightCellNew(k, iCell) = & @@ -168,32 +165,46 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) end do !$omp end do !$omp end parallel - layerThickness = layerThicknessNew + + call ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, layerThickness, layerThickEdge) + call ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, layerThicknessNew, layerThickEdgeNew) !$omp parallel - !$omp do schedule(runtime) private(k, iCell, jCell) + !$omp do schedule(runtime) private(k,kTop,kBot,cell1,cell2) do iEdge = 1, nEdges - iCell = cellsOnEdge(1, iEdge) - jCell = cellsOnEdge(2, iEdge) - - heightEdgeNow(:, iEdge) = bottomDepth(iCell) - heightEdgeNew(:, iEdge) = bottomDepth(iCell) - - !- this is just a quick hack for the edge heights - !- really should use layerThicknessEdge here, but it's unclear if this has - !- been properly initiallised at this point... - !- layerThicknessEdge could involve upwinding, so need to be careful! - - do k = maxLevelEdgeTop(iEdge), minLevelEdgeBot(iEdge), -1 - heightEdgeNow(k, iEdge) = & - 0.5_RKIND * (heightCellNow(k, iCell) + heightCellNow(k, jCell)) - heightEdgeNew(k, iEdge) = & - 0.5_RKIND * (heightCellNew(k, iCell) + heightCellNew(k, jCell)) + + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + kTop = minLevelEdgeBot(iEdge) + kBot = maxLevelEdgeTop(iEdge) + + ! We only need to define the kBot+1 depth for one active cell + ! as the absolute kBot+1 depth has no effect on remapping + if (cell1 <= nCellsAll .and. maxLevelCell(cell1) > 0) then + heightEdgeNow(:, iEdge) = heightCellNow(kBot+1, cell1) + heightEdgeNew(:, iEdge) = heightCellNew(kBot+1, cell1) + elseif (cell2 <= nCellsAll .and. maxLevelCell(cell2) > 0) then + heightEdgeNow(:, iEdge) = heightCellNow(kBot+1, cell2) + heightEdgeNew(:, iEdge) = heightCellNew(kBot+1, cell2) + else + call mpas_log_write('No valid cell for edge remapping', MPAS_LOG_WARN) + heightEdgeNow(:,iEdge) = 0.0_RKIND + cycle + endif + + do k = kBot, kTop, -1 + heightEdgeNow(k, iEdge) = heightEdgeNow(k+1, iEdge) - layerThickEdge(k, iEdge) + heightEdgeNew(k, iEdge) = heightEdgeNew(k+1, iEdge) - layerThickEdgeNew(k, iEdge) end do + end do !$omp end do !$omp end parallel + ! Assign new layerThicknesses + layerThickness = layerThicknessNew + !---------------------------------------------------------------------------- ! ACTUAL REMAPPING FROM HERE @@ -391,14 +402,15 @@ subroutine ocn_vertical_remap_init(err) opts%cell_meth = pqm_method ! PPM method opts%edge_meth = p5e_method ! 3rd-order edge interp. endif - opts%cell_lims = mono_limit ! monotone slope limits TODO: consider config_monotonic + opts%cell_lims = mono_limit ! monotone slope limits TODO: consider +onfig_monotonic ! TODO replace hardcoded bcon_loose with config option !select case (config_vert_remap_bc_upper) !case ('extrapolate') - bc_upper = bcon_loose + bc_upper = bcon_loose !case ('dirichlet') !bc_upper = bcon_value - bc_lower = bcon_loose + bc_lower = bcon_loose !bc_lower = bcon_value end subroutine ocn_vertical_remap_init From b3ec79f6be0834af31d5915190c95d32b9b25f0b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 31 Aug 2021 15:10:05 -0500 Subject: [PATCH 11/22] Changes to config options related to VLR: - Add option for WENO slope limiter - Fixup config name for vertical tracer adv order - Add option for remapping every interval number of timesteps - Update namelist entries - Add checks on config options --- components/mpas-ocean/bld/build-namelist | 9 +++-- .../namelist_defaults_mpaso.xml | 9 +++-- .../namelist_definition_mpaso.xml | 38 +++++++++++++++---- components/mpas-ocean/src/Registry.xml | 16 ++++++-- .../src/shared/mpas_ocn_tracer_advection.F | 9 ++++- .../shared/mpas_ocn_tracer_advection_vert.F | 2 +- .../src/shared/mpas_ocn_vertical_remap.F | 33 ++++++++++------ 7 files changed, 86 insertions(+), 30 deletions(-) diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index 5fee738461f8..64a922e9781d 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -768,11 +768,14 @@ add_default($nl, 'config_land_ice_flux_jenkins_salt_transfer_coefficient'); # Namelist group: advection # ############################# -add_default($nl, 'config_vert_tracer_adv'); -add_default($nl, 'config_vert_tracer_adv_order'); +add_default($nl, 'config_vert_advection_method'); +add_default($nl, 'config_vert_remap_order'); +add_default($nl, 'config_vert_remap_interval'); +add_default($nl, 'config_vert_tracer_adv_flux_order'); add_default($nl, 'config_horiz_tracer_adv_order'); add_default($nl, 'config_coef_3rd_order'); -add_default($nl, 'config_monotonic'); +add_default($nl, 'config_flux_limiter'); +add_default($nl, 'config_remap_limiter'); ############################### # Namelist group: bottom_drag # 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 789acecc6e54..a1c31bf79ab2 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -342,11 +342,14 @@ 8.42e-5 -'stencil' -3 +'flux-form' +3 +0 +3 3 0.25 -.true. +'monotonic' +'monotonic' .true. 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 6e0153ead725..0e2a2ae1b177 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -1582,15 +1582,31 @@ Default: Defined in namelist_defaults.xml - -Method for interpolating tracer values from layer centers to layer edges +Method for advecting tracers, momentum, and thickness vertically. -Valid values: 'spline' and 'stencil' +Valid values: 'flux-form' and 'remap' Default: Defined in namelist_defaults.xml - +Order of remapping method used for momentum and tracer advection + +Valid values: 1, 2, 3 and 5 +Default: Defined in namelist_defaults.xml + + + +Number of timesteps between each remapping. If 0, remapping occurs every timestep + +Valid values: Any integer greater than or equal to 0 +Default: Defined in namelist_defaults.xml + + + Order of polynomial used for tracer reconstruction at layer edges @@ -1614,11 +1630,19 @@ Valid values: any real between 0 and 1 Default: Defined in namelist_defaults.xml - -If .true. then fluxes are limited to produce a monotonic advection scheme +Slope limiter for the flux-form advection scheme. -Valid values: .true. and .false. +Valid values: 'none','monotonic' +Default: Defined in namelist_defaults.xml + + + +Slope limiter for the vertical remap advection scheme. + +Valid values: 'none','monotonic','weno' Default: Defined in namelist_defaults.xml diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index c1833cf2df68..fda2d17e651c 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1008,6 +1008,10 @@ description="Order of remapping method used for momentum and tracer advection" possible_values="1, 2, 3 and 5" /> + - + @@ -1344,7 +1352,7 @@ possible_values=".true. or .false." /> Date: Wed, 15 Jun 2022 19:20:21 -0500 Subject: [PATCH 12/22] Fix merge conflict in ocn_diagnostics --- .../src/shared/mpas_ocn_diagnostics.F | 41 +++++++++++++++++-- 1 file changed, 38 insertions(+), 3 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index 97d5bd35029c..c2051ef03506 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -111,8 +111,9 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, scratchPoo integer, intent(in), optional :: timeLevelIn !< Input: Time level in state logical, intent(in), optional :: full !< Input: Run all computations? - integer :: iEdge, iCell, iTracer - integer :: err, nCells, nEdges, nVertices + integer :: iEdge, iCell, iTracer, k, kmin, kmax + integer :: err, nCells, nEdges, nVertices, numTracers + integer :: timeLevel real (kind=RKIND) :: coef_3rd_order @@ -123,7 +124,6 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, scratchPoo real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers - integer :: timeLevel integer, pointer :: indexTemperature, indexSalinity logical :: full_compute = .true. @@ -240,11 +240,46 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, scratchPoo kineticEnergyCell, tangentialVelocity) if ( config_vert_advection_method == 'remap' ) then + ! Overwrite vertVelocityTop with exact Eulerian solution from remapping call mpas_pool_get_array(statePool, 'layerThicknessLag', & layerThicknessLag, timeLevel) call ocn_diagnostic_solve_vertVel_remap(layerThickness, & layerThicknessLag, dt, vertVelocityTop) + + ! Overwrite vertical tracer budget terms with new vertVelocityTop + if ( config_compute_active_tracer_budgets ) then + numTracers = size(activeTracers, dim=1) + do iTracer = 1, numTracers + !$omp parallel + !$omp do schedule(runtime) private(k,kmin,kmax) + do iCell = 1, nCells + kmin = minLevelCell(iCell) + kmax = maxLevelCell(iCell) + + activeTracerVerticalAdvectionTopFlux(iTracer, :, iCell) = 0.0_RKIND + do k = kmin+1, kmax + activeTracerVerticalAdvectionTopFlux(iTracer, k, iCell) = & + min(0.0_RKIND,vertVelocityTop(k, iCell))* & + activeTracers(iTracer, k-1, iCell) & + + max(0.0_RKIND,vertVelocityTop(k,iCell))* & + activeTracers(iTracer, k, iCell) + enddo + + do k = kmin, kmax-1 + activeTracerVerticalAdvectionTendency(iTracer, k, iCell) = ( & + activeTracerVerticalAdvectionTopFlux(iTracer, k+1, iCell) & + - activeTracerVerticalAdvectionTopFlux(iTracer, k, iCell) & + ) / layerThickness(k, iCell) + enddo + activeTracerVerticalAdvectionTendency(iTracer, kmax, iCell) = & + - activeTracerVerticalAdvectionTopFlux(iTracer, kmax, iCell) & + / layerThickness(kmax, iCell) + enddo + !$omp end do + !$omp end parallel + enddo + endif endif ! From 5a1c849dc7a8658074852538169fd980372d431f Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 3 Feb 2022 15:17:31 -0600 Subject: [PATCH 13/22] Aesthetic changes --- .../src/shared/mpas_ocn_vertical_remap.F | 50 +++++++++---------- 1 file changed, 24 insertions(+), 26 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F index cd0222c40db9..e21173210be5 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -43,11 +43,11 @@ module ocn_vertical_remap ! ! Public parameters ! - + type(rmap_opts) :: opts integer :: bc_upper, bc_lower integer :: itimestepLastRemap - + !-------------------------------------------------------------------- !-------------------------------------------------------------------- @@ -69,7 +69,7 @@ module ocn_vertical_remap !> \author Mark Petersen, Doug Jacobsen, Todd Ringler !> \date September 2011 !> \details -!> This routine remaps state variables from current layerThickness to +!> This routine remaps state variables from current layerThickness to !> a new layerThickness ! !----------------------------------------------------------------------- @@ -97,8 +97,8 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) layerThickEdge, & ! layerThickness at edges (lagrangian grid) layerThickEdgeNew, & ! layerThickness at edges (new target grid) heightCellNow, heightEdgeNow, & ! depth at layer interfaces (lagrangian grid) - heightCellNew, heightEdgeNew ! depth at layer interfaces (new target grid) - + heightCellNew, heightEdgeNew ! depth at layer interfaces (new target grid) + real (kind=RKIND), dimension(:, :), pointer :: layerThickness, restingThickness real (kind=RKIND), dimension(:, :), pointer :: normalVelocity real (kind=RKIND), dimension(:, :), pointer :: highFreqThickness,lowFreqDivergence @@ -118,14 +118,14 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) - + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 2) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 2) if (config_use_freq_filtered_thickness) then call mpas_pool_get_array(statePool, 'highFreqThickness', highFreqThickness, 2) call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergence, 2) end if - + nCells = nCellsAll nEdges = nEdgesAll @@ -137,11 +137,11 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) allocate(heightEdgeNew(nVertLevels + 1, nEdges)) allocate(layerThickEdge(nVertLevels, nEdges)) allocate(layerThickEdgeNew(nVertLevels, nEdges)) - + ! There's currently a silly issue in PPR: it assumes the "coordinates" for - ! the data being remapped increases with index. Here, we have the opposite - ! (z becomes more -ve with k)... + ! the data being remapped increases with index. Here, we have the opposite + ! (z becomes more -ve with k)... ! I'll make some updates to PPR to handle this, but for now, just flip the ! signs to make PPR happy. This shouldn't have any actual impact. @@ -212,17 +212,17 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) nDoFs = 1 ! always 1 for MPAS, it'd be >1 for DG methods (with more DoF per layer) ! how many variables to be remapped concurrently per column - if (associated(highFreqThickness) .and. associated(lowFreqDivergence)) then + if (config_use_freq_filtered_thickness) then nVars = 3 else nVars = 1 endif - + allocate(dataNow(nDoFs, nVars, nVertLevels+1)) allocate(dataNew(nDoFs, nVars, nVertLevels+1)) allocate(bcUpper(nVars)) allocate(bcLower(nVars)) - + ! set boundary conditions bcUpper%bcopt = bc_upper bcLower%bcopt = bc_lower @@ -233,7 +233,7 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) call work%init(nVertLevels + 1, nVars, opts) ! this is internal workspace for PPR ! Remap all edge-centred variables from heightEdgeNow to heightEdgeNew - if (associated(highFreqThickness) .and. associated(lowFreqDivergence)) then + if (config_use_freq_filtered_thickness) then !$omp parallel !$omp do schedule(runtime) private(iEdge, kTop, kBot, nLayers, nLevels, dataNow, dataNew) do iEdge = 1, nEdges @@ -288,7 +288,7 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) call mpas_timer_stop("ppr") normalVelocity (kTop:kBot, iEdge) = dataNew(1, 1, 1:nLayers) - + end do !$omp end do !$omp end parallel @@ -296,7 +296,7 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) deallocate(dataNow,dataNew) deallocate(bcUpper,bcLower) call work%free() - + ! Remap all cell-centred variables from heightCellNow to heightCellNew call mpas_pool_begin_iteration(tracersPool) @@ -304,19 +304,19 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) if ( groupItr % memberType == MPAS_POOL_FIELD ) then call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroup, 2) if ( associated(tracersGroup) ) then - + nTracers = size(tracersGroup, dim=1) - + allocate(dataNow(nDoFs, nTracers, nVertLevels)) allocate(dataNew(nDoFs, nTracers, nVertLevels)) - + call work%init(nVertLevels + 1, nTracers, opts) ! this is internal workspace for PPR allocate(bcUpper(nTracers)) allocate(bcLower(nTracers)) bcUpper%bcopt = bc_upper bcLower%bcopt = bc_lower - + !$omp parallel !$omp do schedule(runtime) private(iCell, iTrac, kTop, kBot, nLayers, nLevels, dataNow, dataNew) do iCell = 1, nCells @@ -345,7 +345,7 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) do iTrac = 1, nTracers tracersGroup(iTrac, kTop:kBot, iCell) = dataNew(1, iTrac, 1:nLayers) end do - + end do !$omp end do !$omp end parallel @@ -353,13 +353,11 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) deallocate(dataNow,dataNew) deallocate(bcUpper,bcLower) call work%free() - + end if end if end do - !- dealloc any local data, etc - deallocate(heightCellNow) deallocate(heightEdgeNow) deallocate(heightCellNew) @@ -381,9 +379,9 @@ end subroutine ocn_remap_vert_state !----------------------------------------------------------------------- subroutine ocn_vertical_remap_init(err) - + integer, intent(out) :: err !< Output: Error flag - + err = 0 if ( config_time_integrator == 'RK4' ) then From f88cef9ea22668ffc4f7b8694087c140fa6aee5d Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 18 Feb 2022 12:00:58 -0600 Subject: [PATCH 14/22] Move zstar computation to regrid module --- .../mpas_ocn_time_integration_split.F | 28 +--- components/mpas-ocean/src/ocean.cmake | 1 + components/mpas-ocean/src/shared/Makefile | 5 +- .../src/shared/mpas_ocn_vertical_regrid.F | 131 ++++++++++++++++++ .../src/shared/mpas_ocn_vertical_remap.F | 18 ++- 5 files changed, 151 insertions(+), 32 deletions(-) create mode 100644 components/mpas-ocean/src/shared/mpas_ocn_vertical_regrid.F diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index bb9de5fe95fc..cb17bba1ec79 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -266,10 +266,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersNew ! Remap variables - real (kind=RKIND) :: totalThickness - real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessRemap real (kind=RKIND), dimension(:,:), pointer :: & - restingThickness, & layerThicknessLagNew ! These are public variables used from the diagnostics module @@ -2445,35 +2442,14 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call mpas_timer_start("vertical remap") - call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) call mpas_pool_get_array(statePool, 'layerThicknessLag', layerThicknessLagNew, 2) - allocate(layerThicknessRemap(nVertLevels, nCellsAll)) - ! Initialize target layerThickness for remapping to - ! its current, Lagrangian location - layerThicknessRemap = layerThicknessNew ! Store the current, Lagrangian location for computation of the ! Eulerian vertical velocity induced by remapping layerThicknessLagNew = layerThicknessNew - ! calc. z*-style expansions - !$omp parallel - !$omp do schedule(runtime) private(k,totalThickness) - do iCell = 1, nCellsAll - totalThickness = 0.0_RKIND - do k = maxLevelCell(iCell), minLevelCell(iCell), -1 - totalThickness = totalThickness + layerThicknessNew(k,iCell) - end do - do k = maxLevelCell(iCell), minLevelCell(iCell), -1 - layerThicknessRemap(k, iCell) = restingThickness(k, iCell) * & - totalThickness / bottomDepth(iCell) - end do - end do - !$omp end do - !$omp end parallel - call ocn_remap_vert_state(block, layerThicknessRemap, err) - - deallocate(layerThicknessRemap) + ! Perform the remapping + call ocn_remap_vert_state(block, err) call mpas_timer_stop("vertical remap") diff --git a/components/mpas-ocean/src/ocean.cmake b/components/mpas-ocean/src/ocean.cmake index 422116ef79ab..6bae4d67342b 100644 --- a/components/mpas-ocean/src/ocean.cmake +++ b/components/mpas-ocean/src/ocean.cmake @@ -54,6 +54,7 @@ list(APPEND RAW_SOURCES core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F core_ocean/shared/mpas_ocn_vel_pressure_grad.F core_ocean/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F + core_ocean/shared/mpas_ocn_vertical_regrid.F core_ocean/shared/mpas_ocn_vertical_remap.F core_ocean/shared/mpas_ocn_vmix.F core_ocean/shared/mpas_ocn_vmix_coefs_redi.F diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index 0a4bed12fe3c..15d6c1fdf77c 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -24,6 +24,7 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_vel_forcing_explicit_bottom_drag.o \ mpas_ocn_vel_pressure_grad.o \ mpas_ocn_vertical_remap.o \ + mpas_ocn_vertical_regrid.o \ mpas_ocn_vmix.o \ mpas_ocn_vmix_coefs_redi.o \ mpas_ocn_vmix_cvmix.o \ @@ -218,7 +219,9 @@ mpas_ocn_vel_pressure_grad.o: mpas_ocn_constants.o mpas_ocn_vel_tidal_potential.o: mpas_ocn_diagnostics_variables.o -mpas_ocn_vertical_remap.o: mpas_ocn_diagnostics.o mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_mesh.o +mpas_ocn_vertical_regrid.o: mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_mesh.o + +mpas_ocn_vertical_remap.o: mpas_ocn_diagnostics.o mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_mesh.o mpas_ocn_vertical_regrid.o clean: $(RM) *.o *.i *.mod *.f90 diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_regrid.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_regrid.F new file mode 100644 index 000000000000..369a7919f31e --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_regrid.F @@ -0,0 +1,131 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vertical_regrid +! +!> \brief MPAS ocean vertical regridding +!> \author Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis +!> \date July 2021 +!> \details +!> This module contains the vertical regridding routine, used for +!> vertical Lagrangian remapping. +! +!----------------------------------------------------------------------- + +module ocn_vertical_regrid + + use mpas_derived_types + use mpas_pool_routines + use mpas_constants + use mpas_dmpar + use mpas_threading + use mpas_vector_reconstruction + use mpas_spline_interpolation + use mpas_timer + use mpas_log + + use ocn_constants + use ocn_config + use ocn_mesh + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_vert_regrid + public :: ocn_vert_regrid_init + + contains + +!*********************************************************************** +! +! routine ocn_vert_regrid +! +!> \brief Regridding for vertical Lagrangian remapping +!> \author Carolyn Begeman +!> \date 2021 +!> \details +!> This routine determines the layerThickness to remap to subject to +!> constraints +! +!----------------------------------------------------------------------- + + subroutine ocn_vert_regrid(restingThickness, layerThicknessLag, & + layerThicknessTarget, err) + + real (kind=RKIND), dimension(:, :), intent(in) :: & + layerThicknessLag, & ! layerThickness after the lagrangian step + restingThickness + real (kind=RKIND), dimension(:, :), intent(out) :: & + layerThicknessTarget ! adjusted target layerThickness for remapping + integer, intent(out) :: err !< Output: Error flag + + integer :: iCell, k, kmin, kmax + + real (kind=RKIND) :: totalThickness + + err = 0 + + ! Calculate z-star layer locations + layerThicknessTarget = 0.0_RKIND + + !$omp parallel + !$omp do schedule(runtime) private(k,totalThickness) + do iCell = 1, nCellsAll + totalThickness = 0.0_RKIND + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 + totalThickness = totalThickness + layerThicknessLag(k,iCell) + end do + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 + layerThicknessTarget(k, iCell) = restingThickness(k, iCell) * & + totalThickness / bottomDepth(iCell) + end do + end do + !$omp end do + !$omp end parallel + + + end subroutine ocn_vert_regrid + + +!*********************************************************************** +! +! routine ocn_vert_regrid_init +! +!> \brief Initializes ocean vertical regridding +!> \author Carolyn Begeman +!> \date July 2021 +!> \details +!> This routine initializes parameters required for vertical Lagrangian +!> regridding +! +!----------------------------------------------------------------------- + + subroutine ocn_vert_regrid_init(err) + + integer, intent(out) :: err !< Output: Error flag + + err = 0 + + end subroutine ocn_vert_regrid_init + +end module ocn_vertical_regrid +! vim: foldmethod=marker diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F index e21173210be5..15bbb6b3f7d9 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -32,6 +32,7 @@ module ocn_vertical_remap use ocn_config use ocn_mesh use ocn_diagnostics + use ocn_vertical_regrid use ppr_1d @@ -65,19 +66,18 @@ module ocn_vertical_remap ! ! ocn_remap_vert_state ! -!> \brief MPAS ocean vertical Lagrangian remapping -!> \author Mark Petersen, Doug Jacobsen, Todd Ringler -!> \date September 2011 +!> \brief MPAS ocean vertical Lagrangian remapping +!> \author Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis +!> \date July 2021 !> \details !> This routine remaps state variables from current layerThickness to !> a new layerThickness ! !----------------------------------------------------------------------- - subroutine ocn_remap_vert_state(block, layerThicknessNew, err) + subroutine ocn_remap_vert_state(block, err) type (block_type), intent(in) :: block - real (kind=RKIND), dimension(:, :), intent(in) :: layerThicknessNew integer, intent(out) :: err integer :: nCells, nEdges @@ -94,6 +94,7 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) real (kind=RKIND), dimension(:, :, :), allocatable :: dataNow, dataNew real (kind=RKIND), dimension(:, :), allocatable :: & + layerThicknessNew, & ! layerThickness (new target grid) layerThickEdge, & ! layerThickness at edges (lagrangian grid) layerThickEdgeNew, & ! layerThickness at edges (new target grid) heightCellNow, heightEdgeNow, & ! depth at layer interfaces (lagrangian grid) @@ -135,9 +136,14 @@ subroutine ocn_remap_vert_state(block, layerThicknessNew, err) allocate(heightCellNew(nVertLevels + 1, nCells)) allocate(heightEdgeNow(nVertLevels + 1, nEdges)) allocate(heightEdgeNew(nVertLevels + 1, nEdges)) + allocate(layerThicknessNew(nVertLevels, nCells + 1)) allocate(layerThickEdge(nVertLevels, nEdges)) allocate(layerThickEdgeNew(nVertLevels, nEdges)) + ! Compute new layer thicknesses based on vertical coordinate choice + ! Currently, zstar coordinate is hard-coded + call ocn_vert_regrid(restingThickness, layerThickness, & + layerThicknessNew, err) ! There's currently a silly issue in PPR: it assumes the "coordinates" for ! the data being remapped increases with index. Here, we have the opposite @@ -384,6 +390,8 @@ subroutine ocn_vertical_remap_init(err) err = 0 + call ocn_vert_regrid_init(err) + if ( config_time_integrator == 'RK4' ) then CALL mpas_log_write('Vertical remap not supported for RK4', & MPAS_LOG_CRIT) From b7e337cc0b281c49c6d3973ae08f55d732d55100 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 18 Feb 2022 15:10:14 -0600 Subject: [PATCH 15/22] Move vertical advection config to new module --- .../src/mode_forward/mpas_ocn_forward_mode.F | 7 +- .../mpas_ocn_time_integration_split.F | 18 ++-- components/mpas-ocean/src/ocean.cmake | 1 + components/mpas-ocean/src/shared/Makefile | 7 +- .../src/shared/mpas_ocn_diagnostics.F | 7 +- .../src/shared/mpas_ocn_vertical_advection.F | 98 +++++++++++++++++++ 6 files changed, 119 insertions(+), 19 deletions(-) create mode 100644 components/mpas-ocean/src/shared/mpas_ocn_vertical_advection.F diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F index 001a881131aa..d007eaedd8c3 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F @@ -83,7 +83,8 @@ module ocn_forward_mode use ocn_vmix use ocn_vmix_gotm, only: ocn_vmix_gotm_finalize - use ocn_vertical_remap + use ocn_vertical_advection + use ocn_vertical_remap, only: ocn_vertical_remap_init use ocn_forcing use ocn_time_varying_forcing @@ -464,7 +465,9 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{ ierr = ior(ierr, err_tmp) call ocn_tendency_init(err_tmp) ierr = ior(ierr,err_tmp) - call ocn_vertical_remap_init(err_tmp) + call ocn_vertical_advection_init(err_tmp) + if ( configVertAdvMethod == vertAdvRemap ) & + call ocn_vertical_remap_init(err_tmp) ierr = ior(ierr,err_tmp) if (ierr /= 0) then diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index cb17bba1ec79..d89847e37bf6 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -40,6 +40,7 @@ module ocn_time_integration_split use ocn_equation_of_state use ocn_vmix + use ocn_vertical_advection use ocn_vertical_remap use ocn_time_average_coupled @@ -84,11 +85,10 @@ module ocn_time_integration_split ! outer timestep iteration integer :: & - configVertAdvMethod, &! type of vertical advection method - neededHalos, &! number of halo levels needed - haloDecrement, &! number of halos to invalidate each cycle - numTSIterations, &! number of outer timestep iterations - nBtrSubcycles ! number of barotropic subcycles + neededHalos, &! number of halo levels needed + haloDecrement, &! number of halos to invalidate each cycle + numTSIterations, &! number of outer timestep iterations + nBtrSubcycles ! number of barotropic subcycles integer :: ssh_sal_on @@ -2438,7 +2438,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call mpas_timer_stop("se implicit vert mix") - if ( configVertAdvMethod == 2 ) then + if ( configVertAdvMethod == vertAdvRemap ) then call mpas_timer_start("vertical remap") @@ -2745,12 +2745,6 @@ subroutine ocn_time_integration_split_init(domain)!{{{ useVelocityCorrection = 0.0_RKIND endif - if ( config_vert_advection_method == 'flux-form' ) then - configVertAdvMethod = 1 - elseif ( config_vert_advection_method == 'remap' ) then - configVertAdvMethod = 2 - endif - !*** Compute some halo needs to reduce the number !*** of halo updates during subcycling neededHalos = 1 + config_n_btr_cor_iter diff --git a/components/mpas-ocean/src/ocean.cmake b/components/mpas-ocean/src/ocean.cmake index 6bae4d67342b..fb1dba45d4ba 100644 --- a/components/mpas-ocean/src/ocean.cmake +++ b/components/mpas-ocean/src/ocean.cmake @@ -54,6 +54,7 @@ list(APPEND RAW_SOURCES core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F core_ocean/shared/mpas_ocn_vel_pressure_grad.F core_ocean/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F + core_ocean/shared/mpas_ocn_vertical_advection.F core_ocean/shared/mpas_ocn_vertical_regrid.F core_ocean/shared/mpas_ocn_vertical_remap.F core_ocean/shared/mpas_ocn_vmix.F diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index 15d6c1fdf77c..7eb16e0b42de 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -71,7 +71,8 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_vel_tidal_potential.o \ mpas_ocn_vel_forcing_topographic_wave_drag.o \ mpas_ocn_transport_tests.o \ - mpas_ocn_vel_self_attraction_loading.o + mpas_ocn_vel_self_attraction_loading.o \ + mpas_ocn_vertical_advection.o all: $(OBJS) @@ -79,7 +80,7 @@ mpas_ocn_init_routines.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_tendency.o: mpas_ocn_high_freq_thickness_hmix_del2.o mpas_ocn_tracer_surface_restoring.o mpas_ocn_thick_surface_flux.o mpas_ocn_tracer_short_wave_absorption.o mpas_ocn_tracer_advection.o mpas_ocn_tracer_hmix.o mpas_ocn_tracer_nonlocalflux.o mpas_ocn_surface_bulk_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_tracer_surface_flux_to_tend.o mpas_ocn_tracer_interior_restoring.o mpas_ocn_tracer_exponential_decay.o mpas_ocn_tracer_ideal_age.o mpas_ocn_tracer_TTD.o mpas_ocn_vmix.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_frazil_forcing.o mpas_ocn_tidal_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_tracer_CFC.o mpas_ocn_diagnostics.o mpas_ocn_wetting_drying.o mpas_ocn_vel_self_attraction_loading.o mpas_ocn_vel_tidal_potential.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_thick_hadv.o mpas_ocn_thick_vadv.o mpas_ocn_vel_hadv_coriolis.o mpas_ocn_vel_pressure_grad.o mpas_ocn_vel_vadv.o mpas_ocn_vel_hmix.o mpas_ocn_vel_forcing.o -mpas_ocn_diagnostics.o: mpas_ocn_thick_ale.o mpas_ocn_equation_of_state.o mpas_ocn_gm.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_surface_land_ice_fluxes.o +mpas_ocn_diagnostics.o: mpas_ocn_thick_ale.o mpas_ocn_equation_of_state.o mpas_ocn_gm.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_vertical_advection.o mpas_ocn_diagnostics_variables.o: mpas_ocn_config.o @@ -219,6 +220,8 @@ mpas_ocn_vel_pressure_grad.o: mpas_ocn_constants.o mpas_ocn_vel_tidal_potential.o: mpas_ocn_diagnostics_variables.o +mpas_ocn_vertical_advection.o: mpas_ocn_config.o + mpas_ocn_vertical_regrid.o: mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_mesh.o mpas_ocn_vertical_remap.o: mpas_ocn_diagnostics.o mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_mesh.o mpas_ocn_vertical_regrid.o diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index c2051ef03506..6b72c56a3a63 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -37,6 +37,7 @@ module ocn_diagnostics use ocn_diagnostics_variables use ocn_mesh use ocn_surface_land_ice_fluxes + use ocn_vertical_advection implicit none private @@ -239,7 +240,7 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, scratchPoo vertGMBolusVelocityTop, relativeVorticityCell, divergence, & kineticEnergyCell, tangentialVelocity) - if ( config_vert_advection_method == 'remap' ) then + if ( configVertAdvMethod == vertAdvRemap ) then ! Overwrite vertVelocityTop with exact Eulerian solution from remapping call mpas_pool_get_array(statePool, 'layerThicknessLag', & @@ -2334,8 +2335,8 @@ subroutine ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, oldLayerT err = 0 - if (config_vert_coord_movement.eq.'impermeable_interfaces' .or. & - config_vert_advection_method.eq.'remap') then + if (config_vert_coord_movement == 'impermeable_interfaces' .or. & + configVertAdvMethod == vertAdvRemap) then vertAleTransportTop=0.0_RKIND !$acc update device(vertAleTransportTop) return diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_advection.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_advection.F new file mode 100644 index 000000000000..a390d6529ef5 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_advection.F @@ -0,0 +1,98 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vertical_advection +! +!> \brief MPAS ocean vertical advection module +!> \author Carolyn Begeman +!> \date February 2022 +!> \details +!> This module contains the initialization for vertical advection +!> schemes. +! +!----------------------------------------------------------------------- + +module ocn_vertical_advection + + use ocn_config + use mpas_log + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + integer, public :: & + configVertAdvMethod ! choice of vertical advection method + integer, parameter, public :: &! supported vertical advection methods + vertAdvFluxForm = 1, &! flux form + vertAdvRemap = 2 ! remapping + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + public :: ocn_vertical_advection_init + +contains + +!*********************************************************************** + +!*********************************************************************** +! +! routine ocn_vertical_advection_init +! +!> \brief Initializes vertical advection scheme properties +!> \author Carolyn Begeman +!> \date February 2022 +!> \details +!> This routine initializes quantities related to +!> the vertical advection scheme. +! +!---------------------------------------------------------------------- + + subroutine ocn_vertical_advection_init(err)!{{{ + + !-------------------------------------------------------------------- + + integer, intent(out) :: err + + + err = 0 + + select case (trim(config_vert_advection_method)) + + case ('flux-form') + + configVertAdvMethod = vertAdvFluxForm + + case ('remap') + + configVertAdvMethod = vertAdvRemap + + case default + + call mpas_log_write( & + 'Invalid choice for config_vert_advection_method. Choices are: flux-form, remap') + err = 1 + + end select + + end subroutine ocn_vertical_advection_init!}}} + +!*********************************************************************** + +end module ocn_vertical_advection + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker From 22493f12c222036107b0cc7a89374a90d412595f Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 3 May 2022 14:58:39 -0500 Subject: [PATCH 16/22] Remove layerThickEdge diagnostics call in remap --- components/mpas-ocean/src/shared/Makefile | 2 +- .../src/shared/mpas_ocn_vertical_remap.F | 41 +++++++++++++++---- 2 files changed, 35 insertions(+), 8 deletions(-) diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index 7eb16e0b42de..81d8235bb707 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -224,7 +224,7 @@ mpas_ocn_vertical_advection.o: mpas_ocn_config.o mpas_ocn_vertical_regrid.o: mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_mesh.o -mpas_ocn_vertical_remap.o: mpas_ocn_diagnostics.o mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_mesh.o mpas_ocn_vertical_regrid.o +mpas_ocn_vertical_remap.o: mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_diagnostics_variables.o mpas_ocn_mesh.o mpas_ocn_vertical_regrid.o clean: $(RM) *.o *.i *.mod *.f90 diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F index 15bbb6b3f7d9..b684725f5c6f 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -30,8 +30,8 @@ module ocn_vertical_remap use ocn_constants use ocn_config + use ocn_diagnostics_variables use ocn_mesh - use ocn_diagnostics use ocn_vertical_regrid use ppr_1d @@ -82,7 +82,8 @@ subroutine ocn_remap_vert_state(block, err) integer :: nCells, nEdges integer :: nLayers, nLevels, nVars, nDoFs, nTracers - integer :: cell1, cell2, iCell, jCell, iEdge, iTrac, k, kTop, kBot + integer :: cell1, cell2, iCell, jCell, iEdge, iTrac + integer :: k, kmin, kmax, kTop, kBot type (mpas_pool_type), pointer :: verticalMeshPool type (mpas_pool_type), pointer :: statePool @@ -95,7 +96,6 @@ subroutine ocn_remap_vert_state(block, err) real (kind=RKIND), dimension(:, :, :), allocatable :: dataNow, dataNew real (kind=RKIND), dimension(:, :), allocatable :: & layerThicknessNew, & ! layerThickness (new target grid) - layerThickEdge, & ! layerThickness at edges (lagrangian grid) layerThickEdgeNew, & ! layerThickness at edges (new target grid) heightCellNow, heightEdgeNow, & ! depth at layer interfaces (lagrangian grid) heightCellNew, heightEdgeNew ! depth at layer interfaces (new target grid) @@ -137,7 +137,6 @@ subroutine ocn_remap_vert_state(block, err) allocate(heightEdgeNow(nVertLevels + 1, nEdges)) allocate(heightEdgeNew(nVertLevels + 1, nEdges)) allocate(layerThicknessNew(nVertLevels, nCells + 1)) - allocate(layerThickEdge(nVertLevels, nEdges)) allocate(layerThickEdgeNew(nVertLevels, nEdges)) ! Compute new layer thicknesses based on vertical coordinate choice @@ -173,8 +172,36 @@ subroutine ocn_remap_vert_state(block, err) !$omp end do !$omp end parallel - call ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, layerThickness, layerThickEdge) - call ocn_diagnostic_solve_layerThicknessEdge(normalVelocity, layerThicknessNew, layerThickEdgeNew) +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(layerThickness, layerThickEdgeMean, & + !$acc minLevelEdgeBot, maxLevelEdgeTop, cellsOnEdge) & + !$acc private(k, kmin, kmax, cell1, cell2) +#else + !$omp parallel + !$omp do schedule(runtime) private(k, kmin, kmax, cell1, cell2) +#endif + do iEdge = 1, nEdgesAll + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1,nVertLevels + ! initialize layerThicknessEdgeMean to avoid divide by + ! zero and NaN problems. + layerThickEdgeNew(k,iEdge) = -1.0e34_RKIND + end do + do k = kmin,kmax + ! central differenced + layerThickEdgeNew(k,iEdge) = 0.5_RKIND * & + (layerThicknessNew(k,cell1) + & + layerThicknessNew(k,cell2)) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif !$omp parallel !$omp do schedule(runtime) private(k,kTop,kBot,cell1,cell2) @@ -201,7 +228,7 @@ subroutine ocn_remap_vert_state(block, err) endif do k = kBot, kTop, -1 - heightEdgeNow(k, iEdge) = heightEdgeNow(k+1, iEdge) - layerThickEdge(k, iEdge) + heightEdgeNow(k, iEdge) = heightEdgeNow(k+1, iEdge) - layerThickEdgeMean(k, iEdge) heightEdgeNew(k, iEdge) = heightEdgeNew(k+1, iEdge) - layerThickEdgeNew(k, iEdge) end do From b3d01b0bbc65de26953f3f64abf6965aeda48b81 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 8 Jun 2022 10:50:37 -0600 Subject: [PATCH 17/22] Fixup omp private in vertical remap --- .../mpas-ocean/src/shared/mpas_ocn_vertical_remap.F | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F index b684725f5c6f..f16f009fe6e6 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -268,7 +268,7 @@ subroutine ocn_remap_vert_state(block, err) ! Remap all edge-centred variables from heightEdgeNow to heightEdgeNew if (config_use_freq_filtered_thickness) then !$omp parallel - !$omp do schedule(runtime) private(iEdge, kTop, kBot, nLayers, nLevels, dataNow, dataNew) + !$omp do schedule(runtime) private(kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) do iEdge = 1, nEdges kTop = minLevelEdgeBot(iEdge) @@ -299,7 +299,7 @@ subroutine ocn_remap_vert_state(block, err) !$omp end parallel else !$omp parallel - !$omp do schedule(runtime) private(iEdge, kTop, kBot, nLayers, nLevels, dataNow, dataNew) + !$omp do schedule(runtime) private(kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) do iEdge = 1, nEdges kTop = minLevelEdgeBot(iEdge) @@ -313,7 +313,7 @@ subroutine ocn_remap_vert_state(block, err) dataNow(1, 1, 1:nLayers) = normalVelocity(kTop:kBot, iEdge) call mpas_timer_start("ppr") - call rmap1d(nLevels, nLevels, nVars ,nDoFs, & + call rmap1d(nLevels, nLevels, nVars, nDoFs, & heightEdgeNow(kTop:kBot+1, iEdge), & heightEdgeNew(kTop:kBot+1, iEdge), & dataNow, dataNew, & @@ -351,7 +351,7 @@ subroutine ocn_remap_vert_state(block, err) bcLower%bcopt = bc_lower !$omp parallel - !$omp do schedule(runtime) private(iCell, iTrac, kTop, kBot, nLayers, nLevels, dataNow, dataNew) + !$omp do schedule(runtime) private(iTrac, kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) do iCell = 1, nCells kTop = minLevelCell(iCell) From 76b5f4f41e43396cf776c3381cb75e4182db9510 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 8 Jun 2022 10:38:32 -0600 Subject: [PATCH 18/22] Remove ppr timer --- components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F | 6 ------ 1 file changed, 6 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F index f16f009fe6e6..8a12d32a28e1 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -283,13 +283,11 @@ subroutine ocn_remap_vert_state(block, err) dataNow(1, 2, 1:nLayers) = highFreqThickness(kTop:kBot, iEdge) dataNow(1, 3, 1:nLayers) = lowFreqDivergence(kTop:kBot, iEdge) - call mpas_timer_start("ppr") call rmap1d(nLevels, nLevels, nVars ,nDoFs, & heightEdgeNow(kTop:kBot+1, iEdge), & heightEdgeNew(kTop:kBot+1, iEdge), & dataNow, dataNew, & bcUpper, bcLower, work, opts) - call mpas_timer_stop("ppr") normalVelocity (kTop:kBot, iEdge) = dataNew(1, 1, 1:nLayers) highFreqThickness(kTop:kBot, iEdge) = dataNew(1, 2, 1:nLayers) @@ -312,13 +310,11 @@ subroutine ocn_remap_vert_state(block, err) dataNow(1, 1, 1:nLayers) = normalVelocity(kTop:kBot, iEdge) - call mpas_timer_start("ppr") call rmap1d(nLevels, nLevels, nVars, nDoFs, & heightEdgeNow(kTop:kBot+1, iEdge), & heightEdgeNew(kTop:kBot+1, iEdge), & dataNow, dataNew, & bcUpper, bcLower, work, opts) - call mpas_timer_stop("ppr") normalVelocity (kTop:kBot, iEdge) = dataNew(1, 1, 1:nLayers) @@ -367,13 +363,11 @@ subroutine ocn_remap_vert_state(block, err) dataNow(1, iTrac, 1:nLayers) = tracersGroup(iTrac, kTop:kBot, iCell) end do - call mpas_timer_start("ppr") call rmap1d(nLevels, nLevels, nTracers,nDoFs, & heightCellNow(kTop:kBot+1, iCell), & heightCellNew(kTop:kBot+1, iCell), & dataNow, dataNew, & bcUpper, bcLower, work, opts) - call mpas_timer_stop("ppr") do iTrac = 1, nTracers tracersGroup(iTrac, kTop:kBot, iCell) = dataNew(1, iTrac, 1:nLayers) From 208f75ad2df5df12ab84effec28b7fbf3b2ef802 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 9 Jun 2022 10:10:39 -0600 Subject: [PATCH 19/22] Aesthetic changes in remap module --- .../src/shared/mpas_ocn_vertical_remap.F | 30 +++++++++---------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F index 8a12d32a28e1..7bde91d922b0 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -90,7 +90,7 @@ subroutine ocn_remap_vert_state(block, err) type (mpas_pool_type), pointer :: tracersPool type (mpas_pool_iterator_type) :: groupItr - type(rmap_work) :: work + type(rmap_work) :: work ! Internal workspace for PPR type(rcon_ends), dimension(:), allocatable :: bcUpper, bcLower real (kind=RKIND), dimension(:, :, :), allocatable :: dataNow, dataNew @@ -144,11 +144,8 @@ subroutine ocn_remap_vert_state(block, err) call ocn_vert_regrid(restingThickness, layerThickness, & layerThicknessNew, err) - ! There's currently a silly issue in PPR: it assumes the "coordinates" for - ! the data being remapped increases with index. Here, we have the opposite - ! (z becomes more -ve with k)... - ! I'll make some updates to PPR to handle this, but for now, just flip the - ! signs to make PPR happy. This shouldn't have any actual impact. + ! Compute the layer interface locations between bottomDepth and -ssh + ! (rather than between -bottomDepth and ssh due to PPR limitations) !$omp parallel !$omp do schedule(runtime) private(k) @@ -186,12 +183,12 @@ subroutine ocn_remap_vert_state(block, err) kmax = maxLevelEdgeTop(iEdge) cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) - do k = 1,nVertLevels + do k = 1, nVertLevels ! initialize layerThicknessEdgeMean to avoid divide by ! zero and NaN problems. layerThickEdgeNew(k,iEdge) = -1.0e34_RKIND end do - do k = kmin,kmax + do k = kmin, kmax ! central differenced layerThickEdgeNew(k,iEdge) = 0.5_RKIND * & (layerThicknessNew(k,cell1) + & @@ -204,7 +201,7 @@ subroutine ocn_remap_vert_state(block, err) #endif !$omp parallel - !$omp do schedule(runtime) private(k,kTop,kBot,cell1,cell2) + !$omp do schedule(runtime) private(k, kTop, kBot, cell1, cell2) do iEdge = 1, nEdges cell1 = cellsOnEdge(1,iEdge) @@ -263,12 +260,13 @@ subroutine ocn_remap_vert_state(block, err) dataNow = 0.0_RKIND dataNew = 0.0_RKIND - call work%init(nVertLevels + 1, nVars, opts) ! this is internal workspace for PPR + call work%init(nVertLevels + 1, nVars, opts) ! Remap all edge-centred variables from heightEdgeNow to heightEdgeNew if (config_use_freq_filtered_thickness) then !$omp parallel - !$omp do schedule(runtime) private(kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) + !$omp do schedule(runtime) & + !$omp private(kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) do iEdge = 1, nEdges kTop = minLevelEdgeBot(iEdge) @@ -297,7 +295,8 @@ subroutine ocn_remap_vert_state(block, err) !$omp end parallel else !$omp parallel - !$omp do schedule(runtime) private(kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) + !$omp do schedule(runtime) & + !$omp private(kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) do iEdge = 1, nEdges kTop = minLevelEdgeBot(iEdge) @@ -339,7 +338,7 @@ subroutine ocn_remap_vert_state(block, err) allocate(dataNow(nDoFs, nTracers, nVertLevels)) allocate(dataNew(nDoFs, nTracers, nVertLevels)) - call work%init(nVertLevels + 1, nTracers, opts) ! this is internal workspace for PPR + call work%init(nVertLevels + 1, nTracers, opts) allocate(bcUpper(nTracers)) allocate(bcLower(nTracers)) @@ -347,7 +346,8 @@ subroutine ocn_remap_vert_state(block, err) bcLower%bcopt = bc_lower !$omp parallel - !$omp do schedule(runtime) private(iTrac, kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) + !$omp do schedule(runtime) & + !$omp private(iTrac, kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) do iCell = 1, nCells kTop = minLevelCell(iCell) @@ -438,7 +438,7 @@ subroutine ocn_vertical_remap_init(err) opts%cell_lims = mono_limit ! monotone slope limits elseif ( config_remap_limiter == 'weno' ) then opts%cell_lims = weno_limit ! WENO slope limits - opts%wall_lims = mono_limit ! WENO slope limits + opts%wall_lims = mono_limit ! monotone slope limits elseif ( config_remap_limiter == 'none' ) then opts%cell_lims = null_limit ! no slope limits else From 847dc1458f89a2f7ac02b5966893d8affc785b97 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 20 May 2022 18:01:58 -0500 Subject: [PATCH 20/22] Add vertical remap design doc --- .../mpas-ocean/docs/design_docs/index.rst | 1 + .../vert-lagrang-remap/vert-lagrang-remap.rst | 582 ++++++++++++++++++ 2 files changed, 583 insertions(+) create mode 100644 components/mpas-ocean/docs/design_docs/vert-lagrang-remap/vert-lagrang-remap.rst diff --git a/components/mpas-ocean/docs/design_docs/index.rst b/components/mpas-ocean/docs/design_docs/index.rst index 703a9d4efabb..190686a4a264 100644 --- a/components/mpas-ocean/docs/design_docs/index.rst +++ b/components/mpas-ocean/docs/design_docs/index.rst @@ -8,3 +8,4 @@ Design document describing new capabilities added to MPAS-Ocean. time-varying-wind prognostic-sediment-transport + vert-lagrang-remap/vert-lagrang-remap diff --git a/components/mpas-ocean/docs/design_docs/vert-lagrang-remap/vert-lagrang-remap.rst b/components/mpas-ocean/docs/design_docs/vert-lagrang-remap/vert-lagrang-remap.rst new file mode 100644 index 000000000000..a2328e75be2c --- /dev/null +++ b/components/mpas-ocean/docs/design_docs/vert-lagrang-remap/vert-lagrang-remap.rst @@ -0,0 +1,582 @@ + +Vertical Lagrangian-remap method +================================ + +date: 2020/12/21 +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + + + +Summary +------- + +We propose to implement the vertical Lagrangian-remap method to support the +choice of hybrid vertical coordinates, which may be required for evolving ice- +shelf cavity geometries. It can be seen as an extension of the Arbitrary +Lagrangian-Eulerian method which is already implemented in MPAS-Ocean, a key +difference being that the solution of diasurface transport is coordinate-free. +This has the advantage of removing the vertical CFL constraint, which is of +particular concern for an evolving ice shelf geometry in terrain-following +coordinate schemes. This development may also reduce numerical diffusion in the +vertical, but we do not require this result for this implementation to be +successful. + +Generalization of the vertical coordinate is not in the scope of this development. +Additionally, the infrastructure for specifying a target grid that differs from +the existing z-star and z-tilde options is reserved for future development. +The minimum criteria for success are not degrading the accuracy of the solution +or computational performance. + + +Requirements +------------ + +Requirement: Ability to perform vertical Lagrangian-remapping +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +The core of the proposed development is the ability to perform vertical +Lagrangian-remapping (VLR). When VLR is performed, the momentum, volume and +tracer conservation equations are first solved in a Lagrangian sense; that is, +layer thicknesses evolve such that there is no vertical transport across layer +interfaces. A target grid is then specified. For this stage of development, +we only require that the algorithm be able to use the existing grid +specifications. However, the implementation should be general enough that +a different grid could be used for remapping. This regridding step is followed +by a conservative remapping of velocity and scalars to the target grid. + +Requirement: Grid is sufficiently conditioned for accuracy and stability +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +A severely deformed grid can lead to issues with the accuracy and stability of +the numerical methods used for the dynamics and thermodynamics. If layers are +allowed to become very thin they can also degrade the accuracy and stability of +the solution. Thus, a requirement of the VLR implementation is that there are +safeguards to prevent layers from becoming too thin or deformed. + +Requirement: Ability to specify remapping method and its options +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +There are several choices to make for the remapping operation to balance +accuracy, stability and computational cost. To the extent that it is feasible, +the user should have the ability to control the order of accuracy of the +remapping method, and other remapping options (e.g., via namelist options). + +Requirement: Support for existing time stepping schemes +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Vertical Lagrangian remapping should be supported (eventually) in both the RK4 +and the split-explicit time stepping schemes currently implemented in MPAS- +Ocean. The scheme will also have a straightforward interface so it can be easily +added to new time stepping schemes. + +Requirement: not climate changing +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +For the same target grid, time-stepping scheme and problem, the solution for +layer thicknesses and key prognostic variables using VLR is within a non- +climate-changing threshold value of the solution when VLR is not performed. +Tests should span from idealized tests of simplified dynamics to full climate +simulations. Acceptable thresholds specified for each test could be similar to +those from changing time stepping schemes. + +Requirement: tracer and volume conservation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2020/12/04 + +Contributors: Darren Engwirda, Xylar Asay-Davis, Carolyn Begeman, Mark Petersen + +The column-wise tracer content should be conserved over the remapping step to a +specified threshold, and using a specified method to integrate the tracer. The +column-integrated volume must remain unchanged over the remapping step, as +expected. + +Requirement: performance +^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2020/12/04 + +Contributors: Darren Engwirda, Xylar Asay-Davis, Carolyn Begeman, Mark Petersen + +Performance is not degraded below an acceptable level. Compute time for +stand-alone ocean, without i/o, on full primitive equations, is expected to +remain very similar - within 5% of compute time when using vertical +advection method but likely less because remapping is a column-wise, in-cache +operation. + +Requirement: no interference +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Remapping does not introduce errors in the computation of prognostic and +diagnostic variables, such as through the use of an Eulerian vertical velocity +that is not consistent with the Lagrangian solution. + +Requirement: modularity +^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2020/12/21 + +Contributors: Darren Engwirda, Xylar Asay-Davis, Carolyn Begeman + +To the degree possible, the code for determining the target grid and performing +remapping should be kept in its own Fortran module(s) for better readability. +These modules should be called by both timestepping routines to maximize code +reuse. + + +Algorithm Design (optional) +--------------------------- + +Algorithm Design: Ability to perform vertical Lagrangian-remapping +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2020/12/15 + +Contributors: Darren Engwirda, Xylar Asay-Davis, Carolyn Begeman + +The conservation of momentum, volume, and tracer equations in MPAS-Ocean ( +`Ringler et al. 2013 `_; +`Petersen et al. 2014 `_) +are: + +.. math:: + + \frac{\partial u_k}{\partial t} + q_k h_k u_k^{normal} + \overline{w^t \delta z^t u} = -\frac{1}{\rho_0} \nabla p_k - \frac{\rho_k g}{\rho_0} \nabla z_k - \nabla K_k + [D_h^u]_k + [D_{\nu}^u]_k + F_k^u + + \frac{\partial h_k}{\partial t} + \nabla \cdot (h_k \mathbf{u}_k) + w_k^t - w_{k+1}^t = 0 + + \frac{\partial (h_k \phi_k)}{\partial t} + \nabla \cdot (h_k \mathbf{u}_k \phi_k) + \overline{\phi}_k^t w_k^t - \overline{\phi}_{k+1}^t w_{k+1}^t = [D_h^{\phi}]_k + [D_v^{\phi}]_k + F_k^{\phi} + +For the Lagrangian step, the vertical velocity through the top of the cell, +:math:`w_k^t`, is set to zero in all of the above equations. Thus, these +equations simplify to: + +.. math:: + + \frac{\partial u_k}{\partial t} + q_k h_k u_k^{\perp} = -\frac{1}{\rho_0} \nabla p_k - \frac{\rho_k g}{\rho_0} \nabla z_k - \nabla K_k + [D_h^u]_k + [D_v^u]_k + F_k^u + + \frac{\partial h_k}{\partial t} + \nabla \cdot (h_k \mathbf{u}_k) = 0 + + \frac{\partial (h_k \phi_k)}{\partial t} + \nabla \cdot (h_k \mathbf{u}_k \phi_k) = [D_h^{\phi}]_k + [D_v^{\phi}]_k + F_k^{\phi} + +The time-stepping algorithm (RK4 or split-explicit) advances the prognostic +variables and layer thickness from their values at time n +:math:`u_k^{n},\phi_k^{n},h_k^{n}`, to their values after the Lagrangian step, +designated by the superscript *lg*, :math:`u_k^{lg},h_k^{lg},\phi_k^{lg}`. + +Note that the vertical mixing terms :math:`D_v^h, D_v^{\phi}` +are retained here. We opt to compute these terms prior to remapping as this +allow for future development in which the dynamics are subcycled relative to +the thermodynamics and remapping is scheduled on the thermodynamic timestep. +This computation of vertical mixing terms prior to remapping is similar to +both MOM6 and HYCOM. We anticipate that there could be a trade-off between (a) +loss of accuracy of vertical mixing terms when their computation precedes +remapping due to grid deformation and (b) loss of accuracy when their +computation follows remapping due to remapping errors in vertical gradients of +prognostic variables. We do not intend to test this at this time. + +The target grid needs to be determined after the solution for prognostic +variables so that the vertical Lagrangian-remapping method is general enough to +be used with coordinate systems that depend on the ocean state (this includes +the z-star coordinate system in which SSH perturbations are vertically +distributed between layers). We do not present an algorithmic design for +regridding to coordinate systems not already supported in MPAS-Ocean, as this +will be the subject of future development. For now, the target grid is based on a +constant set of z-star levels that are specified at initialization. + +For the grid selection step, the target grid, :math:`h_k^{target}`, is +determined, conserving volume: + +.. math:: + + \sum_{k=1}^{kmax}h_k^{target} = \sum_{k=1}^{kmax}h_k^{lg} + + +For scalar remapping, layer thicknesses at the next timestep, +:math:`h_k^{n+1}` are set according to the target grid and scalars are remapped +to the target grid using the remapping operations represented by the function +:math:`G`: + +.. math:: + + h_k^{n+1} = h_k^{target} + + hEdge_k^{n+1} = 0.5 (h_{k,cell1}^{n+1} + h_{k,cell2}^{n+1}) + + u_k^{n+1} = G(u_k^{lg},hEdge_k^{lg},hEdge_k^{n+1}) + + \phi_k^{n+1} = G(\phi_k^{lg},\phi_k^{lg},h_k^{n+1}) + +For velocity remapping, we solve for layer thicknesses at edges after the +lagrangian step and the regridded thickness. In this development, we only plan +to support centered edge layer thicknesses consistent with the centered +advection scheme. There does not appear to be a precedent among ocean models +(HYCOM, MOM6) at this time for using upwinded layer thickness in the remapping +operation. We touch on a few of the implementation challenges with using +upwinded layer thicknesses for remapping in the corresponding implementation +section. If VLR is run with an upwinded thickness flux, the horizontal momentum +flux will not be conserved as :math:`hEdge^{n+1}` will be reassigned to the +upwinded layer thickness (errors will likely increase as horizontal gradients +in layer thickness increase). Otherwise, the remapping operation :math:`G` +conserves volume flux and scalar concentration. + +.. math:: + + hEdge_k^{lg} = 0.5 (h_{k,cell1}^{lg} + h_{k,cell2}^{lg}) + + hEdge_k^{n+1} = 0.5 (h_{k,cell1}^{n+1} + h_{k,cell2}^{n+1}) + + u_k^{n+1} = G(u_k^{lg},hEdge_k^{lg},hEdge_k^{n+1}) + + +.. math:: + + \sum_{k=1}^{kmax} u_k^{n+1} h_k^{n+1} = \sum_{k=1}^{kmax} u_k^{lg} h_k^{lg} + + \sum_{k=1}^{kmax} \phi_k^{n+1} h_k^{n+1} = \sum_{k=1}^{kmax} \phi_k^{lg} h_k^{lg} + +The vertical velocity across layer interfaces may be computed anytime after +regridding. It can be computed as + +.. math:: + + w = - \nabla \cdot (h_k \mathbf{u}_k) - (h_k^{t+1} - h_k^t)/dt + +or + +.. math:: + + w = (h_k^{t+1} - h_k^{lg})/dt + +The choice between the two is discussed in the Implementation section. + + +Implementation +-------------- + +Implementation: Ability to perform vertical Lagrangian-remapping +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Namelist options: + +- To turn VLR on/off: + :code:`advection, config_vert_advection_method = 'remap' or 'flux-form'` + + +Lagrangian step: + +The solution for prognostic variables in RK4 and split-explicit remains +largely the same. The main difference is that the vertical velocity through +the top of layers is set to zero in the routine +:code:`ocn_vert_transport_velocity_top`. This is the same as what is done when +:code:`config_vert_coord_movement` is :code:`impermeable_interfaces`. + + +Grid selection step: + +#. :math:`z_k^{target}`, the depth of the top of the layer, is determined based + on an analytical expression for the grid. + The simplest case is constant z-levels, :math:`z_k^{target} = z_k^{init}`. + Since :math:`z_k^{target}` can be a function of the ocean state (e.g., + :math:`\rho` for isopycnal coordinates, grid selection doesn't occur until + after the solution for prognostic variables. +#. Superimpose SSH perturbations according to one of the existing depth- + dependent functions, :math:`z_k^{target} = z_k^{target} + c(z) \: \eta`. As + in :code:`ocn_ALE_thickness`, layer thicknesses are adjusted from the + seafloor upwards. Currently :code:`layerThicknessTarget` is computed in a + new routine :code:`ocn_vert_regrid` following other instances in which a z- + star grid is recomputed. +#. Apply conditioning steps outlined in the following section. (Not + implemented until coordinates other than z-star are implemented.) +#. In preparation for remapping, compute :math:`z_k^{target}` from + :code:`layerThicknessTarget`. + +All of the grid selection steps will be performed from a separate module called +:code:`ocn_vert_regrid`. This topic is further addressed in section +Implementation: modularity. + + +Remapping step: + +There is a new remapping routine :code:`ocn_remap_vert_state` which updates +all state variables with depth coordinates at :code:`timeLevel=2`. On input, +:code:`layerThickness(tlev=2)` reflects the Lagrangian layer thickness +determined by :code:`ocn_tend_thick`. On output, it is equal to +:code:`layerThicknessNew` as determined by the regridding routine +:code:`ocn_vert_regrid`. + +Members of :code:`statePool` that will be remapped: + +- All members of :code:`tracerPool` unless :code:`activeTracersOnly`, in which + case only the :code:`activeTracers` will be remapped +- :code:`normalVelocity` +- :code:`highFreqThickness` +- :code:`lowFreqDivergence` +- :code:`normalBarotropicVelocity`, only for split-explicit time-stepping + +In preparation for remapping, we compute the scratch variables +:code:`layerThickEdgeNew` as the mean of :code:`layerThicknessNew` at +neighboring cells. We do these locally rather than through a call to +:code:`ocn_diagnostic_solve_layerThicknessEdge`. For the Lagrangian +:code:`layerThickEdgeMean` we use the existing solution from +:code:`ocn_diagnostic_solve_layerThicknessEdge`. This approach requires that +:code:`config_thickness_flux_type` is :code:`'centered'`. At initialization, +we throw an error but do not terminate the run if +:code:`config_thickness_flux_type` is not :code:`'centered'` and VLR is active. + + +A note about difficulties of implementing upwinded :code:`layerThicknessEdge` +for remapping: + +Currently, the PPR library assumes that the total height is the same before and +after remapping (i.e., :code:`sum_k(layerThicknessOld)` equals +:code:`sum_k(layerThicknessTarget)`. Over the course of remapping, the upwind +cell could change for one or more layers and thus the total column height could +change. PPR would have to be carefully adapted to deal with this condition in +order to preserve the total volume flux as well as the vertical distribution of +momentum during remapping. + +An alternative to modifying the remapping library is to use centered edge layer +thicknesses for remapping and correct :code:`uNormal` prior to remapping such +that :code:`uNormalCorr(k) * layerThicknessEdgeCntr(k) = uNormal(k) * layerThicknessEdgeUpwind(k)`. +When edge layer thicknesses are upwinded based on the remapped :code:`uNormal`, +:code:`uNormal` must be corrected again to preserve layerwise fluxes. There +will still be some error in the vertical distribution of volume flux with this +approach. Given the complexity of either of these implementation options, we +leave this issue for future development. + + +After determining the layer thicknesses to remap to, the vertical layer +interface locations are determined, :code:`heightCellNow`, +:code:`heightCellNew`, :code:`heightEdgeNow`, :code:`heightEdgeNew`. +These variables and the Lagrangian state variable are the inputs to the PPR +routine :code:`rmap1d`, and the output is the remapped state variable. + + +Some implementation considerations for PPR: + +- Error-checking in PPR: make consistent with MPAS errors, consider additional + error checks. (Not implemented) + +After remapping, :code:`ocn_diagnostic_solve` is called. This is needed to +compute the density and pressure fields based on the remapped ocean state and +the diagnostic field :code:`vertVelocityTop` which is the vertical velocity +through the top of the layer. This is only used as a diagnostic variable for +computing the MOC streamfunction and tracer budgets. None of the mixing +parameterizations require a vertical velocity (Eulerian or diasurface velocity). + +Note: if `vertVelocityTop` is computed between regridding and remapping then it +can be computed as + +.. code:: + + vertVelocityTop(k) = vertVelocityTop(k+1) - div_hu(k) - + (layerThickness(k,tlev=2) - layerThickness(k,tlev=1))/dt + +In our implementation, `vertVelocityTop` is computed after remapping and +:code:`div_hu` is no longer appropriate as it has been remapped. In this case, +the Lagrangian layer thickness is stored in the state variable +:code:`layerThicknessLag` and then the vertical velocity through the top of the +layer is computed as: + +.. code:: + + layerThicknessLag = layerThickness(tlev=2) + + layerThickness(tlev=2) = layerThicknessTarget + + vertVelocityTop = (layerThickness(tlev=2) - layerThicknessLag)/dt + +The function that performs this computation is +:code:`ocn_diagnostic_solve_vertVel_remap`. + +If :code:`normalGMBolusVelocity` is computed based on the remapped ocean state +then the computation of :code:`vertTransportVelocityTop` and +:code:`vertGMBolusVelocitytop` is unchanged as these fields represent Eulerian +velocities. + +However, the current implementation will not compute +:code:`normalGMBolusVelocity` based on the remapped ocean state before +:code:`ocn_diagnostic_solve` is called. Thus, :code:`vertGMBolusVelocityTop` and +:code:`vertTransportVelocityTop` will be inaccurate. This will only pose an +issue when the number of vertical levels changes during regridding, a +capability which isn't included in this development scope. + +Tracer tendencies that are computed as diagnostics will also be inaccurate +after regridding, as they will not be remapped. Remapping these variables does +not make physical sense without also computing vertical tracer fluxes, which +would be overly burdensome. Analysis members that currently use these diagnostic +variables are :code:`mpas_ocn_layer_volume_weighted_averages` and +:code:`mpas_ocn_mixed_layer_heat_budget`. + +Implementation: Grid is sufficiently conditioned for accuracy and stability +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Since the grid is remapped back to z-star in the current implementation, we +assume that the grid is sufficiently regular. + +For future coordinate choices, consdider applying the following after +determining the target grid: + +#. Assign :math:`h_k^{t+1}` to :math:`h_k^{lg}` if :math:`h_k^{t+1} - h_k^{lg}` + is less than a minimum thickness alteration. This is motivated by accuracy + considerations, as each remapping may introduce errors. (Not implemented) +#. Apply minimum layer thickness criterion. (Not implemented) +#. Smooth layers in space and time. (Not implemented) + +Namelist options: + +- Minimum layer thickness +- Minimum thickness change for remapping to occur + + +Implementation: Ability to specify remapping method and its options +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Namelist options: + +- Order of remapping: + :code:`advection, config_vert_remap_order = 1, 2, 3, or 5` +- Interval number of timesteps for remapping: + :code:`advection, config_vert_remap_interval = integer >= 0` +- Slope limiter: + :code:`advection, config_remap_limiter = 'none', 'monotonic', or 'weno'` + +Implementation: Support for existing time stepping schemes +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Currently, only the split-explicit time integration scheme is supported. We +explored the RK4 implementation but it was not found to be stable. + +Implementation: performance +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +OMP calls were implemented for all cell loops in the remapping routine. Our +main strategy for improving performance is remapping an interval of time steps +rather than every time step. For EC60to30 global configurations, an interval of +3-4 timesteps was found to not significantly degrade the solution. + +Options for improving performance (not implemented): + +- Splitting the scalar and momentum timesteps +- Only remapping when the change in thickness exceeds given threshold +- Optimizing/parallelizing PPR + +Implementation: no interference +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2020/12/15 + +Contributors: Darren Engwirda, Xylar Asay-Davis, Carolyn Begeman + +Ensure that no calculations are made outside of remapping itself using +prognostic (or diagnostic) variables both before and after remapping to avoid +introducing additional errors. + +Look for places in the code where prognostic variables are used at the previous +timestep. + +Implementation: modularity +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Remapping operations (PPR) are performed in a separate routine, +:code:`ocn_vert_remap_state`. + +Target grid levels are determined in a separate routine, +:code:`ocn_vert_regrid`. + + +Testing +------- + +Testing and Validation: Ability to perform vertical Lagrangian-remapping +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Ability to handle strong vertical velocities: + +- Baroclinic channel test case + +Evaluating spurious mixing due to remapping: Compare with and without VLR + +- Internal wave test case +- Baroclinic channel test case + + +Testing and Validation: not climate changing +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Global ocean case :code:`GMPAS-JRA1p4.TL319_EC30to60E2r2` 50 year run comparison +between VLR branch with VLR on and master without remapping: +`Link to analysis `_ + +Testing and Validation: performance +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +For the same global ocean case above we show that remapping increases run time +by 20% for remapping every time step and that it scales almost perfectly with as +the remapping interval increases, i.e., 5% increase for remapping every 4 time +steps. We also show that 90% of the remapping time is spent in the PPR routine +:code:`rmap1d`, indicating that this would be the target for performance +improvements. From 27dd2c0d453b2de4c31cf457ca90b0c94cd5d802 Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Mon, 27 Jun 2022 13:32:29 -0500 Subject: [PATCH 21/22] Bld file updates corresponding to Registry changes, from auto scripts --- components/mpas-ocean/bld/build-namelist-section | 9 ++++++--- .../bld/namelist_files/namelist_definition_mpaso.xml | 4 ++-- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index d7d4efe740f0..0a9d5459111c 100644 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -301,11 +301,14 @@ add_default($nl, 'config_land_ice_flux_jenkins_salt_transfer_coefficient'); # Namelist group: advection # ############################# -add_default($nl, 'config_vert_tracer_adv'); -add_default($nl, 'config_vert_tracer_adv_order'); +add_default($nl, 'config_vert_advection_method'); +add_default($nl, 'config_vert_remap_order'); +add_default($nl, 'config_vert_remap_interval'); +add_default($nl, 'config_vert_tracer_adv_flux_order'); add_default($nl, 'config_horiz_tracer_adv_order'); add_default($nl, 'config_coef_3rd_order'); -add_default($nl, 'config_monotonic'); +add_default($nl, 'config_flux_limiter'); +add_default($nl, 'config_remap_limiter'); ############################### # Namelist group: bottom_drag # 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 0e2a2ae1b177..841de7c4e7fc 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -1608,7 +1608,7 @@ Default: Defined in namelist_defaults.xml -Order of polynomial used for tracer reconstruction at layer edges +Order of polynomial used for tracer reconstruction at layer edges for flux-form method Valid values: 2, 3 and 4 Default: Defined in namelist_defaults.xml @@ -2212,7 +2212,7 @@ Default: Defined in namelist_defaults.xml -Enables a change on tracer monotonicity at the end of the monotonic advection routine. Only used if config_monotonic is set to .true. +Enables a change on tracer monotonicity at the end of the monotonic advection routine. Only used if config_flux_limiter is set to monotonic Valid values: .true. or .false. Default: Defined in namelist_defaults.xml From e30698f6eeaff3f673b3f551a649f61283d4e134 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 27 Jun 2022 13:16:20 -0700 Subject: [PATCH 22/22] Replace PPR submodule --- .gitmodules | 3 +++ components/mpas-ocean/src/ppr | 1 + 2 files changed, 4 insertions(+) create mode 160000 components/mpas-ocean/src/ppr diff --git a/.gitmodules b/.gitmodules index 01c7283aa319..0610f98d3e22 100644 --- a/.gitmodules +++ b/.gitmodules @@ -61,3 +61,6 @@ [submodule "components/mpas-ocean/src/FFTW"] path = components/mpas-ocean/src/FFTW url = https://github.com/knbarton/fftw-3.3.8.git +[submodule "components/mpas-ocean/src/ppr"] + path = components/mpas-ocean/src/ppr + url = https://github.com/dengwirda/PPR.git diff --git a/components/mpas-ocean/src/ppr b/components/mpas-ocean/src/ppr new file mode 160000 index 000000000000..3c2246463d5d --- /dev/null +++ b/components/mpas-ocean/src/ppr @@ -0,0 +1 @@ +Subproject commit 3c2246463d5de14f492139987bd84978b84e86ea