module m_MEF90_EXO
#include "petsc/finclude/petsc.h"
   use m_MEF90_Parameters
   use m_MEF90_Utils
   use m_MEF90_Elements
   use m_MEF90_Ctx
   use m_MEF90_DMPlex

   implicit none(type)
#include "../mef90version.h"

   private
   PetscInt, public                                 :: exo_ver

   public :: MEF90CtxOpenEXO
   public :: MEF90CtxCloseEXO
   public :: MEF90EXOFormat
   public :: MEF90EXODMView
   public :: MEF90EXOVecView
   public :: MEF90EXOVecLoad

contains
#undef __FUNCT__
#define __FUNCT__ "MEF90CtxOpenEXO"
!!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca)
!!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca)
!!!
!!!  MEF90CtxOpenEXO:
!!!
   subroutine MEF90CtxOpenEXO(MEF90Ctx, Viewer, mode, ierr)
      type(MEF90Ctx_Type), target, intent(IN)                   :: MEF90Ctx
      type(tPetscViewer), intent(INOUT)               :: Viewer
      type(ePetscFileMode), intent(IN)                :: mode
      PetscErrorCode, intent(INOUT)                   :: ierr

      integer                                         :: opts

#ifdef PETSC_USE_DEBUG
      opts = EXVRBS + EXDEBG
#else
      opts = 0
#endif
      call exopts(opts, ierr)
      PetscCall(PetscViewerExodusIIOpen(MEF90Ctx%Comm, MEF90Ctx%resultFile, mode, Viewer, ierr))
   end subroutine MEF90CtxOpenEXO

#undef __FUNCT__
#define __FUNCT__ "MEF90CtxCloseEXO"
!!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca)
!!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca)
!!!
!!!  MEF90CtxCloseEXO:
!!!
   subroutine MEF90CtxCloseEXO(Viewer, ierr)
      type(tPetscViewer), intent(INOUT)                :: Viewer
      PetscErrorCode, intent(INOUT)                    :: ierr

      PetscCall(PetscViewerDestroy(Viewer, ierr))
   end subroutine MEF90CtxCloseEXO

#undef __FUNCT__
#define __FUNCT__ "MEF90EXOFormat"
!!! author: Blaise Bourdin (2012-2022, bourdin@lsu.edu)
!!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca)
!!! author: Blaise Bourdin (2024, bourdin@mcmaster.ca)
!!!
!!!  MEF90EXOFormat:
!!!
   subroutine MEF90EXOFormat(Viewer, nameG, nameC, nameV, time, ierr)
      type(tPetscViewer), intent(IN)                        :: Viewer
      character(len=*), dimension(:), intent(IN)            :: nameG, nameC, nameV
      PetscReal, dimension(:), pointer                      :: time
      PetscErrorCode, intent(INOUT)                         :: ierr

      PetscInt                                              :: numCS, numC
      integer                                               :: i, exoid
      PetscInt                                              :: step
      character(len=MXSTLN)                                 :: sJunk
      PetscReal                                             :: rJunk
      logical, dimension(:, :), pointer                     :: truthtable

      if (size(nameV) > 0) then
         PetscCall(PetscViewerExodusIISetNodalVariable(Viewer, size(nameV), ierr))
         do i = 1, size(nameV)
            PetscCall(PetscViewerExodusIISetNodalVariableName(Viewer, i - 1, nameV(i), ierr))
         end do
      end if
      if (size(nameC) > 0) then
         PetscCall(PetscViewerExodusIISetZonalVariable(Viewer, size(nameC), ierr))
         do i = 1, size(nameC)
            PetscCall(PetscViewerExodusIISetZonalVariableName(Viewer, i - 1, nameC(i), ierr))
         end do
      end if

      if (.not. associated(time)) then
         SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_POINTER, "Time value must be allocated prior to calling MEF90EXOFormat")
         stop
      end if

      PetscCall(PetscViewerExodusIIGetId(Viewer, exoid, ierr))
      call exinq(exoid, EX_INQ_ELEM_BLK, numCS, rJunk, sjunk, ierr)

   !! Write truth tables
      numC = size(nameC)
      if (numC > 0) then
         allocate (truthtable(numCS, numC))
         truthtable = .true.
         call expvtt(exoid, numCS, numC, truthtable, ierr)
         deallocate (truthtable)
      end if

      do step = 1, size(time)
         call exptim(exoid, step, time(step), ierr)
      end do
      PetscCall(PetscViewerFlush(Viewer, ierr))
   end subroutine MEF90EXOFormat

#undef __FUNCT__
#define __FUNCT__ "MEF90EXODMView"
!!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca)
!!!
!!!  MEF90EXODMView:
!!!
   subroutine MEF90EXODMView(dm, Viewer, order, ierr)
      type(tPetscViewer), intent(IN)                      :: Viewer
      type(tDM), intent(IN)                               :: dm
      PetscInt, intent(IN)                                :: order
      PetscErrorCode, intent(INOUT)                       :: ierr

      character(len=PETSC_MAX_PATH_LEN)                   :: IOBuffer

      if ((order > 2) .or. (order < 1)) then
         write (IOBuffer, '("Unsupported polynomial order ", I2, " not in [1,2]")') order
         SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, IOBuffer)
      end if
      PetscCall(PetscViewerExodusIISetOrder(Viewer, order, ierr))
      PetscCall(DMView(dm, Viewer, ierr))
   end subroutine MEF90EXODMView

#undef __FUNCT__
#define __FUNCT__ "MEF90EXOVecView"
!!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca)
!!! author: Blaise Bourdin (2023, bourdin@mcmaster.ca)
!!!
!!!  MEF90EXOVecView:
!!!
   subroutine MEF90EXOVecView(v, sf, invSF, Viewer, step, bs, ierr)
      type(tVec), intent(IN)                             :: v
      type(tPetscSF), intent(IN)                         :: sf, invSF
      type(tPetscViewer), intent(IN)                     :: Viewer
      PetscExodusIIInt, intent(IN)                       :: step
      PetscInt, intent(IN)                               :: bs
      PetscErrorCode, intent(INOUT)                      :: ierr

      integer                                            :: exoid
      integer                                            :: offsetN, offsetZ
      type(tVec)                                         :: iov
      character(len=PETSC_MAX_PATH_LEN)                  :: vecname, IOBuffer

      offsetN = -1
      offsetZ = -1
      PetscCall(PetscObjectGetName(v, vecname, ierr))

      PetscCall(MEF90VecCreateIO(iov, bs, sf, ierr))
      PetscCall(PetscObjectSetName(iov, vecname, ierr))
      PetscCall(MEF90VecCopySF(v, iov, sf, ierr))

      PetscCall(PetscViewerExodusIIGetNodalVariableIndex(Viewer, vecname, offsetN, ierr))
      PetscCall(PetscViewerExodusIIGetZonalVariableIndex(Viewer, vecname, offsetZ, ierr))
      PetscCall(PetscViewerExodusIIGetId(Viewer, exoid, ierr))
      if (offsetN >= 0) then
         PetscCall(MEF90EXOVecViewNodal_Private(iov, exoid, step, offsetN + 1, ierr))
      else if (offsetZ >= 0) then
         PetscCall(MEF90EXOVecViewZonal_Private(iov, exoid, step, offsetZ + 1, ierr))
      else
         write (IOBuffer, '("Could not find nodal or zonal variable ", A, " in exodus file. ")') trim(vecname)
         SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_FILE_UNEXPECTED, IOBuffer)
      end if
      PetscCall(VecDestroy(iov, ierr))
   end subroutine MEF90EXOVecView

#undef __FUNCT__
#define __FUNCT__ "MEF90EXOVecLoad"
!!! author: Alexis Marboeuf (2022, marboeua@mcmaster.ca)
!!!
!!!  MEF90EXOVecLoad:
!!!
   subroutine MEF90EXOVecLoad(v, sf, invSF, Viewer, step, bs, ierr)
      type(tVec), intent(INOUT)                          :: v
      type(tPetscSF), intent(IN)                         :: sf, invSF
      type(tPetscViewer), intent(IN)                     :: Viewer
      PetscExodusIIInt, intent(IN)                       :: step
      PetscInt, intent(IN)                               :: bs
      PetscErrorCode, intent(INOUT)                      :: ierr

      integer                                            :: exoid
      integer                                            :: offsetN, offsetZ
      type(tVec)                                         :: iov
      type(tDM)                                          :: locDM, oDM
      type(tVec)                                         :: vGlob
      character(len=PETSC_MAX_PATH_LEN)                  :: vecname, IOBuffer

      offsetN = -1
      offsetZ = -1
      PetscCall(PetscObjectGetName(v, vecname, ierr))

      PetscCall(MEF90VecCreateIO(iov, bs, sf, ierr))
      PetscCall(PetscObjectSetName(iov, vecname, ierr))

      PetscCall(PetscViewerExodusIIGetNodalVariableIndex(Viewer, vecname, offsetN, ierr))
      PetscCall(PetscViewerExodusIIGetZonalVariableIndex(Viewer, vecname, offsetZ, ierr))
      PetscCall(PetscViewerExodusIIGetId(Viewer, exoid, ierr))
      if (offsetN >= 0) then
         PetscCall(MEF90EXOVecLoadNodal_Private(iov, exoid, step, offsetN + 1, ierr))
      else if (offsetZ >= 0) then
         PetscCall(MEF90EXOVecLoadZonal_Private(iov, exoid, step, offsetZ + 1, ierr))
      else
         write (IOBuffer, '("Could not find nodal or zonal variable ", A, " in exodus file. ")') trim(vecname)
         SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_FILE_UNEXPECTED, IOBuffer)
      end if
      PetscCall(MEF90VecCopySF(iov, v, invSF, ierr))
      PetscCall(VecDestroy(iov, ierr))

      !! Make sure that the halo values in the local vector are updated by doing a L2G followed by a G2L
      !! This is probably only needed for nodal vectors, unless we have distributed the mesh with an overlap.
      PetscCall(VecGetDM(V, locDM, ierr))
      PetscCall(DMGetOutputDM(locDM, oDM, ierr))
      PetscCall(DMGetGlobalVector(oDM, vGlob, ierr))
      PetscCall(DMLocalToGlobal(oDM, v, INSERT_ALL_VALUES, vGlob, ierr))
      PetscCall(DMGlobalToLocal(oDM, vGlob, INSERT_ALL_VALUES, v, ierr))
      PetscCall(DMRestoreGlobalVector(oDM, vGlob, ierr))

   end subroutine MEF90EXOVecLoad

#undef __FUNCT__
#define __FUNCT__ "MEF90EXOVecViewNodal_Private"
   subroutine MEF90EXOVecViewNodal_Private(v, exoid, step, offset, ierr)
      integer, intent(IN)                :: exoid
      PetscExodusIIInt, intent(IN)       :: step, offset
      type(tVec), intent(IN)             :: v
      PetscErrorCode, intent(INOUT)      :: ierr

      PetscExodusIIInt                   :: c
      PetscInt                           :: xs, xe, bs
      PetscScalar, dimension(:), pointer :: varray
      type(tVec)                         :: vComp
      type(tIS)                          :: compIS

      PetscCall(VecGetOwnershipRange(v, xs, xe, ierr))
      PetscCall(VecGetBlockSize(v, bs, ierr))

      if (bs == 1) then
         PetscCall(VecGetArrayRead(v, varray, ierr))
         call expnvs(exoid, step, offset, xs + 1, xe - xs, varray, ierr)
         PetscCall(VecRestoreArrayRead(v, varray, ierr))
      else
         PetscCall(ISCreateStride(PETSC_COMM_WORLD, (xe - xs) / bs, xs, bs, compIS, ierr))
         do c = 0, bs - 1
            PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs, ierr))
            PetscCall(VecGetSubVector(v, compIS, vComp, ierr))
            PetscCall(VecGetArrayRead(vComp, varray, ierr))
            call expnvs(exoid, step, offset + c, xs / bs + 1, (xe - xs) / bs, varray, ierr)
            PetscCall(VecRestoreArrayRead(vComp, varray, ierr))
            PetscCall(VecRestoreSubVector(v, compIS, vComp, ierr))
         end do ! c
         PetscCall(ISDestroy(compIS, ierr))
      end if ! bs
   end subroutine MEF90EXOVecViewNodal_Private

#undef __FUNCT__
#define __FUNCT__ "MEF90EXOVecLoadNodal_Private"
   subroutine MEF90EXOVecLoadNodal_Private(v, exoid, step, offset, ierr)
      integer, intent(IN)                :: exoid
      PetscExodusIIInt, intent(IN)       :: step, offset
      type(tVec), intent(INOUT)          :: v
      PetscErrorCode, intent(INOUT)      :: ierr

      PetscExodusIIInt                   :: c
      PetscInt                           :: xs, xe, bs
      PetscScalar, dimension(:), pointer :: varray
      type(tVec)                         :: vComp
      type(tIS)                          :: compIS

      PetscCall(VecGetOwnershipRange(v, xs, xe, ierr))
      PetscCall(VecGetBlockSize(v, bs, ierr))
      if (bs == 1) then
         PetscCall(VecGetArrayRead(v, varray, ierr))
         call exgnnv(exoid, step, offset, xs + 1, xe - xs, varray, ierr)
         PetscCall(VecRestoreArrayRead(v, varray, ierr))
      else
         PetscCall(ISCreateStride(PETSC_COMM_WORLD, (xe - xs) / bs, xs, bs, compIS, ierr))
         do c = 0, bs - 1
            PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs, ierr))
            PetscCall(VecGetSubVector(v, compIS, vComp, ierr))
            PetscCall(VecGetArrayRead(vComp, varray, ierr))
            call exgnnv(exoid, step, offset + c, xs / bs + 1, (xe - xs) / bs, varray, ierr)
            PetscCall(VecRestoreArrayRead(vComp, varray, ierr))
            PetscCall(VecISCopy(v, compIS, SCATTER_FORWARD, vComp, ierr))
            PetscCall(VecRestoreSubVector(v, compIS, vComp, ierr))
         end do
         PetscCall(ISDestroy(compIS, ierr))
      end if
   end subroutine MEF90EXOVecLoadNodal_Private

#undef __FUNCT__
#define __FUNCT__ "MEF90EXOVecViewZonal_Private"
   subroutine MEF90EXOVecViewZonal_Private(v, exoid, step, offset, ierr)
      integer, intent(IN)                :: exoid
      PetscExodusIIInt, intent(IN)       :: step, offset
      type(tVec), intent(IN)             :: v
      PetscErrorCode, intent(INOUT)      :: ierr

      PetscExodusIIInt                   :: c
      PetscInt                           :: xs, xe, bs, numCS, set, csLocalSize, csxs
      PetscScalar, dimension(:), pointer :: varray
      PetscInt, dimension(:), pointer    :: csID, csSize
      type(tVec)                         :: vComp
      type(tIS)                          :: compIS
      character(len=MXSTLN)              :: elemType
      PetscMPIInt                        :: rank

      csxs = 0_ki
      PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
      numCS = exinqi(exoid, EX_INQ_ELEM_BLK)
      allocate (csID(numCS))
      allocate (csSize(numCS))
      call exgebi(exoid, csID, ierr)
      do set = 1, numCS
         call exgelb(exoid, csID(set), elemType, csSize(set), PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr)
      end do
      PetscCall(VecGetOwnershipRange(v, xs, xe, ierr))
      PetscCall(VecGetBlockSize(v, bs, ierr))
      if (bs > 1) then
         PetscCall(ISCreateStride(PETSC_COMM_WORLD, (xe - xs) / bs, xs, bs, compIS, ierr))
      end if
      do set = 1, numCS
         !  range of indices for set setID[set]: csxs:csxs + csSize[set]-1
         !  local slice of zonal values:         xs/bs,xm/bs-1
         !  intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1)
         csLocalSize = max(0, min(xe / bs, csxs + csSize(set)) - max(xs / bs, csxs))
         if (bs == 1) then
            PetscCall(VecGetArrayRead(v, varray, ierr))
            call expevs(exoid, step, offset, csID(set), max(xs - csxs, 0) + 1, csLocalSize, varray, ierr)
            PetscCall(VecRestoreArrayRead(v, varray, ierr))
         else
            do c = 0, bs - 1
               PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs, ierr))
               PetscCall(VecGetSubVector(v, compIS, vComp, ierr))
               PetscCall(VecGetArrayRead(vComp, varray, ierr))
               call expevs(exoid, step, offset + c, csID(set), max(xs / bs - csxs, 0) + 1, csLocalSize, varray(max(0, csxs - xs / bs) + 1:max(0, csxs - xs / bs) + csLocalSize), ierr)
               PetscCall(VecRestoreArrayRead(vComp, varray, ierr))
               PetscCall(VecRestoreSubVector(v, compIS, vComp, ierr))
            end do
         end if
         csxs = csxs + csSize(set)
      end do
      if (bs > 1) then
         PetscCall(ISDestroy(compIS, ierr))
      end if
      deallocate (csID)
      deallocate (csSize)
   end subroutine MEF90EXOVecViewZonal_Private

#undef __FUNCT__
#define __FUNCT__ "MEF90EXOVecLoadZonal_Private"
   subroutine MEF90EXOVecLoadZonal_Private(v, exoid, step, offset, ierr)
      integer, intent(IN)                :: exoid
      PetscExodusIIInt, intent(IN)       :: step, offset
      type(tVec), intent(INOUT)          :: v
      PetscErrorCode, intent(INOUT)      :: ierr

      PetscExodusIIInt                   :: c
      PetscInt                           :: xs, xe, bs, numCS, set, csLocalSize, csxs
      PetscScalar, dimension(:), pointer :: varray
      PetscInt, dimension(:), pointer    :: csID, csSize
      type(tVec)                         :: vComp
      type(tIS)                          :: compIS
      character(len=MXSTLN)              :: elemType
      PetscMPIInt                        :: rank

      csxs = 0_ki
      PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
      numCS = exinqi(exoid, EX_INQ_ELEM_BLK)
      allocate (csID(numCS))
      allocate (csSize(numCS))
      call exgebi(exoid, csID, ierr)
      do set = 1, numCS
         call exgelb(exoid, csID(set), elemType, csSize(set), PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr)
      end do
      PetscCall(VecGetOwnershipRange(v, xs, xe, ierr))
      PetscCall(VecGetBlockSize(v, bs, ierr))
      if (bs > 1) then
         PetscCall(ISCreateStride(PETSC_COMM_WORLD, (xe - xs) / bs, xs, bs, compIS, ierr))
      end if
      do set = 1, numCS
         !  range of indices for set setID[set]: csxs:csxs + csSize[set]-1
         !  local slice of zonal values:         xs/bs,xm/bs-1
         !  intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1)
         csLocalSize = max(0, min(xe / bs, csxs + csSize(set)) - max(xs / bs, csxs))
         if (bs == 1) then
            PetscCall(VecGetArray(v, varray, ierr))
            call exgnev(exoid, step, offset, csID(set), csSize(set), max(xs - csxs, 0) + 1, csLocalSize, varray, ierr)
            PetscCall(VecRestoreArray(v, varray, ierr))
         else
            do c = 0, bs - 1
               PetscCall(ISStrideSetStride(compIS, (xe - xs) / bs, xs + c, bs, ierr))
               PetscCall(VecGetSubVector(v, compIS, vComp, ierr))
               PetscCall(VecGetArray(vComp, varray, ierr))
               call exgnev(exoid, step, offset + c, csID(set), 0_ki, max(xs / bs - csxs, 0) + 1, csLocalSize, varray(max(0, csxs - xs / bs) + 1:max(0, csxs - xs / bs) + csLocalSize), ierr)
               ! the 5th argument of exgnev is unused
               PetscCall(VecRestoreArray(vComp, varray, ierr))
               PetscCall(VecISCopy(v, compIS, SCATTER_FORWARD, vComp, ierr))
               PetscCall(VecRestoreSubVector(v, compIS, vComp, ierr))
            end do
         end if
         csxs = csxs + csSize(set)
      end do
      if (bs > 1) then
         PetscCall(ISDestroy(compIS, ierr))
      end if
      deallocate (csID)
      deallocate (csSize)
   end subroutine MEF90EXOVecLoadZonal_Private

#undef __FUNCT__
#define __FUNCT__ "MEF90EXOVecViewSide_Private"
   subroutine MEF90EXOVecViewSide_Private(v, exoid, step, offset, ierr)
      integer, intent(IN)                :: exoid
      PetscExodusIIInt, intent(IN)       :: step
      PetscInt, intent(IN)               :: offset
      type(tVec), intent(IN)             :: v
      PetscErrorCode, intent(INOUT)      :: ierr

      PetscInt                           :: xs, xe, bs, c, numSS, set, ssLocalSize, ssxs, sscs
      PetscScalar, dimension(:), pointer :: varray
      PetscInt, dimension(:), pointer    :: ssID, ssSize
      type(tVec)                         :: vComp
      type(tIS)                          :: compIS
      PetscMPIInt                        :: rank

      ssxs = 0_ki
      sscs = 0_ki
      PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
      numSS = exinqi(exoid, EX_INQ_SIDE_SETS)
      allocate (ssID(numSS))
      allocate (ssSize(numSS))
      call exgssi(exoid, ssID, ierr)
      do set = 1, numSS
         call exgsp(exoid, ssID(set), ssSize(set), PETSC_NULL_INTEGER, ierr)
      end do
      PetscCall(VecGetOwnershipRange(v, xs, xe, ierr))
      PetscCall(VecGetBlockSize(v, bs, ierr))
      do set = 1, numSS
         !  range of indices for set setID[set]: csxs:csxs + csSize[set]-1
         !  local slice of zonal values:         xs/bs,xm/bs-1
         !  intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1)
         ssLocalSize = max(0, min(xe / bs, ssxs + ssSize(set)) - max(xs / bs, ssxs))
         if (bs == 1) then
            PetscCall(VecGetArrayRead(v, varray, ierr))
            call expssv(exoid, step, offset, ssID(set), ssLocalSize, varray, ierr)
            PetscCall(VecRestoreArrayRead(v, varray, ierr))
         else
            PetscCall(ISCreateStride(PETSC_COMM_WORLD, ssLocalSize, xs + sscs, bs, compIS, ierr))
            do c = 0, bs - 1
               PetscCall(ISStrideSetStride(compIS, ssLocalSize, xs + sscs + c, bs, ierr))
               PetscCall(VecGetSubVector(v, compIS, vComp, ierr))
               PetscCall(VecGetArrayRead(vComp, varray, ierr))
               call exppv(exoid, step, EX_SIDE_SET, offset + c, ssID(set), max(xs / bs - ssxs, 0) + 1, ssLocalSize, varray, ierr)
               PetscCall(VecRestoreArrayRead(vComp, varray, ierr))
               PetscCall(VecRestoreSubVector(v, compIS, vComp, ierr))
            end do
            PetscCall(ISDestroy(compIS, ierr))
         end if
         ssxs = ssxs + ssSize(set)
         sscs = sscs + bs * ssLocalSize
      end do
      deallocate (ssID)
      deallocate (ssSize)
   end subroutine MEF90EXOVecViewSide_Private

#undef __FUNCT__
#define __FUNCT__ "MEF90EXOVecLoadSide_Private"
   subroutine MEF90EXOVecLoadSide_Private(v, exoid, step, offset, ierr)
      integer, intent(IN)                :: exoid
      PetscExodusIIInt, intent(IN)       :: step
      PetscInt, intent(IN)               :: offset
      type(tVec), intent(IN)             :: v
      PetscErrorCode, intent(INOUT)      :: ierr

      PetscInt                           :: xs, xe, bs, c, numSS, set, ssLocalSize, ssxs, sscs
      PetscScalar, dimension(:), pointer :: varray
      PetscInt, dimension(:), pointer    :: ssID, ssSize
      type(tVec)                         :: vComp
      type(tIS)                          :: compIS
      PetscMPIInt                        :: rank

      ssxs = 0_ki
      sscs = 0_ki
      PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
      numSS = exinqi(exoid, EX_INQ_SIDE_SETS)
      allocate (ssID(numSS))
      allocate (ssSize(numSS))
      call exgssi(exoid, ssID, ierr)
      do set = 1, numSS
         call exgsp(exoid, ssID(set), ssSize(set), PETSC_NULL_INTEGER, ierr)
      end do
      PetscCall(VecGetOwnershipRange(v, xs, xe, ierr))
      PetscCall(VecGetBlockSize(v, bs, ierr))
      do set = 1, numSS
         !  range of indices for set setID[set]: csxs:csxs + csSize[set]-1
         !  local slice of zonal values:         xs/bs,xm/bs-1
         !  intersection:                        max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1)
         ssLocalSize = max(0, min(xe / bs, ssxs + ssSize(set)) - max(xs / bs, ssxs))
         if (bs == 1) then
            PetscCall(VecGetArray(v, varray, ierr))
            PetscCall(exgssv(exoid, step, offset, ssID(set), ssLocalSize, varray, ierr))
            PetscCall(VecRestoreArray(v, varray, ierr))
         else
            PetscCall(ISCreateStride(PETSC_COMM_WORLD, ssLocalSize, xs + sscs, bs, compIS, ierr))
            do c = 0, bs - 1
               PetscCall(ISStrideSetStride(compIS, ssLocalSize, xs + sscs + c, bs, ierr))
               PetscCall(VecGetSubVector(v, compIS, vComp, ierr))
               PetscCall(VecGetArray(vComp, varray, ierr))
               call exgpv(exoid, step, EX_SIDE_SET, offset + c, ssID(set), max(xs / bs - ssxs, 0) + 1, ssLocalSize, varray, ierr)
               PetscCall(VecRestoreArray(vComp, varray, ierr))
               PetscCall(VecISCopy(v, compIS, SCATTER_FORWARD, vComp, ierr))
               PetscCall(VecRestoreSubVector(v, compIS, vComp, ierr))
            end do
            PetscCall(ISDestroy(compIS, ierr))
         end if
         ssxs = ssxs + ssSize(set)
         sscs = sscs + bs * ssLocalSize
      end do
      deallocate (ssID)
      deallocate (ssSize)
   end subroutine MEF90EXOVecLoadSide_Private
end module m_MEF90_EXO