Repository: parflow/parflow (language: python)

Issue:
header protection in input_database did not match #ifdef to #define
Pretty obvious, just a typo in which an:
```
#ifdef PROTECTION_HEADER
#define PROTETION
...
#endif
```

didnt' match.  (Note the Hypre and Silo issues were fixed when I went to the github repo instead of the svn one).


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml
index ed11c7d1d..a1ab38b87 100644
--- a/.github/workflows/linux.yml
+++ b/.github/workflows/linux.yml
@@ -46,6 +46,14 @@ jobs:
             backend: "cuda",
             cudaos: 'ubuntu1804'
           }
+        - {
+            name: "Ubuntu 18.04 Kokkos Build",
+            os: ubuntu-18.04,
+            cc: "gcc", cxx: "g++",
+            python: "false",
+            backend: "kokkos",
+            cudaos: 'ubuntu1804'
+          }
         - {
             name: "Ubuntu 16.04",
             os: ubuntu-16.04,
@@ -64,7 +72,7 @@ jobs:
         sudo apt -qq install gfortran libhdf5-openmpi-dev libhdf5-openmpi-103 hdf5-helpers tcl-dev tk-dev
 
     - name: Package Install 18.04 CUDA
-      if: matrix.config.backend == 'cuda' && matrix.config.os == 'ubuntu-18.04'
+      if: (matrix.config.backend == 'cuda' || matrix.config.backend == 'kokkos') && matrix.config.os == 'ubuntu-18.04'
       run: |
         sudo apt-get -qq update
         sudo apt-get -qq install -y software-properties-common
@@ -133,14 +141,14 @@ jobs:
         echo "$HOME/depend/cmake-3.14.0-Linux-x86_64/bin" >> $GITHUB_PATH
 
     - name: NVidia Driver Install
-      if: matrix.config.backend == 'cuda'
+      if: matrix.config.backend == 'cuda' || matrix.config.backend == 'kokkos'
       run: |
         sudo add-apt-repository -y ppa:graphics-drivers
         sudo apt-get update -qq
         sudo apt-get install -y nvidia-driver-440
 
     - name: CUDA Install
-      if: matrix.config.backend == 'cuda'
+      if: matrix.config.backend == 'cuda' || matrix.config.backend == 'kokkos'
       run: |
         CUDA_SHORT=10.1
         CUDA_VER=10.1.105-1
@@ -158,7 +166,7 @@ jobs:
         echo "${CUDA_HOME}/bin" >> $GITHUB_PATH
 
     - name: UCX
-      if: matrix.config.backend == 'cuda'
+      if: matrix.config.backend == 'cuda' || matrix.config.backend == 'kokkos'
       env:
         CACHE_HIT: ${{steps.cache-parflow-dependencies.outputs.cache-hit}}
       run: |
@@ -174,7 +182,7 @@ jobs:
     - name: RMM Install
       env:
         CACHE_HIT: ${{steps.cache-parflow-dependencies.outputs.cache-hit}}
-      if: matrix.config.backend == 'cuda'
+      if: matrix.config.backend == 'cuda' || matrix.config.backend == 'kokkos'
       run: |
         if [[ "$CACHE_HIT" != 'true' ]]; then
           git clone -b branch-0.10 --single-branch --recurse-submodules https://github.com/rapidsai/rmm.git
@@ -182,12 +190,27 @@ jobs:
           cmake . -DCMAKE_INSTALL_PREFIX=$PARFLOW_DEP_DIR
           make -j 2 install
         fi
-        echo "RMM_FLAGS=-DRMM_ROOT=$PARFLOW_DEP_DIR" >> $GITHUB_ENV     
+        echo "RMM_FLAGS=-DRMM_ROOT=$PARFLOW_DEP_DIR" >> $GITHUB_ENV  
+
+    - name: Kokkos Install
+      env:
+        CACHE_HIT: ${{steps.cache-parflow-dependencies.outputs.cache-hit}}
+      if: matrix.config.backend == 'kokkos'
+      run: |
+        if [[ "$CACHE_HIT" != 'true' ]]; then
+          wget -nv https://github.com/kokkos/kokkos/archive/refs/tags/3.3.01.tar.gz
+          tar -xf 3.3.01.tar.gz
+          cd kokkos-3.3.01
+          mkdir build && cd build
+          cmake .. -DCMAKE_CXX_COMPILER=$(pwd)/../bin/nvcc_wrapper -DCMAKE_INSTALL_PREFIX=$PARFLOW_DEP_DIR -DKokkos_ENABLE_CUDA=On -DKokkos_ENABLE_CUDA_LAMBDA=On -DKokkos_ENABLE_CUDA_UVM=On -DKokkos_ARCH_VOLTA70=On  
+          make -j 2 install
+        fi
+        echo "KOKKOS_FLAGS=-DKOKKOS_ROOT=$PARFLOW_DEP_DIR" >> $GITHUB_ENV     
 
     - name: OpenMPI Install for CUDA
       env:
         CACHE_HIT: ${{steps.cache-parflow-dependencies.outputs.cache-hit}}
-      if: matrix.config.backend == 'cuda'
+      if: matrix.config.backend == 'cuda' || matrix.config.backend == 'kokkos'
       run: |
         if [[ "$CACHE_HIT" != 'true' ]]; then
           wget -nv https://download.open-mpi.org/release/open-mpi/v4.0/openmpi-4.0.3.tar.gz
@@ -233,14 +256,14 @@ jobs:
         which cmake
         export OMPI_MCA_rmaps_base_oversubscribe=1
         export OMP_NUM_THREADS=1
-        CC=mpicc CXX=mpicxx F77=mpif77 FC=mpif90 cmake -Bbuild -DCMAKE_BUILD_TYPE=Release -DCMAKE_C_FLAGS="-Wall -Werror -Wno-unused-result -Wno-unused-function" -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_AMPS_LAYER=mpi1 -DMPIEXEC_POSTFLAGS='--oversubscribe' -DPARFLOW_ACCELERATOR_BACKEND=${{ matrix.config.backend }} -DPARFLOW_AMPS_SEQUENTIAL_IO=true -DPARFLOW_HAVE_CLM=ON -DHYPRE_ROOT=$PARFLOW_DEP_DIR -DSILO_ROOT=$PARFLOW_DEP_DIR -DPARFLOW_ENABLE_PYTHON=${{ matrix.config.python }} -DPARFLOW_PYTHON_VIRTUAL_ENV=${{ matrix.config.python }} $RMM_FLAGS -DCMAKE_INSTALL_PREFIX=$PARFLOW_DIR
+        CC=mpicc CXX=mpicxx F77=mpif77 FC=mpif90 cmake -Bbuild -DCMAKE_BUILD_TYPE=Release -DCMAKE_C_FLAGS="-Wall -Werror -Wno-unused-result -Wno-unused-function" -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_AMPS_LAYER=mpi1 -DMPIEXEC_POSTFLAGS='--oversubscribe' -DPARFLOW_ACCELERATOR_BACKEND=${{ matrix.config.backend }} -DPARFLOW_AMPS_SEQUENTIAL_IO=true -DPARFLOW_HAVE_CLM=ON -DHYPRE_ROOT=$PARFLOW_DEP_DIR -DSILO_ROOT=$PARFLOW_DEP_DIR -DPARFLOW_ENABLE_PYTHON=${{ matrix.config.python }} -DPARFLOW_PYTHON_VIRTUAL_ENV=${{ matrix.config.python }} $KOKKOS_FLAGS $RMM_FLAGS -DCMAKE_INSTALL_PREFIX=$PARFLOW_DIR
 
     - name: ParFlow CMake Build
       run: (cd build; make -j 2 install)
 
-    # Can't test with CUDA since have no CUDA hardware on testing nodes, CUDA test is build only
+    # Can't test with GPU since have no GPU hardware on testing nodes, GPU test is build only
     - name: ParFlow CTest
-      if: matrix.config.backend != 'cuda'
+      if: matrix.config.backend != 'cuda' && matrix.config.backend != 'kokkos'
       run: |
         export OMPI_MCA_rmaps_base_oversubscribe=1
         export OMP_NUM_THREADS=1
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1f88cec44..a49dbb3fb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -34,34 +34,54 @@ endforeach ()
 # Set accelerator backend
 #-----------------------------------------------------------------------------
 set(PARFLOW_ACCELERATOR_BACKEND "none" CACHE STRING "Set accelerator backend")
-set_property(CACHE PARFLOW_ACCELERATOR_BACKEND PROPERTY STRINGS none cuda omp)
+set_property(CACHE PARFLOW_ACCELERATOR_BACKEND PROPERTY STRINGS none cuda kokkos omp)
 
 if(PARFLOW_ACCELERATOR_BACKEND STREQUAL "none")
-  set(PARFLOW_ACC_BACKEND 0)
 elseif(PARFLOW_ACCELERATOR_BACKEND STREQUAL "cuda")
-  # Initialize CUDA
-  include(CheckLanguage)
-  check_language(CUDA)
-  if(CMAKE_CUDA_COMPILER)
-    enable_language(CUDA)
-    set(CMAKE_CUDA_STANDARD 14)
-    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -gencode=arch=compute_60,code=sm_60 -gencode=arch=compute_70,code=sm_70 --expt-extended-lambda")
-    if(CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER 11)
-      set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -gencode=arch=compute_80,code=sm_80")
-    endif(CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER 11)
-    set(CMAKE_CUDA_FLAGS_DEBUG "${CMAKE_CUDA_FLAGS_DEBUG} -G --ptxas-options=-v")
-    set(CMAKE_CUDA_FLAGS_RELEASE "-Xptxas ${CMAKE_CUDA_FLAGS_RELEASE}")
-    set(CMAKE_CUDA_HOST_COMPILER mpicxx)
-    message(STATUS "Found CUDA ${CMAKE_CUDA_COMPILER_VERSION} at ${CMAKE_CUDA_COMPILER}")
-    set(PARFLOW_HAVE_CUDA "yes")
-    set(PARFLOW_ACC_BACKEND 1)
-  else(CMAKE_CUDA_COMPILER)
-    message(FATAL_ERROR "ERROR: CUDA installation not found")
-  endif(CMAKE_CUDA_COMPILER)
-
-  # Include RapidsAI Memory Manager (RMM) for pool allocation
-  # @MCB 04/16/2020: This should probably be moved outside of the CUDA block at some point
-  #   so other accelerators have the option to use it as well.
+  set(PARFLOW_HAVE_CUDA "yes")
+  if(NOT DEFINED KOKKOS_ROOT)
+    # Initialize CUDA
+    include(CheckLanguage)
+    check_language(CUDA)
+    if(CMAKE_CUDA_COMPILER)
+      enable_language(CUDA)
+      set(CMAKE_CUDA_STANDARD 14)
+      set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -gencode=arch=compute_60,code=sm_60 -gencode=arch=compute_70,code=sm_70 --expt-extended-lambda")
+      if(CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER 11)
+        set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -gencode=arch=compute_80,code=sm_80")
+      endif(CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER 11)
+      set(CMAKE_CUDA_FLAGS_DEBUG "${CMAKE_CUDA_FLAGS_DEBUG} -G --ptxas-options=-v")
+      set(CMAKE_CUDA_HOST_COMPILER mpicxx)
+      message(STATUS "Found CUDA ${CMAKE_CUDA_COMPILER_VERSION} at ${CMAKE_CUDA_COMPILER}")
+    else(CMAKE_CUDA_COMPILER)
+      message(FATAL_ERROR "ERROR: CUDA installation not found")
+    endif(CMAKE_CUDA_COMPILER)
+  else(NOT DEFINED KOKKOS_ROOT)
+    # Allow mixed use of CUDA and Kokkos for development purposes
+    find_package(Kokkos REQUIRED HINTS ${KOKKOS_ROOT})
+    set(KOKKOS_INCLUDE "${KOKKOS_ROOT}/include")
+    set(PARFLOW_HAVE_KOKKOS "yes")
+  endif(NOT DEFINED KOKKOS_ROOT)
+elseif(PARFLOW_ACCELERATOR_BACKEND STREQUAL "kokkos")
+  if(DEFINED KOKKOS_ROOT)
+    find_package(Kokkos REQUIRED HINTS ${KOKKOS_ROOT})
+    set(KOKKOS_INCLUDE "${KOKKOS_ROOT}/include")
+    set(PARFLOW_HAVE_KOKKOS "yes")
+  else(DEFINED KOKKOS_ROOT)
+    message(FATAL_ERROR "ERROR: KOKKOS_ROOT not specified with PARFLOW_ACCELERATOR_BACKEND=kokkos")
+  endif(DEFINED KOKKOS_ROOT)
+elseif(PARFLOW_ACCELERATOR_BACKEND STREQUAL "omp")
+  message(STATUS "ACCELERATOR: Compiling ParFlow with backend accelerator OpenMP")
+  # Enable C and CXX -fopenmp flag, enable ParFlow defines.
+  set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp")
+  set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
+  set(PARFLOW_HAVE_OMP "yes")
+else()
+  message(FATAL_ERROR "ERROR: Unknown backend type! PARFLOW_ACCELERATOR_BACKEND=${PARFLOW_ACCELERATOR_BACKEND} does not exist!")
+endif()
+
+# Include RapidsAI Memory Manager (RMM) for pool allocation
+if((${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_KOKKOS}))
   if(DEFINED RMM_ROOT)
     find_path(RMM_INCLUDE "rmm" HINTS "${RMM_ROOT}/include")
     find_library(RMM_LIBRARY "rmm" HINTS "${RMM_ROOT}/lib")
@@ -78,20 +98,19 @@ elseif(PARFLOW_ACCELERATOR_BACKEND STREQUAL "cuda")
     else(RMM_INCLUDE AND RMM_LIBRARY)
       message(FATAL_ERROR "ERROR: RMM installation not found")
     endif(RMM_INCLUDE AND RMM_LIBRARY)
-    message(STATUS "ACCELERATOR: Compiling with CUDA and RMM")
+    if(${PARFLOW_HAVE_CUDA}) 
+      message(STATUS "ACCELERATOR: Compiling with CUDA and RMM")
+    else(${PARFLOW_HAVE_CUDA}) 
+      message(STATUS "ACCELERATOR: Compiling with Kokkos and RMM")
+    endif(${PARFLOW_HAVE_CUDA}) 
   else(DEFINED RMM_ROOT)
-    message(STATUS "ACCELERATOR: Compiling with CUDA (RMM not used)")
+    if(${PARFLOW_HAVE_CUDA}) 
+      message(STATUS "ACCELERATOR: Compiling with CUDA (RMM not used)")
+    else(${PARFLOW_HAVE_CUDA}) 
+      message(STATUS "ACCELERATOR: Compiling with Kokkos (RMM not used)")
+    endif(${PARFLOW_HAVE_CUDA}) 
   endif(DEFINED RMM_ROOT)
-elseif(PARFLOW_ACCELERATOR_BACKEND STREQUAL "omp")
-  message(STATUS "ACCELERATOR: Compiling ParFlow with backend accelerator OpenMP")
-  # Enable C and CXX -fopenmp flag, enable ParFlow defines.
-  set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp")
-  set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
-  set(PARFLOW_HAVE_OMP "yes")
-  set(PARFLOW_ACC_BACKEND 2)
-else()
-  message(FATAL_ERROR "ERROR: Unknown backend type! PARFLOW_ACCELERATOR_BACKEND=${PARFLOW_ACCELERATOR_BACKEND} does not exist!")
-endif()
+endif((${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_KOKKOS}))
 
 #-----------------------------------------------------------------------------
 # Version number
diff --git a/README-CUDA.md b/README-CUDA.md
deleted file mode 100644
index cad42999b..000000000
--- a/README-CUDA.md
+++ /dev/null
@@ -1,32 +0,0 @@
-# Building Parflow with GPU acceleration
-
-**WARNING!** Parflow GPU support is still in beta and there may be issues with it.
-
-The GPU acceleration is currently compatible only with NVIDIA GPUs due to CUDA being the only backend option so far. The minimum supported CUDA compute capability for the hardware is 6.0 (NVIDIA Pascal architecture).
-
-Building with CUDA may improve the performance significantly for large problems but is often slower for test cases and small problems due to initialization overhead associated with the GPU use. Currently, the only preconditioner resulting in a good performance is MGSemi. Installation reference can be found in [Ubuntu recipe](/cmake/recipes/ubuntu-18.10-CUDA), [Dockerfile](Dockerfile_CUDA), and [travis.yml](.travis.yml).
-
-
-## CMake
-
-Building with GPU acceleration requires a [CUDA](https://docs.nvidia.com/cuda/cuda-installation-guide-linux/index.html) installation. However, the performance can be further improved by using pool allocation for Unified Memory (requires [RMM v0.10](https://github.com/rapidsai/rmm/tree/branch-0.10) installation) and direct communication between GPUs (requires a CUDA-Aware MPI library).
-The GPU acceleration is activated by specifying *PARFLOW_ACCELERATOR_BACKEND=cuda* option to the CMake, e.g.,
-
-```shell
-cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=cuda
-```
-where *DPARFLOW_AMPS_LAYER=mpi1* leverages GPU-based data packing and unpacking. By default, the packed data is copied to a pinned host staging buffer which is then passed for MPI to avoid special requirements for the MPI library. Direct communication between GPUs (with [GPUDirect P2P/RDMA](https://developer.nvidia.com/gpudirect)) can be activated by specifying an environment variable *PARFLOW_USE_GPUDIRECT=1* during runtime in which case the memory copy between CPU and GPU is avoided and a GPU pointer is passed for MPI, but this requires a CUDA-Aware MPI library (support for Unified Memory is not required because the pointers passed to the MPI library point to pinned GPU memory allocations).
-
-Furthermore, RMM library can be activated by specifying the RMM root directory with *DRMM_ROOT=/path/to/rmm_root* as follows:
-```shell
-cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=cuda -DRMM_ROOT=/path/to/RMM
-```
-## Running Parflow with GPU acceleration
-
-Running Parflow built with GPU support requires that each MPI process has access to a GPU device. Best performance is typically achieved by launching one MPI process per available GPU device. The MPI processes map to the available GPUs by 
-
-```cudaSetDevice(node_local_rank % local_num_devices);```
-
-where node_local_rank and local_num_devices are the node-local rank of the process and the number of GPUs associated with the corresponding node, respectively. Therefore, launching 4 MPI processes per node that has 4 GPUs automatically means that each process uses a different GPU. Launching more processes (than the number of available GPUs) is only supported when using CUDA Multi-Process Service (MPS), but this typically results in reduced performance.
-
-Any existing input script can be run with GPU acceleration; no changes are necessary. However, it is noted that two subsequent runs of the same input script using the same compiled executable are not guaranteed to produce identical results. This is expected behavior due to atomic operations performed by the GPU device (i.e., the order of floating-point operations may change between two subsequent runs).
diff --git a/README-GPU.md b/README-GPU.md
new file mode 100644
index 000000000..533870dcf
--- /dev/null
+++ b/README-GPU.md
@@ -0,0 +1,40 @@
+# Building Parflow with GPU acceleration
+
+**WARNING!** Parflow GPU support is still in beta and there may be issues with it.
+
+The GPU acceleration is currently compatible only with NVIDIA GPUs either using CUDA directly or using CUDA through the Kokkos library. The minimum supported CUDA compute capability for the hardware is 6.0 (NVIDIA Pascal architecture).
+
+Building with CUDA or Kokkos may improve the performance significantly for large problems but is often slower for test cases and small problems due to initialization overhead associated with the GPU use. Currently, the only preconditioner resulting in a good performance is MGSemi. Installation reference can be found in [Ubuntu recipe](/cmake/recipes/ubuntu-18.10-CUDA), [Dockerfile](Dockerfile_CUDA), and [linux.yml](/.github/workflows/linux.yml).
+
+
+## CMake
+
+Building with GPU acceleration requires a [CUDA](https://docs.nvidia.com/cuda/cuda-installation-guide-linux/index.html) installation, and in case Kokkos is used, [Kokkos](https://github.com/kokkos/kokkos) installation. However, the performance can be further improved by using pool allocation for Unified Memory (requires [RMM v0.10](https://github.com/rapidsai/rmm/tree/branch-0.10) installation) and direct communication between GPUs (requires using a CUDA-Aware MPI library).
+The GPU acceleration is activated by specifying either *PARFLOW_ACCELERATOR_BACKEND=cuda* option to the CMake, i.e.,
+
+```shell
+cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=cuda
+```
+or *PARFLOW_ACCELERATOR_BACKEND=kokkos* and *DKokkos_ROOT=/path/to/Kokkos* i.e.,
+```shell
+cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=kokkos -DKokkos_ROOT=/path/to/Kokkos
+```
+where *DPARFLOW_AMPS_LAYER=mpi1* leverages GPU-based data packing and unpacking. By default, the packed data is copied to a host staging buffer which is then passed for MPI to avoid special requirements for the MPI library. Direct communication between GPUs (with [GPUDirect P2P/RDMA](https://developer.nvidia.com/gpudirect)) can be activated by specifying an environment variable *PARFLOW_USE_GPUDIRECT=1* during runtime in which case the memory copy between CPU and GPU is avoided and a GPU pointer is passed for MPI, but this requires a CUDA-Aware MPI library (support for Unified Memory is not required with the native CUDA backend because the pointers passed to the MPI library point to pinned GPU memory allocations, but is required with the Kokkos backend).
+
+Furthermore, RMM library can be activated by specifying the RMM root directory with *DRMM_ROOT=/path/to/rmm_root* as follows:
+```shell
+cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=cuda -DRMM_ROOT=/path/to/RMM
+```
+or
+```shell
+cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=kokkos -DKokkos_ROOT=/path/to/Kokkos -DRMM_ROOT=/path/to/RMM
+```
+## Running Parflow with GPU acceleration
+
+Running Parflow built with GPU support requires that each MPI process has access to a GPU device. Best performance is typically achieved by launching one MPI process per available GPU device. The MPI processes map to the available GPUs by 
+
+```cudaSetDevice(node_local_rank % local_num_devices);```
+
+where node_local_rank and local_num_devices are the node-local rank of the process and the number of GPUs associated with the corresponding node, respectively. Therefore, launching 4 MPI processes per node that has 4 GPUs automatically means that each process uses a different GPU. Launching more processes (than the number of available GPUs) is only supported when using CUDA Multi-Process Service (MPS), but this typically results in reduced performance.
+
+Any existing input script can be run with GPU acceleration; no changes are necessary. However, it is noted that two subsequent runs of the same input script using the same compiled executable are not guaranteed to produce identical results. This is expected behavior due to atomic operations performed by the GPU device (i.e., the order of floating-point operations may change between two subsequent runs).
diff --git a/README.md b/README.md
index 715019fbe..0d3e8014f 100644
--- a/README.md
+++ b/README.md
@@ -153,7 +153,7 @@ https://cmake.org/runningcmake/
 #### Building with the cmake command line
 
 CMake may also be configured from the command line using the cmake
-command. Instructions to build with different accelerator backends are found from separate documents: [CUDA](README-CUDA.md), [OpenMP](README-OPENMP.md). The default will configure a sequential version of ParFlow
+command. Instructions to build with different accelerator backends are found from the following documents: [CUDA, KOKKOS](README-GPU.md), [OpenMP](README-OPENMP.md). The default will configure a sequential version of ParFlow
 using MPI libraries.  CLM is being enabled.
 
 ```shell
diff --git a/cmake/parflow_config.h.in b/cmake/parflow_config.h.in
index 51bff97d6..311e74cef 100644
--- a/cmake/parflow_config.h.in
+++ b/cmake/parflow_config.h.in
@@ -66,6 +66,8 @@
 
 #cmakedefine PARFLOW_HAVE_CUDA
 
+#cmakedefine PARFLOW_HAVE_KOKKOS
+
 /* PARFLOW_HAVE_RMM is not defined here because because it is only set 
    for the compilation of parflow_exe, parflow_lib, and kinsol sources. 
    RMM must be initialized before it is used, so, for example, 
@@ -73,9 +75,4 @@
 
 #cmakedefine PARFLOW_HAVE_OMP
 
-#define PARFLOW_BACKEND_NONE 0
-#define PARFLOW_BACKEND_CUDA 1
-#define PARFLOW_BACKEND_OMP 2
-#define PARFLOW_ACC_BACKEND @PARFLOW_ACC_BACKEND@
-
 #endif // PARFLOW_CONFIG_H
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 7c8bd44f1..2c02c5702 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -8,7 +8,7 @@ if((${PARFLOW_HAVE_HYPRE}) AND (${PARFLOW_HAVE_SILO}))
 endif()
 
 #This test is suspected to be too precision sensitive
-if((${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_OMP}))
+if((${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_KOKKOS}) OR (${PARFLOW_HAVE_OMP}))
   list(REMOVE_ITEM TESTS test_XPlusYPlusZ.tcl)
 endif()
 
diff --git a/pfsimulator/amps/CMakeLists.txt b/pfsimulator/amps/CMakeLists.txt
index 4ead402c4..0575f9198 100644
--- a/pfsimulator/amps/CMakeLists.txt
+++ b/pfsimulator/amps/CMakeLists.txt
@@ -14,12 +14,18 @@ add_library(amps STATIC ${COMMON_SRC_FILES} ${AMPS_SRC_FILES})
 
 if(${PARFLOW_HAVE_CUDA})
   target_include_directories(amps PUBLIC ${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
-  if( ${PARFLOW_HAVE_RMM} )
-    target_link_libraries(amps rmm)
-    target_include_directories(amps PUBLIC ${RMM_INCLUDE})
-  endif( ${PARFLOW_HAVE_RMM} )
 endif(${PARFLOW_HAVE_CUDA})
 
+if( ${PARFLOW_HAVE_RMM} )
+  target_link_libraries(amps rmm)
+  target_include_directories(amps PUBLIC ${RMM_INCLUDE})
+endif( ${PARFLOW_HAVE_RMM} )
+
+if( ${PARFLOW_HAVE_KOKKOS} )
+  target_link_libraries(amps Kokkos::kokkos)
+  target_include_directories(amps PUBLIC ${KOKKOS_INCLUDE})
+endif( ${PARFLOW_HAVE_KOKKOS} )
+
 if (${PARFLOW_HAVE_MPI})
   target_include_directories (amps PUBLIC "${MPI_C_INCLUDE_PATH}")
 endif (${PARFLOW_HAVE_MPI})
diff --git a/pfsimulator/amps/mpi1/CMakeLists.cmake b/pfsimulator/amps/mpi1/CMakeLists.cmake
index 5618f018f..e81959b33 100644
--- a/pfsimulator/amps/mpi1/CMakeLists.cmake
+++ b/pfsimulator/amps/mpi1/CMakeLists.cmake
@@ -1,4 +1,3 @@
-
 set(AMPS_SRC_FILES
   amps_allreduce.c
   amps_bcast.c
@@ -20,6 +19,10 @@ set(AMPS_SRC_FILES
   amps_vector.c
   )
 
-if(${PARFLOW_HAVE_CUDA})
+if((${PARFLOW_HAVE_CUDA}) AND (NOT (${PARFLOW_HAVE_KOKKOS})))
   list(APPEND AMPS_SRC_FILES amps_gpupacking.cu)
-endif(${PARFLOW_HAVE_CUDA})
+endif((${PARFLOW_HAVE_CUDA}) AND (NOT (${PARFLOW_HAVE_KOKKOS})))
+
+if(${PARFLOW_HAVE_KOKKOS})
+  list(APPEND AMPS_SRC_FILES amps_gpupacking.cpp)
+endif(${PARFLOW_HAVE_KOKKOS})
diff --git a/pfsimulator/amps/mpi1/amps.h b/pfsimulator/amps/mpi1/amps.h
index e6642fdfb..c5269e05a 100644
--- a/pfsimulator/amps/mpi1/amps.h
+++ b/pfsimulator/amps/mpi1/amps.h
@@ -48,15 +48,16 @@
 
 #include <stdbool.h>
 #include <stdio.h>
+#include <string.h>
 #include <sys/times.h>
 
 #ifdef PARFLOW_HAVE_CUDA
 #include <cuda.h>
 #include <cuda_runtime.h>
+#endif
 #ifdef PARFLOW_HAVE_RMM
 #include <rmm/rmm_api.h>
 #endif
-#endif
 
 /*
  * Prevent inclusion of mpi C++ bindings in mpi.h includes.
@@ -1033,7 +1034,7 @@ void amps_ReadDouble(amps_File file, double *ptr, int len);
 #define amps_Error(name, type, comment, operation) \
   printf("%s : %s\n", name, comment)
 
-#ifdef PARFLOW_HAVE_CUDA
+#if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS)
 /*--------------------------------------------------------------------------
  * Amps defines with CUDA
  *--------------------------------------------------------------------------*/
@@ -1056,6 +1057,7 @@ void amps_ReadDouble(amps_File file, double *ptr, int len);
  */
 #define AMPS_MPI_NOT_USE_PERSISTENT
 
+#ifdef PARFLOW_HAVE_CUDA
 /*--------------------------------------------------------------------------
  *  GPU error handling macros
  *--------------------------------------------------------------------------*/
@@ -1074,6 +1076,7 @@ static inline void amps_cuda_error(cudaError_t err, const char *file, int line)
 		exit(1);
 	}
 }
+#endif // PARFLOW_HAVE_CUDA
 
 #ifdef PARFLOW_HAVE_RMM
 /**
@@ -1090,12 +1093,32 @@ static inline void amps_rmm_error(rmmError_t err, const char *file, int line) {
 		exit(1);
 	}
 }
-#endif
+#endif // PARFLOW_HAVE_RMM
 
 /*--------------------------------------------------------------------------
- * Define static unified memory allocation routines for CUDA
+ * Define static unified memory allocation routines for devices
  *--------------------------------------------------------------------------*/
 
+/**
+ * @brief Kokkos C wrapper declaration for memory allocation.
+ */
+void* kokkosUVMAlloc(size_t size);
+
+/**
+ * @brief Kokkos C wrapper declaration for memory deallocation.
+ */
+void kokkosUVMFree(void *ptr);
+
+/**
+ * @brief Kokkos C wrapper declaration for memory copy.
+ */
+void kokkosMemCpyUVMToUVM(char *dest, char *src, size_t size);
+
+/**
+ * @brief Kokkos C wrapper declaration for memset.
+ */
+void kokkosMemSetAmps(char *ptr, size_t size);
+
 /**
  * @brief Allocates unified memory
  * 
@@ -1106,13 +1129,15 @@ static inline void amps_rmm_error(rmmError_t err, const char *file, int line) {
  * @param size bytes to be allocated [IN]
  * @return a void pointer to the allocated dataspace
  */
-static inline void *_amps_talloc_cuda(size_t size)
+static inline void *_amps_talloc_device(size_t size)
 {
   void *ptr = NULL;  
   
 #ifdef PARFLOW_HAVE_RMM
   RMM_ERRCHK(rmmAlloc(&ptr,size,0,__FILE__,__LINE__));
-#else
+#elif defined(PARFLOW_HAVE_KOKKOS)
+  ptr = kokkosUVMAlloc(size);
+#elif defined(PARFLOW_HAVE_CUDA)
   CUDA_ERRCHK(cudaMallocManaged((void**)&ptr, size, cudaMemAttachGlobal));
   // CUDA_ERRCHK(cudaHostAlloc((void**)&ptr, size, cudaHostAllocMapped));  
 #endif
@@ -1130,18 +1155,25 @@ static inline void *_amps_talloc_cuda(size_t size)
  * @param size bytes to be allocated [IN]
  * @return a void pointer to the allocated dataspace
  */
-static inline void *_amps_ctalloc_cuda(size_t size)
+static inline void *_amps_ctalloc_device(size_t size)
 {
   void *ptr = NULL;  
 
 #ifdef PARFLOW_HAVE_RMM
   RMM_ERRCHK(rmmAlloc(&ptr,size,0,__FILE__,__LINE__));
-#else
+#elif defined(PARFLOW_HAVE_KOKKOS)
+  ptr = kokkosUVMAlloc(size);
+#elif defined(PARFLOW_HAVE_CUDA)
   CUDA_ERRCHK(cudaMallocManaged((void**)&ptr, size, cudaMemAttachGlobal));
   // CUDA_ERRCHK(cudaHostAlloc((void**)&ptr, size, cudaHostAllocMapped));
 #endif  
-  // memset(ptr, 0, size);
+
+#if defined(PARFLOW_HAVE_CUDA)
   CUDA_ERRCHK(cudaMemset(ptr, 0, size));  
+#else
+  // memset(ptr, 0, size);
+  kokkosMemSetAmps((char*)ptr, size);
+#endif
   
   return ptr;
 }
@@ -1153,32 +1185,34 @@ static inline void *_amps_ctalloc_cuda(size_t size)
  *
  * @param ptr a void pointer to the allocated dataspace [IN]
  */
-static inline void _amps_tfree_cuda(void *ptr)
+static inline void _amps_tfree_device(void *ptr)
 {
 #ifdef PARFLOW_HAVE_RMM
   RMM_ERRCHK(rmmFree(ptr,0,__FILE__,__LINE__));
-#else
+#elif defined(PARFLOW_HAVE_KOKKOS)
+  kokkosUVMFree(ptr);
+#elif defined(PARFLOW_HAVE_CUDA)
   CUDA_ERRCHK(cudaFree(ptr));
   // CUDA_ERRCHK(cudaFreeHost(ptr));
 #endif
 }
 
 /** 
- * Same as \ref amps_TAlloc but allocates managed memory (CUDA required) 
+ * Same as \ref amps_TAlloc but allocates managed memory
  */
-#define amps_TAlloc_managed(type, count) ((count>0) ? (type*)_amps_talloc_cuda((unsigned int)(sizeof(type) * (count))) : NULL)
+#define amps_TAlloc_managed(type, count) ((count>0) ? (type*)_amps_talloc_device((unsigned int)(sizeof(type) * (count))) : NULL)
 
 /** 
- * Same as \ref amps_CTAlloc but allocates managed memory (CUDA required) 
+ * Same as \ref amps_CTAlloc but allocates managed memory
  */
-#define amps_CTAlloc_managed(type, count) ((count) ? (type*)_amps_ctalloc_cuda((unsigned int)(sizeof(type) * (count))) : NULL)
+#define amps_CTAlloc_managed(type, count) ((count) ? (type*)_amps_ctalloc_device((unsigned int)(sizeof(type) * (count))) : NULL)
 
 /** 
- * Same as \ref amps_TFree but deallocates managed memory (CUDA required) 
+ * Same as \ref amps_TFree but deallocates managed memory
  */
-#define amps_TFree_managed(ptr) if (ptr) _amps_tfree_cuda(ptr); else {}
+#define amps_TFree_managed(ptr) if (ptr) _amps_tfree_device(ptr); else {}
 
-#endif // PARFLOW_HAVE_CUDA
+#endif // PARFLOW_HAVE_CUDA || PARFLOW_HAVE_KOKKOS
 
 #include "amps_proto.h"
 
diff --git a/pfsimulator/amps/mpi1/amps_exchange.c b/pfsimulator/amps/mpi1/amps_exchange.c
index 2cecb2439..29fb845f0 100644
--- a/pfsimulator/amps/mpi1/amps_exchange.c
+++ b/pfsimulator/amps/mpi1/amps_exchange.c
@@ -29,7 +29,7 @@
 #include "amps.h"
 
 /* This CUDA stuff could be combined with AMPS_MPI_NOT_USE_PERSISTENT case */
-#ifdef PARFLOW_HAVE_CUDA
+#if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS)
 
 void _amps_wait_exchange(amps_Handle handle)
 {
diff --git a/pfsimulator/amps/mpi1/amps_finalize.c b/pfsimulator/amps/mpi1/amps_finalize.c
index 9b8bbda7c..f89b19917 100644
--- a/pfsimulator/amps/mpi1/amps_finalize.c
+++ b/pfsimulator/amps/mpi1/amps_finalize.c
@@ -70,9 +70,8 @@ int amps_Finalize()
 
     MPI_Finalize();
   }
-#ifdef PARFLOW_HAVE_CUDA
-  amps_gpu_free_bufs();
-  amps_gpu_destroy_streams();
+#if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS)
+  amps_gpu_finalize();
 #endif
 
 #ifdef AMPS_MALLOC_DEBUG
diff --git a/pfsimulator/amps/mpi1/amps_gpupacking.cpp b/pfsimulator/amps/mpi1/amps_gpupacking.cpp
new file mode 100644
index 000000000..e4892a210
--- /dev/null
+++ b/pfsimulator/amps/mpi1/amps_gpupacking.cpp
@@ -0,0 +1,693 @@
+/**********************************************************************
+ *
+ *  Please read the LICENSE file for the GNU Lesser General Public License.
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License (as published
+ *  by the Free Software Foundation) version 2.1 dated February 1999.
+ *
+ *  This program is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the IMPLIED WARRANTY OF
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the terms
+ *  and conditions of the GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Lesser General Public
+ *  License along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+ *  USA
+ ***********************************************************************/
+
+#include <Kokkos_Core.hpp>
+
+extern "C"{
+
+#include <stdlib.h>
+#include <string.h>
+#include "amps.h"
+
+/* Disable GPU packing (useful for debugging) */
+// #define DISABLE_GPU_PACKING
+
+/**
+ * @brief A global GPU staging buffer for amps_exchange
+ * 
+ * A separate instance is allocated for recv and send operations
+ * which are shared between all amps packages. The space is allocated
+ * according to the largest invoice size.
+ */
+ typedef struct _amps_GpuBuffer {
+  /** 
+  * A list of GPU pointers for the packed data for each recv/send invoice. 
+  */
+  char **buf;
+  /** 
+  * A list of CPU pointers for the packed data for each recv/send invoice. 
+  * Only allocated if host staging is used.
+  */
+  char **buf_host;
+  /** 
+  * A list of buf sizes for each invoice (applies to *buf and *host_buf).
+  */
+  int *buf_size;
+  /** 
+  * The number of buffers allocated (applies to *buf and *host_buf).
+  */
+  int num_bufs;
+} amps_GpuBuffer;
+
+#ifdef DISABLE_GPU_PACKING
+/* Dummy definitions if no GPU packing */
+void amps_gpu_finalize(){}
+
+void amps_gpu_sync_streams(int id){
+  (void)id;
+}
+
+int amps_gpupacking(int action, amps_Invoice inv, int inv_num, char **buffer_out, int *size_out){
+  (void)action;
+  (void)inv;
+  (void)inv_num;
+  (void)buffer_out;
+  (void)size_out;
+  return 1;
+}    
+#else
+
+amps_GpuBuffer amps_gpu_recvbuf = 
+{
+  .buf = NULL, 
+  .buf_host = NULL,
+  .buf_size = NULL,
+  .num_bufs = 0
+};
+
+amps_GpuBuffer amps_gpu_sendbuf = 
+{
+  .buf = NULL, 
+  .buf_host = NULL,
+  .buf_size = NULL,
+  .num_bufs = 0
+};
+
+/**
+ * @brief Defines whether host staging is used
+ *
+ * 0: Use page-locked device staging buffer and pass GPU pointer
+ * to MPI library
+ *
+ * 1: Use page-locked host staging buffer and pass CPU pointer
+ * to MPI library (default)
+ */
+static int ENFORCE_HOST_STAGING = 1;
+
+/**
+ * @brief Check for GPUDirect environment variable
+ * 
+ * If PARFLOW_USE_GPUDIRECT=1 environment variable is found,
+ * disable host staging.
+ *
+ */
+void amps_gpu_check_env(){
+  if(getenv("PARFLOW_USE_GPUDIRECT") != NULL){
+    if(atoi(getenv("PARFLOW_USE_GPUDIRECT")) == 1){
+      if(amps_rank == 0 && ENFORCE_HOST_STAGING != 0)
+        printf("Node %d: Using GPUDirect (GPU-Aware MPI required)\n", amps_rank);
+      ENFORCE_HOST_STAGING = 0;
+    }
+  }
+}
+
+/**
+ * @brief Destroy all GPU streams and free allocations
+ */
+void amps_gpu_destroy_streams(){
+}
+
+/**
+ * @brief Synchronizes GPU streams associated with the id
+ *
+ * Kokkos does not support multiple streams (only CUDA backend does).
+ *
+ * @param id The integer id associated with the synchronized streams [IN]
+ */
+void amps_gpu_sync_streams(int id){
+  if(id == 0) Kokkos::fence();
+}
+
+/**
+ * @brief Allocate device memory with Kokkos
+ *
+ * @param size The size of the allocation in bytes
+ */
+void* kokkosDeviceAlloc(size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  return Kokkos::kokkos_malloc<Kokkos::CudaSpace>(size);
+#else
+  return Kokkos::kokkos_malloc(size);
+#endif
+}
+
+/**
+ * @brief Free device memory with Kokkos
+ *
+ * @param ptr Freed pointer
+ */
+void kokkosDeviceFree(void *ptr){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::kokkos_free<Kokkos::CudaSpace>(ptr);
+#else
+  Kokkos::kokkos_free(ptr);
+#endif
+}
+
+/**
+ * @brief Allocate pinned host memory with Kokkos
+ *
+ * @param size The size of the allocation in bytes
+ */
+void* kokkosHostAlloc(size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  return Kokkos::kokkos_malloc<Kokkos::CudaHostPinnedSpace>(size);
+#else
+  return Kokkos::kokkos_malloc<Kokkos::HostSpace>(size);
+#endif
+}
+
+/**
+ * @brief Free pinned host memory with Kokkos
+ *
+ * @param ptr Freed pointer
+ */
+void kokkosHostFree(void *ptr){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::kokkos_free<Kokkos::CudaHostPinnedSpace>(ptr);
+#else
+  Kokkos::kokkos_free<Kokkos::HostSpace>(ptr);
+#endif
+}
+
+/**
+ * @brief Allocate Unified Memory with Kokkos
+ *
+ * @param size The size of the allocation in bytes
+ */
+void* kokkosUVMAlloc(size_t size){
+  return Kokkos::kokkos_malloc(size);
+}
+
+/**
+ * @brief Free Unified Memory with Kokkos
+ *
+ * @param ptr Freed pointer
+ */
+void kokkosUVMFree(void *ptr){
+  Kokkos::kokkos_free(ptr);
+}
+
+/**
+ * @brief Device-device memcopy
+ *
+ * @param dest Destination pointer
+ * @param src Source pointer
+ * @param size Bytes to be copied
+ */
+void kokkosMemCpyDeviceToDevice(char *dest, char *src, size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::View<char*, Kokkos::CudaSpace> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::CudaSpace> src_view(src, size);
+#else
+  Kokkos::View<char*> dest_view(dest, size);
+  Kokkos::View<char*> src_view(src, size);
+#endif
+  Kokkos::deep_copy(dest_view, src_view);
+}
+
+
+/**
+ * @brief Device-host memcopy
+ *
+ * @param dest Destination pointer
+ * @param src Source pointer
+ * @param size Bytes to be copied
+ */
+void kokkosMemCpyDeviceToHost(char *dest, char *src, size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::View<char*, Kokkos::CudaHostPinnedSpace> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::CudaSpace> src_view(src, size);
+#else
+  Kokkos::View<char*, Kokkos::HostSpace> dest_view(dest, size);
+  Kokkos::View<char*> src_view(src, size);
+#endif
+  Kokkos::deep_copy(dest_view, src_view);
+}
+
+/**
+ * @brief Host-device memcopy
+ *
+ * @param dest Destination pointer
+ * @param src Source pointer
+ * @param size Bytes to be copied
+ */
+void kokkosMemCpyHostToDevice(char *dest, char *src, size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::View<char*, Kokkos::CudaSpace> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::CudaHostPinnedSpace> src_view(src, size);
+#else
+  Kokkos::View<char*> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::HostSpace> src_view(src, size);
+#endif
+  Kokkos::deep_copy(dest_view, src_view);
+}
+
+/**
+ * @brief Host-host memcopy
+ *
+ * @param dest Destination pointer
+ * @param src Source pointer
+ * @param size Bytes to be copied
+ */
+void kokkosMemCpyHostToHost(char *dest, char *src, size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::View<char*, Kokkos::CudaHostPinnedSpace> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::CudaHostPinnedSpace> src_view(src, size);
+#else
+  Kokkos::View<char*, Kokkos::HostSpace> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::HostSpace> src_view(src, size);
+#endif
+  Kokkos::deep_copy(dest_view, src_view);
+}
+
+/**
+ * @brief UVM-UVM memcopy
+ *
+ * @param dest Destination pointer
+ * @param src Source pointer
+ * @param size Bytes to be copied
+ */
+void kokkosMemCpyUVMToUVM(char *dest, char *src, size_t size){
+  Kokkos::View<char*> dest_view(dest, size);
+  Kokkos::View<char*> src_view(src, size);
+  Kokkos::deep_copy(dest_view, src_view);
+}
+
+/**
+ * @brief Kokkos memset
+ *
+ * @param ptr Pointer for the data to be set to zero
+ * @param size Bytes to be zeroed
+ */
+void kokkosMemSetAmps(char *ptr, size_t size){
+  /* Loop style initialization */
+  if(size % sizeof(int))
+  {
+    /* Writing 1 byte / thread is slow */
+    Kokkos::parallel_for(size, KOKKOS_LAMBDA(int i){ptr[i] = 0;});
+  }
+  else
+  {
+    /* Writing 4 bytes / thread is fast */
+    Kokkos::parallel_for(size / sizeof(int), KOKKOS_LAMBDA(int i){((int*)ptr)[i] = 0;});
+  }
+  Kokkos::fence(); 
+
+  /* Deep_copy style initialization for char* should be fast in future Kokkos releases */
+  // Kokkos::View<char*> ptr_view(ptr, size);
+  // Kokkos::deep_copy(ptr_view, 0);
+}
+  
+/**
+ * @brief Free the staging buffers associated with gpubuf
+ *
+ * @param gpubuf A pointer to the amps_GpuBuffer [IN]
+ */
+static void _amps_gpubuf_free(amps_GpuBuffer *gpubuf){
+  if(gpubuf->num_bufs > 0){
+    for(int i = 0; i < gpubuf->num_bufs; i++){
+      if(gpubuf->buf_size[i] > 0){
+        kokkosDeviceFree(gpubuf->buf[i]);
+        if(ENFORCE_HOST_STAGING)
+          kokkosHostFree(gpubuf->buf_host[i]);
+      }
+    }
+    free(gpubuf->buf);
+    free(gpubuf->buf_size);
+    if(ENFORCE_HOST_STAGING)
+      free(gpubuf->buf_host);
+  }
+}
+
+/**
+ * @brief Free the global recv and send staging buffers
+ */
+void amps_gpu_free_bufs(){
+  _amps_gpubuf_free(&amps_gpu_recvbuf);
+  _amps_gpubuf_free(&amps_gpu_sendbuf);
+}
+
+/**
+ * @brief Finalize GPU resource usage
+ */
+void amps_gpu_finalize(){
+  amps_gpu_free_bufs();
+  amps_gpu_destroy_streams();
+  if(Kokkos::is_initialized) Kokkos::finalize();
+}
+
+/**
+ * @brief Initialize Kokkos if not initialized
+ */
+void amps_kokkos_initialization(){
+  if(!Kokkos::is_initialized()){
+    Kokkos::InitArguments args;
+    args.ndevices = 1;
+    Kokkos::initialize(args);  
+  }
+}
+
+/**
+ * @brief Allocate/reallocate page-locked GPU/CPU staging buffers.
+ *
+ * Page-locked host memory allocation mirrors the GPU allocation
+ * if host staging is enabled.
+ *
+ * @param gpubuf A pointer to the amps_GpuBuffer [IN]
+ * @param id The id associated with the buffer (typically amps invoice index) [IN]
+ * @param pos The position in bytes up to which the buffer is already filled [IN]
+ * @param size The total required buffer size subtracted by pos [IN]
+ * @return A pointer to the allocated/reallocated GPU buffer
+ */
+static char* _amps_gpubuf_realloc(amps_GpuBuffer *gpubuf, int id, int pos, int size){
+  /* Check if arrays are big enough for id */
+  if(id >= gpubuf->num_bufs){
+    if(gpubuf->num_bufs == 0){
+      amps_gpu_check_env();
+    }
+    if(id != gpubuf->num_bufs){
+      printf("ERROR at %s:%d: Unexpected id\n", __FILE__, __LINE__);
+      return NULL;
+    }
+
+    char **buf_array = (char**)malloc((id + 1) * sizeof(char*));
+    int *size_array = (int*)malloc((id + 1) * sizeof(int));
+    if(gpubuf->num_bufs != 0){
+      memcpy(buf_array, 
+             gpubuf->buf, 
+               gpubuf->num_bufs * sizeof(char*));
+      memcpy(size_array, 
+             gpubuf->buf_size, 
+               gpubuf->num_bufs * sizeof(int));
+      free(gpubuf->buf);
+      free(gpubuf->buf_size);
+    }
+    gpubuf->buf = buf_array;
+    gpubuf->buf_size = size_array;
+    gpubuf->buf_size[id] = 0;
+
+    if(ENFORCE_HOST_STAGING){
+      char **buf_host_array = (char**)malloc((id + 1) * sizeof(char*));
+      if(gpubuf->num_bufs != 0){
+        memcpy(buf_host_array, 
+               gpubuf->buf_host, 
+                 gpubuf->num_bufs * sizeof(char*));
+        free(gpubuf->buf_host);
+      }
+      gpubuf->buf_host = buf_host_array;
+    }
+    gpubuf->num_bufs++;
+  }
+
+  /* Check to see if enough space is allocated for id */
+  int size_total = pos + size;
+  if (gpubuf->buf_size[id] < size_total)
+  {
+    char *newbuf = (char*)kokkosDeviceAlloc(size_total);
+    if(pos != 0){
+      kokkosMemCpyDeviceToDevice(newbuf,
+        gpubuf->buf[id], pos);
+      kokkosDeviceFree(gpubuf->buf[id]);
+    }
+    gpubuf->buf[id] = newbuf;
+    
+    if(ENFORCE_HOST_STAGING){
+      char *newbuf_host = (char*)kokkosHostAlloc(size_total);
+      if(pos != 0){
+        kokkosMemCpyHostToHost(newbuf_host,
+          gpubuf->buf_host[id], pos);
+        kokkosHostFree(gpubuf->buf_host[id]);
+      }
+      gpubuf->buf_host[id] = newbuf_host;
+    }
+    gpubuf->buf_size[id] = size_total;
+  }
+  return gpubuf->buf[id] + pos;
+}
+
+/**
+ * @brief Get recv buffer associated with id to be passed for MPI
+ *
+ * @param id The id of the recv buffer [IN]
+ * @return A pointer to the recv staging buffer
+ */
+static char* _amps_gpu_recvbuf(int id){
+  if(id >= amps_gpu_recvbuf.num_bufs){
+    printf("ERROR at %s:%d: Not enough space is allocated for recvbufs\n", __FILE__, __LINE__);
+    exit(1);
+  }
+  if(ENFORCE_HOST_STAGING)
+    return amps_gpu_recvbuf.buf_host[id];
+  else
+    return amps_gpu_recvbuf.buf[id];
+}
+
+/**
+ * @brief Allocate/reallocate recv GPU/CPU staging buffers
+ *
+ * @param id The id associated with the buffer (typically amps invoice index) [IN]
+ * @param pos The position in bytes up to which the buffer is already filled [IN]
+ * @param size The total required buffer size subtracted by pos [IN]
+ * @return A pointer to the allocated/reallocated GPU recv buffer
+ */
+static char* _amps_gpu_recvbuf_realloc(int id, int pos, int size){
+  return _amps_gpubuf_realloc(&amps_gpu_recvbuf, id, pos, size);
+}
+
+/**
+ * @brief Get send buffer associated with id to be passed for MPI
+ *
+ * @param id The id of the send buffer [IN]
+ * @return A pointer to the send staging buffer
+ */
+static char* _amps_gpu_sendbuf(int id){
+  if(id >= amps_gpu_sendbuf.num_bufs){
+    printf("ERROR at %s:%d: Not enough space is allocated for sendbufs\n", __FILE__, __LINE__);
+    exit(1);
+  }
+  if(ENFORCE_HOST_STAGING)
+    return amps_gpu_sendbuf.buf_host[id];
+  else
+    return amps_gpu_sendbuf.buf[id];
+}
+
+/**
+ * @brief Allocate/reallocate send GPU/CPU staging buffers
+ *
+ * @param id The id associated with the buffer (typically amps invoice index) [IN]
+ * @param pos The position in bytes up to which the buffer is already filled [IN]
+ * @param size The total required buffer size subtracted by pos [IN]
+ * @return A pointer to the allocated/reallocated GPU send buffer
+ */
+static char* _amps_gpu_sendbuf_realloc(int id, int pos, int size){
+  return _amps_gpubuf_realloc(&amps_gpu_sendbuf, id, pos, size);
+}
+
+/**
+* @brief The main amps GPU packing driver function
+*
+* Operates in 3 different modes:
+*
+* Mode 1: AMPS_GETRBUF / AMPS_GETSBUF
+* -determines the size of the invoice message in bytes
+* -allocates/reallocates the recv/send staging buffer if necessary
+* -*buffer_out is set to point to the recv/send staging buffer
+* -*size_out is set to the size of the invoice message in bytes
+*
+* Mode 2: AMPS_PACK (requires calling amps_gpu_sync_streams() before the packed data is accessed)
+* -determines the size of the invoice message in bytes
+* -allocates/reallocates the send staging buffer if necessary
+* -packs data to the send staging buffer parallel using a GPU
+* -*buffer_out is set to point to the send staging buffer
+* -*size_out is set to the size of invoice message in bytes (size of the packed data)
+*
+* Mode 3: AMPS_UNPACK (requires calling amps_gpu_sync_streams() before the unpacked data is accessed)
+* -determines the size of the invoice message in bytes
+* -allocates/reallocates the recv staging buffer if necessary
+* -unpacks data from the recv staging buffer parallel using a GPU
+* -*buffer_out is set to point to the recv staging buffer
+* -*size_out is set to the size of invoice message in bytes (size of the unpacked data)
+*
+* @param action either AMPS_GETRBUF, AMPS_GETSBUF, AMPS_PACK or AMPS_UNPACK [IN]
+* @param inv amps invoice [IN]
+* @param inv_num amps invoice order number [IN]
+* @param buffer_out pointer to the pointer of the staging buffer [OUT]
+* @param size_out pointer to the invoice message size in bytes [OUT]
+* @return error code (line number), 0 if succesful
+*/
+int amps_gpupacking(int action, amps_Invoice inv, int inv_num, char **buffer_out, int *size_out){
+  
+  char *buffer;
+  int pos = 0;
+
+  amps_kokkos_initialization();
+
+  amps_InvoiceEntry *ptr = inv->list;
+  while (ptr != NULL)
+  {
+    if (ptr->len_type != AMPS_INVOICE_POINTER){
+      printf("ERROR at %s:%d: ptr->len_type must be a pointer\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+  
+    if (ptr->stride_type != AMPS_INVOICE_POINTER){
+      printf("ERROR at %s:%d: ptr->stride_type must be a pointer\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+  
+    if (ptr->data_type == AMPS_INVOICE_POINTER){
+      printf("ERROR at %s:%d: Overlayed invoices not supported\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+  
+    /* Make sure amps invoice vector type is supported */
+    if(ptr->type - AMPS_INVOICE_LAST_CTYPE != AMPS_INVOICE_DOUBLE_CTYPE){
+      printf("ERROR at %s:%d: Only \"AMPS_INVOICE_DOUBLE_CTYPE\" is supported\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+  
+    /* Preparations for the kernel launch */
+    int len_x = 1;
+    int len_y = 1;
+    int len_z = 1;  
+    int stride_x = 0;
+    int stride_y = 0;
+    int stride_z = 0;
+  
+    int dim = (ptr->dim_type == AMPS_INVOICE_POINTER) ? *(ptr->ptr_dim) : ptr->dim;
+    switch (dim)
+    {
+      case 3:
+        len_z = ptr->ptr_len[2];
+        stride_z = ptr->ptr_stride[2];
+      case 2:
+        len_y = ptr->ptr_len[1];
+        stride_y = ptr->ptr_stride[1];
+      case 1:
+        len_x = ptr->ptr_len[0];
+        stride_x = ptr->ptr_stride[0];
+        break;
+      default:
+        printf("ERROR at %s:%d: Only dimensions 1 - 3 are supported\n", __FILE__, __LINE__);
+        return __LINE__;
+    }  
+    
+    /* Get data location and its properties */
+    char *data = (char*)ptr->data;
+
+#ifdef PARFLOW_HAVE_CUDA
+    struct cudaPointerAttributes attributes;
+    cudaPointerGetAttributes(&attributes, (void *)data);  
+  
+    /* Check that the data location (not MPI staging buffer) is accessible by the GPU */
+    if(cudaGetLastError() != cudaSuccess || attributes.type < 2){
+      printf("ERROR at %s:%d: The data location (not MPI staging buffer) is not accessible by the GPU(s)\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+#endif
+
+    int size = len_x * len_y * len_z * sizeof(double);
+    if((action == AMPS_GETSBUF) || (action == AMPS_PACK)){
+      buffer = _amps_gpu_sendbuf_realloc(inv_num, pos, size);
+    }
+    else if((action == AMPS_GETRBUF) || (action == AMPS_UNPACK)){
+      buffer = _amps_gpu_recvbuf_realloc(inv_num, pos, size);
+    }
+    else{
+      printf("ERROR at %s:%d: Unknown action argument (val = %d)\n", __FILE__, __LINE__, action);
+      return __LINE__;
+    }
+
+#ifdef PARFLOW_HAVE_CUDA
+    /* Get buffer location and its properties */
+    cudaPointerGetAttributes(&attributes, (void *)buffer);  
+
+    /* Check that the staging buffer is accessible by the GPU */
+    if(cudaGetLastError() != cudaSuccess || attributes.type < 2){
+      printf("ERROR at %s:%d: The MPI staging buffer location is not accessible by the GPU(s)\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+#endif
+
+    /* Run packing or unpacking kernel */
+    using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;
+    MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{len_x, len_y, len_z}});
+
+    if(action == AMPS_PACK){
+      Kokkos::parallel_for(mdpolicy_3d, KOKKOS_LAMBDA(int i, int j, int k)
+      {
+        double *ptr_buf = (double*)buffer;
+        double *ptr_data = (double*)data;
+        *(ptr_buf + k * len_y * len_x + j * len_x + i) = 
+          *(ptr_data + k * (stride_z + (len_y - 1) * stride_y + len_y * (len_x - 1) * stride_x) + 
+            j * (stride_y + (len_x - 1) * stride_x) + i * stride_x);                                          
+      });
+      if(ENFORCE_HOST_STAGING){
+        /* Copy device buffer to host after packing */
+        kokkosMemCpyDeviceToHost(amps_gpu_sendbuf.buf_host[inv_num] + pos,
+                                    amps_gpu_sendbuf.buf[inv_num] + pos, size);
+      }
+      inv->flags |= AMPS_PACKED;
+    }
+    else if(action == AMPS_UNPACK){
+      if(ENFORCE_HOST_STAGING){
+        /* Copy host buffer to device before unpacking */
+        kokkosMemCpyHostToDevice(amps_gpu_recvbuf.buf[inv_num] + pos,
+                                    amps_gpu_recvbuf.buf_host[inv_num] + pos, size);
+      }
+      Kokkos::parallel_for(mdpolicy_3d, KOKKOS_LAMBDA(int i, int j, int k)
+      {
+        double *ptr_buf = (double*)buffer;
+        double *ptr_data = (double*)data;
+        *(ptr_data + k * (stride_z + (len_y - 1) * stride_y + len_y * (len_x - 1) * stride_x) + 
+          j * (stride_y + (len_x - 1) * stride_x) + i * stride_x) = 
+            *(ptr_buf + k * len_y * len_x + j * len_x + i);
+      });
+      inv->flags &= ~AMPS_PACKED;
+    }
+
+    pos += size;  
+    ptr = ptr->next;
+  }
+
+  /* Check that the size is calculated right */
+  // if(pos != amps_sizeof_invoice(amps_CommWorld, inv)){
+    // printf("ERROR at %s:%d: The size does not match the invoice size\n", __FILE__, __LINE__);
+    // return __LINE__;
+  // }
+
+  /* Set the out values here if everything went fine */
+  if((action == AMPS_GETSBUF) || (action == AMPS_PACK)){
+    *buffer_out = _amps_gpu_sendbuf(inv_num);
+  }
+  else if((action == AMPS_GETRBUF) || (action == AMPS_UNPACK)){
+    *buffer_out = _amps_gpu_recvbuf(inv_num);
+  }
+  else{
+    printf("ERROR at %s:%d: Unknown action argument (val = %d)\n", __FILE__, __LINE__, action);
+    return __LINE__;
+  }
+
+  *size_out = pos;
+
+  return 0;
+} 
+#endif // DISABLE_GPU_PACKING
+}   
diff --git a/pfsimulator/amps/mpi1/amps_gpupacking.cu b/pfsimulator/amps/mpi1/amps_gpupacking.cu
index 4ceb14100..5ff8706a8 100644
--- a/pfsimulator/amps/mpi1/amps_gpupacking.cu
+++ b/pfsimulator/amps/mpi1/amps_gpupacking.cu
@@ -106,9 +106,7 @@ typedef struct _amps_GpuStreams {
 
 #ifdef DISABLE_GPU_PACKING
 /* Dummy definitions if no GPU packing */
-void amps_gpu_free_bufs(){}
-
-void amps_gpu_destroy_streams(){}
+void amps_gpu_finalize(){}
 
 void amps_gpu_sync_streams(int id){
   (void)id;
@@ -286,6 +284,14 @@ void amps_gpu_free_bufs(){
   _amps_gpubuf_free(&amps_gpu_sendbuf);
 }
 
+/**
+ * @brief Finalize GPU resource usage
+ */
+ void amps_gpu_finalize(){
+  amps_gpu_free_bufs();
+  amps_gpu_destroy_streams();
+}
+
 /**
  * @brief Allocate/reallocate page-locked GPU/CPU staging buffers.
  *
diff --git a/pfsimulator/amps/mpi1/amps_proto.h b/pfsimulator/amps/mpi1/amps_proto.h
index 6999a9d26..9dcb25a0d 100644
--- a/pfsimulator/amps/mpi1/amps_proto.h
+++ b/pfsimulator/amps/mpi1/amps_proto.h
@@ -34,8 +34,7 @@ void amps_FindPowers(int N, int *log, int *Nnext, int *Nprev);
 amps_File amps_Fopen(char *filename, char *type);
 
 /* amps_gpupacking.cu */
-void amps_gpu_free_bufs();
-void amps_gpu_destroy_streams();
+void amps_gpu_finalize();
 void amps_gpu_sync_streams(int id);
 int amps_gpupacking(int action, amps_Invoice inv, int inv_num, char **buffer_out, int *size_out);
 
@@ -59,7 +58,6 @@ void amps_ScanInt(amps_File file, int *data, int len, int stride);
 void amps_ScanLong(amps_File file, long *data, int len, int stride);
 void amps_ScanFloat(amps_File file, float *data, int len, int stride);
 void amps_ScanDouble(amps_File file, double *data, int len, int stride);
-
 #ifndef CASC_HAVE_BIGENDIAN
 void amps_WriteDouble(amps_File file, double *ptr, int len);
 void amps_ReadDouble(amps_File file, double *ptr, int len);
@@ -67,7 +65,6 @@ void amps_WriteInt(amps_File file, int *ptr, int len);
 void amps_ReadInt(amps_File file, int *ptr, int len);
 #endif
 
-
 /* amps_irecv.c */
 amps_Handle amps_IRecv(amps_Comm comm, int source, amps_Invoice invoice);
 
diff --git a/pfsimulator/amps/oas3/CMakeLists.cmake b/pfsimulator/amps/oas3/CMakeLists.cmake
index 4d6e86f0c..8996c8b86 100644
--- a/pfsimulator/amps/oas3/CMakeLists.cmake
+++ b/pfsimulator/amps/oas3/CMakeLists.cmake
@@ -27,6 +27,10 @@ set(AMPS_SRC_FILES
   send_fld2_clm.F90
   )
 
-if(${PARFLOW_HAVE_CUDA})
+if((${PARFLOW_HAVE_CUDA}) AND (NOT (${PARFLOW_HAVE_KOKKOS})))
   list(APPEND AMPS_SRC_FILES amps_gpupacking.cu)
-endif(${PARFLOW_HAVE_CUDA})
\ No newline at end of file
+endif((${PARFLOW_HAVE_CUDA}) AND (NOT (${PARFLOW_HAVE_KOKKOS})))
+
+if(${PARFLOW_HAVE_KOKKOS})
+  list(APPEND AMPS_SRC_FILES amps_gpupacking.cpp)
+endif(${PARFLOW_HAVE_KOKKOS})
diff --git a/pfsimulator/amps/oas3/amps.h b/pfsimulator/amps/oas3/amps.h
index cec982523..9c676382c 100644
--- a/pfsimulator/amps/oas3/amps.h
+++ b/pfsimulator/amps/oas3/amps.h
@@ -48,15 +48,16 @@
 
 #include <stdbool.h>
 #include <stdio.h>
+#include <string.h>
 #include <sys/times.h>
 
 #ifdef PARFLOW_HAVE_CUDA
 #include <cuda.h>
 #include <cuda_runtime.h>
+#endif
 #ifdef PARFLOW_HAVE_RMM
 #include <rmm/rmm_api.h>
 #endif
-#endif
 
 /*
  * Prevent inclusion of mpi C++ bindings in mpi.h includes.
@@ -1065,7 +1066,7 @@ void amps_ReadDouble(amps_File file, double *ptr, int len);
 #define amps_Error(name, type, comment, operation) \
   printf("%s : %s\n", name, comment)
 
-#ifdef PARFLOW_HAVE_CUDA
+#if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS)
 /*--------------------------------------------------------------------------
  * Amps defines with CUDA
  *--------------------------------------------------------------------------*/
@@ -1088,6 +1089,7 @@ void amps_ReadDouble(amps_File file, double *ptr, int len);
  */
 #define AMPS_MPI_NOT_USE_PERSISTENT
 
+#ifdef PARFLOW_HAVE_CUDA
 /*--------------------------------------------------------------------------
  *  GPU error handling macros
  *--------------------------------------------------------------------------*/
@@ -1106,6 +1108,7 @@ static inline void amps_cuda_error(cudaError_t err, const char *file, int line)
 		exit(1);
 	}
 }
+#endif // PARFLOW_HAVE_CUDA
 
 #ifdef PARFLOW_HAVE_RMM
 /**
@@ -1122,12 +1125,32 @@ static inline void amps_rmm_error(rmmError_t err, const char *file, int line) {
 		exit(1);
 	}
 }
-#endif
+#endif // PARFLOW_HAVE_RMM
 
 /*--------------------------------------------------------------------------
- * Define static unified memory allocation routines for CUDA
+ * Define static unified memory allocation routines for devices
  *--------------------------------------------------------------------------*/
 
+/**
+ * @brief Kokkos C wrapper declaration for memory allocation.
+ */
+void* kokkosUVMAlloc(size_t size);
+
+/**
+ * @brief Kokkos C wrapper declaration for memory deallocation.
+ */
+void kokkosUVMFree(void *ptr);
+
+/**
+ * @brief Kokkos C wrapper declaration for memory copy.
+ */
+void kokkosMemCpyUVMToUVM(char *dest, char *src, size_t size);
+
+/**
+ * @brief Kokkos C wrapper declaration for memset.
+ */
+void kokkosMemSetAmps(char *ptr, size_t size);
+
 /**
  * @brief Allocates unified memory
  * 
@@ -1138,13 +1161,15 @@ static inline void amps_rmm_error(rmmError_t err, const char *file, int line) {
  * @param size bytes to be allocated [IN]
  * @return a void pointer to the allocated dataspace
  */
-static inline void *_amps_talloc_cuda(size_t size)
+static inline void *_amps_talloc_device(size_t size)
 {
   void *ptr = NULL;  
   
 #ifdef PARFLOW_HAVE_RMM
   RMM_ERRCHK(rmmAlloc(&ptr,size,0,__FILE__,__LINE__));
-#else
+#elif defined(PARFLOW_HAVE_KOKKOS)
+  ptr = kokkosUVMAlloc(size);
+#elif defined(PARFLOW_HAVE_CUDA)
   CUDA_ERRCHK(cudaMallocManaged((void**)&ptr, size, cudaMemAttachGlobal));
   // CUDA_ERRCHK(cudaHostAlloc((void**)&ptr, size, cudaHostAllocMapped));  
 #endif
@@ -1162,18 +1187,25 @@ static inline void *_amps_talloc_cuda(size_t size)
  * @param size bytes to be allocated [IN]
  * @return a void pointer to the allocated dataspace
  */
-static inline void *_amps_ctalloc_cuda(size_t size)
+static inline void *_amps_ctalloc_device(size_t size)
 {
   void *ptr = NULL;  
 
 #ifdef PARFLOW_HAVE_RMM
   RMM_ERRCHK(rmmAlloc(&ptr,size,0,__FILE__,__LINE__));
-#else
+#elif defined(PARFLOW_HAVE_KOKKOS)
+  ptr = kokkosUVMAlloc(size);
+#elif defined(PARFLOW_HAVE_CUDA)
   CUDA_ERRCHK(cudaMallocManaged((void**)&ptr, size, cudaMemAttachGlobal));
   // CUDA_ERRCHK(cudaHostAlloc((void**)&ptr, size, cudaHostAllocMapped));
 #endif  
-  // memset(ptr, 0, size);
+
+#if defined(PARFLOW_HAVE_CUDA)
   CUDA_ERRCHK(cudaMemset(ptr, 0, size));  
+#else
+  // memset(ptr, 0, size);
+  kokkosMemSetAmps((char*)ptr, size);
+#endif
   
   return ptr;
 }
@@ -1185,32 +1217,34 @@ static inline void *_amps_ctalloc_cuda(size_t size)
  *
  * @param ptr a void pointer to the allocated dataspace [IN]
  */
-static inline void _amps_tfree_cuda(void *ptr)
+static inline void _amps_tfree_device(void *ptr)
 {
 #ifdef PARFLOW_HAVE_RMM
   RMM_ERRCHK(rmmFree(ptr,0,__FILE__,__LINE__));
-#else
+#elif defined(PARFLOW_HAVE_KOKKOS)
+  kokkosUVMFree(ptr);
+#elif defined(PARFLOW_HAVE_CUDA)
   CUDA_ERRCHK(cudaFree(ptr));
   // CUDA_ERRCHK(cudaFreeHost(ptr));
 #endif
 }
 
 /** 
- * Same as \ref amps_TAlloc but allocates managed memory (CUDA required) 
+ * Same as \ref amps_TAlloc but allocates managed memory
  */
-#define amps_TAlloc_managed(type, count) ((count>0) ? (type*)_amps_talloc_cuda((unsigned int)(sizeof(type) * (count))) : NULL)
+#define amps_TAlloc_managed(type, count) ((count>0) ? (type*)_amps_talloc_device((unsigned int)(sizeof(type) * (count))) : NULL)
 
 /** 
- * Same as \ref amps_CTAlloc but allocates managed memory (CUDA required) 
+ * Same as \ref amps_CTAlloc but allocates managed memory
  */
-#define amps_CTAlloc_managed(type, count) ((count) ? (type*)_amps_ctalloc_cuda((unsigned int)(sizeof(type) * (count))) : NULL)
+#define amps_CTAlloc_managed(type, count) ((count) ? (type*)_amps_ctalloc_device((unsigned int)(sizeof(type) * (count))) : NULL)
 
 /** 
- * Same as \ref amps_TFree but deallocates managed memory (CUDA required) 
+ * Same as \ref amps_TFree but deallocates managed memory
  */
-#define amps_TFree_managed(ptr) if (ptr) _amps_tfree_cuda(ptr); else {}
+#define amps_TFree_managed(ptr) if (ptr) _amps_tfree_device(ptr); else {}
 
-#endif // PARFLOW_HAVE_CUDA
+#endif // PARFLOW_HAVE_CUDA || PARFLOW_HAVE_KOKKOS
 
 #include "amps_proto.h"
 #include "oas3_external.h"
diff --git a/pfsimulator/amps/oas3/amps_exchange.c b/pfsimulator/amps/oas3/amps_exchange.c
index fd27730b7..700daa3c9 100644
--- a/pfsimulator/amps/oas3/amps_exchange.c
+++ b/pfsimulator/amps/oas3/amps_exchange.c
@@ -29,7 +29,7 @@
 #include "amps.h"
 
 /* This CUDA stuff could be combined with AMPS_MPI_NOT_USE_PERSISTENT case */
-#ifdef PARFLOW_HAVE_CUDA
+#if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS)
 
 void _amps_wait_exchange(amps_Handle handle)
 {
diff --git a/pfsimulator/amps/oas3/amps_finalize.c b/pfsimulator/amps/oas3/amps_finalize.c
index 7bf65ea5f..e6b314d24 100644
--- a/pfsimulator/amps/oas3/amps_finalize.c
+++ b/pfsimulator/amps/oas3/amps_finalize.c
@@ -69,9 +69,8 @@ int amps_Finalize()
     MPI_Comm_free(&amps_CommWrite);
     CALL_oas_pfl_finalize(&dummy2_oas3);
   }
-#ifdef PARFLOW_HAVE_CUDA
-  amps_gpu_free_bufs();
-  amps_gpu_destroy_streams();
+#if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS)
+  amps_gpu_finalize();
 #endif
 
 #ifdef AMPS_MALLOC_DEBUG
diff --git a/pfsimulator/amps/oas3/amps_gpupacking.cpp b/pfsimulator/amps/oas3/amps_gpupacking.cpp
new file mode 100644
index 000000000..e4892a210
--- /dev/null
+++ b/pfsimulator/amps/oas3/amps_gpupacking.cpp
@@ -0,0 +1,693 @@
+/**********************************************************************
+ *
+ *  Please read the LICENSE file for the GNU Lesser General Public License.
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License (as published
+ *  by the Free Software Foundation) version 2.1 dated February 1999.
+ *
+ *  This program is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the IMPLIED WARRANTY OF
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the terms
+ *  and conditions of the GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Lesser General Public
+ *  License along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+ *  USA
+ ***********************************************************************/
+
+#include <Kokkos_Core.hpp>
+
+extern "C"{
+
+#include <stdlib.h>
+#include <string.h>
+#include "amps.h"
+
+/* Disable GPU packing (useful for debugging) */
+// #define DISABLE_GPU_PACKING
+
+/**
+ * @brief A global GPU staging buffer for amps_exchange
+ * 
+ * A separate instance is allocated for recv and send operations
+ * which are shared between all amps packages. The space is allocated
+ * according to the largest invoice size.
+ */
+ typedef struct _amps_GpuBuffer {
+  /** 
+  * A list of GPU pointers for the packed data for each recv/send invoice. 
+  */
+  char **buf;
+  /** 
+  * A list of CPU pointers for the packed data for each recv/send invoice. 
+  * Only allocated if host staging is used.
+  */
+  char **buf_host;
+  /** 
+  * A list of buf sizes for each invoice (applies to *buf and *host_buf).
+  */
+  int *buf_size;
+  /** 
+  * The number of buffers allocated (applies to *buf and *host_buf).
+  */
+  int num_bufs;
+} amps_GpuBuffer;
+
+#ifdef DISABLE_GPU_PACKING
+/* Dummy definitions if no GPU packing */
+void amps_gpu_finalize(){}
+
+void amps_gpu_sync_streams(int id){
+  (void)id;
+}
+
+int amps_gpupacking(int action, amps_Invoice inv, int inv_num, char **buffer_out, int *size_out){
+  (void)action;
+  (void)inv;
+  (void)inv_num;
+  (void)buffer_out;
+  (void)size_out;
+  return 1;
+}    
+#else
+
+amps_GpuBuffer amps_gpu_recvbuf = 
+{
+  .buf = NULL, 
+  .buf_host = NULL,
+  .buf_size = NULL,
+  .num_bufs = 0
+};
+
+amps_GpuBuffer amps_gpu_sendbuf = 
+{
+  .buf = NULL, 
+  .buf_host = NULL,
+  .buf_size = NULL,
+  .num_bufs = 0
+};
+
+/**
+ * @brief Defines whether host staging is used
+ *
+ * 0: Use page-locked device staging buffer and pass GPU pointer
+ * to MPI library
+ *
+ * 1: Use page-locked host staging buffer and pass CPU pointer
+ * to MPI library (default)
+ */
+static int ENFORCE_HOST_STAGING = 1;
+
+/**
+ * @brief Check for GPUDirect environment variable
+ * 
+ * If PARFLOW_USE_GPUDIRECT=1 environment variable is found,
+ * disable host staging.
+ *
+ */
+void amps_gpu_check_env(){
+  if(getenv("PARFLOW_USE_GPUDIRECT") != NULL){
+    if(atoi(getenv("PARFLOW_USE_GPUDIRECT")) == 1){
+      if(amps_rank == 0 && ENFORCE_HOST_STAGING != 0)
+        printf("Node %d: Using GPUDirect (GPU-Aware MPI required)\n", amps_rank);
+      ENFORCE_HOST_STAGING = 0;
+    }
+  }
+}
+
+/**
+ * @brief Destroy all GPU streams and free allocations
+ */
+void amps_gpu_destroy_streams(){
+}
+
+/**
+ * @brief Synchronizes GPU streams associated with the id
+ *
+ * Kokkos does not support multiple streams (only CUDA backend does).
+ *
+ * @param id The integer id associated with the synchronized streams [IN]
+ */
+void amps_gpu_sync_streams(int id){
+  if(id == 0) Kokkos::fence();
+}
+
+/**
+ * @brief Allocate device memory with Kokkos
+ *
+ * @param size The size of the allocation in bytes
+ */
+void* kokkosDeviceAlloc(size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  return Kokkos::kokkos_malloc<Kokkos::CudaSpace>(size);
+#else
+  return Kokkos::kokkos_malloc(size);
+#endif
+}
+
+/**
+ * @brief Free device memory with Kokkos
+ *
+ * @param ptr Freed pointer
+ */
+void kokkosDeviceFree(void *ptr){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::kokkos_free<Kokkos::CudaSpace>(ptr);
+#else
+  Kokkos::kokkos_free(ptr);
+#endif
+}
+
+/**
+ * @brief Allocate pinned host memory with Kokkos
+ *
+ * @param size The size of the allocation in bytes
+ */
+void* kokkosHostAlloc(size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  return Kokkos::kokkos_malloc<Kokkos::CudaHostPinnedSpace>(size);
+#else
+  return Kokkos::kokkos_malloc<Kokkos::HostSpace>(size);
+#endif
+}
+
+/**
+ * @brief Free pinned host memory with Kokkos
+ *
+ * @param ptr Freed pointer
+ */
+void kokkosHostFree(void *ptr){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::kokkos_free<Kokkos::CudaHostPinnedSpace>(ptr);
+#else
+  Kokkos::kokkos_free<Kokkos::HostSpace>(ptr);
+#endif
+}
+
+/**
+ * @brief Allocate Unified Memory with Kokkos
+ *
+ * @param size The size of the allocation in bytes
+ */
+void* kokkosUVMAlloc(size_t size){
+  return Kokkos::kokkos_malloc(size);
+}
+
+/**
+ * @brief Free Unified Memory with Kokkos
+ *
+ * @param ptr Freed pointer
+ */
+void kokkosUVMFree(void *ptr){
+  Kokkos::kokkos_free(ptr);
+}
+
+/**
+ * @brief Device-device memcopy
+ *
+ * @param dest Destination pointer
+ * @param src Source pointer
+ * @param size Bytes to be copied
+ */
+void kokkosMemCpyDeviceToDevice(char *dest, char *src, size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::View<char*, Kokkos::CudaSpace> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::CudaSpace> src_view(src, size);
+#else
+  Kokkos::View<char*> dest_view(dest, size);
+  Kokkos::View<char*> src_view(src, size);
+#endif
+  Kokkos::deep_copy(dest_view, src_view);
+}
+
+
+/**
+ * @brief Device-host memcopy
+ *
+ * @param dest Destination pointer
+ * @param src Source pointer
+ * @param size Bytes to be copied
+ */
+void kokkosMemCpyDeviceToHost(char *dest, char *src, size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::View<char*, Kokkos::CudaHostPinnedSpace> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::CudaSpace> src_view(src, size);
+#else
+  Kokkos::View<char*, Kokkos::HostSpace> dest_view(dest, size);
+  Kokkos::View<char*> src_view(src, size);
+#endif
+  Kokkos::deep_copy(dest_view, src_view);
+}
+
+/**
+ * @brief Host-device memcopy
+ *
+ * @param dest Destination pointer
+ * @param src Source pointer
+ * @param size Bytes to be copied
+ */
+void kokkosMemCpyHostToDevice(char *dest, char *src, size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::View<char*, Kokkos::CudaSpace> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::CudaHostPinnedSpace> src_view(src, size);
+#else
+  Kokkos::View<char*> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::HostSpace> src_view(src, size);
+#endif
+  Kokkos::deep_copy(dest_view, src_view);
+}
+
+/**
+ * @brief Host-host memcopy
+ *
+ * @param dest Destination pointer
+ * @param src Source pointer
+ * @param size Bytes to be copied
+ */
+void kokkosMemCpyHostToHost(char *dest, char *src, size_t size){
+#ifdef PARFLOW_HAVE_CUDA
+  Kokkos::View<char*, Kokkos::CudaHostPinnedSpace> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::CudaHostPinnedSpace> src_view(src, size);
+#else
+  Kokkos::View<char*, Kokkos::HostSpace> dest_view(dest, size);
+  Kokkos::View<char*, Kokkos::HostSpace> src_view(src, size);
+#endif
+  Kokkos::deep_copy(dest_view, src_view);
+}
+
+/**
+ * @brief UVM-UVM memcopy
+ *
+ * @param dest Destination pointer
+ * @param src Source pointer
+ * @param size Bytes to be copied
+ */
+void kokkosMemCpyUVMToUVM(char *dest, char *src, size_t size){
+  Kokkos::View<char*> dest_view(dest, size);
+  Kokkos::View<char*> src_view(src, size);
+  Kokkos::deep_copy(dest_view, src_view);
+}
+
+/**
+ * @brief Kokkos memset
+ *
+ * @param ptr Pointer for the data to be set to zero
+ * @param size Bytes to be zeroed
+ */
+void kokkosMemSetAmps(char *ptr, size_t size){
+  /* Loop style initialization */
+  if(size % sizeof(int))
+  {
+    /* Writing 1 byte / thread is slow */
+    Kokkos::parallel_for(size, KOKKOS_LAMBDA(int i){ptr[i] = 0;});
+  }
+  else
+  {
+    /* Writing 4 bytes / thread is fast */
+    Kokkos::parallel_for(size / sizeof(int), KOKKOS_LAMBDA(int i){((int*)ptr)[i] = 0;});
+  }
+  Kokkos::fence(); 
+
+  /* Deep_copy style initialization for char* should be fast in future Kokkos releases */
+  // Kokkos::View<char*> ptr_view(ptr, size);
+  // Kokkos::deep_copy(ptr_view, 0);
+}
+  
+/**
+ * @brief Free the staging buffers associated with gpubuf
+ *
+ * @param gpubuf A pointer to the amps_GpuBuffer [IN]
+ */
+static void _amps_gpubuf_free(amps_GpuBuffer *gpubuf){
+  if(gpubuf->num_bufs > 0){
+    for(int i = 0; i < gpubuf->num_bufs; i++){
+      if(gpubuf->buf_size[i] > 0){
+        kokkosDeviceFree(gpubuf->buf[i]);
+        if(ENFORCE_HOST_STAGING)
+          kokkosHostFree(gpubuf->buf_host[i]);
+      }
+    }
+    free(gpubuf->buf);
+    free(gpubuf->buf_size);
+    if(ENFORCE_HOST_STAGING)
+      free(gpubuf->buf_host);
+  }
+}
+
+/**
+ * @brief Free the global recv and send staging buffers
+ */
+void amps_gpu_free_bufs(){
+  _amps_gpubuf_free(&amps_gpu_recvbuf);
+  _amps_gpubuf_free(&amps_gpu_sendbuf);
+}
+
+/**
+ * @brief Finalize GPU resource usage
+ */
+void amps_gpu_finalize(){
+  amps_gpu_free_bufs();
+  amps_gpu_destroy_streams();
+  if(Kokkos::is_initialized) Kokkos::finalize();
+}
+
+/**
+ * @brief Initialize Kokkos if not initialized
+ */
+void amps_kokkos_initialization(){
+  if(!Kokkos::is_initialized()){
+    Kokkos::InitArguments args;
+    args.ndevices = 1;
+    Kokkos::initialize(args);  
+  }
+}
+
+/**
+ * @brief Allocate/reallocate page-locked GPU/CPU staging buffers.
+ *
+ * Page-locked host memory allocation mirrors the GPU allocation
+ * if host staging is enabled.
+ *
+ * @param gpubuf A pointer to the amps_GpuBuffer [IN]
+ * @param id The id associated with the buffer (typically amps invoice index) [IN]
+ * @param pos The position in bytes up to which the buffer is already filled [IN]
+ * @param size The total required buffer size subtracted by pos [IN]
+ * @return A pointer to the allocated/reallocated GPU buffer
+ */
+static char* _amps_gpubuf_realloc(amps_GpuBuffer *gpubuf, int id, int pos, int size){
+  /* Check if arrays are big enough for id */
+  if(id >= gpubuf->num_bufs){
+    if(gpubuf->num_bufs == 0){
+      amps_gpu_check_env();
+    }
+    if(id != gpubuf->num_bufs){
+      printf("ERROR at %s:%d: Unexpected id\n", __FILE__, __LINE__);
+      return NULL;
+    }
+
+    char **buf_array = (char**)malloc((id + 1) * sizeof(char*));
+    int *size_array = (int*)malloc((id + 1) * sizeof(int));
+    if(gpubuf->num_bufs != 0){
+      memcpy(buf_array, 
+             gpubuf->buf, 
+               gpubuf->num_bufs * sizeof(char*));
+      memcpy(size_array, 
+             gpubuf->buf_size, 
+               gpubuf->num_bufs * sizeof(int));
+      free(gpubuf->buf);
+      free(gpubuf->buf_size);
+    }
+    gpubuf->buf = buf_array;
+    gpubuf->buf_size = size_array;
+    gpubuf->buf_size[id] = 0;
+
+    if(ENFORCE_HOST_STAGING){
+      char **buf_host_array = (char**)malloc((id + 1) * sizeof(char*));
+      if(gpubuf->num_bufs != 0){
+        memcpy(buf_host_array, 
+               gpubuf->buf_host, 
+                 gpubuf->num_bufs * sizeof(char*));
+        free(gpubuf->buf_host);
+      }
+      gpubuf->buf_host = buf_host_array;
+    }
+    gpubuf->num_bufs++;
+  }
+
+  /* Check to see if enough space is allocated for id */
+  int size_total = pos + size;
+  if (gpubuf->buf_size[id] < size_total)
+  {
+    char *newbuf = (char*)kokkosDeviceAlloc(size_total);
+    if(pos != 0){
+      kokkosMemCpyDeviceToDevice(newbuf,
+        gpubuf->buf[id], pos);
+      kokkosDeviceFree(gpubuf->buf[id]);
+    }
+    gpubuf->buf[id] = newbuf;
+    
+    if(ENFORCE_HOST_STAGING){
+      char *newbuf_host = (char*)kokkosHostAlloc(size_total);
+      if(pos != 0){
+        kokkosMemCpyHostToHost(newbuf_host,
+          gpubuf->buf_host[id], pos);
+        kokkosHostFree(gpubuf->buf_host[id]);
+      }
+      gpubuf->buf_host[id] = newbuf_host;
+    }
+    gpubuf->buf_size[id] = size_total;
+  }
+  return gpubuf->buf[id] + pos;
+}
+
+/**
+ * @brief Get recv buffer associated with id to be passed for MPI
+ *
+ * @param id The id of the recv buffer [IN]
+ * @return A pointer to the recv staging buffer
+ */
+static char* _amps_gpu_recvbuf(int id){
+  if(id >= amps_gpu_recvbuf.num_bufs){
+    printf("ERROR at %s:%d: Not enough space is allocated for recvbufs\n", __FILE__, __LINE__);
+    exit(1);
+  }
+  if(ENFORCE_HOST_STAGING)
+    return amps_gpu_recvbuf.buf_host[id];
+  else
+    return amps_gpu_recvbuf.buf[id];
+}
+
+/**
+ * @brief Allocate/reallocate recv GPU/CPU staging buffers
+ *
+ * @param id The id associated with the buffer (typically amps invoice index) [IN]
+ * @param pos The position in bytes up to which the buffer is already filled [IN]
+ * @param size The total required buffer size subtracted by pos [IN]
+ * @return A pointer to the allocated/reallocated GPU recv buffer
+ */
+static char* _amps_gpu_recvbuf_realloc(int id, int pos, int size){
+  return _amps_gpubuf_realloc(&amps_gpu_recvbuf, id, pos, size);
+}
+
+/**
+ * @brief Get send buffer associated with id to be passed for MPI
+ *
+ * @param id The id of the send buffer [IN]
+ * @return A pointer to the send staging buffer
+ */
+static char* _amps_gpu_sendbuf(int id){
+  if(id >= amps_gpu_sendbuf.num_bufs){
+    printf("ERROR at %s:%d: Not enough space is allocated for sendbufs\n", __FILE__, __LINE__);
+    exit(1);
+  }
+  if(ENFORCE_HOST_STAGING)
+    return amps_gpu_sendbuf.buf_host[id];
+  else
+    return amps_gpu_sendbuf.buf[id];
+}
+
+/**
+ * @brief Allocate/reallocate send GPU/CPU staging buffers
+ *
+ * @param id The id associated with the buffer (typically amps invoice index) [IN]
+ * @param pos The position in bytes up to which the buffer is already filled [IN]
+ * @param size The total required buffer size subtracted by pos [IN]
+ * @return A pointer to the allocated/reallocated GPU send buffer
+ */
+static char* _amps_gpu_sendbuf_realloc(int id, int pos, int size){
+  return _amps_gpubuf_realloc(&amps_gpu_sendbuf, id, pos, size);
+}
+
+/**
+* @brief The main amps GPU packing driver function
+*
+* Operates in 3 different modes:
+*
+* Mode 1: AMPS_GETRBUF / AMPS_GETSBUF
+* -determines the size of the invoice message in bytes
+* -allocates/reallocates the recv/send staging buffer if necessary
+* -*buffer_out is set to point to the recv/send staging buffer
+* -*size_out is set to the size of the invoice message in bytes
+*
+* Mode 2: AMPS_PACK (requires calling amps_gpu_sync_streams() before the packed data is accessed)
+* -determines the size of the invoice message in bytes
+* -allocates/reallocates the send staging buffer if necessary
+* -packs data to the send staging buffer parallel using a GPU
+* -*buffer_out is set to point to the send staging buffer
+* -*size_out is set to the size of invoice message in bytes (size of the packed data)
+*
+* Mode 3: AMPS_UNPACK (requires calling amps_gpu_sync_streams() before the unpacked data is accessed)
+* -determines the size of the invoice message in bytes
+* -allocates/reallocates the recv staging buffer if necessary
+* -unpacks data from the recv staging buffer parallel using a GPU
+* -*buffer_out is set to point to the recv staging buffer
+* -*size_out is set to the size of invoice message in bytes (size of the unpacked data)
+*
+* @param action either AMPS_GETRBUF, AMPS_GETSBUF, AMPS_PACK or AMPS_UNPACK [IN]
+* @param inv amps invoice [IN]
+* @param inv_num amps invoice order number [IN]
+* @param buffer_out pointer to the pointer of the staging buffer [OUT]
+* @param size_out pointer to the invoice message size in bytes [OUT]
+* @return error code (line number), 0 if succesful
+*/
+int amps_gpupacking(int action, amps_Invoice inv, int inv_num, char **buffer_out, int *size_out){
+  
+  char *buffer;
+  int pos = 0;
+
+  amps_kokkos_initialization();
+
+  amps_InvoiceEntry *ptr = inv->list;
+  while (ptr != NULL)
+  {
+    if (ptr->len_type != AMPS_INVOICE_POINTER){
+      printf("ERROR at %s:%d: ptr->len_type must be a pointer\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+  
+    if (ptr->stride_type != AMPS_INVOICE_POINTER){
+      printf("ERROR at %s:%d: ptr->stride_type must be a pointer\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+  
+    if (ptr->data_type == AMPS_INVOICE_POINTER){
+      printf("ERROR at %s:%d: Overlayed invoices not supported\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+  
+    /* Make sure amps invoice vector type is supported */
+    if(ptr->type - AMPS_INVOICE_LAST_CTYPE != AMPS_INVOICE_DOUBLE_CTYPE){
+      printf("ERROR at %s:%d: Only \"AMPS_INVOICE_DOUBLE_CTYPE\" is supported\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+  
+    /* Preparations for the kernel launch */
+    int len_x = 1;
+    int len_y = 1;
+    int len_z = 1;  
+    int stride_x = 0;
+    int stride_y = 0;
+    int stride_z = 0;
+  
+    int dim = (ptr->dim_type == AMPS_INVOICE_POINTER) ? *(ptr->ptr_dim) : ptr->dim;
+    switch (dim)
+    {
+      case 3:
+        len_z = ptr->ptr_len[2];
+        stride_z = ptr->ptr_stride[2];
+      case 2:
+        len_y = ptr->ptr_len[1];
+        stride_y = ptr->ptr_stride[1];
+      case 1:
+        len_x = ptr->ptr_len[0];
+        stride_x = ptr->ptr_stride[0];
+        break;
+      default:
+        printf("ERROR at %s:%d: Only dimensions 1 - 3 are supported\n", __FILE__, __LINE__);
+        return __LINE__;
+    }  
+    
+    /* Get data location and its properties */
+    char *data = (char*)ptr->data;
+
+#ifdef PARFLOW_HAVE_CUDA
+    struct cudaPointerAttributes attributes;
+    cudaPointerGetAttributes(&attributes, (void *)data);  
+  
+    /* Check that the data location (not MPI staging buffer) is accessible by the GPU */
+    if(cudaGetLastError() != cudaSuccess || attributes.type < 2){
+      printf("ERROR at %s:%d: The data location (not MPI staging buffer) is not accessible by the GPU(s)\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+#endif
+
+    int size = len_x * len_y * len_z * sizeof(double);
+    if((action == AMPS_GETSBUF) || (action == AMPS_PACK)){
+      buffer = _amps_gpu_sendbuf_realloc(inv_num, pos, size);
+    }
+    else if((action == AMPS_GETRBUF) || (action == AMPS_UNPACK)){
+      buffer = _amps_gpu_recvbuf_realloc(inv_num, pos, size);
+    }
+    else{
+      printf("ERROR at %s:%d: Unknown action argument (val = %d)\n", __FILE__, __LINE__, action);
+      return __LINE__;
+    }
+
+#ifdef PARFLOW_HAVE_CUDA
+    /* Get buffer location and its properties */
+    cudaPointerGetAttributes(&attributes, (void *)buffer);  
+
+    /* Check that the staging buffer is accessible by the GPU */
+    if(cudaGetLastError() != cudaSuccess || attributes.type < 2){
+      printf("ERROR at %s:%d: The MPI staging buffer location is not accessible by the GPU(s)\n", __FILE__, __LINE__);
+      return __LINE__;
+    }
+#endif
+
+    /* Run packing or unpacking kernel */
+    using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;
+    MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{len_x, len_y, len_z}});
+
+    if(action == AMPS_PACK){
+      Kokkos::parallel_for(mdpolicy_3d, KOKKOS_LAMBDA(int i, int j, int k)
+      {
+        double *ptr_buf = (double*)buffer;
+        double *ptr_data = (double*)data;
+        *(ptr_buf + k * len_y * len_x + j * len_x + i) = 
+          *(ptr_data + k * (stride_z + (len_y - 1) * stride_y + len_y * (len_x - 1) * stride_x) + 
+            j * (stride_y + (len_x - 1) * stride_x) + i * stride_x);                                          
+      });
+      if(ENFORCE_HOST_STAGING){
+        /* Copy device buffer to host after packing */
+        kokkosMemCpyDeviceToHost(amps_gpu_sendbuf.buf_host[inv_num] + pos,
+                                    amps_gpu_sendbuf.buf[inv_num] + pos, size);
+      }
+      inv->flags |= AMPS_PACKED;
+    }
+    else if(action == AMPS_UNPACK){
+      if(ENFORCE_HOST_STAGING){
+        /* Copy host buffer to device before unpacking */
+        kokkosMemCpyHostToDevice(amps_gpu_recvbuf.buf[inv_num] + pos,
+                                    amps_gpu_recvbuf.buf_host[inv_num] + pos, size);
+      }
+      Kokkos::parallel_for(mdpolicy_3d, KOKKOS_LAMBDA(int i, int j, int k)
+      {
+        double *ptr_buf = (double*)buffer;
+        double *ptr_data = (double*)data;
+        *(ptr_data + k * (stride_z + (len_y - 1) * stride_y + len_y * (len_x - 1) * stride_x) + 
+          j * (stride_y + (len_x - 1) * stride_x) + i * stride_x) = 
+            *(ptr_buf + k * len_y * len_x + j * len_x + i);
+      });
+      inv->flags &= ~AMPS_PACKED;
+    }
+
+    pos += size;  
+    ptr = ptr->next;
+  }
+
+  /* Check that the size is calculated right */
+  // if(pos != amps_sizeof_invoice(amps_CommWorld, inv)){
+    // printf("ERROR at %s:%d: The size does not match the invoice size\n", __FILE__, __LINE__);
+    // return __LINE__;
+  // }
+
+  /* Set the out values here if everything went fine */
+  if((action == AMPS_GETSBUF) || (action == AMPS_PACK)){
+    *buffer_out = _amps_gpu_sendbuf(inv_num);
+  }
+  else if((action == AMPS_GETRBUF) || (action == AMPS_UNPACK)){
+    *buffer_out = _amps_gpu_recvbuf(inv_num);
+  }
+  else{
+    printf("ERROR at %s:%d: Unknown action argument (val = %d)\n", __FILE__, __LINE__, action);
+    return __LINE__;
+  }
+
+  *size_out = pos;
+
+  return 0;
+} 
+#endif // DISABLE_GPU_PACKING
+}   
diff --git a/pfsimulator/amps/oas3/amps_gpupacking.cu b/pfsimulator/amps/oas3/amps_gpupacking.cu
index 4ceb14100..5ff8706a8 100644
--- a/pfsimulator/amps/oas3/amps_gpupacking.cu
+++ b/pfsimulator/amps/oas3/amps_gpupacking.cu
@@ -106,9 +106,7 @@ typedef struct _amps_GpuStreams {
 
 #ifdef DISABLE_GPU_PACKING
 /* Dummy definitions if no GPU packing */
-void amps_gpu_free_bufs(){}
-
-void amps_gpu_destroy_streams(){}
+void amps_gpu_finalize(){}
 
 void amps_gpu_sync_streams(int id){
   (void)id;
@@ -286,6 +284,14 @@ void amps_gpu_free_bufs(){
   _amps_gpubuf_free(&amps_gpu_sendbuf);
 }
 
+/**
+ * @brief Finalize GPU resource usage
+ */
+ void amps_gpu_finalize(){
+  amps_gpu_free_bufs();
+  amps_gpu_destroy_streams();
+}
+
 /**
  * @brief Allocate/reallocate page-locked GPU/CPU staging buffers.
  *
diff --git a/pfsimulator/amps/oas3/amps_proto.h b/pfsimulator/amps/oas3/amps_proto.h
index ef288249a..776fde962 100644
--- a/pfsimulator/amps/oas3/amps_proto.h
+++ b/pfsimulator/amps/oas3/amps_proto.h
@@ -34,8 +34,7 @@ void amps_FindPowers(int N, int *log, int *Nnext, int *Nprev);
 amps_File amps_Fopen(char *filename, char *type);
 
 /* amps_gpupacking.cu */
-void amps_gpu_free_bufs();
-void amps_gpu_destroy_streams();
+void amps_gpu_finalize();
 void amps_gpu_sync_streams(int id);
 int amps_gpupacking(int action, amps_Invoice inv, int inv_num, char **buffer_out, int *size_out);
 
diff --git a/pfsimulator/amps/smpi/amps_pack.c b/pfsimulator/amps/smpi/amps_pack.c
index 953c6fcac..6b76f323f 100644
--- a/pfsimulator/amps/smpi/amps_pack.c
+++ b/pfsimulator/amps/smpi/amps_pack.c
@@ -345,7 +345,7 @@ amps_Invoice inv;
         {
           if (i == dim - 1)
           {
-            MPI_Type_hvector(ptr->ptr_len[i], 1,
+            MPI_Type_create_hvector(ptr->ptr_len[i], 1,
                              base_size +
                              (ptr->ptr_stride[i] - 1) * element_size,
                              *base_type, &mpi_types[element]);
@@ -356,7 +356,7 @@ amps_Invoice inv;
           }
           else
           {
-            MPI_Type_hvector(ptr->ptr_len[i], 1,
+            MPI_Type_create_hvector(ptr->ptr_len[i], 1,
                              base_size +
                              (ptr->ptr_stride[i] - 1) * element_size,
                              *base_type, new_type);
@@ -380,14 +380,14 @@ amps_Invoice inv;
         break;
     }
 
-    MPI_Address(data, &mpi_displacements[element]);
+    MPI_Get_address(data, &mpi_displacements[element]);
 
     mpi_block_len[element] = 1;
     element++;
     ptr = ptr->next;
   }
 
-  MPI_Type_struct(inv->num,
+  MPI_Type_create_struct(inv->num,
                   mpi_block_len,
                   mpi_displacements,
                   mpi_types,
diff --git a/pfsimulator/kinsol/CMakeLists.txt b/pfsimulator/kinsol/CMakeLists.txt
index 0cf341d60..6602c76cc 100644
--- a/pfsimulator/kinsol/CMakeLists.txt
+++ b/pfsimulator/kinsol/CMakeLists.txt
@@ -5,13 +5,14 @@ add_library(pfkinsol ${SRC_FILES})
 
 if(${PARFLOW_HAVE_CUDA})
   target_include_directories(pfkinsol PUBLIC ${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
-  if( ${PARFLOW_HAVE_RMM} )
-    target_link_libraries(pfkinsol rmm)
-    target_include_directories(pfkinsol PRIVATE ${RMM_INCLUDE})
-    target_compile_definitions(pfkinsol PRIVATE PARFLOW_HAVE_RMM)
-  endif( ${PARFLOW_HAVE_RMM} )
 endif(${PARFLOW_HAVE_CUDA})
 
+if( ${PARFLOW_HAVE_RMM} )
+  target_link_libraries(pfkinsol rmm)
+  target_include_directories(pfkinsol PRIVATE ${RMM_INCLUDE})
+  target_compile_definitions(pfkinsol PRIVATE PARFLOW_HAVE_RMM)
+endif( ${PARFLOW_HAVE_RMM} )
+
 if (${PARFLOW_HAVE_MPI})
   target_include_directories (pfkinsol PUBLIC "${MPI_C_INCLUDE_PATH}")
 endif (${PARFLOW_HAVE_MPI})
diff --git a/pfsimulator/parflow_exe/main.c b/pfsimulator/parflow_exe/main.c
index e7a837b47..867b8eb1a 100644
--- a/pfsimulator/parflow_exe/main.c
+++ b/pfsimulator/parflow_exe/main.c
@@ -61,8 +61,8 @@ using namespace SAMRAI;
 #include <cegdb.h>
 #endif
 
-#if PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
-#include "pf_cudamain.h"
+#if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS)
+#include "pf_devices.h"
 #endif
 
 #ifdef PARFLOW_HAVE_ETRACE
@@ -118,8 +118,6 @@ int main(int argc, char *argv [])
     cegdb(&argc, &argv, amps_Rank(MPI_CommWorld));
 #endif
 
-#if PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
-
 #ifndef NDEBUG
     /*-----------------------------------------------------------------------
     * Wait for debugger if PARFLOW_DEBUG_RANK environment variable is set
@@ -136,6 +134,13 @@ int main(int argc, char *argv [])
     }
 #endif // !NDEBUG
 
+    /*-----------------------------------------------------------------------
+     * Initialize acceleration architectures
+     *-----------------------------------------------------------------------*/
+#if defined(PARFLOW_HAVE_KOKKOS)
+    kokkosInit();
+#elif defined(PARFLOW_HAVE_CUDA)
+
     /*-----------------------------------------------------------------------
     * Check CUDA compute capability, set device, and initialize RMM allocator
     *-----------------------------------------------------------------------*/
@@ -166,7 +171,13 @@ int main(int argc, char *argv [])
         amps_Printf("\nThe minimum required GPU compute capability is 6.0.\n");
         exit(1);
       }
+    }
+#endif // PARFLOW_HAVE_KOKKOS
+
 
+  /*-----------------------------------------------------------------------
+  * Initialize RMM pool allocator
+  *-----------------------------------------------------------------------*/
 #ifdef PARFLOW_HAVE_RMM
       // RMM
       rmmOptions_t rmmOptions;
@@ -175,8 +186,6 @@ int main(int argc, char *argv [])
       rmmOptions.enable_logging = false;
       RMM_ERR(rmmInitialize(&rmmOptions));
 #endif // PARFLOW_HAVE_RMM
-    }
-#endif // PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
 
     wall_clock_time = amps_Clock();
 
@@ -485,11 +494,18 @@ int main(int argc, char *argv [])
   tbox::SAMRAI_MPI::finalize();
 #endif
 
+  /*-----------------------------------------------------------------------
+  * Shutdown Kokkos
+  *-----------------------------------------------------------------------*/
+#ifdef PARFLOW_HAVE_KOKKOS
+  kokkosFinalize();
+#endif
+
   /*-----------------------------------------------------------------------
   * Shutdown RMM pool allocator
   *-----------------------------------------------------------------------*/
-#if (PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA) && defined(PARFLOW_HAVE_RMM)
-    RMM_ERR(rmmFinalize());
+#ifdef PARFLOW_HAVE_RMM
+  RMM_ERR(rmmFinalize());
 #endif
 
   return 0;
diff --git a/pfsimulator/parflow_lib/CMakeLists.txt b/pfsimulator/parflow_lib/CMakeLists.txt
index 76ab84b72..0cebbeb99 100644
--- a/pfsimulator/parflow_lib/CMakeLists.txt
+++ b/pfsimulator/parflow_lib/CMakeLists.txt
@@ -1,6 +1,5 @@
 set (SRC_FILES_CONST advect.F
   advection_godunov.c
-  background.c
   bc_lb.c
   bc_pressure.c
   bc_pressure_package.c
@@ -124,11 +123,10 @@ set (SRC_FILES_CONST advect.F
   write_parflow_silo_pmpio.c
 )
 
-if( (${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_OMP}) )
+if( (${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_KOKKOS}) OR (${PARFLOW_HAVE_OMP}) )
 
   set (SRC_FILES_ARCH_TYPE_0
     clustering.cpp
-    globals.cpp
     grgeometry.cpp
   )
   set (SRC_FILES_ARCH_TYPE_1
@@ -154,7 +152,9 @@ if( (${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_OMP}) )
     vector_utilities.cpp
   )
   set (SRC_FILES_ARCH_TYPE_2
+    background.cpp
     calc_elevations.cpp
+    globals.cpp
     pf_module.cpp
     problem_bc.cpp
     problem_bc_pressure.cpp
@@ -162,6 +162,10 @@ if( (${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_OMP}) )
     region.cpp
   )
 
+if(${PARFLOW_HAVE_KOKKOS})
+  list(APPEND SRC_FILES_ARCH_TYPE_0 pf_kokkos.cpp)
+endif(${PARFLOW_HAVE_KOKKOS})
+
   #  PF_COMP_UNIT_TYPE determines the behavior of the accelerated compilation unit
 	#  ------------------------------------------------------------
   #  CUDA
@@ -173,8 +177,8 @@ if( (${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_OMP}) )
 	#  ------------------------------------------------------------
   #  OpenMP
 	#  ------------------------------------------------------------
-  #  1:     CXX compiler, Unified Memory allocation, Parallel loops on CPU
-	#  2:     CXX compiler, Unified Memory allocation, Sequential loops on CPU
+  #  1:     CXX compiler, Standard heap allocation, Parallel loops on CPU
+	#  2:     CXX compiler, Standard heap allocation, Sequential loops on CPU
 	#  Other: CXX compiler, Standard heap allocation, Sequential loops on CPU
 
   set_source_files_properties(${SRC_FILES_ARCH_TYPE_0} PROPERTIES COMPILE_DEFINITIONS "PF_COMP_UNIT_TYPE=0")
@@ -183,31 +187,37 @@ if( (${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_OMP}) )
 
   set (SRC_FILES_ARCH ${SRC_FILES_ARCH_TYPE_0} ${SRC_FILES_ARCH_TYPE_1} ${SRC_FILES_ARCH_TYPE_2})
 
-  if ( ${PARFLOW_HAVE_CUDA} )
+  if((${PARFLOW_HAVE_CUDA}) AND (NOT (${PARFLOW_HAVE_KOKKOS})))
     set_source_files_properties(${SRC_FILES_ARCH} PROPERTIES LANGUAGE CUDA)
-  else ( ${PARFLOW_HAVE_CUDA} )
+  else((${PARFLOW_HAVE_CUDA}) AND (NOT (${PARFLOW_HAVE_KOKKOS})))
     set_source_files_properties(${SRC_FILES_ARCH} PROPERTIES LANGUAGE CXX)
-  endif ( ${PARFLOW_HAVE_CUDA} )
+  endif((${PARFLOW_HAVE_CUDA}) AND (NOT (${PARFLOW_HAVE_KOKKOS})))
+
 
   add_library(pfsimulator ${SRC_FILES_ARCH} ${SRC_FILES_CONST})
 
   if ( ${PARFLOW_HAVE_CUDA} )
-    set_target_properties(pfsimulator PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
     target_link_libraries(pfsimulator nvToolsExt)
     target_include_directories(pfsimulator PUBLIC "../third_party/cub")
     target_include_directories(pfsimulator PUBLIC ${CMAKE_CUDA_TOOLKIT_INCLUDE_DIRECTORIES})
   endif ( ${PARFLOW_HAVE_CUDA} )
 
+  if( ${PARFLOW_HAVE_KOKKOS} )
+    target_link_libraries(pfsimulator Kokkos::kokkos)
+    target_include_directories(pfsimulator PUBLIC ${KOKKOS_INCLUDE})
+  endif( ${PARFLOW_HAVE_KOKKOS} )
+
   if( ${PARFLOW_HAVE_RMM} )
     target_link_libraries(pfsimulator rmm)
     target_include_directories(pfsimulator PRIVATE ${RMM_INCLUDE})
     target_compile_definitions(pfsimulator PRIVATE PARFLOW_HAVE_RMM)
   endif( ${PARFLOW_HAVE_RMM} )
 
-else( (${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_OMP}) )
+else( (${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_KOKKOS}) OR (${PARFLOW_HAVE_OMP}) )
 
   set (SRC_FILES_ARCH
     axpy.c
+    background.c
     calc_elevations.c
     clustering.c
     copy.c
@@ -239,7 +249,7 @@ else( (${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_OMP}) )
   )
   add_library(pfsimulator ${SRC_FILES_ARCH} ${SRC_FILES_CONST})
 
-endif( (${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_OMP}) )
+endif( (${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_KOKKOS}) OR (${PARFLOW_HAVE_OMP}) )
 
 target_link_libraries(pfsimulator pfkinsol amps cjson ${PARFLOW_ETRACE_LIBRARY})
 target_include_directories(pfsimulator PUBLIC "../third_party/cjson")
diff --git a/pfsimulator/parflow_lib/backend_mapping.h b/pfsimulator/parflow_lib/backend_mapping.h
index b2ea4da12..85acb61b7 100644
--- a/pfsimulator/parflow_lib/backend_mapping.h
+++ b/pfsimulator/parflow_lib/backend_mapping.h
@@ -23,7 +23,7 @@
 
 /* PF_COMP_UNIT_TYPE determines the behavior of the NVCC compilation unit and/or OpenMP loops:
 ------------------------------------------------------------
-   CUDA: PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
+   CUDA
 ------------------------------------------------------------
    1:     NVCC compiler, Unified Memory allocation, Parallel loops on GPUs
    2:     NVCC compiler, Unified Memory allocation, Sequential loops on host
@@ -31,19 +31,31 @@
 
 
 ------------------------------------------------------------
-   OpenMP: PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_OMP
+   OpenMP
 ------------------------------------------------------------
-   1:     CXX compiler, Unified Memory allocation, Parallel loops on CPU
-   2:     CXX compiler, Unified Memory allocation, Sequential loops on CPU
+   1:     CXX compiler, Standard heap allocation, Parallel loops on CPU
+   2:     CXX compiler, Standard heap allocation, Sequential loops on CPU
    Other: CXX compiler, Standard heap allocation, Sequential loops on CPU
 */
 
 /* Include headers depending on the accelerator backend */
-#if (PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA)
+#ifdef PARFLOW_HAVE_KOKKOS
+
+  #define ACC_ID _kokkos
+
+  #include "pf_devices.h"
+
+  #if PF_COMP_UNIT_TYPE == 1
+    #include "pf_kokkosloops.h"
+  #elif PF_COMP_UNIT_TYPE == 2
+    #include "pf_kokkosmalloc.h"
+  #endif
+
+#elif defined(PARFLOW_HAVE_CUDA)
 
   #define ACC_ID _cuda
 
-  #include "pf_cudamain.h"
+  #include "pf_devices.h"
 
   #if PF_COMP_UNIT_TYPE == 1
     #include "pf_cudaloops.h"
@@ -51,7 +63,7 @@
     #include "pf_cudamalloc.h"
   #endif
 
-#elif (PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_OMP)
+#elif defined(PARFLOW_HAVE_OMP)
 
   #define ACC_ID _omp
 
@@ -74,67 +86,67 @@
 
 // Memory management
 
-#if defined(talloc_amps_cuda) || defined(talloc_amps_omp)
+#if defined(talloc_amps_cuda) || defined(talloc_amps_kokkos) || defined(talloc_amps_omp)
   #define talloc_amps CHOOSE_BACKEND(DEFER(talloc_amps), ACC_ID)
 #else
   #define talloc_amps talloc_amps_default
 #endif
 
-#if defined(ctalloc_amps_cuda) || defined(ctalloc_amps_omp)
+#if defined(ctalloc_amps_cuda) || defined(ctalloc_amps_kokkos) || defined(ctalloc_amps_omp)
   #define ctalloc_amps CHOOSE_BACKEND(DEFER(ctalloc_amps), ACC_ID)
 #else
   #define ctalloc_amps ctalloc_amps_default
 #endif
 
-#if defined(tfree_amps_cuda) || defined(tfree_amps_omp)
+#if defined(tfree_amps_cuda) || defined(tfree_amps_kokkos) || defined(tfree_amps_omp)
   #define tfree_amps CHOOSE_BACKEND(DEFER(tfree_amps), ACC_ID)
 #else
   #define tfree_amps tfree_amps_default
 #endif
 
-#if defined(talloc_cuda) || defined(talloc_omp)
+#if defined(talloc_cuda) || defined(talloc_kokkos) || defined(talloc_omp)
   #define talloc CHOOSE_BACKEND(DEFER(talloc), ACC_ID)
 #else
   #define talloc talloc_default
 #endif
 
-#if defined(ctalloc_cuda) || defined(ctalloc_omp)
+#if defined(ctalloc_cuda) || defined(ctalloc_kokkos) || defined(ctalloc_omp)
   #define ctalloc CHOOSE_BACKEND(DEFER(ctalloc), ACC_ID)
 #else
   #define ctalloc ctalloc_default
 #endif
 
-#if defined(tfree_cuda) || defined(tfree_omp)
+#if defined(tfree_cuda) || defined(tfree_kokkos) || defined(tfree_omp)
   #define tfree CHOOSE_BACKEND(DEFER(tfree), ACC_ID)
 #else
   #define tfree tfree_default
 #endif
 
-#if defined(tmemcpy_cuda) || defined(tmemcpy_omp)
+#if defined(tmemcpy_cuda) || defined(tmemcpy_kokkos) || defined(tmemcpy_omp)
   #define tmemcpy CHOOSE_BACKEND(DEFER(tmemcpy), ACC_ID)
 #else
   #define tmemcpy tmemcpy_default
 #endif
 
-#if defined(MemPrefetchDeviceToHost_cuda) || defined(MemPrefetchDeviceToHost_omp)
+#if defined(MemPrefetchDeviceToHost_cuda) || defined(MemPrefetchDeviceToHost_kokkos) || defined(MemPrefetchDeviceToHost_omp)
   #define MemPrefetchDeviceToHost CHOOSE_BACKEND(DEFER(MemPrefetchDeviceToHost), ACC_ID)
 #else
   #define MemPrefetchDeviceToHost MemPrefetchDeviceToHost_default
 #endif
 
-#if defined(MemPrefetchHostToDevice_cuda) || defined(MemPrefetchHostToDevice_omp)
+#if defined(MemPrefetchHostToDevice_cuda) || defined(MemPrefetchHostToDevice_kokkos) || defined(MemPrefetchHostToDevice_omp)
   #define MemPrefetchHostToDevice CHOOSE_BACKEND(DEFER(MemPrefetchHostToDevice), ACC_ID)
 #else
   #define MemPrefetchHostToDevice MemPrefetchHostToDevice_default
 #endif
 
-#if defined(SKIP_PARALLEL_SYNC_cuda) || defined(SKIP_PARALLEL_SYNC_omp)
+#if defined(SKIP_PARALLEL_SYNC_cuda) || defined(SKIP_PARALLEL_SYNC_kokkos) || defined(SKIP_PARALLEL_SYNC_omp)
   #define SKIP_PARALLEL_SYNC CHOOSE_BACKEND(DEFER(SKIP_PARALLEL_SYNC), ACC_ID)
 #else
   #define SKIP_PARALLEL_SYNC SKIP_PARALLEL_SYNC_default
 #endif
 
-#if defined(PARALLEL_SYNC_cuda) || defined(PARALLEL_SYNC_omp)
+#if defined(PARALLEL_SYNC_cuda) || defined(PARALLEL_SYNC_kokkos) || defined(PARALLEL_SYNC_omp)
   #define PARALLEL_SYNC CHOOSE_BACKEND(DEFER(PARALLEL_SYNC), ACC_ID)
 #else
   #define PARALLEL_SYNC PARALLEL_SYNC_default
@@ -142,117 +154,117 @@
 
 // General 
 
-#if defined(PlusEquals_cuda) || defined(PlusEquals_omp)
+#if defined(PlusEquals_cuda) || defined(PlusEquals_kokkos) || defined(PlusEquals_omp)
   #define PlusEquals CHOOSE_BACKEND(DEFER(PlusEquals), ACC_ID)
 #else
   #define PlusEquals PlusEquals_default
 #endif
 
-#if defined(ReduceMax_cuda) || defined(ReduceMax_omp)
+#if defined(ReduceMax_cuda) || defined(ReduceMax_kokkos) || defined(ReduceMax_omp)
   #define ReduceMax CHOOSE_BACKEND(DEFER(ReduceMax), ACC_ID)
 #else
   #define ReduceMax ReduceMax_default
 #endif
 
-#if defined(ReduceMin_cuda) || defined(ReduceMin_omp)
+#if defined(ReduceMin_cuda) || defined(ReduceMin_kokkos) || defined(ReduceMin_omp)
   #define ReduceMin CHOOSE_BACKEND(DEFER(ReduceMin), ACC_ID)
 #else
   #define ReduceMin ReduceMin_default
 #endif
 
-#if defined(ReduceSum_cuda) || defined(ReduceSum_omp)
+#if defined(ReduceSum_cuda) || defined(ReduceSum_kokkos) || defined(ReduceSum_omp)
   #define ReduceSum CHOOSE_BACKEND(DEFER(ReduceSum), ACC_ID)
 #else
   #define ReduceSum ReduceSum_default
 #endif
 
-#if defined(PUSH_NVTX_cuda) || defined(PUSH_NVTX_omp)
-  #define PUSH_NVTX CHOOSE_BACKEND(DEFER(PUSH_NVTX), ACC_ID)
+#if defined(PUSH_NVTX_cuda)
+  #define PUSH_NVTX PUSH_NVTX_cuda
 #else
   #define PUSH_NVTX PUSH_NVTX_default
 #endif
 
-#if defined(POP_NVTX_cuda) || defined(POP_NVTX_omp)
-  #define POP_NVTX CHOOSE_BACKEND(DEFER(POP_NVTX), ACC_ID)
+#if defined(POP_NVTX_cuda)
+  #define POP_NVTX POP_NVTX_cuda
 #else
   #define POP_NVTX POP_NVTX_default
 #endif
 
 // Loops
 
-#if defined(BoxLoopI0_cuda) || defined(BoxLoopI0_omp)
+#if defined(BoxLoopI0_cuda) || defined(BoxLoopI0_kokkos) || defined(BoxLoopI0_omp)
   #define BoxLoopI0 CHOOSE_BACKEND(DEFER(BoxLoopI0), ACC_ID)
 #else
   #define BoxLoopI0 BoxLoopI0_default
 #endif
 
-#if defined(BoxLoopI1_cuda) || defined(BoxLoopI1_omp)
+#if defined(BoxLoopI1_cuda) || defined(BoxLoopI1_kokkos) || defined(BoxLoopI1_omp)
   #define BoxLoopI1 CHOOSE_BACKEND(DEFER(BoxLoopI1), ACC_ID)
 #else
   #define BoxLoopI1 BoxLoopI1_default
 #endif
 
-#if defined(BoxLoopI2_cuda) || defined(BoxLoopI2_omp)
+#if defined(BoxLoopI2_cuda) || defined(BoxLoopI2_kokkos) || defined(BoxLoopI2_omp)
   #define BoxLoopI2 CHOOSE_BACKEND(DEFER(BoxLoopI2), ACC_ID)
 #else
   #define BoxLoopI2 BoxLoopI2_default
 #endif
 
-#if defined(BoxLoopI3_cuda) || defined(BoxLoopI3_omp)
+#if defined(BoxLoopI3_cuda) || defined(BoxLoopI3_kokkos) || defined(BoxLoopI3_omp)
   #define BoxLoopI3 CHOOSE_BACKEND(DEFER(BoxLoopI3), ACC_ID)
 #else
   #define BoxLoopI3 BoxLoopI3_default
 #endif
 
-#if defined(BoxLoopReduceI1_cuda) || defined(BoxLoopReduceI1_omp)
+#if defined(BoxLoopReduceI1_cuda) || defined(BoxLoopReduceI1_kokkos) || defined(BoxLoopReduceI1_omp)
   #define BoxLoopReduceI1 CHOOSE_BACKEND(DEFER(BoxLoopReduceI1), ACC_ID)
 #else
   #define BoxLoopReduceI1 BoxLoopReduceI1_default
 #endif
 
-#if defined(BoxLoopReduceI2_cuda) || defined(BoxLoopReduceI2_omp)
+#if defined(BoxLoopReduceI2_cuda) || defined(BoxLoopReduceI2_kokkos) || defined(BoxLoopReduceI2_omp)
   #define BoxLoopReduceI2 CHOOSE_BACKEND(DEFER(BoxLoopReduceI2), ACC_ID)
 #else
   #define BoxLoopReduceI2 BoxLoopReduceI2_default
 #endif
 
-#if defined(BoxLoopReduceI3_cuda) || defined(BoxLoopReduceI3_omp)
+#if defined(BoxLoopReduceI3_cuda) || defined(BoxLoopReduceI3_kokkos) || defined(BoxLoopReduceI3_omp)
   #define BoxLoopReduceI3 CHOOSE_BACKEND(DEFER(BoxLoopReduceI3), ACC_ID)
 #else
   #define BoxLoopReduceI3 BoxLoopReduceI3_default
 #endif
 
-#if defined(GrGeomInLoopBoxes_cuda) || defined(GrGeomInLoopBoxes_omp)
+#if defined(GrGeomInLoopBoxes_cuda) || defined(GrGeomInLoopBoxes_kokkos) || defined(GrGeomInLoopBoxes_omp)
   #define GrGeomInLoopBoxes CHOOSE_BACKEND(DEFER(GrGeomInLoopBoxes), ACC_ID)
 #else
   #define GrGeomInLoopBoxes GrGeomInLoopBoxes_default
 #endif
 
-#if defined(GrGeomSurfLoopBoxes_cuda) || defined(GrGeomSurfLoopBoxes_omp)
+#if defined(GrGeomSurfLoopBoxes_cuda) || defined(GrGeomSurfLoopBoxes_kokkos) || defined(GrGeomSurfLoopBoxes_omp)
   #define GrGeomSurfLoopBoxes CHOOSE_BACKEND(DEFER(GrGeomSurfLoopBoxes), ACC_ID)
 #else
   #define GrGeomSurfLoopBoxes GrGeomSurfLoopBoxes_default
 #endif
 
-#if defined(GrGeomPatchLoopBoxes_cuda) || defined(GrGeomPatchLoopBoxes_omp)
+#if defined(GrGeomPatchLoopBoxes_cuda) || defined(GrGeomPatchLoopBoxes_kokkos) || defined(GrGeomPatchLoopBoxes_omp)
   #define GrGeomPatchLoopBoxes CHOOSE_BACKEND(DEFER(GrGeomPatchLoopBoxes), ACC_ID)
 #else
   #define GrGeomPatchLoopBoxes GrGeomPatchLoopBoxes_default
 #endif
 
-#if defined(GrGeomPatchLoopBoxesNoFdir_cuda) || defined(GrGeomPatchLoopBoxesNoFdir_omp)
+#if defined(GrGeomPatchLoopBoxesNoFdir_cuda) || defined(GrGeomPatchLoopBoxesNoFdir_kokkos) || defined(GrGeomPatchLoopBoxesNoFdir_omp)
   #define GrGeomPatchLoopBoxesNoFdir CHOOSE_BACKEND(DEFER(GrGeomPatchLoopBoxesNoFdir), ACC_ID)
 #else
   #define GrGeomPatchLoopBoxesNoFdir GrGeomPatchLoopBoxesNoFdir_default
 #endif
 
-#if defined(GrGeomOctreeExteriorNodeLoop_cuda) || defined(GrGeomOctreeExteriorNodeLoop_omp)
+#if defined(GrGeomOctreeExteriorNodeLoop_cuda) || defined(GrGeomOctreeExteriorNodeLoop_kokkos) || defined(GrGeomOctreeExteriorNodeLoop_omp)
   #define GrGeomOctreeExteriorNodeLoop CHOOSE_BACKEND(DEFER(GrGeomOctreeExteriorNodeLoop), ACC_ID)
 #else
   #define GrGeomOctreeExteriorNodeLoop GrGeomOctreeExteriorNodeLoop_default
 #endif
 
-#if defined(GrGeomOutLoop_cuda) || defined(GrGeomOutLoop_omp)
+#if defined(GrGeomOutLoop_cuda) || defined(GrGeomOutLoop_kokkos) || defined(GrGeomOutLoop_omp)
   #define GrGeomOutLoop CHOOSE_BACKEND(DEFER(GrGeomOutLoop), ACC_ID)
 #else
   #define GrGeomOutLoop GrGeomOutLoop_default
diff --git a/pfsimulator/parflow_lib/background.cpp b/pfsimulator/parflow_lib/background.cpp
new file mode 100644
index 000000000..95bd814c8
--- /dev/null
+++ b/pfsimulator/parflow_lib/background.cpp
@@ -0,0 +1,23 @@
+/* PF_COMP_UNIT_TYPE determines the behavior of the accelerated compilation unit
+	 ------------------------------------------------------------
+   CUDA
+	 ------------------------------------------------------------
+	 1:     NVCC compiler, Unified Memory allocation, Parallel loops on GPUs
+	 2:     NVCC compiler, Unified Memory allocation, Sequential loops on host
+	 Other: NVCC compiler, Standard heap allocation, Sequential loops on host
+
+	 ------------------------------------------------------------
+   OpenMP
+	 ------------------------------------------------------------
+     1:     CXX compiler, Unified Memory allocation, Parallel loops on CPU
+	 2:     CXX compiler, Unified Memory allocation, Sequential loops on CPU
+	 Other: CXX compiler, Standard heap allocation, Sequential loops on CPU
+*/
+
+// This definition is now controlled by CMake (./CMakeLists.txt)
+// #define PF_COMP_UNIT_TYPE 1
+
+/* extern "C" is required for the C source files when compiled with CPP compiler */
+extern "C"{
+  #include "background.c"
+}
diff --git a/pfsimulator/parflow_lib/clustering.c b/pfsimulator/parflow_lib/clustering.c
index 6d61acc66..32d8afe3a 100644
--- a/pfsimulator/parflow_lib/clustering.c
+++ b/pfsimulator/parflow_lib/clustering.c
@@ -51,15 +51,8 @@
 
 /**
  * Tag count.
- *
- * The __managed__ keyword makes this variable to automatically migrate
- * between host and device, but unfortunately requires global scope.
- * 
- * @TODO Find a better way to do this, some options:
- *  - Put this variable into heap
- *  - Use parallel reduction loops (currently not compatible with GrGeomLoops) 
  */
-__managed__ static int tag_count;
+static int tag_count;
 
 /**
  * Maximum number of ghost layers.
@@ -328,7 +321,6 @@ void ReduceTags(HistogramBox *histogram_box, Vector *vector, int dim, DoubleTags
           }
         });
       }
-      MemPrefetchDeviceToHost(&tag_count, sizeof(int), 0);
       HistogramBoxAddTags(histogram_box, dim, ic_sb, tag_count);
     }
   }
diff --git a/pfsimulator/parflow_lib/globals.c b/pfsimulator/parflow_lib/globals.c
index add5209ff..6972788c4 100644
--- a/pfsimulator/parflow_lib/globals.c
+++ b/pfsimulator/parflow_lib/globals.c
@@ -88,7 +88,7 @@ void   NewGlobals(char *run_name)
 
 void  FreeGlobals()
 {
-  free(globals_ptr);
+  tfree(globals_ptr);
 }
 
 /*--------------------------------------------------------------------------
@@ -117,39 +117,3 @@ void  LogGlobals()
     CloseLogFile(log_file);
   }
 }
-
-#if PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
-
-/**
- * @brief A struct for constant global data in read-only device memory.
- **/
-__constant__ Globals dev_globals;
-
-/**
- * @brief A struct member for constant global data in read-only device memory.
- **/
-__constant__ Background dev_background;
-
-/**
- * @brief A function to copy constant global data to read-only device memory.
- **/
-void CopyGlobalsToDevice()
-{
-  // Temporary pointers for getting the symbol addresses
-  Globals *tmp_globals_ptr;
-  Background *tmp_background_ptr;
-
-  // Get the symbol addresses from the __constant__ allocations
-  CUDA_ERR(cudaGetSymbolAddress((void**)&tmp_globals_ptr, dev_globals));
-  CUDA_ERR(cudaGetSymbolAddress((void**)&tmp_background_ptr, dev_background));
-
-  // Copy background data from host to device
-  CUDA_ERR(cudaMemcpy(tmp_background_ptr, globals->background, sizeof(Background), cudaMemcpyHostToDevice));
-
-  // Assign dev_background address to dev_globals.background
-  CUDA_ERR(cudaMemcpyToSymbol(dev_globals, &tmp_background_ptr, sizeof(Background*), offsetof(Globals, background), cudaMemcpyHostToDevice));
-
-  // Assign dev_globals address to dev_globals_ptr
-  CUDA_ERR(cudaMemcpyToSymbol(dev_globals_ptr, &tmp_globals_ptr, sizeof(Globals*), 0, cudaMemcpyHostToDevice));
-}
-#endif // PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
diff --git a/pfsimulator/parflow_lib/globals.h b/pfsimulator/parflow_lib/globals.h
index 1637a6076..26d8e7c2a 100644
--- a/pfsimulator/parflow_lib/globals.h
+++ b/pfsimulator/parflow_lib/globals.h
@@ -101,29 +101,7 @@ amps_ThreadLocalDcl(extern Globals  *, globals_ptr);
 amps_ThreadLocalDcl(extern IDB *, input_database);
 #endif
 
-/*--------------------------------------------------------------------------
- * Define __constant__ device pointer for globals struct (CUDA)  
- *--------------------------------------------------------------------------*/
-#if (PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA) && defined(__CUDACC__)
-#ifdef PARFLOW_GLOBALS
-__constant__  Globals *dev_globals_ptr;
-#else
-/* This extern requires CUDA separate compilation, otherwise nvcc compiler 
-   treats *dev_globals_ptr as static variable for each compilation unit   */
-extern __constant__  Globals *dev_globals_ptr;
-#endif // PARFLOW_GLOBALS
-#endif // PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA && __CUDACC__
-
-
-/*--------------------------------------------------------------------------
- * The globals accessor macro depends on compilation trajectory (host/device)
- *--------------------------------------------------------------------------*/
-#ifdef __CUDA_ARCH__
-  #define globals amps_ThreadLocal(dev_globals_ptr)
-#else
-  #define globals amps_ThreadLocal(globals_ptr)
-#endif
-
+#define globals amps_ThreadLocal(globals_ptr)
 
 /*--------------------------------------------------------------------------
  * Accessor macros: Globals
diff --git a/pfsimulator/parflow_lib/grgeometry.c b/pfsimulator/parflow_lib/grgeometry.c
index 7101b3563..59baca0ca 100644
--- a/pfsimulator/parflow_lib/grgeometry.c
+++ b/pfsimulator/parflow_lib/grgeometry.c
@@ -295,7 +295,7 @@ GrGeomSolid   *GrGeomNewSolid(
 
   new_grgeomsolid->interior_boxes = NULL;
 
-#if PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
+#if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS)
   GrGeomSolidCellFlagData(new_grgeomsolid) = NULL; 
   GrGeomSolidCellFlagDataSize(new_grgeomsolid) = 0;
   GrGeomSolidCellFlagInitialized(new_grgeomsolid) = 0;
@@ -361,9 +361,9 @@ void          GrGeomFreeSolid(
     tfree(solid->patch_boxes[f]);
   }
 
-#if PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
-  // Internal _tfree_cuda function is used because unified memory is not active in this comp unit
-  if(GrGeomSolidCellFlagData(solid)) _tfree_cuda(GrGeomSolidCellFlagData(solid));
+#if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS)
+  // Internal _tfree_device function is used because unified memory is not active in this comp unit
+  if(GrGeomSolidCellFlagData(solid)) _tfree_device(GrGeomSolidCellFlagData(solid));
 
   for (int f = 0; f < GrGeomOctreeNumFaces; f++)
   {
@@ -371,7 +371,7 @@ void          GrGeomFreeSolid(
     {
       int *ival = GrGeomSolidCellIval(solid, ipatch, f);
       if(ival)
-        _tfree_cuda(ival);
+        _tfree_device(ival);
     }
     tfree(solid->ival[f]);
   }
diff --git a/pfsimulator/parflow_lib/grgeometry.h b/pfsimulator/parflow_lib/grgeometry.h
index 4c3ebd13d..5a665ce79 100644
--- a/pfsimulator/parflow_lib/grgeometry.h
+++ b/pfsimulator/parflow_lib/grgeometry.h
@@ -59,7 +59,7 @@ typedef struct {
   GrGeomOctree **patches;
   int num_patches;
 
-#if PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
+#if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS)
   /* Cell flags for 3 GrGeomLoops (0: do not evaluate cell, 1 = evaluate cell)
    *  Bit 0: GrGeomInLoop
    *  Bit 1: GrGeomOutLoop
diff --git a/pfsimulator/parflow_lib/logging.c b/pfsimulator/parflow_lib/logging.c
index 32ceec0b2..7e0e2874a 100644
--- a/pfsimulator/parflow_lib/logging.c
+++ b/pfsimulator/parflow_lib/logging.c
@@ -141,11 +141,15 @@ void PrintVersionInfo(FILE *log_file)
   fprintf(log_file, "\tWith F77 flags  : %s\n", FFLAGS);
 #endif
 
-#if (PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA) && !defined(PARFLOW_HAVE_RMM)
+#if defined(PARFLOW_HAVE_KOKKOS) && !defined(PARFLOW_HAVE_RMM)
+  fprintf(log_file, "\tWith acc backend: KOKKOS\n");
+#elif defined(PARFLOW_HAVE_KOKKOS) && defined(PARFLOW_HAVE_RMM)
+  fprintf(log_file, "\tWith acc backend: KOKKOS+RMM\n");
+#elif defined(PARFLOW_HAVE_CUDA) && !defined(PARFLOW_HAVE_RMM)
   fprintf(log_file, "\tWith acc backend: CUDA\n");
-#elif (PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA) && defined(PARFLOW_HAVE_RMM)
+#elif defined(PARFLOW_HAVE_CUDA) && defined(PARFLOW_HAVE_RMM)
   fprintf(log_file, "\tWith acc backend: CUDA+RMM\n");
-#elif (PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_OMP)
+#elif defined(PARFLOW_HAVE_OMP)
   fprintf(log_file, "\tWith acc backend: OMP\n");
 #endif
 
diff --git a/pfsimulator/parflow_lib/overlandflow_eval.c b/pfsimulator/parflow_lib/overlandflow_eval.c
index 5046f34bc..b68e84f3e 100644
--- a/pfsimulator/parflow_lib/overlandflow_eval.c
+++ b/pfsimulator/parflow_lib/overlandflow_eval.c
@@ -38,7 +38,7 @@
 
 #include "parflow.h"
 
-#if PARFLOW_ACC_BACKEND != PARFLOW_BACKEND_CUDA
+#if !defined(PARFLOW_HAVE_CUDA) && !defined(PARFLOW_HAVE_KOKKOS)
 #include "llnlmath.h"
 //#include "llnltyps.h"
 #endif
@@ -105,7 +105,7 @@ void    OverlandFlowEval(
   {
     if (qx_v == NULL || qy_v == NULL)  /* do not return velocity fluxes */
     {
-      ForPatchCellsPerFace(ALL,
+      ForPatchCellsPerFace(BC_ALL,
                            BeforeAllCells(DoNothing),
                            LoopVars(i, j, k, ival, bc_struct, ipatch, sg),
                            Locals(int io, itop, ip, k1, ii, step;
@@ -187,7 +187,7 @@ void    OverlandFlowEval(
     }
     else   /* return velocity fluxes */
     {
-      ForPatchCellsPerFace(ALL,
+      ForPatchCellsPerFace(BC_ALL,
                            BeforeAllCells(DoNothing),
                            LoopVars(i, j, k, ival, bc_struct, ipatch, sg),
                            Locals(int io, itop, ip, k1, step, ii;
@@ -267,7 +267,7 @@ void    OverlandFlowEval(
   {
     if (qx_v == NULL || qy_v == NULL)  /* Do not return derivs of velocity fluxes */
     {
-      ForPatchCellsPerFace(ALL,
+      ForPatchCellsPerFace(BC_ALL,
                            BeforeAllCells(DoNothing),
                            LoopVars(i, j, k, ival, bc_struct, ipatch, sg),
                            Locals(int io, ip;
@@ -316,7 +316,7 @@ void    OverlandFlowEval(
     }
     else   /* return derivs of velocity fluxes */
     {
-      ForPatchCellsPerFace(ALL,
+      ForPatchCellsPerFace(BC_ALL,
                            BeforeAllCells(DoNothing),
                            LoopVars(i, j, k, ival, bc_struct, ipatch, sg),
                            Locals(int io, ip;
diff --git a/pfsimulator/parflow_lib/overlandflow_eval_Kin.c b/pfsimulator/parflow_lib/overlandflow_eval_Kin.c
index 4932e0140..fa78b3070 100644
--- a/pfsimulator/parflow_lib/overlandflow_eval_Kin.c
+++ b/pfsimulator/parflow_lib/overlandflow_eval_Kin.c
@@ -25,7 +25,7 @@
 
 #include "parflow.h"
 
-#if PARFLOW_ACC_BACKEND != PARFLOW_BACKEND_CUDA
+#if !defined(PARFLOW_HAVE_CUDA) && !defined(PARFLOW_HAVE_KOKKOS)
 #include "llnlmath.h"
 #endif
 /*--------------------------------------------------------------------------
@@ -102,7 +102,7 @@ void    OverlandFlowEvalKin(
 
   if (fcn == CALCFCN)
   {
-    ForPatchCellsPerFaceWithGhost(ALL,
+    ForPatchCellsPerFaceWithGhost(BC_ALL,
                                   BeforeAllCells(DoNothing),
                                   LoopVars(i, j, k, ival, bc_struct, ipatch, sg),
                                   Locals(int io, itop, ip, ipp1, ippsy;
@@ -195,7 +195,7 @@ void    OverlandFlowEvalKin(
                                   AfterAllCells(DoNothing)
       );
 
-    ForPatchCellsPerFace(ALL,
+    ForPatchCellsPerFace(BC_ALL,
                          BeforeAllCells(DoNothing),
                          LoopVars(i, j, k, ival, bc_struct, ipatch, sg),
                          Locals(int io;),
@@ -217,7 +217,7 @@ void    OverlandFlowEvalKin(
   }
   else          //fcn = CALCDER calculates the derivs
   {
-    ForPatchCellsPerFaceWithGhost(ALL,
+    ForPatchCellsPerFaceWithGhost(BC_ALL,
                                   BeforeAllCells(DoNothing),
                                   LoopVars(i, j, k, ival, bc_struct, ipatch, sg),
                                   Locals(int io, itop, ip, ipp1, ippsy;
diff --git a/pfsimulator/parflow_lib/overlandflow_eval_diffusive.c b/pfsimulator/parflow_lib/overlandflow_eval_diffusive.c
index 3b6798450..c55edd7b2 100644
--- a/pfsimulator/parflow_lib/overlandflow_eval_diffusive.c
+++ b/pfsimulator/parflow_lib/overlandflow_eval_diffusive.c
@@ -36,7 +36,7 @@
 
 #include "parflow.h"
 
-#if PARFLOW_ACC_BACKEND != PARFLOW_BACKEND_CUDA
+#if !defined(PARFLOW_HAVE_CUDA) && !defined(PARFLOW_HAVE_KOKKOS)
 #include "llnlmath.h"
 //#include "llnltyps.h"
 #endif
@@ -122,7 +122,7 @@ void    OverlandFlowEvalDiff(
 
   if (fcn == CALCFCN)
   {
-    ForPatchCellsPerFaceWithGhost(ALL,
+    ForPatchCellsPerFaceWithGhost(BC_ALL,
                                   BeforeAllCells(DoNothing),
                                   LoopVars(i, j, k, ival, bc_struct, ipatch, sg),
                                   Locals(int io, itop, ip=0, ipp1, ippsy;
@@ -226,7 +226,7 @@ void    OverlandFlowEvalDiff(
                                   AfterAllCells(DoNothing)
       );
 
-    ForPatchCellsPerFace(ALL,
+    ForPatchCellsPerFace(BC_ALL,
                          BeforeAllCells(DoNothing),
                          LoopVars(i, j, k, ival, bc_struct, ipatch, sg),
                          Locals(int io;),
@@ -249,7 +249,7 @@ void    OverlandFlowEvalDiff(
   }
   else          //fcn = CALCDER calculates the derivs of KE KW KN KS wrt to current cell (i,j,k)
   {
-    ForPatchCellsPerFace(ALL,
+    ForPatchCellsPerFace(BC_ALL,
                          BeforeAllCells(DoNothing),
                          LoopVars(i, j, k, ival, bc_struct, ipatch, sg),
                          Locals(int io, itop, ip=0, ipp1, ippsy;
diff --git a/pfsimulator/parflow_lib/parflow_proto.h b/pfsimulator/parflow_lib/parflow_proto.h
index 5f6a3bc08..c7fc543d4 100644
--- a/pfsimulator/parflow_lib/parflow_proto.h
+++ b/pfsimulator/parflow_lib/parflow_proto.h
@@ -1,6 +1,4 @@
 /* Header.c */
-
-
 typedef void (*AdvectionConcentrationInvoke) (ProblemData *problem_data, int phase, int concentration, Vector *old_concentration, Vector *new_concentration, Vector *x_velocity, Vector *y_velocity, Vector *z_velocity, Vector *solid_mass_factor, double time, double deltat, int order);
 typedef PFModule *(*AdvectionConcentrationInitInstanceXtraType) (Problem *problem, Grid *grid, double *temp_data);
 
@@ -200,9 +198,6 @@ void IntersectLineWithTriangle(unsigned int line_direction, double coord_0, doub
 void NewGlobals(char *run_name);
 void FreeGlobals(void);
 void LogGlobals(void);
-#if PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
-  void CopyGlobalsToDevice(void);
-#endif
 
 /* grgeom_list.c */
 ListMember *NewListMember(double value, int normal_component, int triangle_id);
@@ -348,6 +343,10 @@ int KinsolPCSizeOfTempData(void);
 
 typedef void (*L2ErrorNormInvoke) (double time, Vector *pressure, ProblemData *problem_data, double *l2_error_norm);
 
+/* kokkos.cpp */
+void kokkosInit();
+void kokkosFinalize();
+
 /* l2_error_norm.c */
 void L2ErrorNorm(double time, Vector *pressure, ProblemData *problem_data, double *l2_error_norm);
 PFModule *L2ErrorNormInitInstanceXtra(void);
diff --git a/pfsimulator/parflow_lib/pf_cudaloops.h b/pfsimulator/parflow_lib/pf_cudaloops.h
index d30534910..a68a6566f 100644
--- a/pfsimulator/parflow_lib/pf_cudaloops.h
+++ b/pfsimulator/parflow_lib/pf_cudaloops.h
@@ -27,7 +27,7 @@
 /*--------------------------------------------------------------------------
  * Include CUDA headers
  *--------------------------------------------------------------------------*/
-#include "pf_cudamain.h"
+#include "pf_devices.h"
 #include "pf_cudamalloc.h"
 
 extern "C++"{
@@ -491,13 +491,13 @@ DotKernel(LambdaFun loop_fun, const T init_val, T * __restrict__ rslt,
   int flagdata_size = sizeof(char) * (nz * ny * nx);                                \
   if(GrGeomSolidCellFlagDataSize(grgeom) < flagdata_size)                           \
   {                                                                                 \
-    char *flagdata = (char*)_ctalloc_cuda(flagdata_size);                           \
+    char *flagdata = (char*)_ctalloc_device(flagdata_size);                         \
                                                                                     \
     if(GrGeomSolidCellFlagDataSize(grgeom) > 0)                                     \
       CUDA_ERR(cudaMemcpy(flagdata, GrGeomSolidCellFlagData(grgeom),                \
         GrGeomSolidCellFlagDataSize(grgeom), cudaMemcpyDeviceToDevice));            \
                                                                                     \
-    _tfree_cuda(GrGeomSolidCellFlagData(grgeom));                                   \
+    _tfree_device(GrGeomSolidCellFlagData(grgeom));                                 \
     GrGeomSolidCellFlagData(grgeom) = flagdata;                                     \
     GrGeomSolidCellFlagDataSize(grgeom) = flagdata_size;                            \
   }                                                                                 \
@@ -662,7 +662,7 @@ DotKernel(LambdaFun loop_fun, const T init_val, T * __restrict__ rslt,
         loop_body;                                                                  \
       };                                                                            \
                                                                                     \
-    decltype(rslt)*ptr_rslt = (decltype(rslt)*)_talloc_cuda(sizeof(decltype(rslt)));\
+    decltype(rslt)*ptr_rslt = (decltype(rslt)*)_talloc_device(sizeof(decltype(rslt)));\
     MemPrefetchDeviceToHost_cuda(ptr_rslt, sizeof(decltype(rslt)), 0);              \
     *ptr_rslt = rslt;                                                               \
     MemPrefetchHostToDevice_cuda(ptr_rslt, sizeof(decltype(rslt)), 0);              \
@@ -675,7 +675,7 @@ DotKernel(LambdaFun loop_fun, const T init_val, T * __restrict__ rslt,
                                                                                     \
     MemPrefetchDeviceToHost_cuda(ptr_rslt, sizeof(decltype(rslt)), 0);              \
     rslt = *ptr_rslt;                                                               \
-    _tfree_cuda(ptr_rslt);                                                          \
+    _tfree_device(ptr_rslt);                                                        \
   }                                                                                 \
   (void)i;(void)j;(void)k;                                                          \
 }
@@ -715,7 +715,7 @@ DotKernel(LambdaFun loop_fun, const T init_val, T * __restrict__ rslt,
         loop_body;                                                                  \
       };                                                                            \
                                                                                     \
-    decltype(rslt)*ptr_rslt = (decltype(rslt)*)_talloc_cuda(sizeof(decltype(rslt)));\
+    decltype(rslt)*ptr_rslt = (decltype(rslt)*)_talloc_device(sizeof(decltype(rslt)));\
     MemPrefetchDeviceToHost_cuda(ptr_rslt, sizeof(decltype(rslt)), 0);              \
     *ptr_rslt = rslt;                                                               \
     MemPrefetchHostToDevice_cuda(ptr_rslt, sizeof(decltype(rslt)), 0);              \
@@ -728,7 +728,7 @@ DotKernel(LambdaFun loop_fun, const T init_val, T * __restrict__ rslt,
                                                                                     \
     MemPrefetchDeviceToHost_cuda(ptr_rslt, sizeof(decltype(rslt)), 0);              \
     rslt = *ptr_rslt;                                                               \
-    _tfree_cuda(ptr_rslt);                                                          \
+    _tfree_device(ptr_rslt);                                                        \
   }                                                                                 \
   (void)i;(void)j;(void)k;                                                          \
 }
@@ -770,6 +770,7 @@ DotKernel(LambdaFun loop_fun, const T init_val, T * __restrict__ rslt,
         dim3 block, grid;                                                           \
         FindDims(grid, block, PV_nx, PV_ny, PV_nz, 1);                              \
                                                                                     \
+        Globals *globals = ::globals;                                               \
         auto lambda_body =                                                          \
           GPU_LAMBDA(int i, int j, int k)                                           \
           {                                                                         \
@@ -810,6 +811,7 @@ DotKernel(LambdaFun loop_fun, const T init_val, T * __restrict__ rslt,
       dim3 block, grid;                                                             \
       FindDims(grid, block, nx_gpu, ny_gpu, nz_gpu, 1);                             \
                                                                                     \
+      Globals *globals = ::globals;                                                 \
       char *inflag = GrGeomSolidCellFlagData(grgeom);                               \
       auto lambda_body =                                                            \
         GPU_LAMBDA(int i, int j, int k)                                             \
@@ -1040,7 +1042,7 @@ DotKernel(LambdaFun loop_fun, const T init_val, T * __restrict__ rslt,
     if(!(ptr_ival))                                                                 \
     {                                                                               \
       GrGeomSolidCellIval(grgeom, patch_loc, PV_f) =                                \
-        (int*)_talloc_cuda(sizeof(int) * nx_bxs * ny_bxs * nz_bxs);                 \
+        (int*)_talloc_device(sizeof(int) * nx_bxs * ny_bxs * nz_bxs);               \
                                                                                     \
       ptr_ival = GrGeomSolidCellIval(grgeom, patch_loc, PV_f);                      \
       for (int idx = 0; idx < nx_bxs * ny_bxs * nz_bxs; idx++)                      \
diff --git a/pfsimulator/parflow_lib/pf_cudamalloc.h b/pfsimulator/parflow_lib/pf_cudamalloc.h
index 94563324a..ea8060b73 100644
--- a/pfsimulator/parflow_lib/pf_cudamalloc.h
+++ b/pfsimulator/parflow_lib/pf_cudamalloc.h
@@ -24,8 +24,7 @@
 #ifndef PF_CUDAMALLOC_H
 #define PF_CUDAMALLOC_H
 
-#include <string.h>
-#include "pf_cudamain.h"
+#include "pf_devices.h"
 
 /*--------------------------------------------------------------------------
  * Memory management macros for CUDA
@@ -38,12 +37,12 @@
 #define tfree_amps_cuda(ptr) amps_TFree_managed(ptr)
 
 #define talloc_cuda(type, count) \
-  ((count) ? (type*)_talloc_cuda(sizeof(type) * (unsigned int)(count)) : NULL)
+  ((count) ? (type*)_talloc_device(sizeof(type) * (unsigned int)(count)) : NULL)
 
 #define ctalloc_cuda(type, count) \
-  ((count) ? (type*)_ctalloc_cuda(sizeof(type) * (unsigned int)(count)) : NULL)
+  ((count) ? (type*)_ctalloc_device(sizeof(type) * (unsigned int)(count)) : NULL)
 
-#define tfree_cuda(ptr) if (ptr) _tfree_cuda(ptr); else {}
+#define tfree_cuda(ptr) if (ptr) _tfree_device(ptr); else {}
 
 #define tmemcpy_cuda(dest, src, bytes) \
   CUDA_ERR(cudaMemcpy(dest, src, bytes, cudaMemcpyDeviceToDevice))
diff --git a/pfsimulator/parflow_lib/pf_cudamain.h b/pfsimulator/parflow_lib/pf_devices.h
similarity index 82%
rename from pfsimulator/parflow_lib/pf_cudamain.h
rename to pfsimulator/parflow_lib/pf_devices.h
index 3c29782f1..b4264cdb4 100644
--- a/pfsimulator/parflow_lib/pf_cudamain.h
+++ b/pfsimulator/parflow_lib/pf_devices.h
@@ -18,12 +18,15 @@
  ***********************************************************************/
 
 /** @file
- * @brief Contains general CUDA related macros and functions.
+ * @brief Contains general device-related macros and functions.
  */
 
-#ifndef PF_CUDAMAIN_H
-#define PF_CUDAMAIN_H
+#ifndef PF_DEVICES_H
+#define PF_DEVICES_H
 
+#include <string.h>
+
+#ifdef PARFLOW_HAVE_CUDA
 #include <cuda.h>
 #include <cuda_runtime.h>
 #include <stdbool.h>
@@ -48,25 +51,6 @@
   }                                                                                    \
 }
 
-#ifdef PARFLOW_HAVE_RMM
-#include <rmm/rmm_api.h>
-/**
- * @brief RMM error handling.
- * 
- * If error detected, print error message and exit.
- *
- * @param expr RMM error (of type rmmError_t) [IN]
- */
-#define RMM_ERR(expr)                                                                  \
-{                                                                                      \
-  rmmError_t err = expr;                                                               \
-  if (err != RMM_SUCCESS) {                                                            \
-    printf("\n\n%s in %s at line %d\n", rmmGetErrorString(err), __FILE__, __LINE__);   \
-    exit(1);                                                                           \
-  }                                                                                    \
-}
-#endif
-
 /*--------------------------------------------------------------------------
  * CUDA profiling macros
  *--------------------------------------------------------------------------*/
@@ -92,10 +76,51 @@
 /** Stop recording an NVTX range for NSYS if accelerator present. */
 #define POP_NVTX_cuda nvtxRangePop();
 
+#endif // PARFLOW_HAVE_CUDA
+
+#ifdef PARFLOW_HAVE_RMM
+#include <rmm/rmm_api.h>
+/**
+ * @brief RMM error handling.
+ * 
+ * If error detected, print error message and exit.
+ *
+ * @param expr RMM error (of type rmmError_t) [IN]
+ */
+#define RMM_ERR(expr)                                                                  \
+{                                                                                      \
+  rmmError_t err = expr;                                                               \
+  if (err != RMM_SUCCESS) {                                                            \
+    printf("\n\n%s in %s at line %d\n", rmmGetErrorString(err), __FILE__, __LINE__);   \
+    exit(1);                                                                           \
+  }                                                                                    \
+}
+#endif
+
 /*--------------------------------------------------------------------------
- * Define static unified memory allocation routines for CUDA
+ * Define static unified memory allocation routines for device
  *--------------------------------------------------------------------------*/
 
+/**
+ * @brief Kokkos C wrapper declaration for memory allocation.
+ */
+void* kokkosAlloc(size_t size);
+
+/**
+ * @brief Kokkos C wrapper declaration for memory deallocation.
+ */
+void kokkosFree(void *ptr);
+
+/**
+ * @brief Kokkos C wrapper declaration for memory copy.
+ */
+void kokkosMemCpy(char *dest, char *src, size_t size);
+
+/**
+ * @brief Kokkos C wrapper declaration for memset.
+ */
+void kokkosMemSet(char *ptr, size_t size);
+
 /**
  * @brief Allocates unified memory.
  * 
@@ -106,17 +131,20 @@
  * @param size bytes to be allocated [IN]
  * @return a void pointer to the allocated dataspace
  */
-static inline void *_talloc_cuda(size_t size)
+
+static inline void *_talloc_device(size_t size)
 {
   void *ptr = NULL;  
   
 #ifdef PARFLOW_HAVE_RMM
   RMM_ERR(rmmAlloc(&ptr,size,0,__FILE__,__LINE__));
-#else
+#elif defined(PARFLOW_HAVE_KOKKOS)
+  ptr = kokkosAlloc(size);
+#elif defined(PARFLOW_HAVE_CUDA)
   CUDA_ERR(cudaMallocManaged((void**)&ptr, size, cudaMemAttachGlobal));
   // CUDA_ERR(cudaHostAlloc((void**)&ptr, size, cudaHostAllocMapped));  
 #endif
-  
+
   return ptr;
 }
 
@@ -130,37 +158,46 @@ static inline void *_talloc_cuda(size_t size)
  * @param size bytes to be allocated [IN]
  * @return a void pointer to the allocated dataspace
  */
-static inline void *_ctalloc_cuda(size_t size)
+static inline void *_ctalloc_device(size_t size)
 {
   void *ptr = NULL;  
 
 #ifdef PARFLOW_HAVE_RMM
   RMM_ERR(rmmAlloc(&ptr,size,0,__FILE__,__LINE__));
-#else
+#elif defined(PARFLOW_HAVE_KOKKOS)
+  ptr = kokkosAlloc(size);
+#elif defined(PARFLOW_HAVE_CUDA)
   CUDA_ERR(cudaMallocManaged((void**)&ptr, size, cudaMemAttachGlobal));
   // CUDA_ERR(cudaHostAlloc((void**)&ptr, size, cudaHostAllocMapped));
 #endif  
-  // memset(ptr, 0, size);
+
+#if defined(PARFLOW_HAVE_CUDA)
   CUDA_ERR(cudaMemset(ptr, 0, size));  
+#else
+  // memset(ptr, 0, size);
+  kokkosMemSet((char*)ptr, size);
+#endif
   
   return ptr;
 }
 
 /**
- * @brief Frees unified memory allocated with \ref _talloc_cuda or \ref _ctalloc_cuda.
+ * @brief Frees unified memory allocated with \ref _talloc_device or \ref _ctalloc_device.
  * 
  * @note Should not be called directly.
  *
  * @param ptr a void pointer to the allocated dataspace [IN]
  */
-static inline void _tfree_cuda(void *ptr)
+static inline void _tfree_device(void *ptr)
 {
 #ifdef PARFLOW_HAVE_RMM
   RMM_ERR(rmmFree(ptr,0,__FILE__,__LINE__));
-#else
+#elif defined(PARFLOW_HAVE_KOKKOS)
+  kokkosFree(ptr);
+#elif defined(PARFLOW_HAVE_CUDA)
   CUDA_ERR(cudaFree(ptr));
   // CUDA_ERR(cudaFreeHost(ptr));
 #endif
 }
 
-#endif // PF_CUDAMAIN_H
+#endif // PF_DEVICES_H
diff --git a/pfsimulator/parflow_lib/pf_kokkos.cpp b/pfsimulator/parflow_lib/pf_kokkos.cpp
new file mode 100644
index 000000000..58ec0eb0c
--- /dev/null
+++ b/pfsimulator/parflow_lib/pf_kokkos.cpp
@@ -0,0 +1,79 @@
+#include "amps.h"
+#include <Kokkos_Core.hpp>
+
+extern "C"{
+#include "pf_devices.h"
+
+/**
+ * @brief Allocate memory with Kokkos
+ *
+ * @param size The size of the allocation in bytes
+ */
+void* kokkosAlloc(size_t size){
+  return Kokkos::kokkos_malloc(size);
+}
+
+/**
+ * @brief Free memory with Kokkos
+ *
+ * @param ptr Freed pointer
+ */
+void kokkosFree(void *ptr){
+  Kokkos::kokkos_free(ptr);
+}
+
+/**
+ * @brief Kokkos memcopy
+ *
+ * @param dest Destination pointer
+ * @param src Source pointer
+ * @param size Bytes to be copied
+ */
+void kokkosMemCpy(char *dest, char *src, size_t size){
+  Kokkos::View<char*> dest_view(dest, size);
+  Kokkos::View<char*> src_view(src, size);
+  Kokkos::deep_copy(dest_view, src_view);
+}
+
+/**
+ * @brief Kokkos memset
+ *
+ * @param ptr Pointer for the data to be set to zero
+ * @param size Bytes to be zeroed
+ */
+void kokkosMemSet(char *ptr, size_t size){
+  /* Loop style initialization */
+  if(size % sizeof(int))
+  {
+    /* Writing 1 byte / thread is slow */
+    Kokkos::parallel_for(size, KOKKOS_LAMBDA(int i){ptr[i] = 0;});
+  }
+  else
+  {
+    /* Writing 4 bytes / thread is fast */
+    Kokkos::parallel_for(size / sizeof(int), KOKKOS_LAMBDA(int i){((int*)ptr)[i] = 0;});
+  }
+  Kokkos::fence(); 
+
+  /* Deep_copy style initialization for char* should be fast in future Kokkos releases */
+  // Kokkos::View<char*> ptr_view(ptr, size);
+  // Kokkos::deep_copy(ptr_view, 0);
+}
+
+/**
+ * @brief Initialize Kokkos
+ */
+void kokkosInit(){
+  Kokkos::InitArguments args;
+  args.ndevices = 1;
+  Kokkos::initialize(args);    
+}
+
+/**
+ * @brief Finalize Kokkos
+ */
+void kokkosFinalize(){
+  if(Kokkos::is_initialized) Kokkos::finalize();
+}
+
+}
diff --git a/pfsimulator/parflow_lib/pf_kokkosloops.h b/pfsimulator/parflow_lib/pf_kokkosloops.h
new file mode 100644
index 000000000..3777ad79b
--- /dev/null
+++ b/pfsimulator/parflow_lib/pf_kokkosloops.h
@@ -0,0 +1,936 @@
+/**********************************************************************
+ *
+ *  Please read the LICENSE file for the GNU Lesser General Public License.
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License (as published
+ *  by the Free Software Foundation) version 2.1 dated February 1999.
+ *
+ *  This program is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the IMPLIED WARRANTY OF
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the terms
+ *  and conditions of the GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Lesser General Public
+ *  License along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+ *  USA
+ ***********************************************************************/
+
+/* @file
+ * @brief Contains macros, functions, and structs for Kokkos compute kernels.
+ */
+
+#ifndef PF_KOKKOSLOOPS_H
+#define PF_KOKKOSLOOPS_H
+
+/*--------------------------------------------------------------------------
+ * Include headers
+ *--------------------------------------------------------------------------*/
+#include "pf_devices.h"
+#include "pf_kokkosmalloc.h"
+
+extern "C++"{
+
+#include <tuple>
+#include <Kokkos_Core.hpp>
+
+/*--------------------------------------------------------------------------
+ * Helper macros and functions
+ *--------------------------------------------------------------------------*/
+#define RAND48_SEED_0   (0x330e)
+#define RAND48_SEED_1   (0xabcd)
+#define RAND48_SEED_2   (0x1234)
+#define RAND48_MULT_0   (0xe66d)
+#define RAND48_MULT_1   (0xdeec)
+#define RAND48_MULT_2   (0x0005)
+#define RAND48_ADD      (0x000b)
+
+ /** Helper struct for type comparison. @note Not for direct use! */
+template <typename T>
+struct function_traits
+    : public function_traits<decltype(&T::operator())>
+{};
+// For generic types, directly use the result of the signature of its 'operator()'
+
+ /** Helper struct for type comparison. @note Not for direct use! */
+template <typename ClassType, typename ReturnType, typename... Args>
+struct function_traits<ReturnType(ClassType::*)(Args...) const>
+// we specialize for pointers to member function
+{
+    enum { arity = sizeof...(Args) };
+    // arity is the number of arguments.
+
+    typedef ReturnType result_type;
+
+    template <size_t i>
+    struct arg
+    {
+        typedef typename std::tuple_element<i, std::tuple<Args...>>::type type;
+        // the i-th argument is equivalent to the i-th tuple element of a tuple
+        // composed of those arguments.
+    };
+};
+
+/** Device-callable dorand48() function for Kokkos compute kernels. */
+KOKKOS_FORCEINLINE_FUNCTION static void dev_dorand48(unsigned short xseed[3])
+{
+  unsigned long accu;
+  
+  unsigned short _rand48_mult[3] = {
+    RAND48_MULT_0,
+    RAND48_MULT_1,
+    RAND48_MULT_2
+  };
+  unsigned short _rand48_add = RAND48_ADD;
+  unsigned short temp[2];
+
+  accu = (unsigned long)_rand48_mult[0] * (unsigned long)xseed[0] +
+         (unsigned long)_rand48_add;
+  temp[0] = (unsigned short)accu;               /* lower 16 bits */
+  accu >>= sizeof(unsigned short) * 8;
+  accu += (unsigned long)_rand48_mult[0] * (unsigned long)xseed[1] +
+          (unsigned long)_rand48_mult[1] * (unsigned long)xseed[0];
+  temp[1] = (unsigned short)accu;               /* middle 16 bits */
+  accu >>= sizeof(unsigned short) * 8;
+  accu += _rand48_mult[0] * xseed[2] + _rand48_mult[1] * xseed[1] + _rand48_mult[2] * xseed[0];
+  xseed[0] = temp[0];
+  xseed[1] = temp[1];
+  xseed[2] = (unsigned short)accu;
+}
+
+/** Device-callable erand48() function for Kokkos compute kernels. */
+KOKKOS_FORCEINLINE_FUNCTION static double dev_erand48(unsigned short xseed[3])
+{
+  dev_dorand48(xseed);
+  return ldexp((double)xseed[0], -48) +
+         ldexp((double)xseed[1], -32) +
+         ldexp((double)xseed[2], -16);
+}
+
+/** Device-callable drand48() function for Kokkos compute kernels. */
+KOKKOS_FORCEINLINE_FUNCTION static double dev_drand48(void)
+{
+  unsigned short _rand48_seed[3] = {
+    RAND48_SEED_0,
+    RAND48_SEED_1,
+    RAND48_SEED_2
+  };
+  return dev_erand48(_rand48_seed);
+}
+
+/** Device-callable RPowerR() function for Kokkos compute kernels. */
+template <typename T>
+KOKKOS_FORCEINLINE_FUNCTION static T RPowerR(T base, T exponent)
+{
+  if (base <= 0.0)
+    return(0.0);
+
+  return((T)pow((double)base, (double)exponent));
+}
+
+ /** Helper struct for type comparison. @note Not for direct use! */
+struct SkipParallelSync {const int dummy = 0;};
+#define SKIP_PARALLEL_SYNC_kokkos struct SkipParallelSync sync_struct; return sync_struct;
+
+#define PARALLEL_SYNC_kokkos Kokkos::fence(); 
+
+#define PlusEquals_kokkos(a, b) Kokkos::atomic_add(&(a), b)
+
+ /** Helper struct for type comparison. @note Not for direct use! */
+template <typename T>
+struct ReduceMaxType {T value;};
+#define ReduceMax_kokkos(a, b) if(lsum < b) { lsum = b; } struct ReduceMaxType<std::decay<decltype(a)>::type> reduce_struct; return reduce_struct;
+
+ /** Helper struct for type comparison. @note Not for direct use! */
+template <typename T>
+struct ReduceMinType {T value;};
+#define ReduceMin_kokkos(a, b) if(lsum > b) { lsum = b; } struct ReduceMinType<std::decay<decltype(a)>::type> reduce_struct; return reduce_struct;
+
+ /** Helper struct for type comparison. @note Not for direct use! */
+template <typename T>
+struct ReduceSumType {T value;};
+#define ReduceSum_kokkos(a, b) lsum += b; struct ReduceSumType<std::decay<decltype(a)>::type> reduce_struct; return reduce_struct;
+
+
+/** A constant table for fdir (borrowed from OMP backend) */
+static const int FDIR_TABLE[6][3] = {
+  {-1, 0,  0}, // Left
+  {1,  0,  0}, // Right
+  {0, -1,  0}, // Down
+  {0,  1,  0}, // Up
+  {0,  0, -1}, // Back
+  {0,  0,  1}, // Front
+};
+
+/*--------------------------------------------------------------------------
+ * Kokkos loop macro redefinitions
+ *--------------------------------------------------------------------------*/
+
+/**
+ * @brief A macro for checking if cell flag data array must be reallocated.
+ *
+ * @note Not for direct use!
+ *
+ * @param grgeom The geometry details [IN/OUT]
+ * @param nx The size of the first dim [IN]
+ * @param ny The size of the second dim [IN]
+ * @param nz The size of the third dim [IN]
+ */
+#define CheckCellFlagAllocation(grgeom, nx, ny, nz)                                 \
+{                                                                                   \
+  int flagdata_size = sizeof(char) * (nz * ny * nx);                                \
+  if(GrGeomSolidCellFlagDataSize(grgeom) < flagdata_size)                           \
+  {                                                                                 \
+    char *flagdata = (char*)_ctalloc_device(flagdata_size);                         \
+                                                                                    \
+    if(GrGeomSolidCellFlagDataSize(grgeom) > 0)                                     \
+      memcpy(flagdata, GrGeomSolidCellFlagData(grgeom),                             \
+        GrGeomSolidCellFlagDataSize(grgeom));                                       \
+                                                                                    \
+    _tfree_device(GrGeomSolidCellFlagData(grgeom));                                 \
+    GrGeomSolidCellFlagData(grgeom) = flagdata;                                     \
+    GrGeomSolidCellFlagDataSize(grgeom) = flagdata_size;                            \
+  }                                                                                 \
+}
+
+ /** Loop definition for Kokkos. */
+#define BoxLoopI1_kokkos(i, j, k,                                                   \
+  ix, iy, iz, nx, ny, nz,                                                           \
+  i1, nx1, ny1, nz1, sx1, sy1, sz1,                                                 \
+  loop_body)                                                                        \
+{                                                                                   \
+  if(nx > 0 && ny > 0 && nz > 0)                                                    \
+  {                                                                                 \
+    DeclareInc(PV_jinc_1, PV_kinc_1, nx, ny, nz, nx1, ny1, nz1, sx1, sy1, sz1);     \
+                                                                                    \
+    const auto &ref_i1 = i1;                                                        \
+                                                                                    \
+    auto lambda_body =                                                              \
+      KOKKOS_LAMBDA(int i, int j, int k)                                            \
+      {                                                                             \
+        const int i1 = k * PV_kinc_1 + (k * ny + j) * PV_jinc_1                     \
+          + (k * ny * nx + j * nx + i) * sx1 + ref_i1;                              \
+                                                                                    \
+        i += ix;                                                                    \
+        j += iy;                                                                    \
+        k += iz;                                                                    \
+                                                                                    \
+        loop_body;                                                                  \
+      };                                                                            \
+                                                                                    \
+    using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+    MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{nx, ny, nz}});                       \
+    Kokkos::parallel_for(mdpolicy_3d, lambda_body);                                 \
+                                                                                    \
+    typedef function_traits<decltype(lambda_body)> traits;                          \
+    if(!std::is_same<traits::result_type, struct SkipParallelSync>::value)          \
+      Kokkos::fence();                                                              \
+  }                                                                                 \
+  (void)i;(void)j;(void)k;                                                          \
+}
+
+ /** Loop definition for Kokkos. */
+#define BoxLoopI2_kokkos(i, j, k,                                                   \
+  ix, iy, iz, nx, ny, nz,                                                           \
+  i1, nx1, ny1, nz1, sx1, sy1, sz1,                                                 \
+  i2, nx2, ny2, nz2, sx2, sy2, sz2,                                                 \
+  loop_body)                                                                        \
+{                                                                                   \
+  if(nx > 0 && ny > 0 && nz > 0)                                                    \
+  {                                                                                 \
+    DeclareInc(PV_jinc_1, PV_kinc_1, nx, ny, nz, nx1, ny1, nz1, sx1, sy1, sz1);     \
+    DeclareInc(PV_jinc_2, PV_kinc_2, nx, ny, nz, nx2, ny2, nz2, sx2, sy2, sz2);     \
+                                                                                    \
+    const auto &ref_i1 = i1;                                                        \
+    const auto &ref_i2 = i2;                                                        \
+                                                                                    \
+    auto lambda_body =                                                              \
+      KOKKOS_LAMBDA(int i, int j, int k)                                            \
+      {                                                                             \
+        const int i1 = k * PV_kinc_1 + (k * ny + j) * PV_jinc_1                     \
+          + (k * ny * nx + j * nx + i) * sx1 + ref_i1;                              \
+        const int i2 = k * PV_kinc_2 + (k * ny + j) * PV_jinc_2                     \
+          + (k * ny * nx + j * nx + i) * sx2 + ref_i2;                              \
+                                                                                    \
+        i += ix;                                                                    \
+        j += iy;                                                                    \
+        k += iz;                                                                    \
+                                                                                    \
+        loop_body;                                                                  \
+      };                                                                            \
+                                                                                    \
+    using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+    MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{nx, ny, nz}});                       \
+    Kokkos::parallel_for(mdpolicy_3d, lambda_body);                                 \
+                                                                                    \
+    typedef function_traits<decltype(lambda_body)> traits;                          \
+    if(!std::is_same<traits::result_type, struct SkipParallelSync>::value)          \
+      Kokkos::fence();                                                              \
+  }                                                                                 \
+  (void)i;(void)j;(void)k;                                                          \
+}
+
+ /** Loop definition for Kokkos. */
+#define BoxLoopI3_kokkos(i, j, k,                                                   \
+  ix, iy, iz, nx, ny, nz,                                                           \
+  i1, nx1, ny1, nz1, sx1, sy1, sz1,                                                 \
+  i2, nx2, ny2, nz2, sx2, sy2, sz2,                                                 \
+  i3, nx3, ny3, nz3, sx3, sy3, sz3,                                                 \
+  loop_body)                                                                        \
+{                                                                                   \
+  if(nx > 0 && ny > 0 && nz > 0)                                                    \
+  {                                                                                 \
+    DeclareInc(PV_jinc_1, PV_kinc_1, nx, ny, nz, nx1, ny1, nz1, sx1, sy1, sz1);     \
+    DeclareInc(PV_jinc_2, PV_kinc_2, nx, ny, nz, nx2, ny2, nz2, sx2, sy2, sz2);     \
+    DeclareInc(PV_jinc_3, PV_kinc_3, nx, ny, nz, nx3, ny3, nz3, sx3, sy3, sz3);     \
+                                                                                    \
+    const auto &ref_i1 = i1;                                                        \
+    const auto &ref_i2 = i2;                                                        \
+    const auto &ref_i3 = i3;                                                        \
+                                                                                    \
+    auto lambda_body =                                                              \
+      KOKKOS_LAMBDA(int i, int j, int k)                                            \
+      {                                                                             \
+        const int i1 = k * PV_kinc_1 + (k * ny + j) * PV_jinc_1                     \
+          + (k * ny * nx + j * nx + i) * sx1 + ref_i1;                              \
+        const int i2 = k * PV_kinc_2 + (k * ny + j) * PV_jinc_2                     \
+          + (k * ny * nx + j * nx + i) * sx2 + ref_i2;                              \
+        const int i3 = k * PV_kinc_3 + (k * ny + j) * PV_jinc_3                     \
+          + (k * ny * nx + j * nx + i) * sx3 + ref_i3;                              \
+                                                                                    \
+        i += ix;                                                                    \
+        j += iy;                                                                    \
+        k += iz;                                                                    \
+                                                                                    \
+        loop_body;                                                                  \
+      };                                                                            \
+                                                                                    \
+    using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+    MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{nx, ny, nz}});                       \
+    Kokkos::parallel_for(mdpolicy_3d, lambda_body);                                 \
+                                                                                    \
+    typedef function_traits<decltype(lambda_body)> traits;                          \
+    if(!std::is_same<traits::result_type, struct SkipParallelSync>::value)          \
+      Kokkos::fence();                                                              \
+  }                                                                                 \
+  (void)i;(void)j;(void)k;                                                          \
+}
+
+ /** Loop definition for Kokkos. */
+#define BoxLoopReduceI1_kokkos(rslt, i, j, k,                                       \
+  ix, iy, iz, nx, ny, nz,                                                           \
+  i1, nx1, ny1, nz1, sx1, sy1, sz1,                                                 \
+  loop_body)                                                                        \
+{                                                                                   \
+  if(nx > 0 && ny > 0 && nz > 0)                                                    \
+  {                                                                                 \
+    DeclareInc(PV_jinc_1, PV_kinc_1, nx, ny, nz, nx1, ny1, nz1, sx1, sy1, sz1);     \
+                                                                                    \
+    const auto &ref_rslt = rslt;                                                    \
+    const auto &ref_i1 = i1;                                                        \
+                                                                                    \
+    auto lambda_body =                                                              \
+      KOKKOS_LAMBDA(int i, int j, int k, double& lsum)                              \
+      {                                                                             \
+        const int i1 = k * PV_kinc_1 + (k * ny + j) * PV_jinc_1                     \
+          + (k * ny * nx + j * nx + i) * sx1 + ref_i1;                              \
+                                                                                    \
+        i += ix;                                                                    \
+        j += iy;                                                                    \
+        k += iz;                                                                    \
+                                                                                    \
+        loop_body;                                                                  \
+      };                                                                            \
+                                                                                    \
+    auto lambda_outer =                                                             \
+      KOKKOS_LAMBDA(int i, int j, int k, double& lsum)                              \
+      {                                                                             \
+        lambda_body(i, j, k, lsum);                                                 \
+      };                                                                            \
+                                                                                    \
+    using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+    MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{nx, ny, nz}});                       \
+    typedef function_traits<decltype(lambda_body)> traits;                          \
+    if(std::is_same<traits::result_type, struct ReduceSumType<double>>::value)      \
+    {                                                                               \
+      Kokkos::parallel_reduce(mdpolicy_3d, lambda_outer, rslt);                     \
+    }                                                                               \
+    else if(std::is_same<traits::result_type, struct ReduceMaxType<double>>::value) \
+    {                                                                               \
+      Kokkos::parallel_reduce(mdpolicy_3d, lambda_outer, Kokkos::Max<double>(rslt));\
+    }                                                                               \
+    else if(std::is_same<traits::result_type, struct ReduceMinType<double>>::value) \
+    {                                                                               \
+      Kokkos::parallel_reduce(mdpolicy_3d, lambda_outer, Kokkos::Min<double>(rslt));\
+    }                                                                               \
+    else                                                                            \
+    {                                                                               \
+      printf("ERROR at %s:%d: Invalid reduction identifier,                         \
+      likely a problem with a BoxLoopReduce body.", __FILE__, __LINE__);            \
+    }                                                                               \
+    Kokkos::fence();                                                                \
+  }                                                                                 \
+  (void)i;(void)j;(void)k;                                                          \
+}
+
+ /** Loop definition for Kokkos. */
+#define BoxLoopReduceI2_kokkos(rslt, i, j, k,                                       \
+  ix, iy, iz, nx, ny, nz,                                                           \
+  i1, nx1, ny1, nz1, sx1, sy1, sz1,                                                 \
+  i2, nx2, ny2, nz2, sx2, sy2, sz2,                                                 \
+  loop_body)                                                                        \
+{                                                                                   \
+  if(nx > 0 && ny > 0 && nz > 0)                                                    \
+  {                                                                                 \
+    DeclareInc(PV_jinc_1, PV_kinc_1, nx, ny, nz, nx1, ny1, nz1, sx1, sy1, sz1);     \
+    DeclareInc(PV_jinc_2, PV_kinc_2, nx, ny, nz, nx2, ny2, nz2, sx2, sy2, sz2);     \
+                                                                                    \
+    const auto &ref_rslt = rslt;                                                    \
+    const auto &ref_i1 = i1;                                                        \
+    const auto &ref_i2 = i2;                                                        \
+                                                                                    \
+    auto lambda_body =                                                              \
+      KOKKOS_LAMBDA(int i, int j, int k, double& lsum)                              \
+      {                                                                             \
+        const int i1 = k * PV_kinc_1 + (k * ny + j) * PV_jinc_1                     \
+          + (k * ny * nx + j * nx + i) * sx1 + ref_i1;                              \
+        const int i2 = k * PV_kinc_2 + (k * ny + j) * PV_jinc_2                     \
+          + (k * ny * nx + j * nx + i) * sx2 + ref_i2;                              \
+                                                                                    \
+        i += ix;                                                                    \
+        j += iy;                                                                    \
+        k += iz;                                                                    \
+                                                                                    \
+        loop_body;                                                                  \
+      };                                                                            \
+                                                                                    \
+    auto lambda_outer =                                                             \
+      KOKKOS_LAMBDA(int i, int j, int k, double& lsum)                              \
+      {                                                                             \
+        lambda_body(i, j, k, lsum);                                                 \
+      };                                                                            \
+                                                                                    \
+    using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+    MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{nx, ny, nz}});                       \
+    typedef function_traits<decltype(lambda_body)> traits;                          \
+    if(std::is_same<traits::result_type, struct ReduceSumType<double>>::value)      \
+    {                                                                               \
+      Kokkos::parallel_reduce(mdpolicy_3d, lambda_outer, rslt);                     \
+    }                                                                               \
+    else if(std::is_same<traits::result_type, struct ReduceMaxType<double>>::value) \
+    {                                                                               \
+      Kokkos::parallel_reduce(mdpolicy_3d, lambda_outer, Kokkos::Max<double>(rslt));\
+    }                                                                               \
+    else if(std::is_same<traits::result_type, struct ReduceMinType<double>>::value) \
+    {                                                                               \
+      Kokkos::parallel_reduce(mdpolicy_3d, lambda_outer, Kokkos::Min<double>(rslt));\
+    }                                                                               \
+    else                                                                            \
+    {                                                                               \
+      printf("ERROR at %s:%d: Invalid reduction identifier,                         \
+      likely a problem with a BoxLoopReduce body.", __FILE__, __LINE__);            \
+    }                                                                               \
+    Kokkos::fence();                                                                \
+  }                                                                                 \
+  (void)i;(void)j;(void)k;                                                          \
+}
+
+ /** Loop definition for Kokkos. */
+#define GrGeomInLoopBoxes_kokkos(i, j, k,                                           \
+  grgeom, ix, iy, iz, nx, ny, nz, loop_body)                                        \
+{                                                                                   \
+  BoxArray* boxes = GrGeomSolidInteriorBoxes(grgeom);                               \
+  int ix_bxs = BoxArrayMinCell(boxes, 0);                                           \
+  int iy_bxs = BoxArrayMinCell(boxes, 1);                                           \
+  int iz_bxs = BoxArrayMinCell(boxes, 2);                                           \
+                                                                                    \
+  int nx_bxs = BoxArrayMaxCell(boxes, 0) - ix_bxs + 1;                              \
+  int ny_bxs = BoxArrayMaxCell(boxes, 1) - iy_bxs + 1;                              \
+  int nz_bxs = BoxArrayMaxCell(boxes, 2) - iz_bxs + 1;                              \
+                                                                                    \
+  if(!(GrGeomSolidCellFlagInitialized(grgeom) & 1))                                 \
+  {                                                                                 \
+    CheckCellFlagAllocation(grgeom, nx_bxs, ny_bxs, nz_bxs);                        \
+    char *inflag = GrGeomSolidCellFlagData(grgeom);                                 \
+                                                                                    \
+    for (int PV_box = 0; PV_box < BoxArraySize(boxes); PV_box++)                    \
+    {                                                                               \
+      Box box = BoxArrayGetBox(boxes, PV_box);                                      \
+      int PV_ixl = box.lo[0];                                                       \
+      int PV_iyl = box.lo[1];                                                       \
+      int PV_izl = box.lo[2];                                                       \
+      int PV_ixu = box.up[0];                                                       \
+      int PV_iyu = box.up[1];                                                       \
+      int PV_izu = box.up[2];                                                       \
+                                                                                    \
+      if(PV_ixl <= PV_ixu && PV_iyl <= PV_iyu && PV_izl <= PV_izu)                  \
+      {                                                                             \
+        int PV_nx = PV_ixu - PV_ixl + 1;                                            \
+        int PV_ny = PV_iyu - PV_iyl + 1;                                            \
+        int PV_nz = PV_izu - PV_izl + 1;                                            \
+                                                                                    \
+        Globals *globals = ::globals;                                               \
+        auto lambda_body =                                                          \
+          KOKKOS_LAMBDA(int i, int j, int k)                                        \
+          {                                                                         \
+            i += PV_ixl;                                                            \
+            j += PV_iyl;                                                            \
+            k += PV_izl;                                                            \
+                                                                                    \
+            /* Set inflag for all cells in boxes regardless of loop limits */       \
+            inflag[(k - iz_bxs) * ny_bxs * nx_bxs +                                 \
+              (j - iy_bxs) * nx_bxs + (i - ix_bxs)] |= 1;                           \
+                                                                                    \
+            /* Only evaluate loop body if the cell is within loop limits */         \
+            if(i >= ix && j >= iy && k >= iz &&                                     \
+              i < ix + nx && j < iy + ny && k < iz + nz)                            \
+            {                                                                       \
+              loop_body;                                                            \
+            }                                                                       \
+          };                                                                        \
+                                                                                    \
+        using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+        MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{PV_nx, PV_ny, PV_nz}});          \
+        Kokkos::parallel_for(mdpolicy_3d, lambda_body);                             \
+        Kokkos::fence();                                                            \
+      }                                                                             \
+    }                                                                               \
+    GrGeomSolidCellFlagInitialized(grgeom) |= 1;                                    \
+  }                                                                                 \
+  else                                                                              \
+  {                                                                                 \
+    int ixl_gpu = pfmax(ix, ix_bxs);                                                \
+    int iyl_gpu = pfmax(iy, iy_bxs);                                                \
+    int izl_gpu = pfmax(iz, iz_bxs);                                                \
+    int nx_gpu = pfmin((ix + nx - 1), BoxArrayMaxCell(boxes, 0)) - ixl_gpu + 1;     \
+    int ny_gpu = pfmin((iy + ny - 1), BoxArrayMaxCell(boxes, 1)) - iyl_gpu + 1;     \
+    int nz_gpu = pfmin((iz + nz - 1), BoxArrayMaxCell(boxes, 2)) - izl_gpu + 1;     \
+                                                                                    \
+    if(nx_gpu > 0 && ny_gpu > 0 && nz_gpu > 0)                                      \
+    {                                                                               \
+      Globals *globals = ::globals;                                                 \
+      char *inflag = GrGeomSolidCellFlagData(grgeom);                               \
+      auto lambda_body =                                                            \
+        KOKKOS_LAMBDA(int i, int j, int k)                                          \
+        {                                                                           \
+          i += ixl_gpu;                                                             \
+          j += iyl_gpu;                                                             \
+          k += izl_gpu;                                                             \
+          if(inflag[(k - iz_bxs) * ny_bxs * nx_bxs +                                \
+            (j - iy_bxs) * nx_bxs + (i - ix_bxs)] & 1)                              \
+          {                                                                         \
+            loop_body;                                                              \
+          }                                                                         \
+        };                                                                          \
+                                                                                    \
+      using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+      MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{nx_gpu, ny_gpu, nz_gpu}});         \
+      Kokkos::parallel_for(mdpolicy_3d, lambda_body);                               \
+      Kokkos::fence();                                                              \
+    }                                                                               \
+  }                                                                                 \
+  (void)i;(void)j;(void)k;                                                          \
+}
+
+ /** Loop definition for Kokkos. */
+#define GrGeomSurfLoopBoxes_kokkos(i, j, k, fdir, grgeom,                           \
+  ix, iy, iz, nx, ny, nz, loop_body)                                                \
+{                                                                                   \
+  for (int PV_f = 0; PV_f < GrGeomOctreeNumFaces; PV_f++)                           \
+  {                                                                                 \
+    const int *fdir = FDIR_TABLE[PV_f];                                             \
+                                                                                    \
+    BoxArray* boxes = GrGeomSolidSurfaceBoxes(grgeom, PV_f);                        \
+    int ix_bxs = BoxArrayMinCell(boxes, 0);                                         \
+    int iy_bxs = BoxArrayMinCell(boxes, 1);                                         \
+    int iz_bxs = BoxArrayMinCell(boxes, 2);                                         \
+                                                                                    \
+    int nx_bxs = BoxArrayMaxCell(boxes, 0) - ix_bxs + 1;                            \
+    int ny_bxs = BoxArrayMaxCell(boxes, 1) - iy_bxs + 1;                            \
+    int nz_bxs = BoxArrayMaxCell(boxes, 2) - iz_bxs + 1;                            \
+                                                                                    \
+    if(!(GrGeomSolidCellFlagInitialized(grgeom) & (1 << (2 + PV_f))))               \
+    {                                                                               \
+      CheckCellFlagAllocation(grgeom, nx_bxs, ny_bxs, nz_bxs);                      \
+      char *surfflag = GrGeomSolidCellFlagData(grgeom);                             \
+                                                                                    \
+      for (int PV_box = 0; PV_box < BoxArraySize(boxes); PV_box++)                  \
+      {                                                                             \
+        Box box = BoxArrayGetBox(boxes, PV_box);                                    \
+        int PV_ixl = box.lo[0];                                                     \
+        int PV_iyl = box.lo[1];                                                     \
+        int PV_izl = box.lo[2];                                                     \
+        int PV_ixu = box.up[0];                                                     \
+        int PV_iyu = box.up[1];                                                     \
+        int PV_izu = box.up[2];                                                     \
+                                                                                    \
+        if(PV_ixl <= PV_ixu && PV_iyl <= PV_iyu && PV_izl <= PV_izu)                \
+        {                                                                           \
+          int PV_nx = PV_ixu - PV_ixl + 1;                                          \
+          int PV_ny = PV_iyu - PV_iyl + 1;                                          \
+          int PV_nz = PV_izu - PV_izl + 1;                                          \
+                                                                                    \
+          const int _fdir0 = fdir[0];                                               \
+          const int _fdir1 = fdir[1];                                               \
+          const int _fdir2 = fdir[2];                                               \
+                                                                                    \
+          auto lambda_body =                                                        \
+            KOKKOS_LAMBDA(int i, int j, int k)                                      \
+            {                                                                       \
+              i += PV_ixl;                                                          \
+              j += PV_iyl;                                                          \
+              k += PV_izl;                                                          \
+                                                                                    \
+              /* Set surfflag for all cells in boxes regardless of loop limits */   \
+              surfflag[(k - iz_bxs) * ny_bxs * nx_bxs +                             \
+                (j - iy_bxs) * nx_bxs + (i - ix_bxs)] |= (1 << (2 + PV_f));         \
+                                                                                    \
+              /* Only evaluate loop body if the cell is within loop limits */       \
+              if(i >= ix && j >= iy && k >= iz &&                                   \
+                i < ix + nx && j < iy + ny && k < iz + nz)                          \
+              {                                                                     \
+                const int fdir[3] = {_fdir0, _fdir1, _fdir2};                       \
+                loop_body;                                                          \
+                (void)fdir;                                                         \
+              }                                                                     \
+            };                                                                      \
+                                                                                    \
+          using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+          MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{PV_nx, PV_ny, PV_nz}});        \
+          Kokkos::parallel_for(mdpolicy_3d, lambda_body);                           \
+          Kokkos::fence();                                                          \
+        }                                                                           \
+      }                                                                             \
+      GrGeomSolidCellFlagInitialized(grgeom) |= (1 << (2 + PV_f));                  \
+    }                                                                               \
+    else                                                                            \
+    {                                                                               \
+      int ixl_gpu = pfmax(ix, ix_bxs);                                              \
+      int iyl_gpu = pfmax(iy, iy_bxs);                                              \
+      int izl_gpu = pfmax(iz, iz_bxs);                                              \
+      int nx_gpu = pfmin((ix + nx - 1), BoxArrayMaxCell(boxes, 0)) - ixl_gpu + 1;   \
+      int ny_gpu = pfmin((iy + ny - 1), BoxArrayMaxCell(boxes, 1)) - iyl_gpu + 1;   \
+      int nz_gpu = pfmin((iz + nz - 1), BoxArrayMaxCell(boxes, 2)) - izl_gpu + 1;   \
+                                                                                    \
+      if(nx_gpu > 0 && ny_gpu > 0 && nz_gpu > 0)                                    \
+      {                                                                             \
+        char *surfflag = GrGeomSolidCellFlagData(grgeom);                           \
+                                                                                    \
+        const int _fdir0 = fdir[0];                                                 \
+        const int _fdir1 = fdir[1];                                                 \
+        const int _fdir2 = fdir[2];                                                 \
+                                                                                    \
+        auto lambda_body =                                                          \
+          KOKKOS_LAMBDA(int i, int j, int k)                                        \
+          {                                                                         \
+            i += ixl_gpu;                                                           \
+            j += iyl_gpu;                                                           \
+            k += izl_gpu;                                                           \
+                                                                                    \
+            if(surfflag[(k - iz_bxs) * ny_bxs * nx_bxs +                            \
+              (j - iy_bxs) * nx_bxs + (i - ix_bxs)] & (1 << (2 + PV_f)))            \
+            {                                                                       \
+              const int fdir[3] = {_fdir0, _fdir1, _fdir2};                         \
+              loop_body;                                                            \
+              (void)fdir;                                                           \
+            }                                                                       \
+          };                                                                        \
+                                                                                    \
+        using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+        MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{nx_gpu, ny_gpu, nz_gpu}});       \
+        Kokkos::parallel_for(mdpolicy_3d, lambda_body);                             \
+        Kokkos::fence();                                                            \
+      }                                                                             \
+    }                                                                               \
+  }                                                                                 \
+  (void)i;(void)j;(void)k;                                                          \
+}
+
+ /** Loop definition for Kokkos. */
+#define GrGeomPatchLoopBoxes_kokkos(i, j, k, fdir, grgeom, patch_num,               \
+  ix, iy, iz, nx, ny, nz, loop_body)                                                \
+{                                                                                   \
+  for (int PV_f = 0; PV_f < GrGeomOctreeNumFaces; PV_f++)                           \
+  {                                                                                 \
+    const int *fdir = FDIR_TABLE[PV_f];                                             \
+                                                                                    \
+    int n_prev = 0;                                                                 \
+    BoxArray* boxes = GrGeomSolidPatchBoxes(grgeom, patch_num, PV_f);               \
+    for (int PV_box = 0; PV_box < BoxArraySize(boxes); PV_box++)                    \
+    {                                                                               \
+      Box box = BoxArrayGetBox(boxes, PV_box);                                      \
+      /* find octree and region intersection */                                     \
+      int PV_ixl = pfmax(ix, box.lo[0]);                                            \
+      int PV_iyl = pfmax(iy, box.lo[1]);                                            \
+      int PV_izl = pfmax(iz, box.lo[2]);                                            \
+      int PV_ixu = pfmin((ix + nx - 1), box.up[0]);                                 \
+      int PV_iyu = pfmin((iy + ny - 1), box.up[1]);                                 \
+      int PV_izu = pfmin((iz + nz - 1), box.up[2]);                                 \
+                                                                                    \
+      if(PV_ixl <= PV_ixu && PV_iyl <= PV_iyu && PV_izl <= PV_izu)                  \
+      {                                                                             \
+        int PV_diff_x = PV_ixu - PV_ixl;                                            \
+        int PV_diff_y = PV_iyu - PV_iyl;                                            \
+        int PV_diff_z = PV_izu - PV_izl;                                            \
+                                                                                    \
+        int nx = PV_diff_x + 1;                                                     \
+        int ny = PV_diff_y + 1;                                                     \
+        int nz = PV_diff_z + 1;                                                     \
+                                                                                    \
+        const int _fdir0 = fdir[0];                                                 \
+        const int _fdir1 = fdir[1];                                                 \
+        const int _fdir2 = fdir[2];                                                 \
+                                                                                    \
+        auto lambda_body =                                                          \
+          KOKKOS_LAMBDA(int i, int j, int k)                                        \
+          {                                                                         \
+            const int fdir[3] = {_fdir0, _fdir1, _fdir2};                           \
+            int ival = n_prev + k * ny * nx + j * nx + i;                           \
+                                                                                    \
+            i += PV_ixl;                                                            \
+            j += PV_iyl;                                                            \
+            k += PV_izl;                                                            \
+                                                                                    \
+            loop_body;                                                              \
+            (void)fdir;                                                             \
+          };                                                                        \
+        n_prev += nz * ny * nx;                                                     \
+                                                                                    \
+        using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+        MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{nx, ny, nz}});                   \
+        Kokkos::parallel_for(mdpolicy_3d, lambda_body);                             \
+        Kokkos::fence();                                                            \
+      }                                                                             \
+    }                                                                               \
+  }                                                                                 \
+  (void)i;(void)j;(void)k;                                                          \
+}
+
+ /** Loop definition for Kokkos. */
+#define GrGeomPatchLoopBoxesNoFdir_kokkos(i, j, k, grgeom, patch_num, ovrlnd,       \
+  ix, iy, iz, nx, ny, nz, locals, setup,                                            \
+  f_left, f_right, f_down, f_up, f_back, f_front, finalize)                         \
+{                                                                                   \
+  int n_ival = 0;                                                                   \
+  for (int PV_f = 0; PV_f < GrGeomOctreeNumFaces; PV_f++)                           \
+  {                                                                                 \
+    BoxArray* boxes = GrGeomSolidPatchBoxes(grgeom, patch_num, PV_f);               \
+                                                                                    \
+    int ix_bxs = BoxArrayMinCell(boxes, 0);                                         \
+    int iy_bxs = BoxArrayMinCell(boxes, 1);                                         \
+    int iz_bxs = BoxArrayMinCell(boxes, 2);                                         \
+                                                                                    \
+    int nx_bxs = BoxArrayMaxCell(boxes, 0) - ix_bxs + 1;                            \
+    int ny_bxs = BoxArrayMaxCell(boxes, 1) - iy_bxs + 1;                            \
+    int nz_bxs = BoxArrayMaxCell(boxes, 2) - iz_bxs + 1;                            \
+                                                                                    \
+    int patch_loc;                                                                  \
+    if(ovrlnd)                                                                      \
+      patch_loc = GrGeomSolidNumPatches(grgeom) + patch_num;                        \
+    else                                                                            \
+      patch_loc = patch_num;                                                        \
+                                                                                    \
+    int *ptr_ival = GrGeomSolidCellIval(grgeom, patch_loc, PV_f);                   \
+    if(!(ptr_ival))                                                                 \
+    {                                                                               \
+      GrGeomSolidCellIval(grgeom, patch_loc, PV_f) =                                \
+        (int*)_talloc_device(sizeof(int) * nx_bxs * ny_bxs * nz_bxs);               \
+                                                                                    \
+      ptr_ival = GrGeomSolidCellIval(grgeom, patch_loc, PV_f);                      \
+      for (int idx = 0; idx < nx_bxs * ny_bxs * nz_bxs; idx++)                      \
+        ptr_ival[idx] = -1;                                                         \
+                                                                                    \
+      for (int PV_box = 0; PV_box < BoxArraySize(boxes); PV_box++)                  \
+      {                                                                             \
+        Box box = BoxArrayGetBox(boxes, PV_box);                                    \
+        int PV_ixl = pfmax(ix, box.lo[0]);                                          \
+        int PV_iyl = pfmax(iy, box.lo[1]);                                          \
+        int PV_izl = pfmax(iz, box.lo[2]);                                          \
+        int PV_ixu = pfmin((ix + nx - 1), box.up[0]);                               \
+        int PV_iyu = pfmin((iy + ny - 1), box.up[1]);                               \
+        int PV_izu = pfmin((iz + nz - 1), box.up[2]);                               \
+                                                                                    \
+        for (k = PV_izl; k <= PV_izu; k++)                                          \
+          for (j = PV_iyl; j <= PV_iyu; j++)                                        \
+            for (i = PV_ixl; i <= PV_ixu; i++)                                      \
+            {                                                                       \
+              UNPACK(locals);                                                       \
+              setup;                                                                \
+              switch(PV_f)                                                          \
+              {                                                                     \
+                f_left;                                                             \
+                f_right;                                                            \
+                f_down;                                                             \
+                f_up;                                                               \
+                f_back;                                                             \
+                f_front;                                                            \
+              }                                                                     \
+              finalize;                                                             \
+              ptr_ival[(k - iz_bxs) * ny_bxs * nx_bxs + (j - iy_bxs) *              \
+                nx_bxs + (i - ix_bxs)] = n_ival++;                                  \
+            }                                                                       \
+      }                                                                             \
+    }                                                                               \
+    else                                                                            \
+    {                                                                               \
+      int ixl_gpu = pfmax(ix, ix_bxs);                                              \
+      int iyl_gpu = pfmax(iy, iy_bxs);                                              \
+      int izl_gpu = pfmax(iz, iz_bxs);                                              \
+      int nx_gpu = pfmin((ix + nx - 1), BoxArrayMaxCell(boxes, 0)) - ixl_gpu + 1;   \
+      int ny_gpu = pfmin((iy + ny - 1), BoxArrayMaxCell(boxes, 1)) - iyl_gpu + 1;   \
+      int nz_gpu = pfmin((iz + nz - 1), BoxArrayMaxCell(boxes, 2)) - izl_gpu + 1;   \
+                                                                                    \
+      if(nx_gpu > 0 && ny_gpu > 0 && nz_gpu > 0)                                    \
+      {                                                                             \
+        auto lambda_body =                                                          \
+          KOKKOS_LAMBDA(int i, int j, int k)                                        \
+          {                                                                         \
+            i += ixl_gpu;                                                           \
+            j += iyl_gpu;                                                           \
+            k += izl_gpu;                                                           \
+                                                                                    \
+            int ival = ptr_ival[(k - iz_bxs) * ny_bxs * nx_bxs +                    \
+              (j - iy_bxs) * nx_bxs + (i - ix_bxs)];                                \
+            if(ival >= 0)                                                           \
+            {                                                                       \
+              UNPACK(locals);                                                       \
+              setup;                                                                \
+              switch(PV_f)                                                          \
+              {                                                                     \
+                f_left;                                                             \
+                f_right;                                                            \
+                f_down;                                                             \
+                f_up;                                                               \
+                f_back;                                                             \
+                f_front;                                                            \
+              }                                                                     \
+              finalize;                                                             \
+            }                                                                       \
+          };                                                                        \
+                                                                                    \
+        using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+        MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{nx_gpu, ny_gpu, nz_gpu}});       \
+        Kokkos::parallel_for(mdpolicy_3d, lambda_body);                             \
+        Kokkos::fence();                                                            \
+      }                                                                             \
+    }                                                                               \
+  }                                                                                 \
+  (void)i;(void)j;(void)k;                                                          \
+}
+
+ /** Loop definition for Kokkos. */
+#define GrGeomOctreeExteriorNodeLoop_kokkos(i, j, k, node, octree, level,           \
+  ix, iy, iz, nx, ny, nz, val_test, loop_body)                                      \
+{                                                                                   \
+  int PV_i, PV_j, PV_k, PV_l;                                                       \
+  int PV_ixl, PV_iyl, PV_izl, PV_ixu, PV_iyu, PV_izu;                               \
+                                                                                    \
+  PV_i = i;                                                                         \
+  PV_j = j;                                                                         \
+  PV_k = k;                                                                         \
+                                                                                    \
+  GrGeomOctreeExteriorLoop(PV_i, PV_j, PV_k, PV_l, node, octree, level, val_test,   \
+  {                                                                                 \
+    if ((PV_i >= ix) && (PV_i < (ix + nx)) &&                                       \
+        (PV_j >= iy) && (PV_j < (iy + ny)) &&                                       \
+        (PV_k >= iz) && (PV_k < (iz + nz)))                                         \
+    {                                                                               \
+      i = PV_i;                                                                     \
+      j = PV_j;                                                                     \
+      k = PV_k;                                                                     \
+      loop_body;                                                                    \
+    }                                                                               \
+  },                                                                                \
+  {                                                                                 \
+    /* find octree and region intersection */                                       \
+    PV_ixl = pfmax(ix, PV_i);                                                       \
+    PV_iyl = pfmax(iy, PV_j);                                                       \
+    PV_izl = pfmax(iz, PV_k);                                                       \
+    PV_ixu = pfmin((ix + nx), (PV_i + (int)PV_inc));                                \
+    PV_iyu = pfmin((iy + ny), (PV_j + (int)PV_inc));                                \
+    PV_izu = pfmin((iz + nz), (PV_k + (int)PV_inc));                                \
+                                                                                    \
+    if(PV_ixl < PV_ixu && PV_iyl < PV_iyu && PV_izl < PV_izu)                       \
+    {                                                                               \
+      const int PV_diff_x = PV_ixu - PV_ixl;                                        \
+      const int PV_diff_y = PV_iyu - PV_iyl;                                        \
+      const int PV_diff_z = PV_izu - PV_izl;                                        \
+                                                                                    \
+      auto lambda_body =                                                            \
+        KOKKOS_LAMBDA(int i, int j, int k)                                          \
+        {                                                                           \
+          i += PV_ixl;                                                              \
+          j += PV_iyl;                                                              \
+          k += PV_izl;                                                              \
+                                                                                    \
+          loop_body;                                                                \
+        };                                                                          \
+                                                                                    \
+      using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+      MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{PV_diff_x, PV_diff_y, PV_diff_z}});\
+      Kokkos::parallel_for(mdpolicy_3d, lambda_body);                               \
+      Kokkos::fence();                                                              \
+    }                                                                               \
+    i = PV_ixu;                                                                     \
+    j = PV_iyu;                                                                     \
+    k = PV_izu;                                                                     \
+  })                                                                                \
+  (void)i;(void)j;(void)k;                                                          \
+}
+
+ /** Loop definition for Kokkos. */
+#define GrGeomOutLoop_kokkos(i, j, k, grgeom, r,                                    \
+  ix, iy, iz, nx, ny, nz, body)                                                     \
+{                                                                                   \
+  if(nx > 0 && ny > 0 && nz > 0)                                                    \
+  {                                                                                 \
+    if(!(GrGeomSolidCellFlagInitialized(grgeom) & (1 << 1)))                        \
+    {                                                                               \
+      CheckCellFlagAllocation(grgeom, nx, ny, nz);                                  \
+      char *outflag = GrGeomSolidCellFlagData(grgeom);                              \
+                                                                                    \
+      GrGeomOctree  *PV_node;                                                       \
+      double PV_ref = pow(2.0, r);                                                  \
+                                                                                    \
+      i = GrGeomSolidOctreeIX(grgeom) * (int)PV_ref;                                \
+      j = GrGeomSolidOctreeIY(grgeom) * (int)PV_ref;                                \
+      k = GrGeomSolidOctreeIZ(grgeom) * (int)PV_ref;                                \
+      GrGeomOctreeExteriorNodeLoop(i, j, k, PV_node,                                \
+                                   GrGeomSolidData(grgeom),                         \
+                                   GrGeomSolidOctreeBGLevel(grgeom) + r,            \
+                                   ix, iy, iz, nx, ny, nz,                          \
+                                   TRUE,                                            \
+      {                                                                             \
+        body;                                                                       \
+        outflag[(k - iz) * ny * nx + (j - iy) * nx + (i - ix)] |= (1 << 1);         \
+      });                                                                           \
+      GrGeomSolidCellFlagInitialized(grgeom) |= (1 << 1);                           \
+    }                                                                               \
+    else                                                                            \
+    {                                                                               \
+      char *outflag = GrGeomSolidCellFlagData(grgeom);                              \
+      auto lambda_body =                                                            \
+        KOKKOS_LAMBDA(int i, int j, int k)                                          \
+        {                                                                           \
+          i += ix;                                                                  \
+          j += iy;                                                                  \
+          k += iz;                                                                  \
+                                                                                    \
+          if(outflag[(k - iz) * ny * nx + (j - iy) * nx + (i - ix)] & (1 << 1))     \
+          {                                                                         \
+            body;                                                                   \
+          }                                                                         \
+        };                                                                          \
+                                                                                    \
+      using MDPolicyType_3D = typename Kokkos::Experimental::MDRangePolicy<Kokkos::Experimental::Rank<3> >;\
+      MDPolicyType_3D mdpolicy_3d({{0, 0, 0}}, {{nx, ny, nz}});                     \
+      Kokkos::parallel_for(mdpolicy_3d, lambda_body);                               \
+      Kokkos::fence();                                                              \
+    }                                                                               \
+  }                                                                                 \
+  (void)i;(void)j;(void)k;                                                          \
+}
+
+}
+#endif // PF_KOKKOSLOOPS_H
diff --git a/pfsimulator/parflow_lib/pf_kokkosmalloc.h b/pfsimulator/parflow_lib/pf_kokkosmalloc.h
new file mode 100644
index 000000000..71504d7b4
--- /dev/null
+++ b/pfsimulator/parflow_lib/pf_kokkosmalloc.h
@@ -0,0 +1,49 @@
+/**********************************************************************
+ *
+ *  Please read the LICENSE file for the GNU Lesser General Public License.
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License (as published
+ *  by the Free Software Foundation) version 2.1 dated February 1999.
+ *
+ *  This program is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the IMPLIED WARRANTY OF
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the terms
+ *  and conditions of the GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Lesser General Public
+ *  License along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+ *  USA
+ ***********************************************************************/
+
+/** @file
+ * @brief Contains macro redefinitions for unified memory management.
+ */
+
+#ifndef PF_KOKKOSMALLOC_H
+#define PF_KOKKOSMALLOC_H
+
+#include "pf_devices.h"
+
+/*--------------------------------------------------------------------------
+ * Memory management macros for Kokkos
+ *--------------------------------------------------------------------------*/
+
+#define talloc_amps_kokkos(type, count) amps_TAlloc_managed(type, count)
+
+#define ctalloc_amps_kokkos(type, count) amps_CTAlloc_managed(type, count)
+
+#define tfree_amps_kokkos(ptr) amps_TFree_managed(ptr)
+
+#define talloc_kokkos(type, count) \
+  ((count) ? (type*)_talloc_device(sizeof(type) * (unsigned int)(count)) : NULL)
+
+#define ctalloc_kokkos(type, count) \
+  ((count) ? (type*)_ctalloc_device(sizeof(type) * (unsigned int)(count)) : NULL)
+
+#define tfree_kokkos(ptr) if (ptr) _tfree_device(ptr); else {}
+
+#define tmemcpy_kokkos(dest, src, bytes) kokkosMemCpy((char*)dest, (char*)src, (size_t)bytes);
+
+#endif // PF_KOKKOSMALLOC_H
diff --git a/pfsimulator/parflow_lib/pf_module.h b/pfsimulator/parflow_lib/pf_module.h
index dea20e281..bb11efa04 100644
--- a/pfsimulator/parflow_lib/pf_module.h
+++ b/pfsimulator/parflow_lib/pf_module.h
@@ -61,20 +61,6 @@ amps_ThreadLocalDcl(PFModule  *, global_ptr_this_pf_module);
 amps_ThreadLocalDcl(extern PFModule *, global_ptr_this_pf_module);
 #endif
 
-/*--------------------------------------------------------------------------
- * Define __device__ pointer for global_ptr_this_pf_module (CUDA)  
- *--------------------------------------------------------------------------*/
-
-#if (PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA) && defined(__CUDACC__)
-#ifdef PARFLOW_GLOBALS
-__device__ PFModule *dev_global_ptr_this_pf_module;
-#else
-/* This extern requires CUDA separate compilation, otherwise nvcc compiler 
-   treats the pointer as static variable for each compilation unit          */
-extern __device__ PFModule *dev_global_ptr_this_pf_module;
-#endif // PARFLOW_GLOBALS
-#endif // PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA && __CUDACC__
-
 /*--------------------------------------------------------------------------
  * Accessor macros
  *--------------------------------------------------------------------------*/
@@ -82,14 +68,8 @@ extern __device__ PFModule *dev_global_ptr_this_pf_module;
 #define PFModuleInstanceXtra(pf_module)      (pf_module->instance_xtra)
 #define PFModulePublicXtra(pf_module)        (pf_module->public_xtra)
 
-/* These accessor macros depends on compilation trajectory (host/device)    */
-#ifdef __CUDA_ARCH__
-#define dev_global_this_pf_module amps_ThreadLocal(dev_global_ptr_this_pf_module)
-#define ThisPFModule dev_global_this_pf_module
-#else
 #define global_this_pf_module amps_ThreadLocal(global_ptr_this_pf_module)
 #define ThisPFModule global_this_pf_module
-#endif
 
 /*--------------------------------------------------------------------------
  * PFModule interface macros
diff --git a/pfsimulator/parflow_lib/pf_omploops.h b/pfsimulator/parflow_lib/pf_omploops.h
index 0e986f009..f4167eabe 100644
--- a/pfsimulator/parflow_lib/pf_omploops.h
+++ b/pfsimulator/parflow_lib/pf_omploops.h
@@ -25,7 +25,7 @@
  * @brief OpenMP functions and loop definitions
  */
 
-#if PARFLOW_ACC_BACKEND != PARFLOW_BACKEND_OMP
+#ifndef PARFLOW_HAVE_OMP
 
 #define NO_OMP_PARALLEL
 
@@ -711,5 +711,5 @@ static const int FDIR_TABLE[6][3] = {
     }                                                                   \
   }
 
-#endif // PARFLOW_ACC_BACKEND != PARFLOW_BACKEND_OMP
+#endif // PARFLOW_HAVE_OMP
 #endif // _PF_OMPLOOPS_H_
diff --git a/pfsimulator/parflow_lib/problem_bc.h b/pfsimulator/parflow_lib/problem_bc.h
index 7d2ff4480..f05fcdda6 100644
--- a/pfsimulator/parflow_lib/problem_bc.h
+++ b/pfsimulator/parflow_lib/problem_bc.h
@@ -226,7 +226,7 @@ typedef struct {
 /**
  * @brief For use when a BCLoop should execute no matter what the patch type is
  */
-#define ALL -1
+#define BC_ALL -1
 
 /**
  * @brief For use when a statement body is unnecessary in ForPatchCellsPerFace().
@@ -399,7 +399,7 @@ ForPatchCellsPerFace(NotARealBCType,
                              finalize,                        \
                              after_loop)                      \
   {                                                           \
-    if ( ((bctype) == ALL) ||                                 \
+    if ( ((bctype) == BC_ALL) ||                              \
          ((bctype) == _GetCurrentPatch(loopvars)))            \
     {                                                         \
       before_loop;                                            \
@@ -427,7 +427,7 @@ ForPatchCellsPerFace(NotARealBCType,
                                       finalize,                       \
                                       after_loop)                     \
   {                                                                   \
-    if ( ((bctype) == ALL) ||                                         \
+    if ( ((bctype) == BC_ALL) ||                                      \
          ((bctype) == _GetCurrentPatch(loopvars)))                    \
     {                                                                 \
       before_loop;                                                    \
diff --git a/pfsimulator/parflow_lib/problem_bc_phase_saturation.c b/pfsimulator/parflow_lib/problem_bc_phase_saturation.c
index 860b60f3a..6cd6c291a 100644
--- a/pfsimulator/parflow_lib/problem_bc_phase_saturation.c
+++ b/pfsimulator/parflow_lib/problem_bc_phase_saturation.c
@@ -150,7 +150,7 @@ void          BCPhaseSaturation(
         case 0:
         {
           double constant;
-          ForPatchCellsPerFace(ALL,
+          ForPatchCellsPerFace(BC_ALL,
                                BeforeAllCells({
                                    dummy0 = (Type0*)(public_xtra->data[indx + ipatch]);
                                    constant = (dummy0->constant);
@@ -181,7 +181,7 @@ void          BCPhaseSaturation(
         {
           double height, lower, upper, dz2;
 
-          ForPatchCellsPerFace(ALL,
+          ForPatchCellsPerFace(BC_ALL,
                                BeforeAllCells({
                                    dummy1 = (Type1*)(public_xtra->data[indx + ipatch]);
                                    height = (dummy1->height);
@@ -238,7 +238,7 @@ void          BCPhaseSaturation(
           double dx2, dy2, dz2;
           double unitx, unity, line_min, line_length;
 
-          ForPatchCellsPerFace(ALL,
+          ForPatchCellsPerFace(BC_ALL,
                                BeforeAllCells({
                                    dummy2 = (Type2*)(public_xtra->data[indx + ipatch]);
                                    num_points = (dummy2->num_points);
diff --git a/pfsimulator/parflow_lib/problem_bc_pressure.c b/pfsimulator/parflow_lib/problem_bc_pressure.c
index 86f54d53d..7fd1723d7 100644
--- a/pfsimulator/parflow_lib/problem_bc_pressure.c
+++ b/pfsimulator/parflow_lib/problem_bc_pressure.c
@@ -208,7 +208,7 @@ BCStruct    *BCPressure(
 
             dz2 = SubgridDZ(subgrid) * 0.5;
 
-            ForPatchCellsPerFace(ALL,
+            ForPatchCellsPerFace(BC_ALL,
                                  BeforeAllCells(DoNothing),
                                  LoopVars(i, j, k, ival, bc_struct, ipatch, is),
                                  Locals(int ips, iel, iterations, phase;
@@ -417,7 +417,7 @@ BCStruct    *BCPressure(
             line_min = DirEquilPLinearXLower(interval_data) * unitx
                        + DirEquilPLinearYLower(interval_data) * unity;
 
-            ForPatchCellsPerFace(ALL,
+            ForPatchCellsPerFace(BC_ALL,
                                  BeforeAllCells(DoNothing),
                                  LoopVars(i, j, k, ival, bc_struct, ipatch, is),
                                  Locals(int ips, iterations;
@@ -634,7 +634,7 @@ BCStruct    *BCPressure(
 
             area = 0.0;
             patch_values_size = 0;
-            ForPatchCellsPerFace(ALL,
+            ForPatchCellsPerFace(BC_ALL,
                                  BeforeAllCells(DoNothing),
                                  LoopVars(i, j, k, ival, bc_struct, ipatch, is),
                                  Locals(int ips;),
@@ -816,7 +816,7 @@ BCStruct    *BCPressure(
             {
               case 1:  /* p = x */
               {
-                ForPatchCellsPerFace(ALL,
+                ForPatchCellsPerFace(BC_ALL,
                                      BeforeAllCells(DoNothing),
                                      LoopVars(i, j, k, ival, bc_struct, ipatch, is),
                                      Locals(double x;),
@@ -833,7 +833,7 @@ BCStruct    *BCPressure(
 
               case 2:  /* p = x + y + z */
               {
-                ForPatchCellsPerFace(ALL,
+                ForPatchCellsPerFace(BC_ALL,
                                      BeforeAllCells(DoNothing),
                                      LoopVars(i, j, k, ival, bc_struct, ipatch, is),
                                      Locals(int ips;
@@ -858,7 +858,7 @@ BCStruct    *BCPressure(
 
               case 3:  /* p = x^3y^2 + sinxy + 1*/
               {
-                ForPatchCellsPerFace(ALL,
+                ForPatchCellsPerFace(BC_ALL,
                                      BeforeAllCells(DoNothing),
                                      LoopVars(i, j, k, ival, bc_struct, ipatch, is),
                                      Locals(double x, y;),
@@ -881,7 +881,7 @@ BCStruct    *BCPressure(
 
               case 4:  /* p = x^3 y^4 + x^2 + sinxy cosy + 1 */
               {
-                ForPatchCellsPerFace(ALL,
+                ForPatchCellsPerFace(BC_ALL,
                                      BeforeAllCells(DoNothing),
                                      LoopVars(i, j, k, ival, bc_struct, ipatch, is),
                                      Locals(double x, y;),
@@ -904,7 +904,7 @@ BCStruct    *BCPressure(
 
               case 5:  /* p = xyzt +1 */
               {
-                ForPatchCellsPerFace(ALL,
+                ForPatchCellsPerFace(BC_ALL,
                                      BeforeAllCells(DoNothing),
                                      LoopVars(i, j, k, ival, bc_struct, ipatch, is),
                                      Locals(int ips; double x, y, z;),
@@ -928,7 +928,7 @@ BCStruct    *BCPressure(
 
               case 6:  /* p = xyzt +1 */
               {
-                ForPatchCellsPerFace(ALL,
+                ForPatchCellsPerFace(BC_ALL,
                                      BeforeAllCells(DoNothing),
                                      LoopVars(i, j, k, ival, bc_struct, ipatch, is),
                                      Locals(int ips; double x, y, z;),
diff --git a/pfsimulator/parflow_lib/region.c b/pfsimulator/parflow_lib/region.c
index 5ef608d6e..9ff0dbf02 100644
--- a/pfsimulator/parflow_lib/region.c
+++ b/pfsimulator/parflow_lib/region.c
@@ -58,7 +58,6 @@ Subregion  *NewSubregion(
   Subregion *new_subregion;
 
   new_subregion = talloc(Subregion, 1);
-  MemPrefetchDeviceToHost(new_subregion, sizeof(Subregion), 0);
 
   (new_subregion->ix) = ix;
   (new_subregion->iy) = iy;
diff --git a/pfsimulator/parflow_lib/richards_jacobian_eval.c b/pfsimulator/parflow_lib/richards_jacobian_eval.c
index faeadbe9b..36914b6aa 100644
--- a/pfsimulator/parflow_lib/richards_jacobian_eval.c
+++ b/pfsimulator/parflow_lib/richards_jacobian_eval.c
@@ -844,7 +844,7 @@ void    RichardsJacobianEval(
 
       ForBCStructNumPatches(ipatch, bc_struct)
       {
-        ForPatchCellsPerFace(ALL,
+        ForPatchCellsPerFace(BC_ALL,
                              BeforeAllCells(DoNothing),
                              LoopVars(i, j, k, ival, bc_struct, ipatch, is),
                              Locals(int ip, im;
diff --git a/pfsimulator/parflow_lib/solver.c b/pfsimulator/parflow_lib/solver.c
index e897a59ed..824fe23a1 100644
--- a/pfsimulator/parflow_lib/solver.c
+++ b/pfsimulator/parflow_lib/solver.c
@@ -121,13 +121,6 @@ NewSolver()
     NA_FreeNameArray(solver_na);
   }
 
-  /*-----------------------------------------------------------------------
-   * Copy Globals struct to GPU (does not copy all unneccessary data)
-   *-----------------------------------------------------------------------*/
-#if PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
-  CopyGlobalsToDevice();
-#endif
-
   switch (solver)
   {
     case 0:
diff --git a/pfsimulator/parflow_lib/vector.c b/pfsimulator/parflow_lib/vector.c
index 8faced4a7..fc3270677 100644
--- a/pfsimulator/parflow_lib/vector.c
+++ b/pfsimulator/parflow_lib/vector.c
@@ -924,7 +924,7 @@ void    InitVectorRandom(
     BoxLoopI1(i, j, k, ix, iy, iz, nx, ny, nz,
               iv, nx_v, ny_v, nz_v, 1, 1, 1,
     {
-#if PARFLOW_ACC_BACKEND == PARFLOW_BACKEND_CUDA
+#if defined(PARFLOW_HAVE_CUDA) || defined(PARFLOW_HAVE_KOKKOS)
       vp[iv] = dev_drand48();
 #else
       vp[iv] = drand48();


Test patch:
diff --git a/pfsimulator/amps/test/src/CMakeLists.txt b/pfsimulator/amps/test/src/CMakeLists.txt
index 2f74c7344..17c06a414 100644
--- a/pfsimulator/amps/test/src/CMakeLists.txt
+++ b/pfsimulator/amps/test/src/CMakeLists.txt
@@ -37,9 +37,12 @@ if(${PARFLOW_AMPS_LAYER} STREQUAL "cuda")
 endif(${PARFLOW_AMPS_LAYER} STREQUAL "cuda")
 
 # Add GPU test case (same as test15 with Unified Memory)
-if((${PARFLOW_HAVE_CUDA}) AND (${PARFLOW_AMPS_LAYER} STREQUAL "mpi1"))
+if(((${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_KOKKOS})) AND (${PARFLOW_AMPS_LAYER} STREQUAL "mpi1"))
   list(APPEND PARALLEL_TESTS test18)
-endif((${PARFLOW_HAVE_CUDA}) AND (${PARFLOW_AMPS_LAYER} STREQUAL "mpi1"))
+  if(NOT (${PARFLOW_HAVE_CUDA}))
+    list(REMOVE_ITEM PARALLEL_TESTS test15)  
+  endif(NOT (${PARFLOW_HAVE_CUDA}))
+endif(((${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_KOKKOS})) AND (${PARFLOW_AMPS_LAYER} STREQUAL "mpi1"))
 
 set(ALL_TESTS ${SEQUENTIAL_TESTS} ${PARALLEL_TESTS})
 list(REMOVE_DUPLICATES ALL_TESTS)
diff --git a/pfsimulator/amps/test/src/test18.c b/pfsimulator/amps/test/src/test18.c
index c8e5f124d..455dc81ba 100644
--- a/pfsimulator/amps/test/src/test18.c
+++ b/pfsimulator/amps/test/src/test18.c
@@ -87,7 +87,7 @@ char *argv[];
   stride[0] = size + 2;
   stride[1] = (size + 2);
 
-  a = amps_CTAlloc_managed(double, (size + 2) * (size + 2) * (size + 2));
+  a = amps_TAlloc_managed(double, (size + 2) * (size + 2) * (size + 2));
 
   if (me == 0)
   {
diff --git a/test/tcl/CMakeLists.txt b/test/tcl/CMakeLists.txt
index 232bbc64d..c1dfe7a5f 100644
--- a/test/tcl/CMakeLists.txt
+++ b/test/tcl/CMakeLists.txt
@@ -105,10 +105,10 @@ endforeach()
 
 foreach(inputfile ${PARALLEL_3DTOPO_TESTS})
   foreach(processor_topology "1 1 2" "1 2 1" "2 1 1" "2 2 2" "3 3 3" "1 1 4" "1 4 1" "4 1 1")
-    if(((${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_OMP})) AND (${processor_topology} STREQUAL "3 3 3"))
+    if(((${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_KOKKOS}) OR (${PARFLOW_HAVE_OMP})) AND (${processor_topology} STREQUAL "3 3 3"))
       #These tests should pass with CUDA or OMP, but they take a very long time
       continue()
-    endif(((${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_OMP})) AND (${processor_topology} STREQUAL "3 3 3"))
+    endif(((${PARFLOW_HAVE_CUDA}) OR (${PARFLOW_HAVE_KOKKOS}) OR (${PARFLOW_HAVE_OMP})) AND (${processor_topology} STREQUAL "3 3 3"))
     pf_add_parallel_test(${inputfile} ${processor_topology})
   endforeach()
 endforeach()