diff --git a/.DS_Store b/.DS_Store
deleted file mode 100644
index c6ec70fbb..000000000
Binary files a/.DS_Store and /dev/null differ
diff --git a/.gitignore b/.gitignore
index cb3b80d30..7d690806e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -54,3 +54,4 @@ doc/*
 !readme
 !figures/*
 
+.DS_Store
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2daf2a749..3bc788f49 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,102 +1,40 @@
-cmake_minimum_required(VERSION 3.8)
+cmake_minimum_required(VERSION 3.15)
 project(m3dc1 LANGUAGES Fortran C CXX)
 
-# macro(checkSetParam varName isRequired)
-#   if("${${varName}}" STREQUAL "") # if varName is empty
-#     set(${varName} $ENV{${varName}})
-#     if(("${${varName}}" STREQUAL "") AND ${isRequired}) # if ENV{varName} is empty also
-#       message(FATAL_ERROR
-#         "${varName} is not defined or in the environment, please specify ${varName}")
-#     endif()
-#   endif()
-# endmacro(checkSetParam varName)
-
-option(ENABLE_OPENMP "Enable OpenMP support" OFF)
-option(ENABLE_COMPLEX "Build complex version" OFF)
-option(ENABLE_3D "Build 3D version" OFF)
-option(ENABLE_ST "Build stellarator version " OFF)
-option(ENABLE_PARTICLE "Enable kinetic particle module" OFF)
+option(M3DC1_ENABLE_OPENMP "Enable OpenMP support" OFF)
+if(M3DC1_ENABLE_OPENMP)
+  find_package(OpenMP REQUIRED COMPONENTS Fortran)
+endif()
+option(M3DC1_ENABLE_COMPLEX "Build complex version" OFF)
+option(M3DC1_ENABLE_3D "Build 3D version" OFF)
+option(M3DC1_ENABLE_ST "Build stellarator version " OFF)
+option(M3DC1_ENABLE_PARTICLE "Enable kinetic particle module" OFF)
+option(M3DC1_ENABLE_GPU "Use GPU/OpenACC Fortran kernels" OFF)
+option(M3DC1_ENABLE_ADAS "Enable ADAS radiation data support" OFF)
+option(M3DC1_ENABLE_TRILINOS "Enable Trilinos support" OFF)
+option(M3DC1_ENABLE_PSPLINE "Enable PSPLINE support for readgato/read_jsolver" OFF)
+option(M3DC1_DOWNLOAD_ADAS_DATA "Download OpenADAS ADF11 data as part of the default build" OFF)
+set(M3DC1_ADAS_DATA_DIR "${CMAKE_BINARY_DIR}/adas/adf11"
+  CACHE PATH "Path to OpenADAS ADF11 data directory")
 
 cmake_host_system_information(RESULT HOST QUERY HOSTNAME)
 message(STATUS "Configuring to build on: ${HOST}")
 
-# # rhel/centos place x86-64 libs in lib64 instead of lib
-# # this *shouldn't* effect searching the lib directory
-# # first on other distros where the default install location
-# # is lib instead of lib64 (like debian derivatives)
-# set_property(GLOBAL PROPERTY FIND_LIBRARY_USE_LIB64_PATHS TRUE)
-
-# prefer static libraries unless we explicitly specify no to
-# set(CMAKE_FIND_LIBRARY_SUFFIXES ".a" ".so")
-# if(BUILD_SHARED_LIBS)
-#   set(CMAKE_FIND_LIBRARY_SUFFIXES ".so" ".a")
-# endif()
-
-# set the build type logically
- # default is Debug if .git directory exists
- # otherwise default is Release
-# list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
-# include(BuildType)
 
-# # add C++ compiler flags
-# if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Cray")
-#   # todo :  set cray-specific cxx flags
-#   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}")
-# elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "XL")
-#   # todo :  check for regular xl vs bg/q xl compiler
-#   #         currently assuming bg/q
-#   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -qflag=w -glanglvl=extended0x")
-# elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
-#   # todo :  set intel-specific cxx flags
-#   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}")
-# elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
-#   # todo :  version checks to see which c++ standards are available
-#   #         and warn/error if the version is incompatible
-#   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x -Wall -Wextra")
-# elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
-#   # todo :  version checks to see which c++ standards are available
-#   #         and warn/error if the version is incompatible
-#   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++1z -Wall -Wextra")
-# endif()
+list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_LIST_DIR}/cmake"
+                              "${CMAKE_CURRENT_LIST_DIR}/m3dc1_scorec/cmake")
+find_package(PETSc REQUIRED)
 
-# # add C compiler flags
-# if("${CMAKE_C_COMPILER_ID}" STREQUAL "Cray")
-#   # todo :  set cray-specific c flags
-#   set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS}")
-# elseif("${CMAKE_C_COMPILER_ID}" STREQUAL "XL")
-#   # todo :  set xl-specific c flags
-#   set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS}")
-# elseif("${CMAKE_C_COMPILER_ID}" STREQUAL "Intel")
-#   # todo :  set intel-specific c flags
-#   set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS}")
-# elseif("${CMAKE_C_COMPILER_ID}" STREQUAL "GNU")
-#   # todo :  version checks to see which c standards are available
-#   #         and warn/error if the version is incompatible
-#   set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wall -Wextra")
-# elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
-#   # todo :  version checks to see which c++ standards are available
-#   #         and warn/error if the version is incompatible
-#   set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -Wall -Wextra")
-# endif()
+# Automatically set exe linker flags from PETSc pkg-config data.
+# Construct -L<dir> flags from PETSC_LIBRARY_DIRS.
+string(REPLACE ";" " -L" _m3dc1_petsc_libdirs_flags "${PETSC_LIBRARY_DIRS}")
+if(_m3dc1_petsc_libdirs_flags)
+  set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -L${_m3dc1_petsc_libdirs_flags}")
+endif()
+# Append PETSc's extra LDFLAGS (which include rpath from pkg-config).
+if(PETSC_LDFLAGS_OTHER)
+  set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${PETSC_LDFLAGS_OTHER}")
+endif()
 
-# # add Fortran compiler flags
-# if("${CMAKE_Fortran_COMPILER_ID}" STREQUAL "Cray")
-#   # todo : set cray-specific fortran flags
-#   set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}")
-# elseif("${CMAKE_Fortran_COMPILER_ID}" STREQUAL "XL")
-#   # todo : set xl-specific fortran flags
-#   set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}")
-# elseif("${CMAKE_Fortran_COMPILER_ID}" STREQUAL "Intel")
-#   # todo : set intel-specific fortran flags
-#   set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -g noarg_temp_created")
-# elseif("${CMAKE_Fortran_COMPILER_ID}" STREQUAL "GNU")
-#   set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-real-8 -Wall -cpp")
-# elseif("${CMAKE_Fortran_COMPILER_ID}" STREQUAL "Clang")
-#   set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-real-8 -Wall -cpp")
-# endif()
-
-message(STATUS "Configuring m3dc1_scorec")
-add_subdirectory(m3dc1_scorec EXCLUDE_FROM_ALL)
-message(STATUS "Configuring unstructured")
+add_subdirectory(m3dc1_scorec)
 add_subdirectory(unstructured)
-
diff --git a/cmake/DownloadOpenADAS.cmake b/cmake/DownloadOpenADAS.cmake
new file mode 100644
index 000000000..9b966bc27
--- /dev/null
+++ b/cmake/DownloadOpenADAS.cmake
@@ -0,0 +1,82 @@
+# Download OpenADAS ADF11 data and reader sources
+# Called from unstructured/CMakeLists.txt with:
+#   -DOUT_DIR=<path>     : output directory (CMAKE_BINARY_DIR/adas)
+#   -DADF11_DIR=<path>   : where ADF11 data files are placed
+
+set(OA_DATA_URL "https://open.adas.ac.uk/download/adf11")
+set(OA_SOURCE_DIR "${OUT_DIR}/source")
+
+# ---- Download and extract reader source code ----
+message(STATUS "Downloading OpenADAS reader sources...")
+file(DOWNLOAD
+  "https://open.adas.ac.uk/code/xxdata_11.tar.gz"
+  "${OUT_DIR}/xxdata_11.tar.gz"
+  STATUS _dl_status
+  SHOW_PROGRESS
+)
+list(GET _dl_status 0 _dl_code)
+if(NOT _dl_code EQUAL 0)
+  message(FATAL_ERROR "Failed to download xxdata_11.tar.gz")
+endif()
+
+# Extract into a temporary directory, then move .for → .f to source/
+set(_extract_dir "${OUT_DIR}/_extract_tmp")
+file(REMOVE_RECURSE "${_extract_dir}")
+file(MAKE_DIRECTORY "${_extract_dir}")
+file(MAKE_DIRECTORY "${OA_SOURCE_DIR}")
+
+execute_process(
+  COMMAND ${CMAKE_COMMAND} -E tar xzf "${OUT_DIR}/xxdata_11.tar.gz"
+  WORKING_DIRECTORY "${_extract_dir}"
+  RESULT_VARIABLE _tar_res
+)
+if(NOT _tar_res EQUAL 0)
+  message(FATAL_ERROR "Failed to extract xxdata_11.tar.gz")
+endif()
+
+# Find all .for files (recursively, handles possible subdirectories)
+# and copy/rename them to OA_SOURCE_DIR as .f
+file(GLOB_RECURSE _for_files "${_extract_dir}/*.for")
+foreach(_f IN LISTS _for_files)
+  get_filename_component(_name "${_f}" NAME_WE)
+  file(COPY "${_f}" DESTINATION "${OA_SOURCE_DIR}")
+  file(RENAME "${OA_SOURCE_DIR}/${_name}.for" "${OA_SOURCE_DIR}/${_name}.f")
+endforeach()
+
+# Cleanup
+file(REMOVE_RECURSE "${_extract_dir}")
+file(REMOVE "${OUT_DIR}/xxdata_11.tar.gz")
+message(STATUS "OpenADAS reader sources extracted to ${OA_SOURCE_DIR}")
+
+# ---- Download ADF11 data files ----
+set(_adf11_files
+  "scd85/scd85_ar.dat"
+  "scd89/scd89_b.dat"
+  "plt89/plt89_b.dat"
+  "plt89/plt89_ar.dat"
+  "scd96/scd96_he.dat"
+  "scd96/scd96_be.dat"
+  "scd96/scd96_c.dat"
+  "scd96/scd96_ne.dat"
+  "plt96/plt96_he.dat"
+  "plt96/plt96_be.dat"
+  "plt96/plt96_c.dat"
+  "plt96/plt96_ne.dat"
+)
+
+foreach(_f IN LISTS _adf11_files)
+  get_filename_component(_subdir "${_f}" DIRECTORY)
+  file(MAKE_DIRECTORY "${ADF11_DIR}/${_subdir}")
+  message(STATUS "Downloading ${_f}...")
+  file(DOWNLOAD
+    "${OA_DATA_URL}/${_f}"
+    "${ADF11_DIR}/${_f}"
+    STATUS _dl_status
+  )
+  list(GET _dl_status 0 _dl_code)
+  if(NOT _dl_code EQUAL 0)
+    message(FATAL_ERROR "Failed to download ${_f}")
+  endif()
+endforeach()
+
+message(STATUS "OpenADAS download complete")
\ No newline at end of file
diff --git a/m3dc1_scorec/CMakeLists.txt b/m3dc1_scorec/CMakeLists.txt
index b1b717a0d..b3539b6d6 100644
--- a/m3dc1_scorec/CMakeLists.txt
+++ b/m3dc1_scorec/CMakeLists.txt
@@ -1,62 +1,3 @@
-#This is the top M3DC1_SCOREC CMakeList File for the Build
-
-#Setting Version Number, Project Name
-cmake_minimum_required (VERSION 2.8)
-project (m3dc1_scorec)
-enable_language (Fortran)
-
-# make sure that the default is a RELEASE
-if (NOT CMAKE_BUILD_TYPE)
-  set (CMAKE_BUILD_TYPE RELEASE CACHE STRING
-      "Choose the type of build, options are: None Debug Release."
-      FORCE)
-endif (NOT CMAKE_BUILD_TYPE)
-
-# default installation
-#get_filename_component (default_prefix ".." ABSOLUTE)
-#set (CMAKE_INSTALL_PREFIX ${default_prefix} CACHE STRING
-#      "Choose the installation directory; by default it installs in the NORMA directory."
-#      FORCE)
-
-# FFLAGS depend on the compiler
-get_filename_component (Fortran_COMPILER_NAME ${CMAKE_Fortran_COMPILER} NAME)
-message (Fortran_COMPILER_NAME = ${Fortran_COMPILER_NAME})
-
-#unless building shared libs, then select static libs 
-# if both static and shared libs are available 
-set(CMAKE_FIND_LIBRARY_SUFFIXES ".a" ".so") 
-if(BUILD_SHARED_LIBS)
-  set(CMAKE_FIND_LIBRARY_SUFFIXES ".so" ".a")
-endif()
-
-#Settings options for testing
-enable_testing()
-include(CTest)
-#This will be set to ON by the CTest driver script (and only by that)
-option(ENABLE_TESTING "Build for CTest" OFF)
-set(MPIRUN "mpirun"
-    CACHE string 
-    "the mpirun or srun executable")
-set(MPIRUN_PROCFLAG "-np"
-    CACHE string 
-    "the command line flag to give process count to MPIRUN")
-
-#Doxygen generation system
-find_package(Doxygen)
-if(DOXYGEN_FOUND)
-configure_file(${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in
-               ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile @ONLY)
-add_custom_target(doc
-${DOXYGEN_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile
-WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-COMMENT "Generating API documentation with Doxygen" VERBATIM
-)
-endif(DOXYGEN_FOUND)
-
-configure_file(
-    "${CMAKE_CURRENT_SOURCE_DIR}/m3dc1_scorecConfig.cmake.in"
-    "${CMAKE_BINARY_DIR}/m3dc1_scorecConfig.cmake")
-
 #Source and header files
 set(SOURCES
   api/m3dc1_scorec.cc
@@ -97,143 +38,37 @@ set(HEADERS
   #include/ReducedQuinticExplicit.h
   include/slntransferUtil.h
   include/Expression.h
-
   include/ReducedQuintic.h)
 
-set(CMAKE_MODULE_PATH 
-   ${CMAKE_MODULE_PATH} 
-   "${CMAKE_CURRENT_SOURCE_DIR}/cmake/")
-
-include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include)
-include_directories(${CMAKE_CURRENT_SOURCE_DIR}/api)
-set(TARGET_LIB_NAME m3dc1_scorec)
-
-find_package(Scorec QUIET REQUIRED)
-include_directories(${SCOREC_INCLUDE_DIRS})
-set(DEP_LIBS ${DEP_LIBS} ${SCOREC_LIBRARIES} )
 
-if (ENABLE_SIMMETRIX)
-set(ENABLE_PETSC "OFF")
-set(ENABLE_TRILINOS "OFF")
-find_package(Simmetrix QUIET REQUIRED)
-include_directories(${SIMMETRIX_INCLUDE_DIRS})
-set(DEP_LIBS ${DEP_LIBS} ${SIMMETRIX_LIBRARIES} )
-set(TARGET_LIB_NAME ${TARGET_LIB_NAME}_simmetrix)
-find_package(Lapack QUIET REQUIRED)
-set(DEP_LIBS ${DEP_LIBS} ${LAPACK_LIBRARIES} )
-  add_executable(m3dc1_meshgen test/meshgen/m3dc1_meshgen.cc)
-  target_link_libraries(m3dc1_meshgen ${DEP_LIBS})
-  target_link_libraries(m3dc1_meshgen ${TARGET_LIB_NAME})
-  INSTALL(TARGETS m3dc1_meshgen RUNTIME DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
-  add_executable(convert_sim_sms test/meshgen/convert_sim_sms.cc)
-  target_link_libraries(convert_sim_sms ${DEP_LIBS})
-  target_link_libraries(convert_sim_sms ${TARGET_LIB_NAME})
-  INSTALL(TARGETS convert_sim_sms RUNTIME DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
-  add_executable(convert_polar test/convert_polar/main.cc)
-  target_link_libraries(convert_polar ${DEP_LIBS})
-  target_link_libraries(convert_polar ${TARGET_LIB_NAME})
-  INSTALL(TARGETS convert_polar RUNTIME DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
-elseif (ENABLE_TRILINOS)
-find_package(Trilinos QUIET REQUIRED)
-find_package(Boost QUIET REQUIRED)
-include_directories(${TRILINOS_INCLUDE_DIRS})
-set(DEP_LIBS ${DEP_LIBS} ${TRILINOS_LIBRARIES})
-set(DEP_LIBS ${DEP_LIBS} ${BOOST_LIBRARIES})
-set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DM3DC1_TRILINOS")
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DM3DC1_TRILINOS")
-set(TARGET_LIB_NAME ${TARGET_LIB_NAME}_trilinos)
-else()
-  set(ENABLE_PETSC "ON")
-endif()
+find_package(SCOREC CONFIG REQUIRED)
+find_package(LAPACK REQUIRED)
+find_package(MPI REQUIRED COMPONENTS C)
 
-option(ENABLE_ZOLTAN "Enable Zoltan interface [ON|OFF]" ON)
-message(STATUS "ENABLE_ZOLTAN: " ${ENABLE_ZOLTAN})
-
-if (ENABLE_ZOLTAN)
-find_library(ZOLTAN_LIBRARY zoltan)
-if (NOT EXISTS "${ZOLTAN_LIBRARY}")
-  message(FATAL ERROR "ZOLTAN library not found")
-endif()
-else()
-set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DDISABLE_ZOLTAN")
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DDISABLE_ZOLTAN")
-endif()
-
-if (ENABLE_PETSC)
-set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DM3DC1_PETSC")
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DM3DC1_PETSC")
-#find_package(Hdf5 QUIET REQUIRED)
-#set(DEP_LIBS ${DEP_LIBS} ${HDF5_LIBRARIES} )
-find_package(Petsc QUIET REQUIRED)
-include_directories(${PETSC_INCLUDE_DIRS})
-set(DEP_LIBS ${DEP_LIBS} ${PETSC_LIBRARIES} )
-endif()
-
-if(ENABLE_COMPLEX)
-set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DPETSC_USE_COMPLEX")
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DPETSC_USE_COMPLEX")
-set(TARGET_LIB_NAME ${TARGET_LIB_NAME}_complex)
-endif()
-
-# if(ENABLE_OPENMP)
-#   if(${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU")
-#     set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp")
-#     set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
-#   elseif (${CMAKE_CXX_COMPILER_ID} STREQUAL "Intel")
-#     set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -qopenmp")
-#     set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -qopenmp")
-#   else()
-#     message(WARNING "Unexpected compiler type ${CMAKE_CXX_COMPILER_ID}")
-#   endif()
-#   set(TARGET_LIB_NAME ${TARGET_LIB_NAME}_openmp)
-# endif()
-
-add_library(${TARGET_LIB_NAME} ${SOURCES})
-target_link_libraries(${TARGET_LIB_NAME} ${DEP_LIBS})
-
-INSTALL(FILES ${HEADERS} DESTINATION include)
-INSTALL(TARGETS ${TARGET_LIB_NAME}
-      LIBRARY DESTINATION lib
-      ARCHIVE DESTINATION lib)
-
-if (ENABLE_TESTING)
-if (ENABLE_PETSC)
-  add_executable(solver test/solver/main.cc)
-  target_link_libraries(solver ${DEP_LIBS})
-  target_link_libraries(solver ${TARGET_LIB_NAME})
-  INSTALL(TARGETS solver RUNTIME DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
-endif()
-
-if (ENABLE_TRILINOS)
-    add_executable(epetra test/epetra/main.cc)
-    target_link_libraries(epetra ${DEP_LIBS})
-    target_link_libraries(epetra ${TARGET_LIB_NAME})
-    INSTALL(TARGETS epetra RUNTIME DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
-endif()
-add_executable(show_meshcount test/create_smb/show_meshcount.cc)
-target_link_libraries(show_meshcount ${DEP_LIBS})
-target_link_libraries(show_meshcount ${TARGET_LIB_NAME})
-INSTALL(TARGETS show_meshcount RUNTIME DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
-endif()
-
-#add_executable(ghosting_test test/ghosting/main.cc)
-#target_link_libraries(ghosting_test ${DEP_LIBS})
-#target_link_libraries(ghosting_test ${TARGET_LIB_NAME})
-#INSTALL(TARGETS ghosting_test RUNTIME DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
-
-#binary distribution package
-set(CPACK_GENERATOR "TGZ")
-set(CPACK_PACKAGE_VERSION "1.0.1")
-include(CPack)
-
-message("ENV_FLAGS = $ENV{FFLAGS}")
-message("CMAKE_C_FLAGS = ${CMAKE_C_FLAGS}")
-message("CMAKE_CXX_FLAGS = ${CMAKE_CXX_FLAGS}")
-message("CMAKE_Fortran_COMPILER_INIT = ${CMAKE_Fortran_COMPILER_INIT}")
-message("CMAKE_Fortran_COMPILER_FULLPATH =
-${CMAKE_Fortran_COMPILER_FULLPATH}")
-message("CMAKE_Fortran_COMPILER = ${CMAKE_Fortran_COMPILER}")
-message("CMAKE_Fortran_FLAGS = ${CMAKE_Fortran_FLAGS}")
-message("CMAKE_INSTALL_PREFIX = ${CMAKE_INSTALL_PREFIX}")
+add_library(m3dc1_scorec ${SOURCES})
+add_library(m3dc1::scorec ALIAS m3dc1_scorec)
+target_link_libraries(m3dc1_scorec PUBLIC SCOREC::core PETSc::PETSc LAPACK::LAPACK MPI::MPI_C)
+target_include_directories(m3dc1_scorec PUBLIC
+  "${CMAKE_CURRENT_SOURCE_DIR}/api"
+  "${CMAKE_CURRENT_SOURCE_DIR}/include"
+)
+target_compile_definitions(m3dc1_scorec PRIVATE M3DC1_PETSC)
+
+# --- Utility executables (SCOREC tools) ---
+# Utilities that only depend on SCOREC core
+set(M3DC1_UTILITY_CORE_SOURCES
+  utility/check_smb.cc
+  utility/create_smb.cc
+  utility/split_smb.cc
+  utility/vtk_order.cc
+)
 
+foreach(_util_src IN LISTS M3DC1_UTILITY_CORE_SOURCES)
+  get_filename_component(_util_name ${_util_src} NAME_WE)
+  add_executable(${_util_name} ${_util_src})
+  target_link_libraries(${_util_name} PRIVATE SCOREC::core MPI::MPI_C)
+endforeach()
 
+# show_meshcount depends on the m3dc1_scorec library (m3dc1_mesh, m3dc1_model, etc.)
+add_executable(show_meshcount utility/show_meshcount.cc)
+target_link_libraries(show_meshcount PRIVATE m3dc1_scorec)
diff --git a/m3dc1_scorec/cmake/FindPETSc.cmake b/m3dc1_scorec/cmake/FindPETSc.cmake
new file mode 100644
index 000000000..24b1e0f70
--- /dev/null
+++ b/m3dc1_scorec/cmake/FindPETSc.cmake
@@ -0,0 +1,54 @@
+find_package(PkgConfig REQUIRED QUIET)
+if(DEFINED PETSC_DIR)
+  if(DEFINED PETSC_ARCH)
+    set(ENV{PKG_CONFIG_PATH} "${PETSC_DIR}/${PETSC_ARCH}/lib/pkgconfig:$ENV{PKG_CONFIG_PATH}")
+  else()
+    set(ENV{PKG_CONFIG_PATH} "${PETSC_DIR}/lib/pkgconfig:$ENV{PKG_CONFIG_PATH}")
+  endif()
+endif()
+
+# we give an internal name _petsc
+# so we can fill up the PETSC_VARIABLE based
+# on static or not
+pkg_check_modules(_petsc PETSc QUIET)
+
+if(_petsc_FOUND AND _petsc_VERSION)
+  set(PETSC_VERSION ${_petsc_VERSION})
+endif()
+
+# note, there are a number of additional properties that
+# can be extracted / set on a target. We just do the basic
+# set for now. See: https://cmake.org/cmake/help/latest/module/FindPkgConfig.html
+if(PETSC_LINK_STATIC)
+  set(PETSC_LIBRARIES ${_petsc_STATIC_LIBRARIES})
+  set(PETSC_INCLUDE_DIRS ${_petsc_STATIC_INCLUDE_DIRS})
+  set(PETSC_LIBRARY_DIRS ${_petsc_STATIC_LIBRARY_DIRS})
+  set(PETSC_LDFLAGS ${_petsc_STATIC_LDFLAGS})
+  set(PETSC_LDFLAGS_OTHER ${_petsc_STATIC_LDFLAGS_OTHER})
+elseif(_petsc_FOUND)
+  set(PETSC_LIBRARIES ${_petsc_LIBRARIES})
+  set(PETSC_INCLUDE_DIRS ${_petsc_INCLUDE_DIRS})
+  set(PETSC_LIBRARY_DIRS ${_petsc_LIBRARY_DIRS})
+  set(PETSC_LDFLAGS ${_petsc_LDFLAGS})
+  set(PETSC_LDFLAGS_OTHER ${_petsc_LDFLAGS_OTHER})
+endif()
+
+if(NOT TARGET PETSc::PETSc)
+  add_library(PETSc::PETSc INTERFACE IMPORTED GLOBAL)
+  set_target_properties(PETSc::PETSc PROPERTIES INTERFACE_LINK_LIBRARIES "${PETSC_LIBRARIES}"
+                                     INTERFACE_INCLUDE_DIRECTORIES "${PETSC_INCLUDE_DIRS}"
+                                     INTERFACE_LINK_DIRECTORIES "${PETSC_LIBRARY_DIRS}"
+                                     INTERFACE_LINK_OPTIONS "${PETSC_LDFLAGS_OTHER}")
+endif()
+
+
+
+include(FindPackageHandleStandardArgs)
+# TODO consider adding version check logic
+find_package_handle_standard_args(PETSc
+  REQUIRED_VARS PETSC_LIBRARIES PETSC_INCLUDE_DIRS
+  VERSION_VAR PETSC_VERSION)
+
+#message(STATUS "STATIC_LIBRARIES=${_petsc_STATIC_LIBRARIES}")
+#message(STATUS "STATIC_LDFLAGS=${_petsc_STATIC_LDFLAGS}")
+#message(STATUS "STATIC_LDFLAGS_OTHER=${_petsc_STATIC_LDFLAGS_OTHER}")
\ No newline at end of file
diff --git a/m3dc1_scorec/cmake/FindScorec.cmake b/m3dc1_scorec/cmake/FindScorec.cmake
deleted file mode 100644
index 9ac49bc75..000000000
--- a/m3dc1_scorec/cmake/FindScorec.cmake
+++ /dev/null
@@ -1,114 +0,0 @@
-# - Try to find SCOREC PUMI libraries
-# Once done this will define
-#  SCOREC_FOUND - System has SCOREC
-#  SCOREC_INCLUDE_DIRS - The SCOREC include directories
-#  SCOREC_LIBRARIES - The libraries needed to use SCOREC
-#  SCOREC_DEFINITIONS - Compiler switches required for using SCOREC
-#
-# This implementation assumes a SCOREC install has the following structure
-# VERSION/
-#         include/*.h
-#         lib/*.a
-
-macro(scorecLibCheck libs isRequired)
-  foreach(lib ${libs}) 
-    unset(scoreclib CACHE)
-    find_library(scoreclib "${lib}" PATHS ${SCOREC_LIB_DIR})
-    if(scoreclib MATCHES "^scoreclib-NOTFOUND$")
-      if(${isRequired})
-        message(FATAL_ERROR "SCOREC library ${lib} not found in ${SCOREC_LIB_DIR}")
-      else()
-        message("SCOREC library ${lib} not found in ${SCOREC_LIB_DIR}")
-      endif()
-    else()
-      set("SCOREC_${lib}_FOUND" TRUE CACHE INTERNAL "SCOREC library present")
-      set(SCOREC_LIBS ${SCOREC_LIBS} ${scoreclib})
-    endif()
-  endforeach()
-endmacro(scorecLibCheck)
-
-if(CMAKE_CXX_COMPILER_ID MATCHES "PGI")
-  set(ENABLE_ZOLTAN OFF)
-endif()
- 
-find_library(METIS_LIBRARY metis)
-if (NOT EXISTS "${METIS_LIBRARY}")
-  message(FATAL ERROR "METIS library not found")
-endif()
-
-set(SCOREC_LIBS "")
-if(ENABLE_SIMMETRIX)
-set(SCOREC_LIB_NAMES
-  ma
-  mds
-  apf_sim
-  apf_zoltan
-  apf
-  parma
-  gmi
-  gmi_sim
-  mth
-  pcu
-  sam
-  spr
-  crv
-  lion
-  ph 
-  )
-else()
-set(SCOREC_LIB_NAMES
-crv
-ph
-sam
-spr
-ma
-apf_zoltan
-parma
-mds
-apf
-lion
-mth
-gmi
-pcu
-  )
-endif()
-scorecLibCheck("${SCOREC_LIB_NAMES}" TRUE)
-
-find_path(SCOREC_INCLUDE_DIR 
-  NAMES apf.h PCU.h ma.h 
-  PATHS ${SCOREC_INCLUDE_DIR})
-if(NOT EXISTS "${SCOREC_INCLUDE_DIR}")
-  message(FATAL_ERROR "SCOREC include dir not found")
-endif()
-
-string(REGEX REPLACE 
-  "/include$" "" 
-  SCOREC_INSTALL_DIR
-  "${SCOREC_INCLUDE_DIR}")
-
-set(SCOREC_LIBRARIES ${SCOREC_LIBS} ${ZOLTAN_LIBRARY} ${PARMETIS_LIBRARY} ${METIS_LIBRARY})
-set(SCOREC_INCLUDE_DIRS ${SCOREC_INCLUDE_DIR})
-
-include(FindPackageHandleStandardArgs)
-# handle the QUIETLY and REQUIRED arguments and set PARMETIS_FOUND to TRUE
-# if all listed variables are TRUE
-find_package_handle_standard_args(SCOREC  DEFAULT_MSG
-                                  SCOREC_LIBS SCOREC_INCLUDE_DIR)
-
-mark_as_advanced(SCOREC_INCLUDE_DIR SCOREC_LIBS ${ZOLTAN_LIBRARY} ${PARMETIS_LIBRARY} ${METIS_LIBRARY})
-
-set(SCOREC_LINK_LIBS "")
-foreach(lib ${SCOREC_LIB_NAMES})
-  set(SCOREC_LINK_LIBS "${SCOREC_LINK_LIBS} -l${lib}")
-endforeach()
-
-#pkgconfig  
-set(prefix "${SCOREC_INSTALL_DIR}")
-set(includedir "${SCOREC_INCLUDE_DIR}")
-configure_file(
-  "cmake/libScorec.pc.in"
-  "${CMAKE_BINARY_DIR}/libScorec.pc"
-  @ONLY)
-
-INSTALL(FILES "${CMAKE_BINARY_DIR}/libScorec.pc" DESTINATION lib/pkgconfig)
-
diff --git a/m3dc1_scorec/include/ReducedQuinticExplicit.h b/m3dc1_scorec/include/ReducedQuinticExplicit.h
new file mode 100644
index 000000000..78ca38db1
--- /dev/null
+++ b/m3dc1_scorec/include/ReducedQuinticExplicit.h
@@ -0,0 +1,32 @@
+/****************************************************************************** 
+
+  (c) 2005-2017 Scientific Computation Research Center, 
+      Rensselaer Polytechnic Institute. All rights reserved.
+  
+  This work is open source software, licensed under the terms of the
+  BSD license as described in the LICENSE file in the top-level directory.
+ 
+*******************************************************************************/
+#ifndef REDUCEDQUINTICEXPLICIT_H
+#define REDUCEDQUINTICEXPLICIT_H
+#include "ReducedQuintic.h"
+// see Barnhill 1981
+class ReducedQuinticExplicit : public ReducedQuintic
+{
+public:
+  virtual void setCoord(double coord[3][2]);
+  // dofs in glb coord
+  virtual void setDofs(double dofs_p[18]);
+  virtual void print();
+  virtual void eval_l( double coord[2], double res[6] );
+  virtual void eval_g( double coord[2], double res[6] );
+  virtual void loc2glb( double coord[2]);
+  virtual void glb2loc( double coord[2]);
+  virtual void loc2glbDofs(double dofs_p[18]);
+  virtual void glb2locDofs(double dofs_p[18]);
+  virtual double getJacobi();
+private:
+  double dofs[18];
+}; 
+
+#endif
diff --git a/m3dc1_scorec/src/ReducedQuinticExplicit.cc b/m3dc1_scorec/src/ReducedQuinticExplicit.cc
new file mode 100644
index 000000000..c66b2cde9
--- /dev/null
+++ b/m3dc1_scorec/src/ReducedQuinticExplicit.cc
@@ -0,0 +1,86 @@
+/****************************************************************************** 
+
+  (c) 2005-2017 Scientific Computation Research Center, 
+      Rensselaer Polytechnic Institute. All rights reserved.
+  
+  This work is open source software, licensed under the terms of the
+  BSD license as described in the LICENSE file in the top-level directory.
+ 
+*******************************************************************************/
+#include "ReducedQuinticExplicit.h"
+extern "C"
+{
+  void set_element_coord_i_ ( int*, double*);
+  void eval_field_i_( double *, double*, double*, double*);
+  void loc_xyz_2_glb_i_ ( double*, double*);
+  void glb_xyz_2_loc_i_ ( double*, double*);
+  void glb_dof_2_loc_i_ ( double*, double*);
+  void loc_dof_2_glb_i_ ( double*, double*);
+  void get_jacobi_i_( double *);
+}
+double ReducedQuinticExplicit:: getJacobi()
+{
+  double jac;
+  get_jacobi_i_(&jac);
+  return jac;
+}
+void ReducedQuinticExplicit:: setCoord(double coord[3][2])
+{
+  double coord_f[2][3];
+  // convert to fortran convention
+  for(int i=0; i<3; i++)
+  {
+    for( int j=0; j<2; j++)
+      coord_f[j][i]=coord[i][j];
+  }
+  int dummy=0;
+  set_element_coord_i_(&dummy,(double*)coord_f);
+}
+void ReducedQuinticExplicit:: setDofs(double dofs_p[18])
+{
+  for(int i=0;i<18; i++)
+    dofs[i]=dofs_p[i];
+  glb2locDofs(dofs);
+}
+void ReducedQuinticExplicit:: print()
+{
+}
+void ReducedQuinticExplicit:: eval_g( double coord[2], double res[6] )
+{
+  double coord_loc[2]={coord[0],coord[1]};
+  glb2loc(coord_loc);
+  eval_field_i_(coord_loc,coord_loc+1,dofs,res);
+}
+void ReducedQuinticExplicit:: eval_l( double coord[2], double res[6] )
+{
+  double coord_loc[2]={coord[0],coord[1]};
+  eval_field_i_(coord_loc,coord_loc+1,dofs,res);
+}
+void ReducedQuinticExplicit:: loc2glb( double coord[2])
+{
+  double coord_glb[2];
+  loc_xyz_2_glb_i_(coord,coord_glb);
+  coord[0]=coord_glb[0];
+  coord[1]=coord_glb[1];
+}
+void ReducedQuinticExplicit:: glb2loc( double coord[2])
+{
+  double coord_loc[2];
+  glb_xyz_2_loc_i_(coord,coord_loc);
+  coord[0]=coord_loc[0];
+  coord[1]=coord_loc[1];
+}
+void ReducedQuinticExplicit:: loc2glbDofs(double dofs_p[18])
+{
+  double dofs_glb[18];
+  loc_dof_2_glb_i_(dofs_p,dofs_glb);
+  for(int i=0; i<18; i++)
+    dofs_p[i]=dofs_glb[i];
+}
+void ReducedQuinticExplicit:: glb2locDofs(double dofs_p[18])
+{
+  double dofs_loc[18];
+  glb_dof_2_loc_i_(dofs_p,dofs_loc);
+  for(int i=0; i<18; i++)
+    dofs_p[i]=dofs_loc[i];
+}
diff --git a/m3dc1_scorec/utility/check_smb.cc b/m3dc1_scorec/utility/check_smb.cc
new file mode 100644
index 000000000..371f243ff
--- /dev/null
+++ b/m3dc1_scorec/utility/check_smb.cc
@@ -0,0 +1,128 @@
+#include <gmi_mesh.h>
+#include <gmi_null.h>
+#include <apf.h>
+#include <apfMesh2.h>
+#include <apfMDS.h>
+#include <PCU.h>
+#include <parma.h>
+#ifdef HAVE_SIMMETRIX
+#include <gmi_sim.h>
+#include <SimUtil.h>
+#endif
+#include <apfZoltan.h>
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+
+namespace {
+
+const char* meshFile = 0;
+
+void freeMesh(apf::Mesh* m)
+{
+  m->destroyNative();
+  apf::destroyMesh(m);
+}
+
+void getConfig(int argc, char** argv)
+{
+  if ( argc < 2 ) {
+    if ( !PCU_Comm_Self() )
+      printf("Usage: %s <mesh>\n", argv[0]);
+    MPI_Finalize();
+    exit(EXIT_FAILURE);
+  }
+}
+
+void print_info (apf::Mesh2* m)
+{
+  if (!PCU_Comm_Self()) std::cout<<"\n===== mesh size =====\n";
+
+  int* local_entity_count = new int[4*PCU_Comm_Peers()];
+  int* own_entity_count = new int[4*PCU_Comm_Peers()];
+
+  for (int i=0; i<4*PCU_Comm_Peers();++i)
+    local_entity_count[i]=own_entity_count[i]=0;
+
+  apf::MeshEntity* e;
+  int self = PCU_Comm_Self();
+
+  for (int d=0; d<4;++d)
+  {
+    local_entity_count[4*self+d] = m->count(d);
+    apf::MeshIterator* it = m->begin(d); 
+    while ((e = m->iterate(it)))
+    {
+      if (m->getOwner(e)==self)
+        ++own_entity_count[4*PCU_Comm_Self()+d];
+    }
+    m->end(it);
+  }
+
+  int* global_local_entity_count = new int[4*PCU_Comm_Peers()];
+  int* global_own_entity_count = new int[4*PCU_Comm_Peers()];
+
+  MPI_Allreduce(local_entity_count, global_local_entity_count, 4*PCU_Comm_Peers(),
+                MPI_INT, MPI_SUM, PCU_Get_Comm());
+
+  MPI_Allreduce(own_entity_count, global_own_entity_count, 4*PCU_Comm_Peers(),
+                MPI_INT, MPI_SUM, PCU_Get_Comm());
+
+  if (!PCU_Comm_Self())
+  {
+    int* global_entity_count = new int[4];
+    global_entity_count[0]=global_entity_count[1]=global_entity_count[2]=global_entity_count[3]=0;
+    for (int d=0; d<4;++d)
+    {
+      for (int p=0; p<PCU_Comm_Peers();++p)
+        global_entity_count[d] += global_own_entity_count[p*4+d];
+    }
+
+    for (int p=0; p<PCU_Comm_Peers(); ++p)
+      std::cout<<"(p"<<p<<") # local ent: v "<<global_local_entity_count[p*4]
+        <<", e "<<global_local_entity_count[p*4+1]
+        <<", f "<<global_local_entity_count[p*4+2]
+        <<", r "<<global_local_entity_count[p*4+3]<<"\n";
+    std::cout<<"\n";
+    for (int p=0; p<PCU_Comm_Peers(); ++p)
+      if (global_own_entity_count[p*4])
+        std::cout<<"(p"<<p<<") # own ent: v "<<global_own_entity_count[p*4]
+          <<", e "<<global_own_entity_count[p*4+1]
+          <<", f "<<global_own_entity_count[p*4+2]
+          <<", r "<<global_own_entity_count[p*4+3]<<"\n";
+    std::cout<<"\n";
+
+    std::cout<<"# global ent: v "<<global_entity_count[0]<<", e "<<global_entity_count[1]
+              <<", f "<<global_entity_count[2]<<", r "<<global_entity_count[3]<<"\n\n";
+
+    delete [] global_entity_count;
+  }
+
+  delete [] local_entity_count;
+  delete [] global_local_entity_count;
+  delete [] own_entity_count;
+  delete [] global_own_entity_count;
+}
+}
+
+int main(int argc, char** argv)
+{
+  MPI_Init(&argc,&argv);
+  PCU_Comm_Init();
+  double t0 = PCU_Time();
+
+  gmi_register_null();
+  getConfig(argc,argv);
+  gmi_model* g = gmi_load(".null");
+  apf::Mesh2* m = 0;
+  m = apf::loadMdsMesh(g, argv[1]);
+  apf::verify(m, false);
+  //Parma_PrintPtnStats(m, "");
+  print_info(m);
+
+  if (!PCU_Comm_Self())
+    std::cout<<"check_smb: completed\n";
+  freeMesh(m);
+  PCU_Comm_Free();
+  MPI_Finalize();
+}
diff --git a/m3dc1_scorec/utility/create_smb.cc b/m3dc1_scorec/utility/create_smb.cc
new file mode 100644
index 000000000..aa9b72674
--- /dev/null
+++ b/m3dc1_scorec/utility/create_smb.cc
@@ -0,0 +1,179 @@
+/****************************************************************************** 
+
+  (c) 2005-2016 Scientific Computation Research Center, 
+      Rensselaer Polytechnic Institute. All rights reserved.
+  
+  This work is open source software, licensed under the terms of the
+  BSD license as described in the LICENSE file in the top-level directory.
+ 
+*******************************************************************************/
+#include <ma.h>
+#include <apf.h>
+#include <gmi_analytic.h>
+#include <apfMDS.h>
+#include <apfShape.h>
+#include <PCU.h>
+#include <parma.h>
+#include <apfZoltan.h>
+#include <assert.h>
+#include <stdlib.h>
+#include <iostream>
+
+double  a_param = 0;
+double  b_param = 0;
+double  c_param = 0;
+double  d_param = 0;
+double  e_param = 0;
+double meshSize=0.1;
+void getCenter(double x[3])
+{
+   x[0]=a_param;
+   x[1]=d_param;
+   x[2]=0.;
+}
+void edgeFunction(double const p[2], double x[3], void*)
+{
+  double phi = p[0];
+  x[0] = a_param + b_param*(cos(phi + c_param*sin(phi)));
+  x[1] = d_param + e_param*sin(phi);
+  x[2] = 0;
+}
+
+void snap2Model(apf::Mesh2* mesh, gmi_model* model)
+{
+  for(int i=0; i<5; i++)
+  {
+    apf::MeshEntity* vt = getMdsEntity(mesh, 0, i);
+    assert(vt);
+    gmi_ent* gent= (gmi_ent*)(mesh->toModel(vt));
+    apf::Vector3 param(0,0,0);
+    int gType = gmi_dim(model,gent);
+    double xyz[3];
+    if(gType==1)
+    {
+      mesh->getParam(vt,param);
+      double pa[]={param[0],0,0};
+      edgeFunction(pa, xyz, 0);
+    }
+    else
+    {
+      assert(gType==2);
+      getCenter(xyz);
+    }
+    mesh->setPoint(vt, 0, xyz);
+  }
+}
+void faceFunction(double const p[2], double x[3], void*)
+{
+  (void)p;
+  (void)x;
+}
+
+gmi_model* makeModel(char* modelfile)
+{
+  FILE* fp=fopen(modelfile,"r");
+  if(!fp)
+  {
+    std::cout<<"SCOREC ERROR: fail to open the file \""<<modelfile<<"\"\n";
+    return NULL;
+  }
+  fscanf(fp, "%lf %lf %lf %lf %lf\n", &a_param, &b_param, &c_param, &d_param, &e_param);
+  fclose(fp);
+  gmi_model* model = gmi_make_analytic();
+  int edgePeriodic = 1;
+  double edgeRange[2] = {0, 2 * apf::pi};
+  gmi_add_analytic(model, 1, 1, edgeFunction, &edgePeriodic, &edgeRange, 0);
+  int facePeriodic[2] = {0, 0};
+  double faceRanges[2][2] = {{0,0},{0,0}};
+  gmi_add_analytic(model, 2, 1, faceFunction, facePeriodic, faceRanges, 0);
+  return model;
+}
+
+class Size : public ma:: IsotropicFunction
+{
+  public:
+    Size(ma::Mesh* m){}
+    virtual double getValue(ma::Entity* vert) {return meshSize;};
+};
+
+static void testIndexing(apf::Mesh2* m)
+{
+  for (int d = 0; d <= m->getDimension(); ++d) {
+    apf::MeshIterator* it = m->begin(d);
+    int i = 0;
+    apf::MeshEntity* e;
+    while ((e = m->iterate(it))) {
+      assert( apf::getMdsIndex(m, e) == i );
+      assert( apf::getMdsEntity(m, d, i) == e );
+      ++i;
+    }
+    m->end(it);
+  }
+}
+
+static void fusionAdapt(apf::Mesh2* m)
+{
+  Size sf(m);
+  for(int i=0; i<10; i++)
+  {
+    std::cout<<" iter "<<i<<std::endl;
+    #ifdef OLDMA
+  	ma::Input* in = ma::configure(m, &sf);
+    #else
+	ma::Input* in = ma::makeAdvanced(ma::configure(m,&sf));
+    #endif
+    in->maximumIterations = 1;
+    in->shouldRunPreZoltan = true;
+    in->shouldRunMidParma = true;
+    in->shouldRunPostParma = true;
+    ma::adapt(in);
+  }
+}
+
+int main( int argc, char* argv[])
+{
+  MPI_Init(&argc,&argv);
+  PCU_Comm_Init();
+
+  if (argc!=3)
+  {
+    std::cout<<"SCOREC ERROR: wrong input arguments. \nUsage: ./create_smb model_file relative_mesh_size"
+             <<"\n - model_file shall contain 5 double: a_param, b_param, c_param, d_param, e_param \n"
+             <<" - given relative_mesh_size, mesh edge length is about relative_mesh_size * the longest edge of the bounding box of the model\n";
+    PCU_Comm_Free();
+    MPI_Finalize();    
+    return 1;
+  }
+
+  char* modelfile = argv[1];
+  meshSize =atof(argv[2]);
+
+  // create model
+  gmi_model* model = makeModel(modelfile);
+  if (!model)
+  {
+    PCU_Comm_Free();
+    MPI_Finalize();    
+    return 1;
+  }
+  // generate and refine mesh
+  apf::Mesh2* mesh = apf::loadMdsMesh(model, "seed.smb");
+  snap2Model(mesh, model);
+  mesh->verify();
+  testIndexing(mesh);
+  fusionAdapt(mesh);
+  // export to files
+  char meshfile[256];
+  sprintf(meshfile,"%s.smb",modelfile); 
+  mesh->writeNative(meshfile);
+  apf::writeVtkFiles("Vtk-Mesh-File", mesh);
+  std::cout<<"\n";
+  mesh->verify();
+  std::cout<<"\n";
+  // clean and finalize
+  mesh->destroyNative();
+  apf::destroyMesh(mesh);
+  PCU_Comm_Free();
+  MPI_Finalize();
+  return 0;
+}
diff --git a/m3dc1_scorec/utility/show_meshcount.cc b/m3dc1_scorec/utility/show_meshcount.cc
new file mode 100644
index 000000000..51a51b2b2
--- /dev/null
+++ b/m3dc1_scorec/utility/show_meshcount.cc
@@ -0,0 +1,78 @@
+/****************************************************************************** 
+
+  (c) 2005-2016 Scientific Computation Research Center, 
+      Rensselaer Polytechnic Institute. All rights reserved.
+  
+  This work is open source software, licensed under the terms of the
+  BSD license as described in the LICENSE file in the top-level directory.
+ 
+*******************************************************************************/
+#include <apf.h>
+#include <PCU.h>
+#include <stdlib.h>
+#include "m3dc1_scorec.h"
+#include "m3dc1_mesh.h"
+#include "apfMesh.h"
+
+char* meshFile = 0;
+char* modelFile = 0;
+int num_plane=1;
+
+void getConfig(int argc, char** argv)
+{
+  if ( argc < 4 ) 
+  {
+    if ( !PCU_Comm_Self() )
+      printf("Usage: %s <in-model> <in-mesh.smb> <#plane>\n", argv[0]);
+    m3dc1_scorec_finalize();
+    MPI_Finalize();
+    return;
+  }
+  modelFile = argv[1];
+  meshFile = argv[2];
+  if (atoi(argv[3])>1) num_plane=atoi(argv[3]);
+}
+
+int main(int argc, char** argv)
+{
+  MPI_Init(&argc,&argv);
+  m3dc1_scorec_init();
+
+  getConfig(argc,argv);
+
+  if (num_plane>1) // construct 3D
+  {
+    if (PCU_Comm_Peers()%num_plane!=0) 
+    {
+      if ( !PCU_Comm_Self() )
+        printf("ERROR: invalid #planes\n");
+      m3dc1_scorec_finalize();
+      MPI_Finalize();
+      return 1;
+    }
+
+    m3dc1_model_setnumplane (&num_plane);
+  }
+
+
+  m3dc1_model_load(modelFile);
+  m3dc1_mesh_load(meshFile);
+//  m3dc1_mesh::instance()->print(0);
+
+  if (num_plane>1)
+  {
+    int zero=0;
+    m3dc1_mesh_build3d(&zero, &zero, &zero);
+  }
+
+ // printStats(m3dc1_mesh::instance()->mesh);
+
+//Is there any way to tell how many vertices per plane after partition is applied from just the mesh file without going to actually M3DC1 code? And how many vertices per partition as well, without counting the ghost vertices which belong to other rank?
+  m3dc1_mesh::instance()->print(1);
+
+  // finalize
+  m3dc1_scorec_finalize();
+  MPI_Finalize();
+  return 0;
+}
+
diff --git a/m3dc1_scorec/utility/split_smb.cc b/m3dc1_scorec/utility/split_smb.cc
new file mode 100644
index 000000000..41275f108
--- /dev/null
+++ b/m3dc1_scorec/utility/split_smb.cc
@@ -0,0 +1,124 @@
+#include <gmi_mesh.h>
+#include <gmi_null.h>
+#include <apf.h>
+#include <apfMesh2.h>
+#include <apfMDS.h>
+#include <PCU.h>
+#include <parma.h>
+#ifdef HAVE_SIMMETRIX
+#include <gmi_sim.h>
+#include <SimUtil.h>
+#endif
+#include <apfZoltan.h>
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+
+namespace {
+
+const char* modelFile = 0;
+const char* meshFile = 0;
+const char* outFile = 0;
+int partitionFactor = 1;
+
+void freeMesh(apf::Mesh* m)
+{
+  m->destroyNative();
+  apf::destroyMesh(m);
+}
+
+apf::Migration* getPlan(apf::Mesh* m)
+{
+  apf::Splitter* splitter = apf::makeZoltanSplitter(
+      m, apf::GRAPH, apf::PARTITION, false);
+  apf::MeshTag* weights = Parma_WeighByMemory(m);
+  apf::Migration* plan = splitter->split(weights, 1.05, partitionFactor);
+  apf::removeTagFromDimension(m, weights, m->getDimension());
+  m->destroyTag(weights);
+  delete splitter;
+  return plan;
+}
+
+void switchToOriginals()
+{
+  int self = PCU_Comm_Self();
+  int groupRank = self / partitionFactor;
+  int group = self % partitionFactor;
+  MPI_Comm groupComm;
+  MPI_Comm_split(MPI_COMM_WORLD, group, groupRank, &groupComm);
+  PCU_Switch_Comm(groupComm);
+}
+
+void switchToAll()
+{
+  MPI_Comm prevComm = PCU_Get_Comm();
+  PCU_Switch_Comm(MPI_COMM_WORLD);
+  MPI_Comm_free(&prevComm);
+  PCU_Barrier();
+}
+
+void getConfig(int argc, char** argv)
+{
+  if ( argc < 4 ) {
+    if ( !PCU_Comm_Self() )
+      printf("Usage: %s <mesh> <outMesh> <factor>\n"
+               "Increase the part count of mesh from inPartCount to #rank (inPartCount*factor).\n",
+               argv[0]);
+    MPI_Finalize();
+    exit(EXIT_FAILURE);
+  }
+  if (argc == 5)
+  {
+    meshFile = argv[2];
+    outFile = argv[3];
+    partitionFactor = atoi(argv[4]);
+  }
+  if (argc == 4)
+  {
+    meshFile = argv[1];
+    outFile = argv[2];
+    partitionFactor = atoi(argv[3]);
+  }
+  assert(partitionFactor <= PCU_Comm_Peers());
+}
+
+}
+
+int main(int argc, char** argv)
+{
+  MPI_Init(&argc,&argv);
+  PCU_Comm_Init();
+  double t0 = PCU_Time();
+#ifdef HAVE_SIMMETRIX
+  SimUtil_start();
+  Sim_readLicenseFile(0);
+  gmi_sim_start();
+  gmi_register_sim();
+#endif
+  gmi_register_null();
+  getConfig(argc,argv);
+  bool isOriginal = ((PCU_Comm_Self() % partitionFactor) == 0);
+  gmi_model* g = gmi_load(".null");
+  apf::Mesh2* m = 0;
+  apf::Migration* plan = 0;
+  switchToOriginals();
+  if (isOriginal) {
+    m = apf::loadMdsMesh(g, meshFile);
+    plan = getPlan(m);
+  }
+  switchToAll();
+  m = apf::repeatMdsMesh(m, g, plan, partitionFactor);
+  if (!PCU_Comm_Self())
+    std::cout<<"[PUMI INFO] mesh splitted into "<<PCU_Comm_Peers()<<" parts in "
+             <<PCU_Time()-t0<<"(sec)\n";
+  Parma_PrintPtnStats(m, "");
+  m->writeNative(outFile);
+  freeMesh(m);
+#ifdef HAVE_SIMMETRIX
+  gmi_sim_stop();
+  Sim_unregisterAllKeys();
+  SimUtil_stop();
+#endif
+  PCU_Comm_Free();
+  MPI_Finalize();
+}
diff --git a/m3dc1_scorec/utility/vtk_order.cc b/m3dc1_scorec/utility/vtk_order.cc
new file mode 100644
index 000000000..3e2a33e8c
--- /dev/null
+++ b/m3dc1_scorec/utility/vtk_order.cc
@@ -0,0 +1,180 @@
+/****************************************************************************** 
+
+  (c) 2005-2016 Scientific Computation Research Center, 
+      Rensselaer Polytechnic Institute. All rights reserved.
+  
+  This work is open source software, licensed under the terms of the
+  BSD license as described in the LICENSE file in the top-level directory.
+ 
+*******************************************************************************/
+#include <apf.h>
+#include <apfMesh.h>
+#include <apfMesh2.h>
+#include <apfMDS.h>
+#include <deque>
+#include <set>
+#include <PCU.h>
+#include <assert.h>
+#include <stdlib.h>
+#include <iostream>
+#include <apfVector.h>
+#include <apfNumbering.h>
+#include <apfShape.h>
+#include <gmi.h>
+#include <gmi_null.h>
+
+// seol -- reorder input mesh before conversion
+//         starting vtx: a vtx with min Y
+struct Queue {
+  bool has(apf::MeshEntity* e) { return h.count(e); }
+  void push(apf::MeshEntity* e)
+  {
+    q.push_back(e);
+    h.insert(e);
+  }
+  void pushVector(std::vector<apf::MeshEntity*> const& l)
+  {
+    for (size_t i = 0; i < l.size(); ++i)
+      push(l[i]);
+  }
+  apf::MeshEntity* pop()
+  {
+    apf::MeshEntity* e;
+    e = q.front();
+    q.pop_front();
+    h.erase(e);
+    return e;
+  }
+  bool empty() { return q.empty(); }
+  std::deque<apf::MeshEntity*> q;
+  std::set<apf::MeshEntity*> h;
+};
+
+
+apf::MeshEntity* findFirst(apf::Mesh* m)
+{
+  apf::MeshEntity* v;
+  apf::MeshEntity* best;
+  apf::MeshIterator* it = m->begin(0);
+  best = m->iterate(it);
+  apf::Vector3 coord;
+  m->getPoint(best, 0, coord);
+  double min_Y=coord[1];
+  while ((v = m->iterate(it)))
+  {
+    m->getPoint(v, 0, coord);  
+    if (min_Y > coord[1])
+    {
+      best = v;
+      min_Y = coord[1];
+    }
+  }
+  m->end(it);
+  return best;
+}
+
+bool visited(Queue& q, apf::Numbering* nn, apf::MeshEntity* e)
+{
+  return apf::isNumbered(nn, e, 0, 0) || q.has(e);
+}
+
+bool hasNode(apf::Mesh* m, apf::MeshEntity* e)
+{
+  if (m->getShape()->countNodesOn(m->getType(e))>0) 
+    return true;
+  return false;
+}
+
+int main( int argc, char* argv[])
+{
+  MPI_Init(&argc,&argv);
+  PCU_Comm_Init();
+
+  if (argc!=2)
+  {
+    std::cout<<"SCOREC ERROR: wrong input arguments. \n"
+             <<"  Usage: ./reorder inmesh.smb\n";
+    PCU_Comm_Free();
+    MPI_Finalize();    
+    return 1;
+  }
+  gmi_register_null();
+  gmi_model* nullModel = gmi_load(".null");
+
+  apf::Mesh2* m = apf::loadMdsMesh(nullModel,argv[1]);
+  int mesh_dim = m->getDimension();
+
+  apf::Numbering* n_org=apf::numberOwnedDimension(m, "elem", mesh_dim);
+  writeVtkFiles("org-mesh.vtk", m);
+  std::cout << "[INFO] vtk files for original ordering saved in \"org-mesh.vtk\"\n";
+  destroyNumbering(n_org);
+
+  // reorder mesh
+  apf::Numbering* nn = apf::createNumbering(m, "node", apf::getConstant(0), 1);
+  apf::Numbering* en = apf::createNumbering(m, "elem", apf::getConstant(mesh_dim), 1);
+
+  Queue q;
+  q.push(findFirst(m));
+
+  // node and element number starts from 0
+  int labelnode = 0;
+  int labelelem = 0;
+
+  std::vector<apf::MeshEntity*> node_arr;
+  std::vector<apf::MeshEntity*> elem_arr;
+
+  node_arr.resize(m->count(0)+1);
+  elem_arr.resize(m->count(mesh_dim)+1);
+
+  apf::MeshEntity* otherVtx;
+  apf::MeshEntity* edge;
+  apf::MeshEntity* elem;
+
+  while (!q.empty()) 
+  {
+    apf::MeshEntity* vtx = q.pop();
+    if (!apf::isNumbered(nn, vtx, 0, 0))
+    {
+      node_arr[labelnode] = vtx;
+      apf::number(nn, vtx, 0, 0, labelnode);
+
+      ++labelnode;
+    }
+
+    std::vector<apf::MeshEntity*> entities;
+    apf::Adjacent edges;
+    m->getAdjacent(vtx,1, edges);
+    for (size_t i = 0; i < edges.getSize(); ++i) 
+    {
+      edge = edges[i];
+      apf::Adjacent adjacent;
+      m->getAdjacent(edge, mesh_dim, adjacent);      
+      for (size_t j = 0; j < adjacent.getSize(); ++j) 
+      {
+        elem = adjacent[j];
+        if (!apf::isNumbered(en, elem, 0, 0))
+        {
+          elem_arr[labelelem] = elem;
+          apf::number(en, elem, 0, 0, labelelem);
+          ++labelelem;
+        }
+      }
+      otherVtx = apf::getEdgeVertOppositeVert(m, edge, vtx);
+      if (!visited(q, nn, otherVtx))
+        entities.push_back(otherVtx);
+    }
+    q.pushVector(entities);
+  } // while
+  destroyNumbering(nn);
+
+  writeVtkFiles("reordered-mesh.vtk", m);
+  std::cout << "[INFO] vtk files for adjacency-based ordering saved in \"reordered-mesh.vtk\"\n";
+  destroyNumbering(en);
+
+  m->destroyNative();
+  apf::destroyMesh(m);
+
+  PCU_Comm_Free();
+  MPI_Finalize();
+  return 0;
+}
diff --git a/skeleton/m3dc1_skeleton_big_spec_acc/main.F90 b/skeleton/m3dc1_skeleton_big_spec_acc/main.F90
index 257f47a46..6532fd3f8 100644
--- a/skeleton/m3dc1_skeleton_big_spec_acc/main.F90
+++ b/skeleton/m3dc1_skeleton_big_spec_acc/main.F90
@@ -1,4 +1,5 @@
 Program m3dc1_skeleton
+  use mpi
   use mesh
   use m3dc1_data
   use matdef
@@ -7,7 +8,6 @@ Program m3dc1_skeleton
 
   implicit none
 
-  include 'mpif.h'
 
   integer :: ierr
   real :: tstart, tend
@@ -151,12 +151,12 @@ end subroutine rotation
   ! populates the *tri arrays
   !============================================================
   subroutine tridef
+    use mpi
     use math
     use mesh
 
     implicit none
 
-    include 'mpif.h'
   
     type(element_data) :: d
     integer :: itri, i, j, k, ii, jj, numelms, numnodes, ndofs
diff --git a/skeleton/m3dc1_skeleton_big_spec_acc/main_acc.F90 b/skeleton/m3dc1_skeleton_big_spec_acc/main_acc.F90
index 7a4ddfc2e..aba4be579 100644
--- a/skeleton/m3dc1_skeleton_big_spec_acc/main_acc.F90
+++ b/skeleton/m3dc1_skeleton_big_spec_acc/main_acc.F90
@@ -1,4 +1,5 @@
 Program m3dc1_skeleton
+  use mpi
   use mesh
   use m3dc1_data
   use matdef
@@ -8,7 +9,6 @@ Program m3dc1_skeleton
 
   implicit none
 
-  include 'mpif.h'
 
   integer :: ierr
   real :: tstart, tend
@@ -176,12 +176,12 @@ end subroutine rotation
   ! populates the *tri arrays
   !============================================================
   subroutine tridef
+    use mpi
     use math
     use mesh
 
     implicit none
 
-    include 'mpif.h'
   
     type(element_data) :: d
     integer :: itri, i, j, k, ii, jj, numelms, numnodes, ndofs
diff --git a/skeleton/m3dc1_skeleton_big_spec_acc/main_omp.F90 b/skeleton/m3dc1_skeleton_big_spec_acc/main_omp.F90
index e796e8681..661b6b4ec 100644
--- a/skeleton/m3dc1_skeleton_big_spec_acc/main_omp.F90
+++ b/skeleton/m3dc1_skeleton_big_spec_acc/main_omp.F90
@@ -1,4 +1,5 @@
 Program m3dc1_skeleton
+  use mpi
   use mesh
   use m3dc1_data
   use matdef
@@ -7,7 +8,6 @@ Program m3dc1_skeleton
 
   implicit none
 
-  include 'mpif.h'
 
   integer :: ierr
   real :: tstart, tend
@@ -151,12 +151,12 @@ end subroutine rotation
   ! populates the *tri arrays
   !============================================================
   subroutine tridef
+    use mpi
     use math
     use mesh
 
     implicit none
 
-    include 'mpif.h'
   
     type(element_data) :: d
     integer :: itri, i, j, k, ii, jj, numelms, numnodes, ndofs
diff --git a/skeleton/m3dc1_skeleton_big_spec_acc/mesh.F90 b/skeleton/m3dc1_skeleton_big_spec_acc/mesh.F90
index 0411d57a8..a4daa975c 100644
--- a/skeleton/m3dc1_skeleton_big_spec_acc/mesh.F90
+++ b/skeleton/m3dc1_skeleton_big_spec_acc/mesh.F90
@@ -127,11 +127,11 @@ subroutine create_rectangular_mesh(m,n, nplanes,inodes, ielms, nodes, elms)
   end subroutine create_rectangular_mesh
 
   subroutine load_mesh(isize)
+    use mpi
     use math
 
     implicit none
 
-    include "mpif.h"
 
     integer :: isize, new_mesh_size
     integer :: num_global_elements
diff --git a/skeleton/m3dc1_skeleton_big_spec_acc/mesh_acc.F90 b/skeleton/m3dc1_skeleton_big_spec_acc/mesh_acc.F90
index 03e4ff1e8..75d5b10b2 100644
--- a/skeleton/m3dc1_skeleton_big_spec_acc/mesh_acc.F90
+++ b/skeleton/m3dc1_skeleton_big_spec_acc/mesh_acc.F90
@@ -132,11 +132,11 @@ subroutine create_rectangular_mesh(m,n, nplanes,inodes, ielms, nodes, elms)
   end subroutine create_rectangular_mesh
 
   subroutine load_mesh(isize)
+    use mpi
     use math
 
     implicit none
 
-    include "mpif.h"
 
     integer :: isize, new_mesh_size
     integer :: num_global_elements
diff --git a/skeleton/m3dc1_skeleton_big_spec_omp/main.F90 b/skeleton/m3dc1_skeleton_big_spec_omp/main.F90
index d49e9787d..133a0eb5e 100644
--- a/skeleton/m3dc1_skeleton_big_spec_omp/main.F90
+++ b/skeleton/m3dc1_skeleton_big_spec_omp/main.F90
@@ -1,4 +1,5 @@
 Program m3dc1_skeleton
+  use mpi
   use mesh
   use m3dc1_data
   use matdef
@@ -7,7 +8,6 @@ Program m3dc1_skeleton
 
   implicit none
 
-  include 'mpif.h'
 
   integer :: ierr
   real :: tstart, tend
@@ -146,12 +146,12 @@ end subroutine rotation
   ! populates the *tri arrays
   !============================================================
   subroutine tridef
+    use mpi
     use math
     use mesh
 
     implicit none
 
-    include 'mpif.h'
   
     type(element_data) :: d
     integer :: itri, i, j, k, ii, jj, numelms, numnodes, ndofs
diff --git a/skeleton/m3dc1_skeleton_big_spec_omp/mesh.F90 b/skeleton/m3dc1_skeleton_big_spec_omp/mesh.F90
index c2f8aa99f..849ea2f9a 100644
--- a/skeleton/m3dc1_skeleton_big_spec_omp/mesh.F90
+++ b/skeleton/m3dc1_skeleton_big_spec_omp/mesh.F90
@@ -127,11 +127,11 @@ subroutine create_rectangular_mesh(m,n, nplanes,inodes, ielms, nodes, elms)
   end subroutine create_rectangular_mesh
 
   subroutine load_mesh
+    use mpi
     use math
 
     implicit none
 
-    include "mpif.h"
 
     integer :: num_global_elements
     type(node_type), allocatable :: global_node(:)
diff --git a/unstructured/CMakeLists.txt b/unstructured/CMakeLists.txt
index ac86dd2ad..28ee00258 100644
--- a/unstructured/CMakeLists.txt
+++ b/unstructured/CMakeLists.txt
@@ -1,487 +1,379 @@
-cmake_minimum_required(VERSION 3.12)
-
-# set(CMAKE_DISABLE_SOURCE_CHANGES ON)
-# set(CMAKE_DISABLE_IN_SOURCE_BUILD ON)
-
-project(m3dc1_unstructured LANGUAGES Fortran C CXX)
-
-
-if("${PROJECT_SOURCE_DIR}" STREQUAL "${PROJECT_BINARY_DIR}")
-   message(SEND_ERROR "In-source builds are not allowed.")
-endif("${PROJECT_SOURCE_DIR}" STREQUAL "${PROJECT_BINARY_DIR}")
-
-# enable_language(Fortran)
-
-list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
-
-if(NOT CMAKE_BUILD_TYPE)
-  set(CMAKE_BUILD_TYPE Release)
+set(M3DC1_MESH_MOD "scorec_mesh_mod" CACHE STRING
+  "Mesh module name")
+set(M3DC1_VECTOR_MOD "scorec_vector_mod" CACHE STRING
+  "Vector module name")
+set(M3DC1_MATRIX_MOD "scorec_matrix_mod" CACHE STRING
+  "Matrix module name")
+set(M3DC1_VECTOR_TYPE "scorec_vector" CACHE STRING
+  "Vector derived-type name")
+set(M3DC1_MATRIX_TYPE "scorec_matrix" CACHE STRING
+  "Matrix derived-type name")
+
+message(STATUS "M3DC1 mesh   backend : ${M3DC1_MESH_MOD}")
+message(STATUS "M3DC1 vector backend : ${M3DC1_VECTOR_MOD}  (type = ${M3DC1_VECTOR_TYPE})")
+message(STATUS "M3DC1 matrix backend : ${M3DC1_MATRIX_MOD}  (type = ${M3DC1_MATRIX_TYPE})")
+
+
+set(M3DC1_LIBRARY_SOURCES
+  adapt.f90
+  auxiliary_fields.f90
+  biharmonic.f90
+  bootstrap.f90
+  boundary.f90
+  coil_sets.f90
+  coils.f90
+  control.f90
+  diagnostics.f90
+  electric_field.f90
+  electrostatic_potential.f90
+  element.f90
+  error_estimate.f90
+  field.f90
+  fit_magnetics.f90
+  fftw_fortran.c
+  geometry.f90
+  get_pc_skip_count.f90
+  gsl_wrapper.c
+  gradshafranov.f90
+  gyroviscosity.f90
+  harned_mikic.f90
+  hdf5_output.f90
+  hypervisc.f90
+  init_3ddiffusion.f90
+  init_3dwave.f90
+  init_basicj.f90
+  init_basicq.f90
+  init_circle.f90
+  init_common.f90
+  init_conds.f90
+  init_cyl.f90
+  init_dskbal.f90
+  init_eigen.f90
+  init_eqdsk.f90
+  init_force_free.f90
+  init_frs.f90
+  init_ftz.f90
+  init_gem.f90
+  init_gmode.f90
+  init_intkink.f90
+  init_jsolver.f90
+  read_jsolver.f90
+  readaeqdsk.f90
+  init_kstar.f90
+  init_lz.f90
+  init_mri.f90
+  init_rotating_cylinder.f90
+  init_rwm.f90
+  init_solovev.f90
+  init_strauss.f90
+  init_taylor.f90
+  init_tilt.f90
+  init_vmec.f90
+  init_wave.f90
+  input.f90
+  interpolate.f90
+  iterdb.f90
+  kprad_m3dc1.f90
+  kprad.f90
+  ludef_t.f90
+  m3dc1_nint.f90
+  M3Dmodules.f90
+  math.f90
+  metricterms_new.f90
+  model.f90
+  neutral_beam.f90
+  newvar.f90
+  newpar.f90
+  nintegrate_mod.f90
+  output.f90
+  parallel_heat_flux.f90
+  parse_solver_options.c
+  particle.f90
+  pellet.f90
+  physical_mesh.f90
+  polar.f90
+  radiation.f90
+  random.c
+  read_ascii.f90
+  read_dskbal.f90
+  read_gyro.f90
+  read_neo.f90
+  read_namelist.cpp
+  read_schaffer_field.f90
+  read_vmec.f90
+  readgeqdsk.f90
+  region.f90
+  resistive_wall.f90
+  restart_hdf5.f90
+  rmp.f90
+  runaway_advection.f90
+  runaway.f90
+  scorec_adapt.f90
+  signal_handler.f90
+  spline.f90
+  temperature_plots.f90
+  time_step_split.f90
+  time_step_unsplit.f90
+  time_step.f90
+  transport.f90
+  two_fluid.f90
+  wall.f90
+  mesh_mod.F90
+  vector_mod.F90
+  matrix_mod.F90
+)
+
+if(M3DC1_MESH_MOD STREQUAL "scorec_mesh_mod")
+  list(APPEND M3DC1_LIBRARY_SOURCES scorec_mesh.f90)
+  set(M3DC1_MESH_DEFINE M3DC1_MESH_IS_SCOREC)
+elseif(M3DC1_MESH_MOD STREQUAL "basic_mesh_mod")
+  list(APPEND M3DC1_LIBRARY_SOURCES basic_mesh.f90)
+  set(M3DC1_MESH_DEFINE M3DC1_MESH_IS_BASIC)
+else()
+  message(FATAL_ERROR "Unknown M3DC1_MESH_MOD=${M3DC1_MESH_MOD}")
 endif()
 
-# set(CMAKE_FIND_LIBRARY_SUFFIXES ".a" ".so")
-# if(BUILD_SHARED_LIBS)
-  # set(CMAKE_FIND_LIBRARY_SUFFIXES ".so" ".a")
-# endif()
-message(STATUS "Found PETSc ${CMAKE_FIND_LIBRARY_SUFFIXES}")
-
-# set(CMAKE_CXX_FLAGS "-Wall -Wextra")
-# set(CMAKE_CXX_FLAGS_DEBUG "-g")
-# set(CMAKE_CXX_FLAGS_RELEASE "-O3")
-
-set(CMAKE_Fortran_PREPROCESS "ON")
-
-# option(BUILD_SHARED_LIBS "Build using shared libraries" OFF)
-
-set(SOURCES
-    scorec_mesh.f90
-    scorec_vector.f90
-    scorec_matrix.f90
-    scorec_adapt.f90
-    PETScInterface.cpp
-    physical_mesh.f90
-    read_vmec.f90
-    fftw_fortran.c
-    read_namelist.cpp
-    gsl_wrapper.c
-    # random.c
-    spline.f90
-    kprad.f90
-    math.f90
-    read_ascii.f90
-    interpolate.f90
-    control.f90
-    iterdb.f90
-    read_gyro.f90
-    read_neo.f90
-    region.f90
-    radiation.f90
-    element.f90
-    field.f90
-    nintegrate_mod.f90
-    M3Dmodules.f90
-    resistive_wall.f90
-    m3dc1_nint.f90
-    boundary.f90
-    gyroviscosity.f90
-    bootstrap.f90
-    metricterms_new.f90
-    # metricterms_new_gpu.f90
-    two_fluid.f90
-    harned_mikic.f90
-    biharmonic.f90
-    electric_field.f90
-    pellet.f90
-    parallel_heat_flux.f90
-    temperature_plots.f90
-    electrostatic_potential.f90
-    newvar.f90
-    runaway.f90
-    diagnostics.f90
-    read_schaffer_field.f90
-    neutral_beam.f90
-    coils.f90
-    coil_sets.f90
-    model.f90
-    fit_magnetics.f90
-    wall.f90
-    init_common.f90
-    gradshafranov.f90
-    rmp.f90
-    kprad_m3dc1.f90
-    readgeqdsk.f90
-    read_jsolver.f90
-    read_dskbal.f90
-    init_cyl.f90
-    init_rwm.f90
-    init_solovev.f90
-    init_circle.f90
-    init_basicj.f90
-    init_tilt.f90
-    init_taylor.f90
-    init_force_free.f90
-    init_gem.f90
-    init_wave.f90
-    init_gmode.f90
-    init_strauss.f90
-    init_mri.f90
-    init_rotating_cylinder.f90
-    init_eqdsk.f90
-    init_dskbal.f90
-    init_jsolver.f90
-    init_3dwave.f90
-    init_3ddiffusion.f90
-    init_frs.f90
-    init_ftz.f90
-    init_eigen.f90
-    init_intkink.f90
-    init_lz.f90
-    init_kstar.f90
-    init_basicq.f90
-    init_vmec.f90
-    transport.f90
-    auxiliary_fields.f90
-    runaway_advection.f90
-    time_step_split.f90
-    time_step_unsplit.f90
-    hdf5_output.f90
-    particle.f90
-    time_step.f90
-    output.f90
-    restart_hdf5.f90
-    error_estimate.f90
-    adapt.f90
-    input.f90
-    ludef_t.f90
-    # ludef_t_gpu.f90
-    init_conds.f90
-    # get_pc_skip_count.f90
-    hypervisc.f90
-    geometry.f90
-    newpar.f90
-    )
-
-set(BIN_NAME m3dc1)
-
-# find_package(PETSc)
-# if(PETSC_FOUND AND USEPETSC)
-#   set(PETSC_SOURCES
-#       basic_mesh.f90
-#       petsc_vector.f90
-#       petsc_matrix.f90)
-#   set(SOURCES "${SOURCES};${PETSC_SOURCES}")
-#   list(APPEND DEFS -DUSEPETSC -Dmesh_mod=basic_mesh_mod -Dvector_mod=petsc_vector_mod -Dmatrix_mod=petsc_matrix_mod)
-#   set(POSTFIX "${POSTFIX}-petsc")
-# elseif(PETSC_FOUND)
-#   include(${CMAKE_SOURCE_DIR}/cmake/petsc.cmake)
-#   petsc_get_version()
-#   message(STATUS "Petsc version detected as : ${PETSC_VERSION}")
-#   list(APPEND DEFS -DPETSC_FORTRAN_PETSCTRUTH_INT -DPETSC_VERSION=${PETSC_VERSION})
-# endif()
-find_package(PkgConfig REQUIRED)
-# list(APPEND PKG_CONFIG_EXECUTABLE "--static")
-
-# pkg_check_modules (PETSC REQUIRED PETSc)
-pkg_search_module(PETSC PETSc IMPORTED_TARGET REQUIRED)
-message(STATUS "Found PETSc ${PETSC_VERSION}")
-message(STATUS "Found PETSc ${CMAKE_Fortran_COMPILER_ID}")
-
-# foreach(ITR ${PETSC_LINK_LIBRARIES})
-#     if(ITR MATCHES "liblapack\.a")
-#       set(LAPACK_LIBRARY ${ITR})
-#     endif()
-# endforeach()
-# list(FILTER PETSC_LINK_LIBRARIES EXCLUDE REGEX "liblapack\.a")
-# foreach(ITR ${PETSC_LINK_LIBRARIES})
-#     if(ITR MATCHES "libblas\.a")
-#       list(FIND PETSC_LINK_LIBRARIES ${ITR} IDX)
-#     endif()
-# endforeach()
-# list(INSERT PETSC_LINK_LIBRARIES ${IDX} ${LAPACK_LIBRARY})
-# message(STATUS "Found PETSc ${PETSC_STATIC_LIBRARIES}")
-set(PETSC_LINK_LIBRARIES "")
- # message(STATUS "Found PETSc ${PETSC_LINK_LIBRARIES}")
-message(STATUS "Found PETSc ${PETSC_STATIC_LDFLAGS}")
-# list(FIND PETSC_STATIC_LDFLAGS "-lmetis" _index)
-# math(EXPR _index "${_index}+1")
-# message(STATUS "Found PETSc ${_index}")
-# list(SUBLIST PETSC_STATIC_LDFLAGS 0 ${_index} PETSC_STATIC_LDFLAGS_MODIFIED)
-# message(STATUS "Found PETSc ${PETSC_STATIC_LDFLAGS_MODIFIED}")
+if(M3DC1_VECTOR_MOD STREQUAL "scorec_vector_mod")
+  list(APPEND M3DC1_LIBRARY_SOURCES scorec_vector.f90)
+  set(M3DC1_VECTOR_DEFINE M3DC1_VECTOR_IS_SCOREC)
+elseif(M3DC1_VECTOR_MOD STREQUAL "petsc_vector_mod")
+  list(APPEND M3DC1_LIBRARY_SOURCES petsc_vector.f90)
+  set(M3DC1_VECTOR_DEFINE M3DC1_VECTOR_IS_PETSC)
+else()
+  message(FATAL_ERROR "Unknown M3DC1_VECTOR_MOD=${M3DC1_VECTOR_MOD}")
+endif()
 
-# if(ENABLE_STATIC)
-#   set(STATIC_LIST lapack blas)
-# endif()
-foreach (flag IN LISTS PETSC_STATIC_LDFLAGS)
-  if (flag MATCHES "^-L(.*)")
-    list(APPEND _search_paths ${CMAKE_MATCH_1})
-    continue()
-  endif()
-  if (flag MATCHES "^-l(.*)")
-    set(_pkg_search "${CMAKE_MATCH_1}")
-  endif()
+if(M3DC1_MATRIX_MOD STREQUAL "scorec_matrix_mod")
+  list(APPEND M3DC1_LIBRARY_SOURCES scorec_matrix.f90)
+  set(M3DC1_MATRIX_DEFINE M3DC1_MATRIX_IS_SCOREC)
+elseif(M3DC1_MATRIX_MOD STREQUAL "petsc_matrix_mod")
+  list(APPEND M3DC1_LIBRARY_SOURCES petsc_matrix.f90)
+  set(M3DC1_MATRIX_DEFINE M3DC1_MATRIX_IS_PETSC)
+else()
+  message(FATAL_ERROR "Unknown M3DC1_MATRIX_MOD=${M3DC1_MATRIX_MOD}")
+endif()
 
-  # if (${_pkg_search} IN_LIST STATIC_LIST)
-  #   message(STATUS "Found PETSc ${_pkg_search}")
-  #   set(_pkg_search "lib${_pkg_search}.a")
-  # endif()
+if(M3DC1_ENABLE_GPU)
+  list(REMOVE_ITEM M3DC1_SOURCES ludef_t.f90 metricterms_new.f90)
+  list(APPEND M3DC1_SOURCES ludef_t_gpu.f90 metricterms_new_gpu.f90)
+endif()
 
-  if(_search_paths)
-      # Firstly search in -L paths
-      find_library(pkgcfg_lib_PETSC_${_pkg_search}
-                   NAMES ${_pkg_search}
-                   HINTS ${_search_paths} NO_DEFAULT_PATH)
-  endif()
-  find_library(pkgcfg_lib_PETSC_${_pkg_search}
-               NAMES ${_pkg_search})
-  if(pkgcfg_lib_PETSC_${_pkg_search})
-    list(APPEND PETSC_LINK_LIBRARIES "${pkgcfg_lib_PETSC_${_pkg_search}}")
-  else()
-    list(APPEND PETSC_LINK_LIBRARIES ${_pkg_search})
-  endif()
+if(NOT DEFINED M3DC1_ADAS_DATA_DIR)
+  set(M3DC1_ADAS_DATA_DIR "${CMAKE_BINARY_DIR}/adas/adf11"
+    CACHE PATH "Path to OpenADAS ADF11 data directory")
+endif()
+set(M3DC1_OPENADAS_SOURCE_DIR "${CMAKE_BINARY_DIR}/adas/source")
+set(M3DC1_OPENADAS_SOURCE_NAMES
+  xxrptn
+  xfelem
+  i4unit
+  i4fctn
+  xxword
+  xxcase
+  xxslen
+  xxdata_11
+)
+foreach(_openadas_source_name IN LISTS M3DC1_OPENADAS_SOURCE_NAMES)
+  list(APPEND M3DC1_OPENADAS_SOURCES
+    "${M3DC1_OPENADAS_SOURCE_DIR}/${_openadas_source_name}.f")
 endforeach()
-# message(STATUS "Found PETSc ${PETSC_LIBRARIES}")
-message(STATUS "Found PETSc ${PETSC_LINK_LIBRARIES}")
-# list(TRANSFORM PETSC_LINK_LIBRARIES REPLACE "libpthread\.a" "libpthread.so")
-# list(FILTER PETSC_LINK_LIBRARIES EXCLUDE REGEX "libpthread")
-list(FILTER PETSC_LINK_LIBRARIES EXCLUDE REGEX "libcuda\.so")
-# list(FILTER PETSC_LINK_LIBRARIES EXCLUDE REGEX "libsci_cray\.a")
-# list(REMOVE_DUPLICATES PETSC_LINK_LIBRARIES)
-set_property(TARGET PkgConfig::PETSC PROPERTY
-  INTERFACE_LINK_LIBRARIES "${PETSC_LINK_LIBRARIES}")
-
-pkg_search_module(FFTW fftw3 IMPORTED_TARGET REQUIRED)
-
-pkg_search_module(GSL gsl IMPORTED_TARGET REQUIRED)
-# list(FILTER GSL_LINK_LIBRARIES EXCLUDE REGEX "libm")
-# list(TRANSFORM GSL_LINK_LIBRARIES REPLACE "libm\.a" "libm.so")
-message(STATUS "Found gsl ${GSL_LINK_LIBRARIES}")
-# set_property(TARGET PkgConfig::GSL PROPERTY
-#   INTERFACE_LINK_LIBRARIES "${GSL_LINK_LIBRARIES}")
-
-#find_package(TRILINOS)
-# if(TRILINOS_FOUND AND USETRILINOS)
-#   list(APPEND DEFS -DM3DC1_TRILINOS)
-#   set(POSTFIX "${POSTFIX}-trilinos")
-# endif()
-
-find_package(Scorec)
-# if(SCOREC_FOUND AND USESCOREC)
-#   if(TRILINOS_FOUND AND USETRILINOS)
-#     set(SCOREC_SOURCES
-#         scorec_mesh.f90
-#         scorec_vector.f90
-#         scorec_matrix.f90)
-#   endif()
-#   if(PETSC_FOUND AND USESCOREC)
-#     set(SCOREC_SOURCES
-#         scorec_mesh.f90
-#         scorec_vector.f90
-#         scorec_matrix.f90
-#         PETScInterface.cpp)
-#   endif()
-#   set(SOURCES "${SOURCES};${SCOREC_SOURCES}")
-#   list(APPEND DEFS -DUSESCOREC -Dmesh_mod=scorec_mesh_mod -Dvector_mod=scorec_vector_mod -Dmatrix_mod=scorec_matrix_mod -Dmatrix_type=scorec_matrix -Dvector_type=scorec_vector)
-#   if(USERW)
-#     list(APPEND DEFS -DUSERW)
-#     set(POSTFIX "${POSTFIX}-rw")
-#   else()
-#     list(APPEND DEFS -Dglobalinsertval=insertval -Dglobalentdofs=entdofs)
-#   endif()
-# endif()
-
-# if(MATRIX_DUMP)
-#   list(APPEND DEFS -DMATRIX_DUMP)
-# endif()
-
-find_package(Zoltan)
-#find_package(PARMETIS)
-#find_package(METIS)
-
-find_package(ZLIB)
-message(STATUS "ZLIB found: ${ZLIB_LIBRARIES}")
-find_package(HDF5 REQUIRED COMPONENTS Fortran HL)
-message(STATUS "HDF5 found: ${HDF5_HL_LIBRARIES}")
-if (HDF5_HL_LIBRARIES MATCHES "libhdf5_hl_fortran\.a")
-  set(HDF5_USE_STATIC_LIBRARIES TRUE)
-  find_package(HDF5 REQUIRED COMPONENTS Fortran HL)
-  message(STATUS "HDF5 found: ${HDF5_HL_LIBRARIES}")
+set_source_files_properties(${M3DC1_OPENADAS_SOURCES} PROPERTIES GENERATED TRUE)
+set(M3DC1_OPENADAS_DATA_FILES
+  "${M3DC1_ADAS_DATA_DIR}/scd85/scd85_ar.dat"
+  "${M3DC1_ADAS_DATA_DIR}/scd89/scd89_b.dat"
+  "${M3DC1_ADAS_DATA_DIR}/plt89/plt89_b.dat"
+  "${M3DC1_ADAS_DATA_DIR}/plt89/plt89_ar.dat"
+  "${M3DC1_ADAS_DATA_DIR}/scd96/scd96_he.dat"
+  "${M3DC1_ADAS_DATA_DIR}/scd96/scd96_be.dat"
+  "${M3DC1_ADAS_DATA_DIR}/scd96/scd96_c.dat"
+  "${M3DC1_ADAS_DATA_DIR}/scd96/scd96_ne.dat"
+  "${M3DC1_ADAS_DATA_DIR}/plt96/plt96_he.dat"
+  "${M3DC1_ADAS_DATA_DIR}/plt96/plt96_be.dat"
+  "${M3DC1_ADAS_DATA_DIR}/plt96/plt96_c.dat"
+  "${M3DC1_ADAS_DATA_DIR}/plt96/plt96_ne.dat"
+)
+set(M3DC1_OPENADAS_STAMP "${CMAKE_BINARY_DIR}/adas/.openadas.stamp")
+add_custom_command(
+  OUTPUT
+    "${M3DC1_OPENADAS_STAMP}"
+    ${M3DC1_OPENADAS_SOURCES}
+    ${M3DC1_OPENADAS_DATA_FILES}
+  COMMAND "${CMAKE_COMMAND}"
+    "-DOUT_DIR=${CMAKE_BINARY_DIR}/adas"
+    "-DADF11_DIR=${M3DC1_ADAS_DATA_DIR}"
+    -P "${PROJECT_SOURCE_DIR}/cmake/DownloadOpenADAS.cmake"
+  COMMAND "${CMAKE_COMMAND}" -E touch "${M3DC1_OPENADAS_STAMP}"
+  COMMENT "Downloading OpenADAS ADF11 data and reader sources"
+  VERBATIM)
+set(_download_open_adas_all)
+if(M3DC1_DOWNLOAD_ADAS_DATA)
+  set(_download_open_adas_all ALL)
 endif()
-if (HDF5_FOUND)
-    message(STATUS "HDF5 found: ${HDF5_VERSION}")
-else ()
-    message(FATAL_ERROR "HDF5 not found.")
-endif ()
-
-pkg_search_module(NETCDF netcdf IMPORTED_TARGET REQUIRED)
-message(STATUS "Found NetCDF ${NETCDF_VERSION}")
-pkg_search_module(NETCDF-Fortran netcdf-fortran IMPORTED_TARGET REQUIRED)
-message(STATUS "Found NetCDF-Fortran ${NETCDF-Fortran_VERSION}")
-
-# find_package(netCDF-Fortran REQUIRED)
-# if (netCDF_FOUND)
-#     message(STATUS "netCDF found: ${netCDF_VERSION}")
-# else ()
-#     message(FATAL_ERROR "netCDF not found.")
-# endif ()
+add_custom_target(download_open_adas ${_download_open_adas_all}
+  DEPENDS "${M3DC1_OPENADAS_STAMP}")
 
-#find_package(ZLIB)
-# find_package(GSL REQUIRED)
-# find_package(FFTW REQUIRED)
+if(M3DC1_ENABLE_ADAS)
+  list(APPEND M3DC1_LIBRARY_SOURCES adas_m3dc1.f90)
+  add_library(openadas_reader STATIC ${M3DC1_OPENADAS_SOURCES})
+  add_dependencies(openadas_reader download_open_adas)
+endif()
 
 find_package(MPI REQUIRED)
-get_property(MPI_LINK_LIBRARIES TARGET MPI::MPI_Fortran PROPERTY INTERFACE_LINK_LIBRARIES)
-list(FILTER MPI_LINK_LIBRARIES EXCLUDE REGEX "libcuda\.so")
-set_property(TARGET MPI::MPI_Fortran PROPERTY
-  INTERFACE_LINK_LIBRARIES "${MPI_LINK_LIBRARIES}")
+find_package(HDF5 COMPONENTS Fortran HL REQUIRED)
+find_package(GSL REQUIRED)
+find_package(LAPACK REQUIRED)
+find_path(FFTW3_INCLUDE_DIR NAMES fftw3.h REQUIRED)
+find_library(FFTW3_LIBRARY NAMES fftw3 REQUIRED)
+
+if(M3DC1_ENABLE_ST)
+  find_library(NETCDF_C_LIBRARY NAMES netcdf
+    HINTS ${NetCDF_ROOT}
+    PATH_SUFFIXES lib lib64
+    REQUIRED)
+endif()
 
-set(M3DC1_SCOREC_LIB_NAME m3dc1_scorec)
-if(ENABLE_3D)
-  set(BIN_NAME "${BIN_NAME}_3d")
-  list(APPEND DEFS -DUSE3D -Dvectype=real -DMAX_PTS=125)
-  if(ENABLE_ST)
-    set(BIN_NAME "${BIN_NAME}_st")
-    list(APPEND DEFS -DUSEST)
-    # if(PAR)
-    #   set(POSTFIX "${POSTFIX}-pix")
-    # else()
-    #   set(POSTFIX "${POSTFIX}-complex")
-    # endif()
-  else()
-    list(APPEND DEFS -Dvectype=real)
-  endif()
-  # set(POSTFIX "${POSTFIX}-125")
+if(PETSC_VERSION MATCHES "^([0-9]+)\\.([0-9]+)")
+  math(EXPR M3DC1_PETSC_VERSION_CODE "${CMAKE_MATCH_1} * 10 + ${CMAKE_MATCH_2}")
 else()
-  set(BIN_NAME "${BIN_NAME}_2d")
-  list(APPEND DEFS -DMAX_PTS=25)
-  # set(POSTFIX "${POSTFIX}-25")
-  if(ENABLE_COMPLEX)
-    set(BIN_NAME "${BIN_NAME}_complex")
-    list(APPEND DEFS -Dvectype=complex -DUSECOMPLEX)
-    set(M3DC1_SCOREC_LIB_NAME ${M3DC1_SCOREC_LIB_NAME}_complex)
-    # if(PAR)
-    #   set(POSTFIX "${POSTFIX}-pix")
-    # else()
-    #   set(POSTFIX "${POSTFIX}-complex")
-    # endif()
-  else()
-    list(APPEND DEFS -Dvectype=real)
-  endif()
+  set(M3DC1_PETSC_VERSION_CODE 0)
 endif()
-if(ENABLE_PARTICLE)
-  list(APPEND DEFS -DUSEPARTICLES)
-endif()
-
-# if(USEOMP)
-#   set(POSTFIX "${POSTFIX}-omp")
-# endif()
-
-execute_process(COMMAND date OUTPUT_VARIABLE DATE)
-string(REGEX REPLACE "\n$" "" DATE "${DATE}")
-message(STATUS "Build date : ${DATE}")
-execute_process(COMMAND cat ${CMAKE_CURRENT_SOURCE_DIR}/release_version OUTPUT_VARIABLE VERSION)
-string(REGEX REPLACE "\n$" "" VERSION "${VERSION}")
-message(STATUS "Release version : ${VERSION}")
-list(APPEND DEFS -DUSESCOREC -Dmesh_mod=scorec_mesh_mod -Dvector_mod=scorec_vector_mod -Dmatrix_mod=scorec_matrix_mod -Dmatrix_type=scorec_matrix -Dvector_type=scorec_vector -DPETSC_VERSION=990 -DUSEBLAS -DDATE_BUILT='${DATE}' -DRELEASE_VERSION='${VERSION}')
-
-# if(NOT ARCH)
-#   if(M3DC1_ARCH)
-#     set(ARCH ${M3DC1_ARCH})
-#   else()
-#     set(ARCH ${HOST})
-#   endif()
-# endif()
-# set(POSTFIX "${ARCH}${POSTFIX}")
-
-# include(cmake/${ARCH}.cmake)
-
-#set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-none")
-
-#todo : build depency list as we go instead of adding all at the end (trilinos, etc)
-
-# if(COM)
-#   set(BIN_NAME "${BIN_NAME}_complex")
-# endif()
-
-set_source_files_properties(
-  scorec_vector.f90
-  PROPERTIES
-  OBJECT_DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/${BIN_NAME}.dir/scorec_mesh.f90.o
+string(TIMESTAMP M3DC1_DATE_BUILT "%a %b %d %H:%M:%S %Y")
+file(READ "${CMAKE_CURRENT_SOURCE_DIR}/release_version" M3DC1_RELEASE_VERSION)
+string(STRIP "${M3DC1_RELEASE_VERSION}" M3DC1_RELEASE_VERSION)
+
+add_library(m3dc1_library ${M3DC1_LIBRARY_SOURCES})
+
+target_include_directories(m3dc1_library PRIVATE ${FFTW3_INCLUDE_DIR})
+target_link_libraries(m3dc1_library PUBLIC
+  m3dc1::scorec
+  MPI::MPI_Fortran
+  PETSc::PETSc
+  hdf5::hdf5_fortran
+  hdf5::hdf5_hl_fortran
+  GSL::gsl
+  LAPACK::LAPACK
+  ${FFTW3_LIBRARY}
+)
+if(M3DC1_ENABLE_ST)
+  target_link_libraries(m3dc1_library PUBLIC
+    ${NETCDF_C_LIBRARY}
   )
-
-set_source_files_properties(
-  scorec_matrix.f90
-  PROPERTIES
-  OBJECT_DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/${BIN_NAME}.dir/scorec_vector.f90.o
-  )
-
-set_source_files_properties(
-  field.f90
-  PROPERTIES
-  OBJECT_DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/${BIN_NAME}.dir/scorec_vector.f90.o
-  OBJECT_DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/${BIN_NAME}.dir/scorec_matrix.f90.o
-  )
-
-set_source_files_properties(
-  M3Dmodules.f90
-  PROPERTIES
-  OBJECT_DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/${BIN_NAME}.dir/scorec_mesh.f90.o
-  )
-
-set_source_files_properties(
-  nintegrate_mod.f90
-  PROPERTIES
-  OBJECT_DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/${BIN_NAME}.dir/scorec_mesh.f90.o
-  )
-
-if (CMAKE_Fortran_COMPILER_ID STREQUAL "Cray")
-  set(CMAKE_Fortran_FLAGS_RELEASE "-O2")
-  set_source_files_properties( scorec_mesh.f90 PROPERTIES COMPILE_FLAGS -O0 )
-  set_source_files_properties( transport.f90 PROPERTIES COMPILE_FLAGS -O0 )
-  set_source_files_properties( particle.f90 PROPERTIES COMPILE_FLAGS -O1 )
-  set_source_files_properties( ludef_t.f90 PROPERTIES COMPILE_FLAGS -O1 )
-  set_source_files_properties( ludef_t_gpu.f90 PROPERTIES COMPILE_FLAGS -O1 )
 endif()
-add_executable(${BIN_NAME} ${SOURCES})
-# set_target_properties(${BIN_NAME} PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/${BIN_NAME}_mods)
-set_property(TARGET ${BIN_NAME} PROPERTY LINKER_LANGUAGE Fortran)
-target_include_directories(${BIN_NAME} PRIVATE ${SCOREC_INCLUDE_DIRS})
-target_link_libraries(${BIN_NAME} PRIVATE PkgConfig::NETCDF-Fortran PkgConfig::NETCDF
-  ${HDF5_HL_LIBRARIES} ${HDF5_LIBRARIES}
- # netCDF::netcdff
-  ${M3DC1_SCOREC_LIB_NAME} ${SCOREC_LIBRARIES} ${ZOLTAN_LIBRARIES}
-  PkgConfig::PETSC PkgConfig::FFTW PkgConfig::GSL MPI::MPI_Fortran stdc++)
-# target_link_directories(m3dc1_2d ${PETSC_STATIC_LIBRARY_DIRECTORIES})
-target_compile_definitions(${BIN_NAME} PRIVATE ${DEFS})
-if (CMAKE_Fortran_COMPILER_ID STREQUAL "Cray")
-  target_compile_options(${BIN_NAME} PRIVATE $<$<COMPILE_LANGUAGE:Fortran>:-s real64>)
-elseif (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
-  target_compile_options(${BIN_NAME} PRIVATE $<$<COMPILE_LANGUAGE:Fortran>:-fdefault-real-8>)
-  if (CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER_EQUAL 10)
-  target_compile_options(${BIN_NAME} PRIVATE $<$<COMPILE_LANGUAGE:Fortran>:-fallow-argument-mismatch>)
-  endif()
-  target_compile_options(${BIN_NAME} PRIVATE $<$<COMPILE_LANGUAGE:Fortran>:-ffree-line-length-512>)
+if(M3DC1_ENABLE_ADAS)
+  target_link_libraries(m3dc1_library PRIVATE openadas_reader)
+endif()
+target_compile_definitions(m3dc1_library PRIVATE
+  USESCOREC
+  ${M3DC1_MESH_DEFINE}
+  ${M3DC1_VECTOR_DEFINE}
+  ${M3DC1_MATRIX_DEFINE}
+  vector_type=${M3DC1_VECTOR_TYPE}
+  matrix_type=${M3DC1_MATRIX_TYPE}
+  PETSC_VERSION=${M3DC1_PETSC_VERSION_CODE}
+  DATE_BUILT="${M3DC1_DATE_BUILT}"
+  RELEASE_VERSION="${M3DC1_RELEASE_VERSION}"
+)
+if(M3DC1_ENABLE_3D)
+  target_compile_definitions(m3dc1_library PRIVATE USE3D MAX_PTS=125)
 else()
-  target_compile_options(${BIN_NAME} PRIVATE $<$<COMPILE_LANGUAGE:Fortran>:-r8>)
+  target_compile_definitions(m3dc1_library PRIVATE MAX_PTS=25)
 endif()
-install(TARGETS ${BIN_NAME} DESTINATION ".")
-# get_property(TTT TARGET m3dc1_3d PROPERTY
-#   INTERFACE_LINK_LIBRARIES)
-list(FILTER CMAKE_C_IMPLICIT_LINK_LIBRARIES EXCLUDE REGEX "^c$")
-list(FILTER CMAKE_CXX_IMPLICIT_LINK_LIBRARIES EXCLUDE REGEX "^c$")
-message(STATUS ${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES})
-message(STATUS ${CMAKE_C_IMPLICIT_LINK_LIBRARIES})
-message(STATUS ${CMAKE_CXX_IMPLICIT_LINK_LIBRARIES})
-# set(CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES "")
-# set(CMAKE_C_IMPLICIT_LINK_LIBRARIES "")
-# set(CMAKE_CXX_IMPLICIT_LINK_LIBRARIES "")
-
-if(ENABLE_OPENMP)
-    find_package(OpenMP REQUIRED)
-    # Add OpenMP flags for C and Fortran
-    # if(OpenMP_C_FOUND)
-    #   target_link_libraries(${BIN_NAME} PRIVATE OpenMP::OpenMP_C)
-    # endif()
-
-    if(OpenMP_Fortran_FOUND)
-      target_link_libraries(${BIN_NAME} PRIVATE OpenMP::OpenMP_Fortran)
-    endif()
-
-    message(STATUS "OpenMP enabled for C and Fortran")
+if(M3DC1_ENABLE_ST)
+  target_compile_definitions(m3dc1_library PRIVATE USEST)
+endif()
+if(M3DC1_ENABLE_COMPLEX)
+  target_compile_definitions(m3dc1_library PRIVATE vectype=complex USECOMPLEX)
+else()
+  target_compile_definitions(m3dc1_library PRIVATE vectype=real)
+endif()
+if(M3DC1_ENABLE_OPENMP)
+  target_link_libraries(m3dc1_library PUBLIC OpenMP::OpenMP_Fortran)
+endif()
+if(M3DC1_ENABLE_PARTICLE)
+  target_compile_definitions(m3dc1_library PRIVATE USEPARTICLES)
+endif()
+if(M3DC1_ENABLE_ADAS)
+  target_compile_definitions(m3dc1_library PRIVATE USEADAS)
+endif()
+set_target_properties(m3dc1_library PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+
+target_compile_options(m3dc1_library PUBLIC
+  # GNU
+  $<$<COMPILE_LANG_AND_ID:Fortran,GNU>:-cpp>
+  $<$<COMPILE_LANG_AND_ID:Fortran,GNU>:-fdefault-real-8>
+  $<$<COMPILE_LANG_AND_ID:Fortran,GNU>:-fdefault-double-8>
+  $<$<COMPILE_LANG_AND_ID:Fortran,GNU>:-fallow-argument-mismatch>
+  $<$<COMPILE_LANG_AND_ID:Fortran,GNU>:-ffree-line-length-none>
+
+  # Intel Classic
+  $<$<COMPILE_LANG_AND_ID:Fortran,Intel>:-fpp>
+  $<$<COMPILE_LANG_AND_ID:Fortran,Intel>:-r8>
+
+  # IntelLLVM (ifx)
+  $<$<COMPILE_LANG_AND_ID:Fortran,IntelLLVM>:-fpp>
+  $<$<COMPILE_LANG_AND_ID:Fortran,IntelLLVM>:-r8>
+
+  # LLVM Flang
+  $<$<COMPILE_LANG_AND_ID:Fortran,LLVMFlang>:-cpp>
+  $<$<COMPILE_LANG_AND_ID:Fortran,LLVMFlang>:-fdefault-real-8>
+  $<$<COMPILE_LANG_AND_ID:Fortran,LLVMFlang>:-fdefault-double-8>
+
+  # Cray Compiling Environment
+  $<$<COMPILE_LANG_AND_ID:Fortran,Cray>:-eZ>
+  $<$<COMPILE_LANG_AND_ID:Fortran,Cray>:-s>
+  $<$<COMPILE_LANG_AND_ID:Fortran,Cray>:real64>
+  $<$<COMPILE_LANG_AND_ID:Fortran,Cray>:-ef>
+)
+
+# --- Build main executable name per Makefile convention ---
+set(M3DC1_MAIN_NAME m3dc1)
+if(M3DC1_ENABLE_3D)
+  string(APPEND M3DC1_MAIN_NAME "_3d")
 else()
-    message(STATUS "OpenMP not enabled")
+  string(APPEND M3DC1_MAIN_NAME "_2d")
+endif()
+if(M3DC1_ENABLE_TRILINOS)
+  string(APPEND M3DC1_MAIN_NAME "_trilinos")
+endif()
+if(M3DC1_ENABLE_COMPLEX)
+  string(APPEND M3DC1_MAIN_NAME "_complex")
+endif()
+if(M3DC1_ENABLE_OPENMP)
+  string(APPEND M3DC1_MAIN_NAME "_omp")
+endif()
+if(M3DC1_ENABLE_ST)
+  string(APPEND M3DC1_MAIN_NAME "_st")
 endif()
+if(M3DC1_ENABLE_PARTICLE)
+  string(APPEND M3DC1_MAIN_NAME "_pic")
+endif()
+
+add_executable(${M3DC1_MAIN_NAME} $<TARGET_OBJECTS:m3dc1_library>)
+target_link_libraries(${M3DC1_MAIN_NAME} PRIVATE m3dc1_library)
 
-#find_package(PSPLINE)
-#add_executable(readgato polar.f90 readgato.f90)
-#target_link_libraries(readgato ${PSPLINE_LIBRARIES})
+add_custom_target(m3dc1 ALL
+  COMMAND ${CMAKE_COMMAND} -E create_symlink ${M3DC1_MAIN_NAME} m3dc1
+  DEPENDS ${M3DC1_MAIN_NAME}
+  WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+  COMMENT "Creating symlink: m3dc1 -> ${M3DC1_MAIN_NAME}")
 
-#add_executable(read_jsolver polar.f90 read_jsolver_exec.f90)
-#target_link_libraries(read_jsolver ${PSPLINE_LIBRARIES})
+add_executable(a2cc a2cc.f90)
+target_link_libraries(a2cc PUBLIC m3dc1_library)
 
-#add_executable(a2cc readaeqdsk.f90 a2cc.f90)
-#target_link_libraries(a2cc )
-#target_compile_definitions(a2cc PUBLIC ${DEFS})
+add_executable(mackenbach_profiles mackenbach_profiles.f90)
 
-#add_executable(kprad kprad.f90 kprad_test.f90)
+# add_executable(kprad kprad.f90 kprad_test.f90)
 
+add_executable(test_region test_region.f90)
+target_link_libraries(test_region PUBLIC m3dc1_library)
 
 
+if(M3DC1_ENABLE_PSPLINE)
+  find_package(PSPLINE REQUIRED)
+  add_executable(readgato readgato.f90)
+  target_link_libraries(readgato PUBLIC m3dc1_library PSPLINE::PSPLINE)
+
+  add_executable(read_jsolver read_jsolver.f90 polar.f90)
+  target_compile_definitions(read_jsolver PRIVATE READ_JSOLVER)
+  target_link_libraries(read_jsolver PUBLIC PSPLINE::PSPLINE)
+endif()
diff --git a/unstructured/M3Dmodules.f90 b/unstructured/M3Dmodules.f90
index 8c7e0c958..4745107cf 100644
--- a/unstructured/M3Dmodules.f90
+++ b/unstructured/M3Dmodules.f90
@@ -821,12 +821,7 @@ end module m3dc1_omp
 
 !cj velocity equation profiling
 module m3dc1_vel_prof
-#if PETSC_VERSION >= 39
-#include <petsc/finclude/petscksp.h>
       use petscksp
-#else
-#include <petsc/finclude/petscsysdef.h>
-#endif
       implicit none
-      PetscLogStage  stageA,stageS
+      integer :: stageA, stageS
 end module m3dc1_vel_prof
diff --git a/unstructured/basic_mesh.f90 b/unstructured/basic_mesh.f90
index 77619454c..ea3dec8ae 100644
--- a/unstructured/basic_mesh.f90
+++ b/unstructured/basic_mesh.f90
@@ -287,33 +287,36 @@ end subroutine create_rectangular_mesh
 
   subroutine load_mesh
 
+    use petscmat
+    use petscis
+    use petscao
+    use petscsys
     use math
 
     implicit none
 
-#include "finclude/petsc.h"
 !#ifndef PETSC_31
-!#include "finclude/petscvec.h"
-!#include "finclude/petscmat.h"
-!#include "finclude/petscis.h"
+!#include <petsc/finclude/petscvec.h>
+!#include <petsc/finclude/petscmat.h>
+!#include <petsc/finclude/petscis.h>
 !#endif
-#include "finclude/petscis.h90"
+!#include <petsc/finclude/petscis.h>
 
     integer :: num_global_elements
     type(node_type), allocatable :: global_node(:)
     type(element_type), allocatable :: global_elm(:)
-    Mat :: connectivity, adjacency
-    MatPartitioning :: partitioning
+    type(tMat) :: connectivity, adjacency
+    type(tMatPartitioning) :: partitioning
     integer :: i, j, k, ierr, rank, size, itri
     integer :: min_node, max_node
     integer, allocatable :: local_id(:)
     logical, allocatable :: is_ghost(:)
-    PetscInt, pointer :: idx(:)
+    integer, pointer :: idx(:)
     integer, allocatable :: nodes_per_proc(:), nodes_per_proc_local(:)
     integer :: new_id(1)
-    IS :: global_numbering, node_distribution
-    AO :: ordering
-    PetscViewer :: pv
+    type(tIS) :: global_numbering, node_distribution
+    type(tAO) :: ordering
+    type(tPetscViewer) :: pv
 
     bb(1) = 0.185
     bb(2) = -1.45
@@ -386,7 +389,7 @@ subroutine load_mesh
        call MatDestroy(adjacency, ierr)
        call ISPartitioningToNumbering(node_distribution,global_numbering,ierr)
        
-       call AOCreateBasicIS(global_numbering, PETSC_NULL, ordering, ierr)
+       call AOCreateBasicIS(global_numbering, PETSC_NULL_IS, ordering, ierr)
        
        call PetscViewerASCIIOpen(PETSC_COMM_WORLD,'node_distribution',pv,ierr)
        call ISView(node_distribution, pv, ierr)
@@ -405,11 +408,11 @@ subroutine load_mesh
        allocate(nodes_per_proc_local(size), nodes_per_proc(size))
        nodes_per_proc_local = 0
 
-       call ISGetIndicesF90(node_distribution, idx, ierr)
+       call ISGetIndices(node_distribution, idx, ierr)
        do i=1, num_local_nodes
           nodes_per_proc_local(idx(i)+1) = nodes_per_proc_local(idx(i)+1) + 1
        end do
-       call ISRestoreIndicesF90(node_distribution, idx, ierr)
+       call ISRestoreIndices(node_distribution, idx, ierr)
        
        call mpi_allreduce(nodes_per_proc_local, nodes_per_proc, size, &
             MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr)
diff --git a/unstructured/cmake/FindParmetis.cmake b/unstructured/cmake/FindParmetis.cmake
deleted file mode 100644
index 0224617d9..000000000
--- a/unstructured/cmake/FindParmetis.cmake
+++ /dev/null
@@ -1,46 +0,0 @@
-# - Try to find parmetis
-# Once done this will define
-#  PARMETIS_FOUND - System has PARMETIS
-#  PARMETIS_INCLUDE_DIRS - The PARMETIS include directories
-#  PARMETIS_LIBRARIES - The libraries needed to use PARMETIS
-#  PARMETIS_DEFINITIONS - Compiler switches required for using PARMETIS
-
-set(PARMETIS_PREFIX "" CACHE STRING "ParMETIS install directory")
-
-find_path(PARMETIS_INCLUDE_DIR parmetis.h PATHS "${PARMETIS_PREFIX}/include")
-if(NOT EXISTS "${PARMETIS_INCLUDE_DIR}")
-  message(FATAL_ERROR "parmetis include dir not found")
-endif()
-
-find_path(METIS_INCLUDE_DIR metis.h PATHS "${PARMETIS_PREFIX}/include")
-if(NOT EXISTS "${METIS_INCLUDE_DIR}")
-  message(FATAL_ERROR "metis include dir not found")
-endif()
-
-
-find_library(PARMETIS_LIBRARY parmetis PATHS "${PARMETIS_PREFIX}/lib")
-if(NOT EXISTS "${PARMETIS_LIBRARY}")
-  message(FATAL_ERROR "parmetis library not found")
-endif()
-
-find_library(METIS_LIBRARY metis PATHS "${PARMETIS_PREFIX}/lib")
-if(NOT EXISTS "${METIS_LIBRARY}")
-  message(FATAL_ERROR "metis library not found")
-endif()
-
-find_library(GK_LIBRARY GKlib PATHS "${PARMETIS_PREFIX}/lib")
-if(EXISTS "${GK_LIBRARY}")
-	set(PARMETIS_LIBRARIES ${PARMETIS_LIBRARY} ${METIS_LIBRARY} ${GK_LIBRARY})
-else()
-	set(PARMETIS_LIBRARIES ${PARMETIS_LIBRARY} ${METIS_LIBRARY})
-endif()
-
-set(PARMETIS_INCLUDE_DIRS ${PARMETIS_INCLUDE_DIR} ${METIS_INCLUDE_DIR})
-
-include(FindPackageHandleStandardArgs)
-# handle the QUIETLY and REQUIRED arguments and set PARMETIS_FOUND to TRUE
-# if all listed variables are TRUE
-find_package_handle_standard_args(Parmetis  DEFAULT_MSG
-    PARMETIS_LIBRARY METIS_LIBRARY PARMETIS_INCLUDE_DIR)
-
-mark_as_advanced(PARMETIS_INCLUDE_DIR PARMETIS_LIBRARY METIS_LIBRARY)
diff --git a/unstructured/cmake/FindScorec.cmake b/unstructured/cmake/FindScorec.cmake
deleted file mode 100644
index 4f6e8c45d..000000000
--- a/unstructured/cmake/FindScorec.cmake
+++ /dev/null
@@ -1,106 +0,0 @@
-# - Try to find SCOREC PUMI libraries
-# Once done this will define
-#  SCOREC_FOUND - System has SCOREC
-#  SCOREC_INCLUDE_DIRS - The SCOREC include directories
-#  SCOREC_LIBRARIES - The libraries needed to use SCOREC
-#  SCOREC_DEFINITIONS - Compiler switches required for using SCOREC
-#
-# This implementation assumes a SCOREC install has the following structure
-# VERSION/
-#         include/*.h
-#         lib/*.a
-
-macro(scorecLibCheck libs isRequired)
-  foreach(lib ${libs}) 
-    unset(scoreclib CACHE)
-    find_library(scoreclib "${lib}" PATHS ${SCOREC_LIB_DIR})
-    if(scoreclib MATCHES "^scoreclib-NOTFOUND$")
-      if(${isRequired})
-        message(FATAL_ERROR "SCOREC library ${lib} not found in ${SCOREC_LIB_DIR}")
-      else()
-        message("SCOREC library ${lib} not found in ${SCOREC_LIB_DIR}")
-      endif()
-    else()
-      set("SCOREC_${lib}_FOUND" TRUE CACHE INTERNAL "SCOREC library present")
-      set(SCOREC_LIBS ${SCOREC_LIBS} ${scoreclib})
-    endif()
-  endforeach()
-endmacro(scorecLibCheck)
-
-if(CMAKE_CXX_COMPILER_ID MATCHES "PGI")
-  set(ENABLE_ZOLTAN OFF)
-endif()
- 
-find_library(METIS_LIBRARY metis)
-if (NOT EXISTS "${METIS_LIBRARY}")
-  message(FATAL ERROR "METIS library not found")
-endif()
-
-set(SCOREC_LIBS "")
-if(ENABLE_SIMMETRIX)
-set(SCOREC_LIB_NAMES
-  pumi
-  ma
-  mds
-  apf_sim
-  apf_zoltan
-  apf
-  parma
-  gmi
-  gmi_sim
-  mth
-  pcu
-  sam
-  spr
-  crv
-  lion
-  ph 
-  )
-else()
-set(SCOREC_LIB_NAMES
-pumi
-crv
-ph
-sam
-spr
-ma
-apf_zoltan
-parma
-mds
-apf
-lion
-mth
-gmi
-pcu
-  )
-endif()
-scorecLibCheck("${SCOREC_LIB_NAMES}" TRUE)
-
-find_path(SCOREC_INCLUDE_DIR 
-  NAMES apf.h PCU.h ma.h 
-  PATHS ${SCOREC_INCLUDE_DIR})
-if(NOT EXISTS "${SCOREC_INCLUDE_DIR}")
-  message(FATAL_ERROR "SCOREC include dir not found")
-endif()
-
-string(REGEX REPLACE 
-  "/include$" "" 
-  SCOREC_INSTALL_DIR
-  "${SCOREC_INCLUDE_DIR}")
-
-set(SCOREC_LIBRARIES ${SCOREC_LIBS} ${ZOLTAN_LIBRARY} ${PARMETIS_LIBRARY} ${METIS_LIBRARY})
-set(SCOREC_INCLUDE_DIRS ${SCOREC_INCLUDE_DIR})
-
-include(FindPackageHandleStandardArgs)
-# handle the QUIETLY and REQUIRED arguments and set PARMETIS_FOUND to TRUE
-# if all listed variables are TRUE
-find_package_handle_standard_args(Scorec  DEFAULT_MSG
-                                  SCOREC_LIBS SCOREC_INCLUDE_DIR)
-
-mark_as_advanced(SCOREC_INCLUDE_DIR SCOREC_LIBS ${ZOLTAN_LIBRARY} ${PARMETIS_LIBRARY} ${METIS_LIBRARY})
-
-set(SCOREC_LINK_LIBS "")
-foreach(lib ${SCOREC_LIB_NAMES})
-  set(SCOREC_LINK_LIBS "${SCOREC_LINK_LIBS} -l${lib}")
-endforeach()
-
diff --git a/unstructured/cmake/FindZoltan.cmake b/unstructured/cmake/FindZoltan.cmake
deleted file mode 100644
index 24d8ad654..000000000
--- a/unstructured/cmake/FindZoltan.cmake
+++ /dev/null
@@ -1,31 +0,0 @@
-# - Try to find zoltan
-# Once done this will define
-#  ZOLTAN_FOUND - System has ZOLTAN
-#  ZOLTAN_INCLUDE_DIRS - The ZOLTAN include directories
-#  ZOLTAN_LIBRARIES - The libraries needed to use ZOLTAN
-#  ZOLTAN_DEFINITIONS - Compiler switches required for using ZOLTAN
-
-set(ZOLTAN_PREFIX "${ZOLTAN_PREFIX_DEFAULT}" CACHE STRING "Zoltan install directory")
-if(ZOLTAN_PREFIX)
-  message(STATUS "ZOLTAN_PREFIX ${ZOLTAN_PREFIX}")
-endif()
-
-find_path(ZOLTAN_INCLUDE_DIR zoltan.h PATHS "${ZOLTAN_PREFIX}/include")
-
-find_library(ZOLTAN_LIBRARY zoltan PATHS "${ZOLTAN_PREFIX}/lib")
-
-set(ZOLTAN_LIBRARIES ${ZOLTAN_LIBRARY} )
-set(ZOLTAN_INCLUDE_DIRS ${ZOLTAN_INCLUDE_DIR} )
-
-find_package(Parmetis MODULE REQUIRED)
-
-include(FindPackageHandleStandardArgs)
-# handle the QUIETLY and REQUIRED arguments and set ZOLTAN_FOUND to TRUE
-# if all listed variables are TRUE
-find_package_handle_standard_args(
-    Zoltan
-    DEFAULT_MSG
-    ZOLTAN_LIBRARY ZOLTAN_INCLUDE_DIR
-)
-
-mark_as_advanced(ZOLTAN_INCLUDE_DIR ZOLTAN_LIBRARY )
diff --git a/unstructured/coils.f90 b/unstructured/coils.f90
index 917f1dd42..267039539 100644
--- a/unstructured/coils.f90
+++ b/unstructured/coils.f90
@@ -12,16 +12,15 @@ module coils
   !
   ! reads coil and current data from file
   !======================================================
- subroutine load_coils(xc, zc, ic, numcoils, coil_filename, current_filename, &
-      coil_mask, filaments)
-   use math
-   use read_ascii
+	 subroutine load_coils(xc, zc, ic, numcoils, coil_filename, current_filename, &
+	      coil_mask, filaments)
+	   use math
+	   use read_ascii
+	   use mpi
 
-   implicit none
+	   implicit none
 
-   include 'mpif.h'
-
-   real, intent(out), dimension(maxfilaments) :: xc, zc  ! coordinates of each coil
+	   real, intent(out), dimension(maxfilaments) :: xc, zc  ! coordinates of each coil
    complex, intent(out), dimension(maxfilaments) :: ic   ! current in each coil
    integer, intent(out) :: numcoils                  ! number of coils read
    character*(*) :: coil_filename, current_filename  ! input files
@@ -160,15 +159,13 @@ end subroutine load_coils
  subroutine field_from_coils(xc, zc, ic, nc, f, ipole, ierr)
    use basic
    use field
-   use mesh_mod
-   use m3dc1_nint
-   use newvar_mod
-
-   implicit none
+	   use mesh_mod
+	   use m3dc1_nint
+	   use newvar_mod
 
-   include 'mpif.h'
+	   implicit none
 
-   real, intent(in), dimension(nc) :: xc, zc   ! array of coil positions
+	   real, intent(in), dimension(nc) :: xc, zc   ! array of coil positions
    complex, intent(in), dimension(nc) :: ic    ! array of coil currents
    integer, intent(in) :: nc                   ! number of coils
    type(field_type), intent(inout) :: f        ! poloidal flux field
diff --git a/unstructured/diagnostics.f90 b/unstructured/diagnostics.f90
index c45a335a3..07f7a22a3 100644
--- a/unstructured/diagnostics.f90
+++ b/unstructured/diagnostics.f90
@@ -7,6 +7,8 @@
 
 module diagnostics
 
+  use mpi
+
   implicit none
 
   real :: tflux0
@@ -129,7 +131,6 @@ subroutine distribute_timings
 
     implicit none
 
-    include 'mpif.h'
     integer :: ier
     integer, parameter :: num_scalars = 14
     real, dimension(num_scalars) :: vin, vout
@@ -296,8 +297,6 @@ subroutine distribute_scalars()
 
     implicit none
 
-    include 'mpif.h'
-
     integer, parameter :: num_scalars = 82
     integer :: ier
     double precision, dimension(num_scalars) :: temp, temp2
@@ -509,8 +508,6 @@ subroutine evaluate(x,phi,z,ans,fin,itri,ierr)
 
   implicit none
 
-  include 'mpif.h'
-
   integer, intent(inout) :: itri
   real, intent(in) :: x, phi, z
   type(field_type), intent(in) :: fin
@@ -725,8 +722,6 @@ subroutine calculate_scalars()
 
   implicit none
  
-  include 'mpif.h'
-
   integer :: itri, numelms, def_fields, ier
   integer :: is_edge(3)  ! is inode on boundary
   real :: n(2,3),tpifac,tpirzero, t0
@@ -1205,8 +1200,6 @@ subroutine calculate_Lor_vol()
 
   implicit none
  
-  include 'mpif.h'
-
   integer :: itri, numelms, ier
   integer :: is_edge(3)  ! is inode on boundary
   real :: tpifac,tpirzero
@@ -1365,8 +1358,6 @@ subroutine magaxis(xguess,zguess,psi,psim,imethod,ier)
 
   implicit none
 
-  include 'mpif.h'
-
   real, intent(inout) :: xguess, zguess
   type(field_type), intent(in) :: psi
   integer, intent(in) :: imethod
@@ -1576,8 +1567,6 @@ subroutine te_max(xguess,zguess,te,tem,imethod,ier)
 
   implicit none
 
-  include 'mpif.h'
-
   real, intent(inout) :: xguess, zguess
   type(field_type), intent(in) :: te
   integer, intent(in) :: imethod
@@ -1777,8 +1766,6 @@ subroutine te_max2(xguess,zguess,te,tem,imethod,ier)
 
   implicit none
 
-  include 'mpif.h'
-
   real, intent(inout) :: xguess, zguess
   type(field_type), intent(in) :: te
   integer, intent(in) :: imethod
@@ -1836,8 +1823,6 @@ subroutine te_max3(xguess,zguess,te,tem,imethod,ier)
 
   implicit none
 
-  include 'mpif.h'
-
   real, intent(inout) :: xguess, zguess
   type(field_type), intent(in) :: te
   real, intent(out) :: tem
@@ -1971,8 +1956,6 @@ subroutine te_max4(te,tem,ilin,ier)
 
   implicit none
 
-  include 'mpif.h'
-
   type(field_type), intent(in) :: te
   integer, intent(in) :: ilin
   real, intent(out) :: tem
@@ -2038,8 +2021,6 @@ subroutine lcfs(psi, test_wall, findx)
 
   implicit none
 
-  include 'mpif.h'
-
   type(field_type), intent(in) :: psi
   logical, intent(in), optional :: test_wall
   logical, intent(in), optional :: findx
@@ -2368,7 +2349,6 @@ subroutine calculate_ke()
   use math
 
   implicit none
-  include 'mpif.h'
   integer :: itri, numelms, def_fields
   real :: ke_N, ketotal, fac
   integer :: ier, k, l, numnodes, N, icounter_t
@@ -2735,7 +2715,6 @@ subroutine calculate_bh()
   use boundary_conditions
   use math
   implicit none
-  include 'mpif.h'
   integer :: itri, numelms, def_fields
   real:: bh_N, bhtotal, fac
   integer :: ier, k, l, numnodes, N, icounter_t
@@ -3216,8 +3195,6 @@ subroutine te_max_dev(xguess,zguess,te,tem,imethod,ier)
 
   implicit none
 
-  include 'mpif.h'
-
   real, intent(inout) :: xguess, zguess
   type(field_type), intent(in) :: te
   real, intent(out) :: tem
@@ -3490,8 +3467,6 @@ subroutine phi_int(x,z,ans,fin,itri,ierr)
 
   implicit none
 
-  include 'mpif.h'
-
   integer, intent(inout) :: itri
   real, intent(in) :: x, z
   type(field_type), intent(in) :: fin
diff --git a/unstructured/error_estimate.f90 b/unstructured/error_estimate.f90
index 8f1f25b87..5f7fe7fa8 100644
--- a/unstructured/error_estimate.f90
+++ b/unstructured/error_estimate.f90
@@ -1,6 +1,6 @@
 module error_estimate
+  use mpi
   use m3dc1_nint
-#include "mpif.h"
   integer, parameter :: EOP_1 = 1
   integer, parameter :: EOP_DX = 2
   integer, parameter :: EOP_DY = 3
diff --git a/unstructured/gradshafranov.f90 b/unstructured/gradshafranov.f90
index ed13376ed..3742e8578 100644
--- a/unstructured/gradshafranov.f90
+++ b/unstructured/gradshafranov.f90
@@ -126,6 +126,7 @@ subroutine gradshafranov_per()
 end subroutine gradshafranov_per
 
 subroutine coil_feedback(itnum)
+  use mpi
   use basic
   use arrays
   use coils
@@ -133,7 +134,6 @@ subroutine coil_feedback(itnum)
 
   implicit none
 
-  include 'mpif.h'
 
   integer, intent(in) :: itnum
 
@@ -1093,6 +1093,7 @@ end subroutine extend_pressure
 
 !============================================================
 subroutine gradshafranov_solve
+  use mpi
 
   use math
   use mesh_mod
@@ -1107,7 +1108,6 @@ subroutine gradshafranov_solve
   use model
 
   implicit none
-  include 'mpif.h'
 
   type(field_type) :: b1vecini_vec, b2vecini_vec
   type(field_type) :: b3vecini_vec, b4vecini_vec
@@ -1814,6 +1814,7 @@ subroutine gradshafranov_solve
 end subroutine gradshafranov_solve
 
 subroutine calculate_error(error, error2, psinew)
+  use mpi
   use basic
   use field
   use boundary_conditions
@@ -1821,7 +1822,6 @@ subroutine calculate_error(error, error2, psinew)
 
   implicit none
 
-  include 'mpif.h'
 
   real, intent(out) :: error, error2
   type(field_type), intent(in) :: psinew
@@ -1883,6 +1883,7 @@ end subroutine calculate_error
 ! constrain solution to have the specified current, etc.
 !============================================================
 subroutine calculate_gamma(g2, g3, g4)
+  use mpi
   use basic
   use mesh_mod
   use arrays
@@ -1891,7 +1892,6 @@ subroutine calculate_gamma(g2, g3, g4)
 
   implicit none
 
-  include 'mpif.h'
 
   real, intent(out) :: g2, g3, g4
 
@@ -2009,6 +2009,7 @@ end subroutine calculate_gamma
 ! sets jout_i =  <mu_i | -val*delta(R-x)delta(Z-z)> 
 ! ===========================================================
 subroutine deltafun(x,z,val,jout)
+  use mpi
 
   use mesh_mod
   use basic
@@ -2019,7 +2020,6 @@ subroutine deltafun(x,z,val,jout)
 
   implicit none
 
-  include 'mpif.h'
 
   type(element_data) :: d
   real, intent(in) :: x, z, val
@@ -2080,6 +2080,7 @@ end subroutine deltafun
 
 ! ===========================================================
 subroutine gaussianfun(x,z,val,denom,jout)
+  use mpi
 
   use mesh_mod
   use basic
@@ -2090,7 +2091,6 @@ subroutine gaussianfun(x,z,val,denom,jout)
 
   implicit none
 
-  include 'mpif.h'
 
   real, intent(in) :: x, z, val, denom
   type(field_type), intent(inout) :: jout
@@ -2371,6 +2371,7 @@ end subroutine fundef
 
 
 subroutine fundef2(error)
+  use mpi
 
   use basic
   use mesh_mod
@@ -2382,7 +2383,6 @@ subroutine fundef2(error)
 
   implicit none
 
-  include 'mpif.h'
 
   real, intent(out) :: error
 
@@ -2878,12 +2878,12 @@ end subroutine write_profile
  ! calculates the fractional error in the GS solution
  !=============================================================
  subroutine calculate_gs_error(error)
+   use mpi
    use basic
    use m3dc1_nint
    
    implicit none
    
-   include 'mpif.h'
    
    real, intent(out) :: error
    
diff --git a/unstructured/hdf5_output.f90 b/unstructured/hdf5_output.f90
index f20824d23..0b53643e3 100644
--- a/unstructured/hdf5_output.f90
+++ b/unstructured/hdf5_output.f90
@@ -30,12 +30,12 @@ end subroutine hdf5_flush
   ! hdf5_initialize
   ! ===============
   subroutine hdf5_initialize(restart, error)
+    use mpi
     use hdf5
     use basic
 
     implicit none
 
-    include 'mpif.h'
 
     logical, intent(in) :: restart     ! if true, do not overwrite file
     integer, intent(out) :: error
diff --git a/unstructured/init_common.f90 b/unstructured/init_common.f90
index bbce4f563..14cf48f29 100644
--- a/unstructured/init_common.f90
+++ b/unstructured/init_common.f90
@@ -669,6 +669,7 @@ end subroutine kinetic_eq
 #endif     
 
 subroutine nre_eq
+  use mpi
   use basic
   use arrays
   use diagnostics
@@ -679,7 +680,6 @@ subroutine nre_eq
   use pellet
 
   implicit none
-  include 'mpif.h'
 
   type(field_type) :: nre_vec
   integer :: itri, numelms, def_fields
diff --git a/unstructured/input.f90 b/unstructured/input.f90
index 7279b0859..806dee06a 100644
--- a/unstructured/input.f90
+++ b/unstructured/input.f90
@@ -97,12 +97,12 @@ end subroutine add_group
 ! reads input namelist
 !=========================
 subroutine input
+  use mpi
   use basic
 
   implicit none
 
   integer :: ierr
-#include "mpif.h"
 
   call set_defaults
 
@@ -1349,16 +1349,8 @@ subroutine validate_input
   use resistive_wall
   use kprad_m3dc1
 
-#if PETSC_VERSION >= 38
   use petsc
   implicit none
-#elif PETSC_VERSION >= 36
-  implicit none
-#include "petsc/finclude/petsc.h"
-#else
-  implicit none
-#include "finclude/petsc.h"
-#endif
 
   integer :: ier,i
   real :: de
@@ -1871,4 +1863,3 @@ subroutine validate_input
   end if
 
 end subroutine validate_input
-
diff --git a/unstructured/iterdb.f90 b/unstructured/iterdb.f90
index 9a93645eb..38260d88e 100644
--- a/unstructured/iterdb.f90
+++ b/unstructured/iterdb.f90
@@ -28,10 +28,10 @@ module iterdb
 contains
 
   subroutine load_iterdb(filename, ierr)
+    use mpi
 
     implicit none
 
-    include 'mpif.h'
 
     character(len=*), intent(in) :: filename
     integer, intent(out) :: ierr
diff --git a/unstructured/kprad.f90 b/unstructured/kprad.f90
index b3528ca96..d0071f05b 100644
--- a/unstructured/kprad.f90
+++ b/unstructured/kprad.f90
@@ -83,9 +83,9 @@ subroutine kprad_allocate(Z)
   end subroutine kprad_allocate
 
   subroutine kprad_rebase_dt()
+    use mpi
     implicit none
 
-    include 'mpif.h'
 
     integer :: ier
     real :: temp
diff --git a/unstructured/kprad_m3dc1.f90 b/unstructured/kprad_m3dc1.f90
index d14a0bfa4..f4752cc34 100644
--- a/unstructured/kprad_m3dc1.f90
+++ b/unstructured/kprad_m3dc1.f90
@@ -274,11 +274,10 @@ subroutine kprad_advect(dti)
     use boundary_conditions
     use metricterms_new
     use sparse
+    use mpi
 
     implicit none
 
-    include 'mpif.h'
-
     real, intent(in) :: dti
     type(matrix_type) :: nmat_lhs, nmat_rhs
     type(field_type) :: rhs
@@ -730,11 +729,10 @@ subroutine read_lp_source(filename, ierr)
     use pellet
     use newvar_mod
     use math
+    use mpi
 
     implicit none
 
-    include 'mpif.h'
-
     integer :: mpierr
     integer, intent(out) :: ierr
     character(len=*), intent(in) :: filename
@@ -873,11 +871,10 @@ subroutine deltafuns(n,x,phi,z,m,val,jout, ier)
   use field
   use m3dc1_nint
   use math
+  use mpi
 
   implicit none
 
-  include 'mpif.h'
-
   integer, intent(in) :: n, m
   real, intent(in), dimension(n) :: x, phi, z
   real, intent(in), dimension(n,m) :: val
diff --git a/unstructured/kprad_test.f90 b/unstructured/kprad_test.f90
index 062026da4..41ed3924b 100644
--- a/unstructured/kprad_test.f90
+++ b/unstructured/kprad_test.f90
@@ -4,8 +4,9 @@ program main
   implicit none
   integer :: zimp, i, j, ntimemax,ntime,nzones, ierr
   real :: dts, facimp, time
-  real, allocatable, dimension(:) :: ne,te,pbrem,dw_brem
-  real, allocatable, dimension(:,:) :: imp_rad,dw_rad,nz
+  real, allocatable, dimension(:) :: ne, te, den, ti, p, pbrem, dw_brem
+  real, allocatable, dimension(:,:) :: imp_rad, dw_rad, dw_ion, dw_reck, &
+       dw_recp, nz, source
   
   open(unit=10,file='rad.txt')
 
@@ -13,18 +14,31 @@ program main
   ZIMP = 10   ! NEON
   !      ZIMP = 18   ! Argon
   
+  ikprad = 1
+  ikprad_min_option = 1
+  ikprad_max_dt = 0
+  kprad_max_dt = -1.
+  ikprad_evolve_internal = 0
+
+  call kprad_allocate(zimp)
   call kprad_atomic_data_sub(zimp, ierr)
   dts = 1.e-7
   ntimemax = 100
   nzones = 1   ! one zone
   j = 1
 
-  allocate(ne(nzones),te(nzones))
+  allocate(ne(nzones), te(nzones), den(nzones), ti(nzones), p(nzones))
   ne(j) = 1.e14   ! density in cm-3
   te(j) = 1000.   ! temp in ev
+  den(j) = ne(j)
+  ti(j) = te(j)
+  p(j) = ne(j)*te(j)*1.6022e-12
   allocate(nz(nzones,0:zimp))
   allocate(pbrem(nzones), imp_rad(nzones,zimp+1), &
-       dw_brem(nzones), dw_rad(nzones,zimp+1))
+       dw_brem(nzones), dw_rad(nzones,zimp+1), &
+       dw_ion(nzones,0:zimp), dw_reck(nzones,0:zimp), &
+       dw_recp(nzones,0:zimp), source(nzones,0:zimp))
+  source = 0.
 
   !     initialize
   facimp = .01
@@ -38,8 +52,9 @@ program main
   do ntime=1, ntimemax
      time = time + dts
 
-     call kprad_advance_densities(dts, nzones, zimp, ne, te, nz, &
-          dw_rad, dw_brem)
+     call kprad_advance_densities(dts, nzones, zimp, p, ne, te, den, ti, &
+          nz, dw_rad, dw_brem, dw_ion, dw_reck, &
+          dw_recp, source)
 
      ! calculate power from radiated energy
      pbrem = dw_brem / dts
@@ -54,7 +69,8 @@ program main
 
   call kprad_deallocate()
 
-  deallocate(nz, pbrem, imp_rad, dw_brem, dw_rad)
+  deallocate(ne, te, den, ti, p, nz, pbrem, imp_rad, dw_brem, dw_rad, &
+       dw_ion, dw_reck, dw_recp, source)
 
   
   stop
diff --git a/unstructured/matrix_mod.F90 b/unstructured/matrix_mod.F90
new file mode 100644
index 000000000..8933c67c3
--- /dev/null
+++ b/unstructured/matrix_mod.F90
@@ -0,0 +1,14 @@
+! matrix_mod.F90 — thin wrapper that selects the matrix backend at compile time.
+! CMake defines exactly one of the following:
+!   M3DC1_MATRIX_IS_SCOREC  → use scorec_matrix_mod
+!   M3DC1_MATRIX_IS_PETSC   → use petsc_matrix_mod
+
+module matrix_mod
+#ifdef M3DC1_MATRIX_IS_SCOREC
+  use scorec_matrix_mod
+#elif defined(M3DC1_MATRIX_IS_PETSC)
+  use petsc_matrix_mod
+#else
+  use scorec_matrix_mod
+#endif
+end module matrix_mod
diff --git a/unstructured/mesh_mod.F90 b/unstructured/mesh_mod.F90
new file mode 100644
index 000000000..55b61522c
--- /dev/null
+++ b/unstructured/mesh_mod.F90
@@ -0,0 +1,14 @@
+! mesh_mod.F90 — thin wrapper that selects the mesh backend at compile time.
+! CMake defines exactly one of the following:
+!   M3DC1_MESH_IS_SCOREC   → use scorec_mesh_mod
+!   M3DC1_MESH_IS_BASIC    → use basic_mesh_mod
+
+module mesh_mod
+#ifdef M3DC1_MESH_IS_SCOREC
+  use scorec_mesh_mod
+#elif defined(M3DC1_MESH_IS_BASIC)
+  use basic_mesh_mod
+#else
+  use scorec_mesh_mod
+#endif
+end module mesh_mod
diff --git a/unstructured/newpar.f90 b/unstructured/newpar.f90
index d5f1976df..ecc2dc37e 100644
--- a/unstructured/newpar.f90
+++ b/unstructured/newpar.f90
@@ -34,16 +34,8 @@ Program Reducedquintic
   use openacc
 #endif
 
-#if PETSC_VERSION >= 38
   use petsc
   implicit none
-#elif PETSC_VERSION >= 36
-  implicit none
-#include "petsc/finclude/petsc.h"
-#else
-  implicit none
-#include "finclude/petsc.h"
-#endif
 
   integer :: ier, i, adapt_flag
   real :: tstart, tend, dtsave, t_solve, t_compute
@@ -582,16 +574,8 @@ subroutine safestop(iarg)
   use particles
   use resistive_wall
 
-#if PETSC_VERSION >= 38
   use petsc
   implicit none
-#elif PETSC_VERSION >= 36
-  implicit none
-#include "petsc/finclude/petsc.h"
-#else
-  implicit none
-#include "finclude/petsc.h"
-#endif
       
   integer, intent(in) :: iarg
   integer :: ier
@@ -1030,13 +1014,13 @@ end subroutine rotation
   ! populates the *tri arrays
   !============================================================
   subroutine tridef
+    use mpi
     use basic
     use math
     use mesh_mod
 
     implicit none
 
-    include 'mpif.h'
   
     type(element_data) :: d
     integer :: itri, i, j, k, ii, jj, numelms, numnodes, ndofs, ierr
@@ -1599,12 +1583,12 @@ subroutine calculate_qdfac(itri, z)
 end subroutine calculate_qdfac
 
 subroutine print_normal_curv()
+  use mpi
 !  use mpi
   use basic
   use mesh_mod
   implicit none
 
-    include 'mpif.h'
 
   integer :: ierr, i, icounter_t, numnodes
   integer :: izone, izonedim
diff --git a/unstructured/output.f90 b/unstructured/output.f90
index 3e6557ef0..fb1a2e158 100644
--- a/unstructured/output.f90
+++ b/unstructured/output.f90
@@ -59,6 +59,7 @@ end subroutine finalize_output
   !
   ! ======================================================================
   subroutine marker
+    use mpi
     use basic
     use hdf5_output
     use diagnostics
@@ -66,7 +67,6 @@ subroutine marker
 
     implicit none
 
-    include 'mpif.h'
 
     integer :: ier,i
     call mark_fields(0);
@@ -80,6 +80,7 @@ end subroutine marker
   ! writes output and restart files
   ! ======================================================================
   subroutine output
+    use mpi
     use basic
     use hdf5_output
     use diagnostics
@@ -89,7 +90,6 @@ subroutine output
 
     implicit none
 
-    include 'mpif.h'
 
     integer :: ier
     real :: tstart, tend, diff, gamma_std, gamma_mean
@@ -665,13 +665,13 @@ end subroutine hdf5_write_timings
 ! hdf5_write_time_slice
 ! =====================
 subroutine hdf5_write_time_slice(equilibrium, error)
+  use mpi
   use hdf5
   use hdf5_output
   use basic
 
   implicit none
 
-  include 'mpif.h'
   
   integer, intent(out) :: error
   integer, intent(in) :: equilibrium
diff --git a/unstructured/particle.f90 b/unstructured/particle.f90
index 2041b5352..b220537c1 100644
--- a/unstructured/particle.f90
+++ b/unstructured/particle.f90
@@ -136,9 +136,9 @@ module particles
 ! Note: any changes to the "particle" user-defined datatype must be reflected
 !       in the definitions of pnvars, pblklen, ptyps, and pdspls below.
 subroutine define_mpi_particle(ierr)
+  use mpi
    implicit none
 
-   include 'mpif.h'
 
    integer, intent(out) :: ierr
    integer, parameter :: pnvars = 15
@@ -189,9 +189,9 @@ subroutine define_mpi_particle(ierr)
 end subroutine define_mpi_particle
 
 subroutine define_mpi_elfield(ierr)
+  use mpi
    implicit none
 
-   include 'mpif.h'
 
    integer, intent(out) :: ierr
    integer, parameter :: pnvars = 24
@@ -264,13 +264,13 @@ end subroutine define_mpi_elfield
 
 !---------------------------------------------------------------------------
 subroutine particle_test
+  use mpi
    use basic
    use diagnostics
    use auxiliary_fields
    use m3dc1_nint
    use arrays
    implicit none
-   include 'mpif.h'
 
    integer, parameter :: trunit = 120
    real :: tstart, tend
@@ -382,6 +382,7 @@ real function tri_area(x1, z1, x2, z2, x3, z3)
 end function tri_area
 !---------------------------------------------------------------------------
 subroutine init_particles(lrestart, ierr)
+  use mpi
    use basic
    use arrays
    use hdf5_output
@@ -392,7 +393,6 @@ subroutine init_particles(lrestart, ierr)
    use boundary_conditions
 
    implicit none
-   include 'mpif.h'
 
    logical, intent(in) :: lrestart
    integer, intent(out) :: ierr
@@ -1914,13 +1914,13 @@ subroutine fdot(x, v, w, dxdt, dvdt, dwdt, dEpdt, itri, kel, f00, ierr, sps, B00
 end subroutine fdot
 
 subroutine particle_scaleback(scalefac)
+  use mpi
    use basic
    use field
    use arrays
    use mesh_mod
    implicit none
 
-   include 'mpif.h'
 
    vectype, intent(in) :: scalefac
    integer :: ipart, ielm, ierr
@@ -1951,10 +1951,10 @@ subroutine particle_scaleback(scalefac)
 end subroutine particle_scaleback
 !---------------------------------------------------------------------------
 subroutine delete_particle(exchange)
+  use mpi
    use basic
    use diagnostics
    implicit none
-   include 'mpif.h'
 
    logical, intent(in) :: exchange
    !integer, intent(in) :: ipart
@@ -2029,11 +2029,11 @@ end subroutine delete_particle
 
 !---------------------------------------------------------------------------
 subroutine particle_step(pdt)
+  use mpi
    use basic
    use diagnostics
    use auxiliary_fields
    implicit none
-   include 'mpif.h'
 
    real, intent(in) :: pdt
 
@@ -2093,11 +2093,11 @@ subroutine particle_step(pdt)
 end subroutine particle_step
 !---------------------------------------------------------------------------
 subroutine update_particle_pressure
+  use mpi
    use basic
    use arrays
    use diagnostics
    implicit none
-   include 'mpif.h'
 
    real    :: tstart, tend
    integer :: ierr
@@ -2127,9 +2127,9 @@ subroutine update_particle_pressure
 end subroutine update_particle_pressure
 !---------------------------------------------------------------------------
 subroutine finalize_particles
+  use mpi
    use arrays
    implicit none
-   include 'mpif.h'
 
    !integer :: nelms, ielm
    integer :: ierr
@@ -2631,11 +2631,11 @@ subroutine update_geom_terms_st(gh, fh, ic2)
 end subroutine update_geom_terms_st
 !---------------------------------------------------------------------------
 subroutine get_field_coefs(eq)
+  use mpi
    use arrays
    use basic
    use auxiliary_fields
    implicit none
-   include 'mpif.h'
 
    !type(elfield), intent(out) :: fh  !Field handle
    integer, intent(in) :: eq
@@ -3274,12 +3274,12 @@ subroutine evalf0(x, vpar, vperp, fh, gh, sps, f0, gradcoef, df0de, df0dxi)
 end subroutine evalf0
 !---------------------------------------------------------------------------
 subroutine particle_pressure_rhs
+  use mpi
    use basic
    use arrays
    use math
    use m3dc1_nint
    implicit none
-   include 'mpif.h'
    intrinsic matmul
 
    real, dimension(dofs_per_element, coeffs_per_element) :: cl
@@ -3612,11 +3612,11 @@ end subroutine solve_pi_tensor
 !---------------------------------------------------------------------------
 ! Dump particle data for current timeslice using parallel HDF5.
 subroutine hdf5_write_particles(ierr)
+  use mpi
    use basic
    use hdf5_output
    implicit none
 
-   include 'mpif.h'
 
    integer, intent(out) :: ierr
 
@@ -3825,11 +3825,11 @@ end subroutine hdf5_write_particles
 !---------------------------------------------------------------------------
 !Read stored HDF5 particle data in parallel.
 subroutine hdf5_read_particles(filename, ierr)
+  use mpi
    use basic
    use diagnostics
    use hdf5_output
    implicit none
-   include 'mpif.h'
 
    character(len=*), intent(in) :: filename
    integer, intent(out) :: ierr
diff --git a/unstructured/particle_com.f90 b/unstructured/particle_com.f90
index 126f274f0..0268ff99e 100644
--- a/unstructured/particle_com.f90
+++ b/unstructured/particle_com.f90
@@ -136,9 +136,9 @@ module particles
 ! Note: any changes to the "particle" user-defined datatype must be reflected
 !       in the definitions of pnvars, pblklen, ptyps, and pdspls below.
 subroutine define_mpi_particle(ierr)
+  use mpi
    implicit none
 
-   include 'mpif.h'
 
    integer, intent(out) :: ierr
    integer, parameter :: pnvars = 15
@@ -189,9 +189,9 @@ subroutine define_mpi_particle(ierr)
 end subroutine define_mpi_particle
 
 subroutine define_mpi_elfield(ierr)
+  use mpi
    implicit none
 
-   include 'mpif.h'
 
    integer, intent(out) :: ierr
    integer, parameter :: pnvars = 24
@@ -264,13 +264,13 @@ end subroutine define_mpi_elfield
 
 !---------------------------------------------------------------------------
 subroutine particle_test
+  use mpi
    use basic
    use diagnostics
    use auxiliary_fields
    use m3dc1_nint
    use arrays
    implicit none
-   include 'mpif.h'
 
    integer, parameter :: trunit = 120
    real :: tstart, tend
@@ -376,6 +376,7 @@ real function tri_area(x1, z1, x2, z2, x3, z3)
 end function tri_area
 !---------------------------------------------------------------------------
 subroutine init_particles(lrestart, ierr)
+  use mpi
    use basic
    use arrays
    use hdf5_output
@@ -386,7 +387,6 @@ subroutine init_particles(lrestart, ierr)
    use boundary_conditions
 
    implicit none
-   include 'mpif.h'
 
    logical, intent(in) :: lrestart
    integer, intent(out) :: ierr
@@ -1879,13 +1879,13 @@ subroutine fdot(x, v, w, dxdt, dvdt, dwdt, dEpdt, itri, kel, f00, ierr, sps, B00
 end subroutine fdot
 
 subroutine particle_scaleback(scalefac)
+  use mpi
    use basic
    use field
    use arrays
    use mesh_mod
    implicit none
 
-   include 'mpif.h'
 
    vectype, intent(in) :: scalefac
    integer :: ipart, ielm, ierr
@@ -1916,10 +1916,10 @@ subroutine particle_scaleback(scalefac)
 end subroutine particle_scaleback
 !---------------------------------------------------------------------------
 subroutine delete_particle(exchange)
+  use mpi
    use basic
    use diagnostics
    implicit none
-   include 'mpif.h'
 
    logical, intent(in) :: exchange
    !integer, intent(in) :: ipart
@@ -1994,11 +1994,11 @@ end subroutine delete_particle
 
 !---------------------------------------------------------------------------
 subroutine particle_step(pdt)
+  use mpi
    use basic
    use diagnostics
    use auxiliary_fields
    implicit none
-   include 'mpif.h'
 
    real, intent(in) :: pdt
 
@@ -2057,11 +2057,11 @@ subroutine particle_step(pdt)
 end subroutine particle_step
 !---------------------------------------------------------------------------
 subroutine update_particle_pressure
+  use mpi
    use basic
    use arrays
    use diagnostics
    implicit none
-   include 'mpif.h'
 
    real    :: tstart, tend
    integer :: ierr
@@ -2091,9 +2091,9 @@ subroutine update_particle_pressure
 end subroutine update_particle_pressure
 !---------------------------------------------------------------------------
 subroutine finalize_particles
+  use mpi
    use arrays
    implicit none
-   include 'mpif.h'
 
    !integer :: nelms, ielm
    integer :: ierr
@@ -2595,11 +2595,11 @@ subroutine update_geom_terms_st(gh, fh, ic2)
 end subroutine update_geom_terms_st
 !---------------------------------------------------------------------------
 subroutine get_field_coefs(eq)
+  use mpi
    use arrays
    use basic
    use auxiliary_fields
    implicit none
-   include 'mpif.h'
 
    !type(elfield), intent(out) :: fh  !Field handle
    integer, intent(in) :: eq
@@ -3301,12 +3301,12 @@ subroutine evalf0(x, vpar, vperp, fh, gh, sps, B0, f0, df0dpsi, df0de, df0dxi)
 end subroutine evalf0
 !---------------------------------------------------------------------------
 subroutine particle_pressure_rhs
+  use mpi
    use basic
    use arrays
    use math
    use m3dc1_nint
    implicit none
-   include 'mpif.h'
    intrinsic matmul
 
    real, dimension(dofs_per_element, coeffs_per_element) :: cl
@@ -3628,11 +3628,11 @@ end subroutine solve_pi_tensor
 !---------------------------------------------------------------------------
 ! Dump particle data for current timeslice using parallel HDF5.
 subroutine hdf5_write_particles(ierr)
+  use mpi
    use basic
    use hdf5_output
    implicit none
 
-   include 'mpif.h'
 
    integer, intent(out) :: ierr
 
@@ -3846,11 +3846,11 @@ end subroutine hdf5_write_particles
 !---------------------------------------------------------------------------
 !Read stored HDF5 particle data in parallel.
 subroutine hdf5_read_particles(filename, ierr)
+  use mpi
    use basic
    use diagnostics
    use hdf5_output
    implicit none
-   include 'mpif.h'
 
    character(len=*), intent(in) :: filename
    integer, intent(out) :: ierr
diff --git a/unstructured/pellet.f90 b/unstructured/pellet.f90
index bd3aa3b77..6bc1dd497 100644
--- a/unstructured/pellet.f90
+++ b/unstructured/pellet.f90
@@ -306,12 +306,12 @@ end function pellet_distribution
   ! Keep track if a pellet has left the plasma domain and should no longer
   ! deposit material
   subroutine pellet_domain
+    use mpi
 
     use basic
 
     implicit none
 
-    include 'mpif.h'
 
     integer :: j
     integer :: itri, iz, izone, ier
diff --git a/unstructured/petsc_matrix.f90 b/unstructured/petsc_matrix.f90
index b0bd32c41..bac3ed3de 100644
--- a/unstructured/petsc_matrix.f90
+++ b/unstructured/petsc_matrix.f90
@@ -4,8 +4,8 @@ module petsc_matrix_mod
 
   implicit none
 
-#include "finclude/petscmatdef.h"
-#include "finclude/petsckspdef.h"
+!#include <petsc/finclude/petscmatdef.h>
+!#include <petsc/finclude/petsckspdef.h>
 
   type :: petsc_matrix
      Mat :: data
@@ -107,9 +107,9 @@ subroutine petsc_matrix_create(mat, m, n, icomplex, lhs)
     use mesh_mod
     use vector_mod
     implicit none
-#include "finclude/petsc.h"
+!#include <petsc/finclude/petsc.h>
 !#ifndef PETSC_31
-!#include "finclude/petscmat.h"
+!#include <petsc/finclude/petscmat.h>
 !#endif
 
     type(petsc_matrix) :: mat
@@ -247,7 +247,7 @@ subroutine petsc_matrix_insert_real(mat,val,i,j,iop)
     use vector_mod
     use mesh_mod
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
 
     type(petsc_matrix), intent(in) :: mat
     real, intent(in) :: val
@@ -282,7 +282,7 @@ end subroutine petsc_matrix_insert_real
   subroutine petsc_matrix_insert_complex(mat,val,i,j,iop)
     use vector_mod
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
 
     type(petsc_matrix), intent(in) :: mat
     complex, intent(in) :: val
@@ -334,7 +334,7 @@ end subroutine petsc_matrix_solve
   !====================================================================
   subroutine petsc_matrix_finalize(mat)
     implicit none
-#include "finclude/petscmat.h"
+!#include <petsc/finclude/petscmat.h>
     type(petsc_matrix) :: mat
     integer :: ierr
 
@@ -349,7 +349,7 @@ end subroutine petsc_matrix_finalize
   !====================================================================
   subroutine petsc_matrix_flush(mat)
     implicit none
-#include "finclude/petscmat.h"
+!#include <petsc/finclude/petscmat.h>
     type(petsc_matrix) :: mat
     integer :: ierr
 
@@ -411,7 +411,7 @@ end subroutine petsc_matrix_get_node_indices
   subroutine petsc_matrix_insert_block(mat, itri, m, n, val, iop)
     use mesh_mod
     implicit none
-!#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
 
     type(matrix_type) :: mat
     integer, intent(in) :: itri, m, n, iop
@@ -447,7 +447,7 @@ subroutine identity_row(mat, irow)
     use mesh_mod
     use vector_mod
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
 
     type(petsc_matrix) :: mat
     integer, intent(in) :: irow
@@ -464,7 +464,7 @@ end subroutine identity_row
   subroutine set_row_vals(mat, irow, ncols, icols, vals)
     use mesh_mod
     use vector_mod
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
 
     type(petsc_matrix) :: mat
     integer, intent(in) :: irow, ncols
@@ -485,7 +485,7 @@ end subroutine set_row_vals
 
   subroutine petsc_matrix_write(m, file)
     implicit none
-#include "finclude/petsc.h"    
+!#include <petsc/finclude/petsc.h>    
     type(matrix_type), intent(in) :: m
     character(len=*) :: file
 
diff --git a/unstructured/petsc_vector.f90 b/unstructured/petsc_vector.f90
index 9df101234..532aad74b 100644
--- a/unstructured/petsc_vector.f90
+++ b/unstructured/petsc_vector.f90
@@ -2,7 +2,7 @@ module petsc_vector_mod
 
   implicit none
 
-#include "finclude/petscvecdef.h"
+!#include <petsc/finclude/petscvecdef.h>
 
   type :: petsc_vector
      Vec :: vec
@@ -132,7 +132,7 @@ end function dof_index
   subroutine petsc_vector_create(v,n)
     use mesh_mod
     implicit none
-#include "finclude/petsc.h"
+!#include <petsc/finclude/petsc.h>
 
     type(petsc_vector), intent(inout) :: v
     integer, intent(in) :: n
@@ -183,7 +183,7 @@ end subroutine petsc_vector_destroy
   subroutine petsc_vector_copy(vout,vin)
     use mesh_mod
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
 
     type(petsc_vector), intent(inout) :: vout    
     type(petsc_vector), intent(in) :: vin
@@ -223,7 +223,7 @@ end subroutine petsc_vector_copy
   !======================================================================
   subroutine petsc_vector_const_real(v,s)
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
 
     type(petsc_vector), intent(inout) :: v
     real, intent(in) :: s
@@ -244,7 +244,7 @@ end subroutine petsc_vector_const_real
   subroutine petsc_vector_add(vout,vin)
     use mesh_mod
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
     type(petsc_vector), intent(inout) :: vout    
     type(petsc_vector), intent(in) :: vin
     
@@ -284,7 +284,7 @@ end subroutine petsc_vector_add
   !======================================================================
   subroutine petsc_vector_multiply_real(v, s)
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
     type(petsc_vector), intent(inout) :: v
     real, intent(in) :: s
     integer :: ierr
@@ -319,7 +319,7 @@ end subroutine petsc_vector_multiply_complex
   !======================================================================
   subroutine petsc_vector_insert_real(v, i, s, iop)
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
     
     type(petsc_vector), intent(inout) :: v
     real, intent(in) :: s
@@ -346,7 +346,7 @@ end subroutine petsc_vector_insert_real
 #ifdef USECOMPLEX
   subroutine petsc_vector_insert_complex(v, i, s, iop)
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
     
     type(petsc_vector), intent(inout) :: v
     complex, intent(in) :: s
@@ -374,7 +374,7 @@ end subroutine petsc_vector_insert_complex
 
   subroutine petsc_vector_sum_shared(v)
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
 
     type(petsc_vector), intent(inout) :: v
     integer :: ierr
@@ -391,7 +391,7 @@ end subroutine petsc_vector_sum_shared
 
   subroutine petsc_vector_finalize(v)
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
 
     type(petsc_vector), intent(inout) :: v
     integer :: ierr
@@ -471,7 +471,7 @@ end subroutine petsc_vector_insert_block
 
   logical function petsc_vector_is_nan(v)
     implicit none
-#include "finclude/petsc.h"
+!#include <petsc/finclude/petsc.h>
     type(petsc_vector), intent(in) :: v
     integer :: ierr
     PetscScalar :: y(1)
@@ -489,7 +489,7 @@ subroutine petsc_vector_get_node_data_real(v, iplace, inode, data, rotate)
     use element
     use mesh_mod
     implicit none
-#include "finclude/petsc.h"
+!#include <petsc/finclude/petsc.h>
 
     type(petsc_vector), intent(in) :: v
     integer, intent(in) :: inode, iplace
@@ -546,7 +546,7 @@ subroutine petsc_vector_get_node_data_complex(v, iplace, inode, data, rotate)
     use element
     use mesh_mod
     implicit none
-#include "finclude/petsc.h"
+!#include <petsc/finclude/petsc.h>
     type(petsc_vector), intent(in) :: v
     integer, intent(in) :: inode, iplace
     complex, intent(out), dimension(dofs_per_node) :: data
@@ -600,7 +600,7 @@ subroutine petsc_vector_set_node_data_real(v, iplace, inode, data, rotate)
     use element
     use mesh_mod
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
     type(petsc_vector), intent(inout) :: v
     integer, intent(in) :: inode, iplace
     real, intent(in), dimension(dofs_per_node) :: data
@@ -649,7 +649,7 @@ subroutine petsc_vector_set_node_data_complex(v, iplace, inode, data, rotate)
     use element
     use mesh_mod
     implicit none
-#include "finclude/petscvec.h"
+!#include <petsc/finclude/petscvec.h>
     type(petsc_vector), intent(inout) :: v
     integer, intent(in) :: inode, iplace
     complex, intent(in), dimension(dofs_per_node) :: data
@@ -738,7 +738,7 @@ end function global_dof_id
 
   subroutine petsc_vector_write(v, file)
     implicit none
-#include "finclude/petsc.h"    
+!#include <petsc/finclude/petsc.h>    
     type(vector_type), intent(in) :: v
     character(len=*) :: file
 
diff --git a/unstructured/read_ascii.f90 b/unstructured/read_ascii.f90
index 7b817343f..97929ac48 100644
--- a/unstructured/read_ascii.f90
+++ b/unstructured/read_ascii.f90
@@ -53,9 +53,9 @@ end subroutine get_token
   ! xrow = column to read (default = 1)
   !======================================================================
   subroutine read_ascii_column_real(filename, x, n, skip, icol, read_until)
+    use mpi
     implicit none
 
-    include 'mpif.h'
 
     character(len=*), intent(in) :: filename
     real, allocatable :: x(:)
diff --git a/unstructured/read_gyro.f90 b/unstructured/read_gyro.f90
index e16dee56b..cfe15066b 100644
--- a/unstructured/read_gyro.f90
+++ b/unstructured/read_gyro.f90
@@ -37,10 +37,10 @@ module gyro
 contains
 
   subroutine load_gyro(ierr)
+    use mpi
 
     implicit none
 
-    include 'mpif.h'
 
     integer, intent(out) :: ierr
 
diff --git a/unstructured/read_neo.f90 b/unstructured/read_neo.f90
index 2771b6a9f..f2a94e40c 100644
--- a/unstructured/read_neo.f90
+++ b/unstructured/read_neo.f90
@@ -69,9 +69,9 @@ subroutine unload_neo()
   end subroutine unload_neo
 
   subroutine read_neo_expnorm(ierr)
+    use mpi
     implicit none
 
-    include 'mpif.h'
 
     integer, intent(out) :: ierr
     integer, parameter :: ifile=22
@@ -113,9 +113,9 @@ subroutine read_neo_expnorm(ierr)
   end subroutine read_neo_expnorm
 
   subroutine read_neo_grid(ierr)
+    use mpi
     implicit none
 
-    include 'mpif.h'
 
     integer, intent(out) :: ierr
     integer, parameter :: ifile=22
@@ -185,9 +185,9 @@ end subroutine read_neo_grid
 
 
   subroutine read_neo_vel_fourier(ierr)
+    use mpi
     implicit none
 
-    include 'mpif.h'
 
     integer, intent(out) :: ierr
     integer, parameter :: ifile=23
diff --git a/unstructured/read_schaffer_field.f90 b/unstructured/read_schaffer_field.f90
index ba2c49f15..07ec53c5d 100644
--- a/unstructured/read_schaffer_field.f90
+++ b/unstructured/read_schaffer_field.f90
@@ -15,10 +15,10 @@ module read_schaffer_field
 contains
 
   subroutine load_schaffer_field(sf, filename, isamp, isamp_pol, ierr)
+    use mpi
     use math
     implicit none
 
-    include 'mpif.h'
 
     type(schaffer_field), intent(inout) :: sf
     character(len=*), intent(in) :: filename
@@ -369,9 +369,9 @@ end subroutine get_external_field_ft
 
   subroutine get_external_field(sf,r1,phi1,z1,br_out,bphi_out,bz_out,p_out, &
        npts)
+         use mpi
     implicit none
 
-    include 'mpif.h'
 
     type(schaffer_field), intent(in) :: sf
     integer, intent(in) :: npts
diff --git a/unstructured/readaeqdsk.f90 b/unstructured/readaeqdsk.f90
index 01ec6405e..a2c52179d 100644
--- a/unstructured/readaeqdsk.f90
+++ b/unstructured/readaeqdsk.f90
@@ -71,11 +71,11 @@ module eqdsk_a
 contains
 
 subroutine load_eqdsk_a(filename)
+  use mpi
 
   implicit none
 
 #ifdef USEMPI
-  include 'mpif.h'
 #endif
 
   character(len=*), intent(in) :: filename
diff --git a/unstructured/readgeqdsk.f90 b/unstructured/readgeqdsk.f90
index 9f69dcb8e..4b586ecf3 100644
--- a/unstructured/readgeqdsk.f90
+++ b/unstructured/readgeqdsk.f90
@@ -16,9 +16,9 @@ module eqdsk
 contains
 
 subroutine load_eqdsk(ierr)
+  use mpi
   implicit none
 
-  include 'mpif.h'
 
   integer, intent(out) :: ierr
 
diff --git a/unstructured/region.f90 b/unstructured/region.f90
index af066ad13..b59eccbf5 100644
--- a/unstructured/region.f90
+++ b/unstructured/region.f90
@@ -66,9 +66,9 @@ subroutine destroy_plane(p)
     end subroutine destroy_plane
 
     subroutine create_region_from_file(r, filename, ierr)
+      use mpi
       implicit none
 
-      include 'mpif.h'
 
       type(region_type), intent(inout) :: r
       character(len=*), intent(in) :: filename
diff --git a/unstructured/restart_hdf5.f90 b/unstructured/restart_hdf5.f90
index ebc3b9579..176d612d5 100644
--- a/unstructured/restart_hdf5.f90
+++ b/unstructured/restart_hdf5.f90
@@ -9,6 +9,7 @@ module restart_hdf5
 contains
   
   subroutine rdrestart_hdf5()
+    use mpi
     use basic
     use hdf5_output
     use hdf5
@@ -20,7 +21,6 @@ subroutine rdrestart_hdf5()
 
     implicit none
 
-    include 'mpif.h'
 
     integer :: error
     integer(HID_T) :: root_id, scalar_group_id, time_id, eq_time_id, pel_group_id, mesh_id
diff --git a/unstructured/rmp.f90 b/unstructured/rmp.f90
index d42d9887a..10ca2071a 100644
--- a/unstructured/rmp.f90
+++ b/unstructured/rmp.f90
@@ -409,6 +409,7 @@ end subroutine rmp_field
 
 
 subroutine calculate_external_fields(ilin)
+  use mpi
   use basic
   use math
   use mesh_mod
@@ -422,7 +423,6 @@ subroutine calculate_external_fields(ilin)
 
   implicit none
 
-  include 'mpif.h'
 
   type(matrix_type) :: br_mat, bf_mat
   type(vector_type) :: psi_vec, bz_vec, p_vec, bf_vec
@@ -849,6 +849,7 @@ subroutine boundary_rmp(rhs, mat)
 end subroutine boundary_rmp
 
   subroutine calculate_field_from_j(ilin)
+    use mpi
     use basic
     use math
     use mesh_mod
@@ -861,7 +862,6 @@ subroutine calculate_field_from_j(ilin)
 
     implicit none
 
-    include 'mpif.h'
 
     type(matrix_type) :: jx_mat
     type(vector_type) :: psi_vec
diff --git a/unstructured/runaway_advection.f90 b/unstructured/runaway_advection.f90
index fcb3b9b5d..e7cff9850 100644
--- a/unstructured/runaway_advection.f90
+++ b/unstructured/runaway_advection.f90
@@ -98,9 +98,9 @@ module runaway_advection
 ! Note: any changes to the "particle" user-defined datatype must be reflected
 !       in the definitions of pnvars, pblklen, ptyps, and pdspls below.
 subroutine define_mpi_particle(ierr)
+  use mpi
    implicit none
 
-   include 'mpif.h'
 
    integer, intent(out) :: ierr
    integer, parameter :: pnvars = 10
@@ -144,9 +144,9 @@ subroutine define_mpi_particle(ierr)
 end subroutine define_mpi_particle
 
 subroutine define_mpi_elfield(ierr)
+  use mpi
    implicit none
 
-   include 'mpif.h'
 
    integer, intent(out) :: ierr
    integer, parameter :: pnvars = 12
@@ -194,11 +194,11 @@ subroutine define_mpi_elfield(ierr)
 end subroutine define_mpi_elfield
 !---------------------------------------------------------------------------
 subroutine runaway_advection_initialize
+  use mpi
    use basic
    use diagnostics
    use auxiliary_fields
    implicit none
-   include 'mpif.h'
 
    integer, parameter :: trunit = 120
    real :: tstart, tend
@@ -248,13 +248,13 @@ real function tri_area(x1, z1, x2, z2, x3, z3)
 end function tri_area
 !---------------------------------------------------------------------------
 subroutine init_particles(lrestart, ierr)
+  use mpi
    use basic
    use arrays
    use m3dc1_nint
    use gradshafranov
    use read_ascii
    implicit none
-   include 'mpif.h'
 
    logical, intent(in) :: lrestart
    integer, intent(out) :: ierr
@@ -678,10 +678,10 @@ subroutine init_particles(lrestart, ierr)
 end subroutine init_particles
 
 subroutine advance_particles(tinc)
+  use mpi
    use basic  !For MPI variables
    use omp_lib
    implicit none
-   include 'mpif.h'
 
    real, intent(in) :: tinc  !Time increment for particle advance
    type(elfield), dimension(nneighbors + 1) :: elcoefs
@@ -1091,11 +1091,11 @@ subroutine fdot(x, v, w, dxdt, dvdt, dwdt, dEpdt, itri, ierr)
 end subroutine fdot
 !---------------------------------------------------------------------------
 subroutine runaway_advection_step(pdt)
+  use mpi
    use basic
    use diagnostics
    use auxiliary_fields
    implicit none
-   include 'mpif.h'
 
    real, intent(in) :: pdt
 
@@ -1127,11 +1127,11 @@ subroutine runaway_advection_step(pdt)
 end subroutine runaway_advection_step
 !---------------------------------------------------------------------------
 subroutine update_particle_pressure
+  use mpi
    use basic
    use arrays
    use diagnostics
    implicit none
-   include 'mpif.h'
 
    real    :: tstart, tend
    integer :: ierr
@@ -1540,11 +1540,11 @@ subroutine get_geom_terms(x, ielm, gh, ic2, ierr)
 end subroutine get_geom_terms
 !---------------------------------------------------------------------------
 subroutine get_field_coefs(eq)
+  use mpi
    use arrays
    use basic
    use auxiliary_fields
    implicit none
-   include 'mpif.h'
 
    !type(elfield), intent(out) :: fh  !Field handle
    integer, intent(in) :: eq
diff --git a/unstructured/scorec_matrix.f90 b/unstructured/scorec_matrix.f90
index 53be5efdb..f6eb8f1e4 100644
--- a/unstructured/scorec_matrix.f90
+++ b/unstructured/scorec_matrix.f90
@@ -1,6 +1,7 @@
 module scorec_matrix_mod
 
   use element
+  use mpi
 
   implicit none
 
@@ -187,16 +188,8 @@ subroutine scorec_matrix_create(mat, m, n, icomplex, lhs)
 
 #ifndef M3DC1_TRILINOS
 
-#if PETSC_VERSION >= 38
     use petsc
     implicit none
-#elif PETSC_VERSION >= 36
-    implicit none
-#include "petsc/finclude/petsc.h"
-#else
-    implicit none
-#include "finclude/petsc.h"
-#endif
 
 #endif
 
@@ -444,16 +437,8 @@ end subroutine scorec_matrix_insert_global_complex
   subroutine scorec_matrix_solve(mat, v, ierr)
     use vector_mod
 
-#if PETSC_VERSION >= 38
     use petsc
     implicit none
-#elif PETSC_VERSION >= 36
-    implicit none
-#include "petsc/finclude/petsc.h"
-#else
-    implicit none
-#include "finclude/petsc.h"
-#endif
     
     type(scorec_matrix), intent(in) :: mat
     type(vector_type), intent(inout) :: v
@@ -498,16 +483,8 @@ end subroutine scorec_matrix_solve
 
   subroutine scorec_matrix_solve_with_guess(mat, v, x_guess, ierr)
    use vector_mod
-#if PETSC_VERSION >= 38
    use petsc
    implicit none
-#elif PETSC_VERSION >= 36
-   implicit none
-#include "petsc/finclude/petsc.h"
-#else
-   implicit none
-#include "finclude/petsc.h"
-#endif
    
    type(scorec_matrix), intent(in) :: mat
    type(vector_type), intent(inout) :: v, x_guess
diff --git a/unstructured/scorec_mesh.f90 b/unstructured/scorec_mesh.f90
index 6bf808a4d..2d18d51e5 100644
--- a/unstructured/scorec_mesh.f90
+++ b/unstructured/scorec_mesh.f90
@@ -67,12 +67,12 @@ module scorec_mesh_mod
 contains
 
   subroutine load_mesh
+    use mpi
     use math
     use read_ascii
     implicit none
 
     integer :: myrank, maxrank, ier
-    include 'mpif.h'
 #ifdef USE3D
     real :: angle, beta
     integer :: i,procs_per_plane, full_group, plane_group
@@ -220,8 +220,8 @@ subroutine load_mesh
   end subroutine load_mesh
 
   subroutine update_nodes_owned
+    use mpi
     integer :: numnodes, inode2, inode, ier, myrank, own_process
-    include 'mpif.h'
     call MPI_Comm_rank(MPI_COMM_WORLD,myrank,ier)
     if(allocated(nodes_owned)) deallocate(nodes_owned)
     allocate(nodes_owned(owned_nodes()))
@@ -262,10 +262,10 @@ subroutine print_node_data()
   end subroutine print_node_data
 
   subroutine get_global_dims
+    use mpi
 
     implicit none
 
-    include 'mpif.h'
 
     integer :: nelms, error
 
diff --git a/unstructured/struct2vac.f90 b/unstructured/struct2vac.f90
index a0002542e..79ec8d26b 100644
--- a/unstructured/struct2vac.f90
+++ b/unstructured/struct2vac.f90
@@ -1,10 +1,10 @@
 program struct2vac
+  use mpi
 
   use mesh_mod
 
   implicit none
 
-  include 'mpif.h'
 
   integer :: inode, ier, num_global_nodes, num_local_nodes
   integer :: i, j, ilast, ifirst, numnodes
diff --git a/unstructured/test_region.f90 b/unstructured/test_region.f90
index 20f19d7a1..ba2ee84b3 100644
--- a/unstructured/test_region.f90
+++ b/unstructured/test_region.f90
@@ -1,10 +1,12 @@
 Program test_region
+  use mpi
   use region
 
   implicit none
 
   type(region_type) :: reg
   integer :: ierr
+  character(len=256) :: boundary_file
   
   integer, parameter :: n = 5
   real, dimension(n) :: r, phi, z
@@ -18,13 +20,20 @@ Program test_region
   r(5) = -2.;  z(5) = -12.
   phi = 0.
 
-  call create_region_from_file(reg, "outer_boundary.pts", ierr)
+  call MPI_Init(ierr)
+
+  boundary_file = "unstructured/templates/ITER/rw2_adapt/outer_boundary.pts"
+  if(command_argument_count().ge.1) call get_command_argument(1, boundary_file)
+
+  call create_region_from_file(reg, trim(boundary_file), ierr)
   
   do i=1, n
      in = point_in_region(reg, r(i), phi(i), z(i))
-     write(*,'("Is (", G, ", ", G," in region? ", L)') r(i), z(i), in
+     write(*,'("Is (", G0, ", ", G0, " in region? ", L1)') r(i), z(i), in
   end do
 
   call destroy_region(reg)
 
+  call MPI_Finalize(ierr)
+
 end Program test_region
diff --git a/unstructured/time_step.f90 b/unstructured/time_step.f90
index 470eba0b1..3be287f16 100644
--- a/unstructured/time_step.f90
+++ b/unstructured/time_step.f90
@@ -326,13 +326,13 @@ end subroutine scaleback
 
 
 subroutine variable_timestep
+  use mpi
 
   use basic
   use arrays
   use diagnostics
 
   implicit none
-  include 'mpif.h'
   integer :: ierr
 
 #ifdef USE3D
diff --git a/unstructured/transport.f90 b/unstructured/transport.f90
index 9d4da8503..4b0ed5465 100644
--- a/unstructured/transport.f90
+++ b/unstructured/transport.f90
@@ -1577,6 +1577,7 @@ end function bs_func
 ! define_transport_coefficients
 ! =============================
 subroutine define_transport_coefficients()
+  use mpi
 
   use basic
   use arrays
@@ -1590,7 +1591,6 @@ subroutine define_transport_coefficients()
 
   implicit none
 
-  include 'mpif.h'
 
   integer :: itri, izone, izone_index
   integer :: numelms, def_fields,ier
diff --git a/unstructured/vector_mod.F90 b/unstructured/vector_mod.F90
new file mode 100644
index 000000000..8aa9ec07b
--- /dev/null
+++ b/unstructured/vector_mod.F90
@@ -0,0 +1,14 @@
+! vector_mod.F90 — thin wrapper that selects the vector backend at compile time.
+! CMake defines exactly one of the following:
+!   M3DC1_VECTOR_IS_SCOREC  → use scorec_vector_mod
+!   M3DC1_VECTOR_IS_PETSC   → use petsc_vector_mod
+
+module vector_mod
+#ifdef M3DC1_VECTOR_IS_SCOREC
+  use scorec_vector_mod
+#elif defined(M3DC1_VECTOR_IS_PETSC)
+  use petsc_vector_mod
+#else
+  use scorec_vector_mod
+#endif
+end module vector_mod