! Copyright 2026 - David Minton
! This file is part of Swiftest.
! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License 
! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
! You should have received a copy of the GNU General Public License along with Swiftest. 
! If not, see: https://www.gnu.org/licenses. 

submodule (swiftest) s_swiftest_io
   use symba
   use ringmoons
   use netcdf
contains

   module subroutine swiftest_io_compact_output(self, param, timer)
      !! author: David Minton
      !!
      !! Generates the terminal output displayed when display_style is set to COMPACT. This is used by the Python driver to 
      !! make nice-looking progress reports.
      implicit none

      interface fmt
         !! author: David Minton
         !!
         !! Formats a pair of variables and corresponding values for the compact display output. Generic interface for different 
         !! variable types to format.
         procedure :: fmt_I4B, fmt_I8B, fmt_DP
      end interface

      ! Arguments
      class(swiftest_nbody_system), intent(in) :: self  
         !! Swiftest nbody system object   
      class(swiftest_parameters),   intent(in) :: param 
         !! Input colleciton of user-defined parameters
      class(*),                     intent(in) :: timer 
         !! Object used for computing elapsed wall time  (must be unlimited polymorphic because the walltimer module requires base)
      ! Internals
      character(len=:), allocatable :: formatted_output

      select type(timer)
      class is (walltimer)
         formatted_output = fmt("ILOOP",param%iloop) // fmt("T",self%t) // fmt("NPL",self%pl%nbody) // fmt("NTP",self%tp%nbody) 
         select type(pl => self%pl)
         class is (symba_pl)
            formatted_output = formatted_output // fmt("NPLM",pl%nplm)
         end select
         if (param%lenergy) then
            formatted_output = formatted_output // fmt("LTOTERR",self%L_total_error) // fmt("ETOTERR",self%te_error) &
                             // fmt("MTOTERR",self%Mtot_error) // fmt("KEOERR",self%ke_orbit_error) // fmt("PEERR",self%pe_error) &
                             // fmt("EORBERR",self%E_orbit_error) // fmt("EUNTRERR",self%E_untracked_error) &
                             // fmt("LESCERR",self%L_escape_error) // fmt("MESCERR",self%Mescape_error)
            if (param%lclose) formatted_output = formatted_output // fmt("ECOLLERR",self%Ecoll_error)
            if (param%lrotation) formatted_output = formatted_output // fmt("KEROTERR",self%ke_rot_error) &
                             // fmt("LROTERR",self%L_rot_error) 
         end if

         if (.not. timer%main_is_started) then ! This is the start of a new run
            formatted_output =  formatted_output // fmt("WT",0.0_DP) // fmt("IWT",0.0_DP) // fmt("WTPS",0.0_DP) 
         else
            formatted_output = formatted_output // fmt("WT",timer%wall_main) // fmt("IWT",timer%wall_step) &
                                               // fmt("WTPS",timer%wall_per_substep)
         end if
         write(*,*) formatted_output
      end select
      return

      contains

         function fmt_I4B(varname,val) result(pair_string)
            implicit none
            ! Arguments
            character(*), intent(in) :: varname !! The variable name of the pair
            integer(I4B), intent(in) :: val !! A 4-byte integer value
            ! Result
            character(len=:), allocatable :: pair_string
            ! Internals
            character(len=24) :: str_value
      
            write(str_value,*) val
            pair_string = trim(adjustl(varname)) // " " // trim(adjustl(str_value)) // ";"

            return 
         end function fmt_I4B

         function fmt_I8B(varname, val) result(pair_string)
            implicit none
            ! Arguments
            character(*), intent(in) :: varname !! The variable name of the pair
            integer(I8B), intent(in) :: val     !! An 8-byte integer value
            ! Result
            character(len=:), allocatable :: pair_string
            ! Internals
            character(len=24) :: str_value
      
            write(str_value,*) val
            pair_string = trim(adjustl(varname)) // " " // trim(adjustl(str_value)) // ";"

            return 
         end function fmt_I8B

         function fmt_DP(varname, val) result(pair_string)
            implicit none
            ! Arguments
            character(*), intent(in) :: varname !! The variable name of the pair
            real(DP),     intent(in) :: val     !! A double precision floating point value
            ! Result
            character(len=:), allocatable :: pair_string
            ! Internals
            character(len=24) :: str_value
      
            write(str_value,'(ES24.16)') val
            pair_string = trim(adjustl(varname)) // " " // trim(adjustl(str_value)) // ";"

            return 
         end function fmt_DP

   end subroutine swiftest_io_compact_output


   module subroutine swiftest_io_conservation_report(self, param, lterminal)
      !! author: The Purdue Swiftest Team -  David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
      !!
      !! Reports the current state of energy, mass, and angular momentum conservation in a run
      implicit none
      ! Arguments
      class(swiftest_nbody_system), intent(inout) :: self      !! Swiftest nbody system object
      class(swiftest_parameters),   intent(inout) :: param     !! Input colleciton of user-defined parameters
      logical,                      intent(in)    :: lterminal !! Indicates whether to output information to the terminal screen
      ! Internals
      real(DP), dimension(NDIM)       :: L_total_now,  L_orbit_now,  L_rot_now
      real(DP)                        :: ke_orbit_now,  ke_rot_now,  pe_now,  E_orbit_now, be_now, be_cb_now, be_cb_orig, te_now
      real(DP)                        :: GMtot_now
      character(len=*), parameter     :: EGYTERMFMT = '(" DL/L0 = ", ES12.5, "; DE/|E0| = ", ES12.5, "; DM/M0 = ", ES12.5)'

      associate(nbody_system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, display_unit => param%display_unit)
         call nbody_system%get_energy_and_momentum(param) 
         ke_orbit_now = nbody_system%ke_orbit
         ke_rot_now = nbody_system%ke_rot
         pe_now = nbody_system%pe
         be_now = nbody_system%be
         be_cb_now = nbody_system%be_cb
         te_now = nbody_system%te
         L_orbit_now(:) = nbody_system%L_orbit(:)
         L_rot_now(:) = nbody_system%L_rot(:)
         E_orbit_now = ke_orbit_now + pe_now
         L_total_now(:) = nbody_system%L_total(:) + nbody_system%L_escape(:)
         GMtot_now = nbody_system%GMtot + nbody_system%GMescape 

         if (param%lfirstenergy) then
            nbody_system%ke_orbit_orig = ke_orbit_now
            nbody_system%ke_rot_orig = ke_rot_now
            nbody_system%pe_orig = pe_now
            nbody_system%be_orig = be_now
            nbody_system%te_orig = te_now
            nbody_system%E_orbit_orig = E_orbit_now
            nbody_system%GMtot_orig = GMtot_now
            nbody_system%L_orbit_orig(:) = L_orbit_now(:)
            nbody_system%L_rot_orig(:) = L_rot_now(:)
            nbody_system%L_total_orig(:) = L_total_now(:)
            param%lfirstenergy = .false.
         end if

         if (.not.param%lfirstenergy) then 
            nbody_system%ke_orbit_error = (ke_orbit_now - nbody_system%ke_orbit_orig) / abs(nbody_system%te_orig)
            nbody_system%ke_rot_error = (ke_rot_now - nbody_system%ke_rot_orig) / abs(nbody_system%te_orig)
            nbody_system%pe_error = (pe_now - nbody_system%pe_orig) / abs(nbody_system%te_orig)

            be_cb_orig = -(3 * cb%GM0**2 / param%GU) / (5 * cb%R0)
            nbody_system%be_error = (be_now - nbody_system%be_orig) / abs(nbody_system%te_orig) + (be_cb_now - be_cb_orig) & 
                                    / abs(nbody_system%te_orig)

            if (abs(nbody_system%E_orbit_orig) < 10*tiny(1.0_DP)) then
               nbody_system%E_orbit_error = 0.0_DP
            else
               nbody_system%E_orbit_error = (E_orbit_now - nbody_system%E_orbit_orig) / abs(nbody_system%E_orbit_orig)
            end if
            nbody_system%Ecoll_error = nbody_system%E_collisions / abs(nbody_system%te_orig)
            nbody_system%E_untracked_error = nbody_system%E_untracked / abs(nbody_system%te_orig)
            nbody_system%te_error = (nbody_system%te - nbody_system%te_orig - nbody_system%E_collisions - nbody_system%E_untracked)&
                                   / abs(nbody_system%te_orig) + (be_cb_now - be_cb_orig) / abs(nbody_system%te_orig)

            nbody_system%L_orbit_error = norm2(L_orbit_now(:) - nbody_system%L_orbit_orig(:)) / norm2(nbody_system%L_total_orig(:))
            nbody_system%L_rot_error = norm2(L_rot_now(:) - nbody_system%L_rot_orig(:)) / norm2(nbody_system%L_total_orig(:))
            nbody_system%L_escape_error = norm2(nbody_system%L_escape(:)) / norm2(nbody_system%L_total_orig(:))
            nbody_system%L_total_error = norm2(L_total_now(:) - nbody_system%L_total_orig(:)) / norm2(nbody_system%L_total_orig(:))

            nbody_system%Mescape_error = nbody_system%GMescape / nbody_system%GMtot_orig
#ifdef COARRAY
   if (this_image() == 1 .or. param%log_output) then
#endif 
            if (lterminal) then
               write(display_unit, EGYTERMFMT) nbody_system%L_total_error, nbody_system%te_error, nbody_system%Mtot_error
               if (param%log_output) flush(display_unit)
            end if

#ifdef COARRAY
   end if ! (this_image() == 1) then
#endif 

            if (abs(nbody_system%Mtot_error) > 100 * epsilon(nbody_system%Mtot_error)) then
               write(*,*) "Severe error! Mass not conserved! Halting!"
               ! Save the frame of data to the bin file in the slot just after the present one for diagnostics
               call base_util_exit(FAILURE,param%display_unit)
            end if
         end if
      end associate

      return
   end subroutine swiftest_io_conservation_report


   module subroutine swiftest_io_display_run_information(self, param, integration_timer, phase)
      !! author:  David A. Minton
      !!
      !! Displays helpful information while a run is executing. The format of the output depends on user parameters
      implicit none
      ! Arguments
      class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object
      class(swiftest_parameters),   intent(inout) :: param !! Current run configuration parameters 
      type(walltimer),              intent(inout) :: integration_timer !! Object used for computing elapsed wall time
      character(len=*), optional,   intent(in)    :: phase !! One of "first" or "last" 
      ! Internals
      integer(I4B)                                   :: phase_val   
      real(DP)                                       :: tfrac             !! Fraction of total simulation time completed
      type(progress_bar), save                       :: pbar              !! Object used to print out a progress bar
      character(len=64)                              :: pbarmessage
#ifdef COARRAY
      character(*), parameter :: co_statusfmt = '("Image: ",I4, "; Time = ", ES12.5, "; fraction done = ", F6.3, ' // & 
                                             '"; Number of active pl, tp = ", I6, ", ", I6)'
      character(*), parameter :: co_symbastatfmt = '("Image: ",I4, "; Image: Time = ", ES12.5, "; fraction done = ", F6.3, ' // &
                                                '"; Number of active pl, plm, tp = ", I6, ", ", I6, ", ", I6)'
      character(*), parameter :: co_ringmoonstatfmt = '("Image: ",I4, "; Image: Time = ", ES12.5, "; fraction done = ", F6.3, ' // &
                                                '"; Number of active pl, plm, tp = ", I6, ", ", I6, ", ", I6, ", ", I6)'
#endif
      character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & 
                                             '"; Number of active pl, tp = ", I6, ", ", I6)'
      character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // &
                                                '"; Number of active pl, plm, tp = ", I6, ", ", I6, ", ", I6)'
      character(*), parameter :: ringmoonstatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // &
                                                '"; Number of active pl, plm, tp, seeds = ", I6, ", ", I6, ", ", I6, ", ", I6)'
      character(*), parameter :: pbarfmt = '("Time = ", ES12.5," of ",ES12.5)'

! The following will syncronize the images so that they report in order, and only write to file one at at ime


      phase_val = 1
      if (present(phase)) then
         if (phase == "first") then
            phase_val = 0
         else if (phase == "last") then
            phase_val = -1
         end if
      end if

      tfrac = (self%t - param%t0) / (param%tstop - param%t0)

#ifdef COARRAY
         if (this_image() == 1 .or. param%log_output) then
#endif
         if (phase_val == 0) then
            if (param%lrestart) then
               write(param%display_unit, *) " *************** Swiftest " // VERSION // " restart " // &
                                               trim(adjustl(param%integrator)) // " *************** "
            else
               write(param%display_unit, *) " *************** Swiftest " // VERSION // " start " // &
                                              trim(adjustl(param%integrator)) // " *************** "
            end if
            if (param%display_style == "PROGRESS") then
               call pbar%reset(param%nloops)
            else if (param%display_style == "COMPACT") then
               write(*,*) "SWIFTEST START " // trim(adjustl(param%integrator))
            end if
         end if
#ifdef COARRAY
         end if !(this_image() == 1)
#endif

      if (param%display_style == "PROGRESS") then
#ifdef COARRAY
         if (this_image() == 1) then
#endif
            write(pbarmessage,fmt=pbarfmt) self%t, param%tstop
            call pbar%update(param%iloop,message=pbarmessage)
#ifdef COARRAY
         end if !(this_image() == 1)
#endif
      else if (param%display_style == "COMPACT") then
         call self%compact_output(param,integration_timer)
      end if

      select type(nbody_system => self)
      class is (ringmoons_nbody_system)
#ifdef COARRAY
         if (param%lcoarray) then
            write(param%display_unit, co_ringmoonstatfmt) this_image(),nbody_system%t, tfrac, nbody_system%pl%nbody, nbody_system%pl%nplm, nbody_system%tp%nbody, nbody_system%seed%nbody
         else
#endif
            write(param%display_unit, ringmoonstatfmt) nbody_system%t, tfrac, nbody_system%pl%nbody, nbody_system%pl%nplm, nbody_system%tp%nbody, nbody_system%seed%nbody
#ifdef COARRAY
         end if
#endif
      class is (symba_nbody_system)
#ifdef COARRAY
         if (param%lcoarray) then
            write(param%display_unit, co_symbastatfmt) this_image(),nbody_system%t, tfrac, nbody_system%pl%nbody, nbody_system%pl%nplm, nbody_system%tp%nbody
         else
#endif
            write(param%display_unit, symbastatfmt) nbody_system%t, tfrac, nbody_system%pl%nbody, nbody_system%pl%nplm, nbody_system%tp%nbody
#ifdef COARRAY
         end if
#endif
      class default
#ifdef COARRAY
         if (param%lcoarray) then
            write(param%display_unit, co_statusfmt) this_image(),nbody_system%t, tfrac, nbody_system%pl%nbody, nbody_system%tp%nbody
         else
#endif
            write(param%display_unit, statusfmt) nbody_system%t, tfrac, nbody_system%pl%nbody, nbody_system%tp%nbody
#ifdef COARRAY
         end if
#endif
      end select

#ifdef COARRAY
      if (this_image() == num_images() .or. param%log_output) then
#endif
         if (phase_val == -1) then
            write(param%display_unit, *)" *************** Swiftest stop " // trim(adjustl(param%integrator)) // " *************** "
            if (param%display_style == "COMPACT") write(*,*) "SWIFTEST STOP" // trim(adjustl(param%integrator))
            if (param%log_output) close(param%display_unit)
         else
            if (param%log_output) flush(param%display_unit)
         end if

#ifdef COARRAY
      end if ! this_image() == num_images()
#endif

      return
   end subroutine swiftest_io_display_run_information


   module subroutine swiftest_io_dump_param(self, param_file_name)
      !! author: David A. Minton
      !!
      !! Dump integration parameters to file
      !!
      !! Adapted from David E. Kaufmann's Swifter routine io_dump_param.f90
      !! Adapted from Martin Duncan's Swift routine io_dump_param.f
      implicit none
      ! Arguments
      class(swiftest_parameters),intent(in) :: self    !! Output collection of parameters
      character(len=*),          intent(in) :: param_file_name !! Parameter input file name (i.e. param.in)
      ! Internals
      character(STRMAX)        :: errmsg !! Error message in UDIO procedure
      integer(I4B)             :: ierr

      open(unit = LUN, file = param_file_name, status='replace', form = 'FORMATTED', err = 667, iomsg = errmsg)
      !! todo: Currently this procedure does not work in user-defined derived-type input mode 
      !!    due to compiler incompatabilities
      !write(LUN,'(DT)') param
      call self%writer(LUN, iotype = "none", v_list = [0], iostat = ierr, iomsg = errmsg)
      if (ierr == 0) then
         close(LUN, err = 667, iomsg = errmsg)
         return
      end if

      667 continue
      write(*,*) "Error opening parameter dump file " // trim(adjustl(errmsg))
      call base_util_exit(FAILURE,self%display_unit)
   end subroutine swiftest_io_dump_param


   module subroutine swiftest_io_dump_system(self, param, system_history)
      !! author: David A. Minton
      !!
      !! Dumps the state of the nbody_system to files in case the simulation is interrupted.
      !! As a safety mechanism, there are two dump files that are written in alternating order
      !! so that if a dump file gets corrupted during writing, the user can restart from the older one.
      implicit none
      ! Arguments
      class(swiftest_nbody_system), intent(inout) :: self  !! Swiftest nbody_system object
      class(swiftest_parameters),   intent(inout) :: param !! Current run configuration parameters 
      class(swiftest_storage),      intent(inout) :: system_history    !! Stores the system history between output dumps
      ! Internals
      class(swiftest_parameters), allocatable :: param_restart !! Local parameters variable used to parameters change input file 
                                                               !! names to dump file-specific values without changing the 
                                                               !! user-defined values
      character(len=:), allocatable :: param_file_name
      character(len=STRMAX) :: time_text
   
      ! Dump the encounter history if necessary
      if (param%lenc_save_trajectory &
         .or. param%lenc_save_closest &
         .and. allocated(self%encounter_history)) &
            call self%encounter_history%dump(param)
      if (allocated(self%collision_history)) call self%collision_history%dump(param)

      ! Dump the nbody_system history to file
      call system_history%dump(param)

#ifdef COARRAY
      if (this_image() == 1) then
#endif 
         allocate(param_restart, source=param)
         param_restart%in_form  = "XV"
         param_restart%out_stat = 'APPEND'
         param_restart%in_type = "NETCDF_DOUBLE"
         param_restart%nc_in = param%outfile
         param_restart%lrestart = .true.
         param_restart%tstart = self%t
         param_file_name    = trim(adjustl(PARAM_RESTART_FILE))
         call param_restart%dump(param_file_name)
         write(time_text,'(I0.20)') param%iloop
         param_file_name = "param." // trim(adjustl(time_text)) // ".in"
         call param_restart%dump(param_file_name)
#ifdef COARRAY
      end if ! (this_image() == 1) 
#endif 
      return
   end subroutine swiftest_io_dump_system


   module subroutine swiftest_io_dump_storage(self, param)
      !! author: David A. Minton
      !!
      !! Dumps the time history of the simulation to file. Each time it writes a frame to file, it deallocates the nbody_system
      !! object from inside. It will only dump frames with systems that are allocated, so this can be called at the end of
      !! a simulation for cases when the number of saved frames is not equal to the dump cadence (for instance, if the dump
      !! cadence is not divisible by the total number of loops).
      implicit none
      ! Arguments
      class(swiftest_storage),    intent(inout) :: self   !! Swiftest simulation history storage object
      class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters 
      ! Internals
      integer(I4B) :: i
#ifdef COARRAY
      type(walltimer) :: iotimer
#endif

      if (self%iframe == 0) return
      call self%make_index_map()
      associate(nc => self%nc)
#ifdef COARRAY
         sync all
         write(param%display_unit,*) "File output started"
         call iotimer%start()
         critical
#endif
         call nc%open(param)
         do i = 1, self%iframe
            ! Writing files is more efficient if we write out the common frames from each image before going to the next frame
            if (allocated(self%frame(i)%item)) then
               select type(nbody_system => self%frame(i)%item)
               class is (swiftest_nbody_system)
                  call nbody_system%write_frame(nc, param)
               end select
               deallocate(self%frame(i)%item)
            end if
         end do
         call nc%close()
#ifdef COARRAY  
         end critical
         call iotimer%stop()
         sync all
         call iotimer%report(message="File output :", unit=param%display_unit)
         flush(param%display_unit)
#endif
      end associate

      call self%reset()
      return
   end subroutine swiftest_io_dump_storage


   module subroutine swiftest_io_get_args(integrator, param_file_name, display_style, from_cli)
      !! author: David A. Minton
      !!
      !! Reads in the name of the parameter file from command line arguments. 
      implicit none
      ! Arguments
      character(len=:), intent(inout), allocatable :: integrator      !! Symbolic code of the requested integrator  
      character(len=:), intent(inout), allocatable :: param_file_name !! Name of the input parameters file
      character(len=:), intent(inout), allocatable :: display_style   !! Style of the output display 
                                                                      !! {"CLASSIC", "COMPACT", "PROGRESS", "QUIET"}). 
                                                                      !! Default is "CLASSIC"
      logical,          intent(in)                 :: from_cli        !! If true, get command-line arguments. Otherwise, use the 
                                                                      !! values of the input variables
      ! Internals
      character(len=STRMAX), dimension(:), allocatable :: arg
      character(len=STRMAX), dimension(3)              :: tmp_arg
      integer(I4B), dimension(:), allocatable :: ierr
      integer :: i,narg

      if (from_cli) then
         narg = command_argument_count() 
         if (narg > 0) then
            allocate(arg(narg),ierr(narg))
            do i = 1,narg
               call get_command_argument(i, arg(i), status = ierr(i))
            end do
            if (any(ierr /= 0)) call base_util_exit(USAGE)
         else
            call base_util_exit(USAGE)
         end if
      else
         narg = 0
         if (len(integrator) > 0) then
            narg = narg + 1
            tmp_arg(narg) = integrator
         end if
         if (len(param_file_name) > 0) then
            narg = narg + 1
            tmp_arg(narg) = param_file_name
         end if
         if (len(display_style) > 0) then
            narg = narg + 1
            tmp_arg(narg) = display_style
         end if
         allocate(arg(narg))
         arg(1:narg) = tmp_arg(1:narg)
         deallocate(integrator, param_file_name, display_style)
      end if
   
      if (narg == 1) then
         if (arg(1) == '-v' .or. arg(1) == '--version') then
            call swiftest_util_version() 
         else if (arg(1) == '-h' .or. arg(1) == '--help') then
            call base_util_exit(HELP)
         else
            call base_util_exit(USAGE)
         end if
      else if (narg >= 2) then
         call swiftest_io_toupper(arg(1))
         select case(arg(1))
         case('BS')
            integrator = INT_BS
         case('HELIO')
            integrator = INT_HELIO
         case('RA15')
            integrator = INT_RA15
         case('TU4')
            integrator = INT_TU4
         case('WHM')
            integrator = INT_WHM
         case('RMVS')
            integrator = INT_RMVS
         case('SYMBA')
            integrator = INT_SYMBA
         case('RINGMOONS')
            integrator = INT_RINGMOONS
         case default
            integrator = UNKNOWN_INTEGRATOR
            write(*,*) trim(adjustl(arg(1))) // ' is not a valid integrator.'
            call base_util_exit(USAGE)
         end select
         param_file_name = trim(adjustl(arg(2)))
      end if

      if (narg == 2) then
         display_style = "PROGRESS"
      else if (narg == 3) then
         call swiftest_io_toupper(arg(3))
         display_style = trim(adjustl(arg(3)))
      else
         call base_util_exit(USAGE)
      end if

      return
   end subroutine swiftest_io_get_args


   module function swiftest_io_get_token(buffer, ifirst, ilast, ierr) result(token)
      !! author: David A. Minton
      !!
      !! Retrieves a character token from an input string. Here a token is defined as any set of contiguous non-blank characters not
      !! beginning with or containing "!". If "!" is present, any remaining part of the buffer including the "!" is ignored
      !!
      !! Adapted from David E. Kaufmann's Swifter routine swiftest_io_get_token.f90
      implicit none
      ! Arguments
      character(len=*), intent(in)    :: buffer         !! Input string buffer
      integer(I4B),     intent(inout) :: ifirst         !! Index of the buffer at which to start the search for a token
      integer(I4B),     intent(out)   :: ilast          !! Index of the buffer at the end of the returned token
      integer(I4B),     intent(out)   :: ierr           !! Error code
      ! Result
      character(len=:), allocatable   :: token          !! Returned token string
      ! Internals
      integer(I4B) :: i,ilength
   
      ilength = len(buffer)
   
      if (ifirst > ilength) then
          ilast = ifirst
          ierr = -1 !! Bad input
          token = ''
          return
      end if
      do i = ifirst, ilength
          if (buffer(i:i) /= ' ') exit
      end do
      if ((i > ilength) .or. (buffer(i:i) == '!')) then
          ifirst = i
          ilast = i
          ierr = -2 !! No valid token
          token = ''
          return
      end if
      ifirst = i
      do i = ifirst, ilength
          if ((buffer(i:i) == ' ') .or. (buffer(i:i) == '!')) exit
      end do
      ilast = i - 1
      ierr = 0
   
      token = buffer(ifirst:ilast)

      return
   end function swiftest_io_get_token


   module subroutine swiftest_io_log_one_message(file, message)
      !! author: David A. Minton
      !!
      !! Writes a single message to a log file
      implicit none
      ! Arguments
      character(len=*), intent(in) :: file   !! Name of file to log
      character(len=*), intent(in) :: message
      ! Internals
      character(STRMAX) :: errmsg

      open(unit=LUN, file=trim(adjustl(file)), status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg)
      write(LUN, *) trim(adjustl(message)) 
      close(LUN)

      return
      667 continue
      write(*,*) "Error writing message to log file: " // trim(adjustl(errmsg))
   end subroutine swiftest_io_log_one_message


   module subroutine swiftest_io_log_start(param, file, header)
      !! author: David A. Minton
      !!
      !! Checks to see if a log file needs to be created if this is a new run, or appended if this is a restarted run
      implicit none
      ! Arguments
      class(swiftest_parameters), intent(in) :: param  !! Current Swiftest run configuration parameters
      character(len=*),           intent(in) :: file   !! Name of file to log
      character(len=*),           intent(in) :: header !! Header to print at top of log file
      ! Internals
      character(STRMAX) :: errmsg
      logical           :: fileExists

      inquire(file=trim(adjustl(file)), exist=fileExists)
      if (.not.param%lrestart .or. .not.fileExists) then
         open(unit=LUN, file=file, status="REPLACE", err = 667, iomsg = errmsg)
         write(LUN, *, err = 667, iomsg = errmsg) trim(adjustl(header))
      end if
      close(LUN)

      return

      667 continue
      write(*,*) "Error writing log file: " // trim(adjustl(errmsg))
   end subroutine swiftest_io_log_start


   module subroutine swiftest_io_netcdf_flush(self, param)
      !! author: David A. Minton
      !!
      !! Flushes the current buffer to disk by closing and re-opening the file.
      !!    
      implicit none
      ! Arguments
      class(swiftest_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset
      class(swiftest_parameters),        intent(inout) :: param !! Current run configuration parameters 

      call self%close()
      call self%open(param,readonly=.false.)

      return
   end subroutine swiftest_io_netcdf_flush


   module subroutine swiftest_io_netcdf_get_t0_values_system(self, nc, param) 
      !! author: David A. Minton
      !!
      !! Gets the t0 values of various parameters such as energy and momentum
      !!
      implicit none
      ! Arguments
      class(swiftest_nbody_system),      intent(inout) :: self
      class(swiftest_netcdf_parameters), intent(inout) :: nc     !! Parameters used to identify a particular NetCDF dataset
      class(swiftest_parameters),        intent(inout) :: param
      ! Internals
      integer(I4B)                              :: itmax, idmax, tslot, status
      real(DP), dimension(:), allocatable       :: vals
      logical, dimension(:), allocatable        :: plmask, tpmask
      real(DP), dimension(1)                    :: rtemp
      real(DP), dimension(NDIM)                 :: rot0, Ip0, L
      real(DP) :: mass0

      associate (cb => self%cb)
         call nc%open(param, readonly=.true.)
         call nc%find_tslot(param%t0, tslot)
         call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%time_dimid, len=itmax), "netcdf_io_get_t0_values_system time_dimid")
         call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), "netcdf_io_get_t0_values_system name_dimid")
         allocate(vals(idmax))
         call netcdf_io_check( nf90_get_var(nc%id, nc%time_varid, rtemp, start=[tslot], count=[1]), &
                              "netcdf_io_get_t0_values_system time_varid" )

         if (param%lenergy) then
            call netcdf_io_check( nf90_get_var(nc%id, nc%KE_orb_varid, rtemp, start=[tslot], count=[1]), & 
                                  "netcdf_io_get_t0_values_system KE_orb_varid" )
            self%ke_orbit_orig = rtemp(1)

            call netcdf_io_check( nf90_get_var(nc%id, nc%KE_rot_varid, rtemp, start=[tslot], count=[1]), &
                                 "netcdf_io_get_t0_values_system KE_rot_varid" )
            self%ke_rot_orig = rtemp(1)

            call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, rtemp, start=[tslot], count=[1]), &
                                  "netcdf_io_get_t0_values_system PE_varid" )
            self%pe_orig = rtemp(1)

            call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, rtemp, start=[tslot], count=[1]), &
                                  "netcdf_io_get_t0_values_system BE_varid" )
            self%be_orig = rtemp(1)
            
            call netcdf_io_check( nf90_get_var(nc%id, nc%TE_varid, rtemp, start=[tslot], count=[1]), &
                                  "netcdf_io_get_t0_values_system TE_varid" )
            self%te_orig = rtemp(1)

            self%E_orbit_orig = self%ke_orbit_orig + self%pe_orig

            call netcdf_io_check( nf90_get_var(nc%id, nc%L_orbit_varid, self%L_orbit_orig(:), start=[1,tslot], count=[NDIM,1]), &
                                  "netcdf_io_get_t0_values_system L_orbit_varid" )
            call netcdf_io_check( nf90_get_var(nc%id, nc%L_rot_varid, self%L_rot_orig(:), start=[1,tslot], count=[NDIM,1]), &
                                  "netcdf_io_get_t0_values_system L_rot_varid" )

            self%L_total_orig(:) = self%L_orbit_orig(:) + self%L_rot_orig(:) 

            call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, vals, start=[1,tslot], count=[idmax,1]), &
                                  "netcdf_io_get_t0_values_system Gmass_varid" )
            call nc%get_valid_masks(plmask,tpmask)
            self%GMtot_orig = vals(1) + sum(vals(1:idmax), plmask(:)) 

            cb%GM0 = vals(1)
            cb%dGM = cb%Gmass - cb%GM0


            status = nf90_inq_varid(nc%id, nc%mass_varname, nc%mass_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%mass_varid,  vals, start=[1,tslot], count=[idmax,1]), &
                                    "netcdf_io_get_t0_values_system mass_varid" )
               mass0 = vals(1)
            else
               mass0 = cb%GM0 / param%GU
            end if


            call netcdf_io_check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1,tslot], count=[1,1]), &
                                  "netcdf_io_get_t0_values_system radius_varid" )
            cb%R0 = rtemp(1) 

            if (param%lrotation) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%rot_varid, rot0, start=[1,1,tslot], count=[NDIM,1,1]), &
                                     "netcdf_io_get_t0_values_system rot_varid" )
               call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid, Ip0, start=[1,1,tslot], count=[NDIM,1,1]), &
                                     "netcdf_io_get_t0_values_system Ip_varid" )
               cb%L0(:) = Ip0(3) * mass0 * cb%R0**2 * rot0(:) * DEG2RAD
               L(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:) * DEG2RAD
               cb%dL(:) = L(:) - cb%L0

               status = nf90_inq_varid(nc%id, nc%j2rp2_varname, nc%j2rp2_varid)
               if (status == NF90_NOERR) then
                  call netcdf_io_check( nf90_get_var(nc%id, nc%j2rp2_varid, cb%j2rp2, start=[tslot]), &
                                        "netcdf_io_get_t0_values_system nf90_getvar j2rp2_varid"  )
               else 
                  cb%j2rp2 = 0.0_DP
               end if
      
               status = nf90_inq_varid(nc%id, nc%j4rp4_varname, nc%j4rp4_varid)   
               if (status == NF90_NOERR) then      
                  call netcdf_io_check( nf90_get_var(nc%id, nc%j4rp4_varid, cb%j4rp4, start=[tslot]), &
                                        "netcdf_io_get_t0_values_system nf90_getvar j4rp4_varid"  )
               else 
                  cb%j4rp4 = 0.0_DP
               end if

            end if

            ! Retrieve the current bookkeeping variables
            call nc%find_tslot(self%t, tslot)
            call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%L_escape(:),  start=[1,tslot], count=[NDIM,1]), &
                                  "netcdf_io_get_t0_values_system L_escape_varid" )
            call netcdf_io_check( nf90_get_var(nc%id, nc%GMescape_varid,    self%GMescape,    start=[tslot]), &
                                  "netcdf_io_get_t0_values_system GMescape_varid" )
            call netcdf_io_check( nf90_get_var(nc%id, nc%E_collisions_varid, self%E_collisions, start=[tslot]), &
                                  "netcdf_io_get_t0_values_system E_collisions_varid" )
            call netcdf_io_check( nf90_get_var(nc%id, nc%E_untracked_varid,  self%E_untracked,  start=[tslot]), &
                                   "netcdf_io_get_t0_values_system E_untracked_varid" )

         end if

         deallocate(vals)
         call nc%close()
      end associate
      
      return
   end subroutine swiftest_io_netcdf_get_t0_values_system


   module subroutine swiftest_io_netcdf_initialize_output(self, param)
      !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton
      !!
      !! Initialize a NetCDF file nbody_system and define all variables.
      use, intrinsic :: ieee_arithmetic
      implicit none
      ! Arguments
      class(swiftest_netcdf_parameters), intent(inout) :: self 
         !! Parameters used to for writing a NetCDF dataset to file
      class(swiftest_parameters),        intent(in)    :: param 
         !! Current run configuration parameters 
      ! Internals
      integer(I4B) :: nvar, varid, vartype
      real(DP) :: dfill
      real(SP) :: sfill
      integer(I4B), parameter :: NO_FILL = 0
      logical :: fileExists
      character(len=STRMAX) :: errmsg
      integer(I4B) :: i, ndims

      associate(nc => self)

         dfill = ieee_value(dfill, IEEE_QUIET_NAN)
         sfill = ieee_value(sfill, IEEE_QUIET_NAN)

         select case (param%out_type)
         case("NETCDF_FLOAT")
            nc%out_type = NF90_FLOAT
         case("NETCDF_DOUBLE")
            nc%out_type = NF90_DOUBLE
         case default
            write(*,*) trim(adjustl(param%out_type)), " is an invalid OUT_TYPE"
         end select

         ! Check if the file exists, and if it does, delete it
         inquire(file=nc%file_name, exist=fileExists)
         if (fileExists) then
            open(unit=LUN, file=nc%file_name, status="old", err=667, iomsg=errmsg)
            close(unit=LUN, status="delete")
         end if

         ! Create the file
         call netcdf_io_check( nf90_create(nc%file_name, NF90_NETCDF4, nc%id), "netcdf_io_initialize_output nf90_create" )
         nc%lfile_is_open = .true.

         ! Dimensions
         call netcdf_io_check( nf90_def_dim(nc%id, nc%time_dimname, NF90_UNLIMITED, nc%time_dimid), &
                               "netcdf_io_initialize_output nf90_def_dim time_dimid" ) ! Simulation time dimension
         call netcdf_io_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), &
                               "netcdf_io_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension
         call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, NF90_UNLIMITED, nc%name_dimid), &
                               "netcdf_io_initialize_output nf90_def_dim name_dimid" ) ! dimension to store particle id numbers
         call netcdf_io_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), &
                               "netcdf_io_initialize_output nf90_def_dim str_dimid"  )  ! Dimension for string variables 
        

         ! Dimension coordinates
         call netcdf_io_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, nc%time_dimid, nc%time_varid), &
                              "netcdf_io_initialize_output nf90_def_var time_varid"  )
         call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), &
                               "netcdf_io_initialize_output nf90_def_var space_varid"  )
         call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), &
                               "netcdf_io_initialize_output nf90_def_var name_varid"  )

         ! Variables
         call netcdf_io_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), &
                               "netcdf_io_initialize_output nf90_def_var id_varid"  )
         call netcdf_io_check( nf90_def_var(nc%id, nc%status_varname, NF90_INT, [nc%name_dimid, nc%time_dimid], nc%status_varid), &
                               "netcdf_io_initialize_output nf90_def_var status_varid"  )
         call netcdf_io_check( nf90_def_var(nc%id, nc%npl_varname, NF90_INT, nc%time_dimid, nc%npl_varid), &
                               "netcdf_io_initialize_output nf90_def_var npl_varid"  )
         call netcdf_io_check( nf90_def_var(nc%id, nc%ntp_varname, NF90_INT, nc%time_dimid, nc%ntp_varid), &
                               "netcdf_io_initialize_output nf90_def_var ntp_varid"  )
         if (param%lmtiny_pl) call netcdf_io_check( nf90_def_var(nc%id, nc%nplm_varname, NF90_INT, nc%time_dimid, nc%nplm_varid), &
                               "netcdf_io_initialize_output nf90_def_var nplm_varid"  )
         call netcdf_io_check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%ptype_varid), &
                               "netcdf_io_initialize_output nf90_def_var ptype_varid"  )

         if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then
            call netcdf_io_check( nf90_def_var(nc%id, nc%rh_varname,  nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], &
                                  nc%rh_varid), "netcdf_io_initialize_output nf90_def_var rh_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%vh_varname,  nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], &
                                 nc%vh_varid), "netcdf_io_initialize_output nf90_def_var vh_varid"  )

            !! When GR is enabled, we need to save the pseudovelocity vectors in addition to the true heliocentric velocity vectors,
            !! otherwise we cannnot expect bit-identical runs from restarted runs with GR enabled due to floating point errors 
            !! during the conversion.
            if (param%lgr) then
               call netcdf_io_check( nf90_def_var(nc%id, nc%gr_pseudo_vh_varname, nc%out_type, &
                                     [nc%space_dimid, nc%name_dimid, nc%time_dimid], nc%gr_pseudo_vh_varid), &
                                     "netcdf_io_initialize_output nf90_def_var gr_psuedo_vh_varid"  )
               nc%lpseudo_vel_exists = .true.
            end if

         end if
      
         if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then
            call netcdf_io_check( nf90_def_var(nc%id, nc%a_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%a_varid), &
                                  "netcdf_io_initialize_output nf90_def_var a_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%e_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%e_varid), &
                                  "netcdf_io_initialize_output nf90_def_var e_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%inc_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%inc_varid), &
                                  "netcdf_io_initialize_output nf90_def_var inc_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%capom_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], &
                                               nc%capom_varid), &
                                  "netcdf_io_initialize_output nf90_def_var capom_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%omega_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], &
                                               nc%omega_varid), &
                                  "netcdf_io_initialize_output nf90_def_var omega_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%capm_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], &
                                               nc%capm_varid), &
                                  "netcdf_io_initialize_output nf90_def_var capm_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%varpi_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], &
                                               nc%varpi_varid), &
                                  "netcdf_io_initialize_output nf90_def_var varpi_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%lam_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%lam_varid), &
                                  "netcdf_io_initialize_output nf90_def_var lam_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%f_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%f_varid), &
                                  "netcdf_io_initialize_output nf90_def_var f_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%cape_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%cape_varid),&
                                  "netcdf_io_initialize_output nf90_def_var cape_varid"  )
         end if

         call netcdf_io_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), &
                                  "netcdf_io_initialize_output nf90_def_var Gmass_varid"  )
         call netcdf_io_check( nf90_def_var(nc%id, nc%mass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%mass_varid), &
                                  "netcdf_io_initialize_output nf90_def_var mass_varid"  )
         call netcdf_io_check( nf90_def_var(nc%id, nc%rhill_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%rhill_varid), &
                                  "netcdf_io_initialize_output nf90_def_var rhill_varid"  )

         if (param%lclose) then
            call netcdf_io_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], &
                                               nc%radius_varid), &
                                  "netcdf_io_initialize_output nf90_def_var radius_varid"  )

            call netcdf_io_check( nf90_def_var(nc%id, nc%origin_time_varname, nc%out_type, nc%name_dimid, nc%origin_time_varid), &
                                  "netcdf_io_initialize_output nf90_def_var origin_time_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%origin_type_varname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], &
                                    nc%origin_type_varid), &
                                  "netcdf_io_initialize_output nf90_create"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%origin_rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], &
                                               nc%origin_rh_varid), &
                                  "netcdf_io_initialize_output nf90_def_var origin_rh_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%origin_vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], &
                                               nc%origin_vh_varid), &
                                  "netcdf_io_initialize_output nf90_def_var origin_vh_varid"  )

            call netcdf_io_check( nf90_def_var(nc%id, nc%collision_id_dimname, NF90_INT, nc%name_dimid, nc%collision_id_varid), &
                                  "netcdf_io_initialize_output nf90_def_var collision_id_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%discard_time_varname, nc%out_type, nc%name_dimid, nc%discard_time_varid), &
                                  "netcdf_io_initialize_output nf90_def_var discard_time_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%discard_rh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], &
                                               nc%discard_rh_varid), &
                                  "netcdf_io_initialize_output nf90_def_var discard_rh_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%discard_vh_varname, nc%out_type, [nc%space_dimid, nc%name_dimid], &
                                               nc%discard_vh_varid), &
                                  "netcdf_io_initialize_output nf90_def_var discard_vh_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%discard_body_id_varname, NF90_INT, nc%name_dimid, &
                                               nc%discard_body_id_varid), &
                                  "netcdf_io_initialize_output nf90_def_var discard_body_id_varid"  )
         end if

         if (param%lrotation) then
            call netcdf_io_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], &
                                  nc%Ip_varid), "netcdf_io_initialize_output nf90_def_var Ip_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type, [nc%space_dimid, nc%name_dimid, nc%time_dimid], &
                                   nc%rot_varid), "netcdf_io_initialize_output nf90_def_var rot_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%rotphase_varname, nc%out_type, nc%time_dimid, nc%rotphase_varid), &
                                  "netcdf_io_initialize_output nf90_def_var rotphase_varid"  )


            if (nc%lc_lm_exists) then
               call netcdf_io_check( nf90_def_dim(nc%id, nc%sign_dimname, 2, nc%sign_dimid), &
                                    "swiftest_io_netcdf_open nf90_def_dim sign_dimid")
               call netcdf_io_check( nf90_def_dim(nc%id, nc%l_dimname, nc%l_dim_max, nc%l_dimid), &
                                    "swiftest_io_netcdf_open nf90_def_dim l_dimid")
               call netcdf_io_check( nf90_def_dim(nc%id, nc%m_dimname, nc%m_dim_max, nc%m_dimid), &
                                    "swiftest_io_netcdf_open nf90_def_dim m_dimid")

               call netcdf_io_check( nf90_def_var(nc%id, nc%sign_dimname, nc%out_type, nc%sign_dimid, nc%sign_varid), &
                                    "swiftest_io_netcdf_open nf90_def_var sign_varid")
               call netcdf_io_check( nf90_def_var(nc%id, nc%l_dimname, nc%out_type, nc%l_dimid, nc%l_varid), &
                                    "swiftest_io_netcdf_open nf90_def_var l_varid")
               call netcdf_io_check( nf90_def_var(nc%id, nc%m_dimname, nc%out_type, nc%m_dimid, nc%m_varid), &
                                    "swiftest_io_netcdf_open nf90_def_var m_varid")
               call netcdf_io_check( nf90_def_var(nc%id, nc%c_lm_varname, nc%out_type, [nc%m_dimid, nc%l_dimid, nc%sign_dimid], &
                                     nc%c_lm_varid), "netcdf_io_initialize_output nf90_def_var c_lm_varid" )
              
            else
               call netcdf_io_check( nf90_def_var(nc%id, nc%j2rp2_varname, nc%out_type, nc%time_dimid, nc%j2rp2_varid), &
                                    "netcdf_io_initialize_output nf90_def_var j2rp2_varid"  )
               call netcdf_io_check( nf90_def_var(nc%id, nc%j4rp4_varname, nc%out_type, nc%time_dimid, nc%j4rp4_varid), &
                                    "netcdf_io_initialize_output nf90_def_var j4rp4_varid"  )
            end if
                        
         end if

         if (param%lyarkovsky) then
            call netcdf_io_check( nf90_def_var(nc%id, nc%albedo_varname, nc%out_type, nc%name_dimid, nc%albedo_varid), &
                                  "netcdf_io_initialize_output nf90_def_var albedo_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%emissivity_varname, nc%out_type, nc%name_dimid, nc%emissivity_varid), &
                                  "netcdf_io_initialize_output nf90_def_var emissivity_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%rot_k_varname, nc%out_type, nc%name_dimid, nc%rot_k_varid), &
                                  "netcdf_io_initialize_output nf90_def_var rot_k_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%gamma_varname, nc%out_type, nc%name_dimid, nc%gamma_varid), &
                                  "netcdf_io_initialize_output nf90_def_var gamma_varid"  )
         end if

         if (param%ltides) then
            call netcdf_io_check( nf90_def_var(nc%id, nc%k2_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%k2_varid), &
                                "netcdf_io_initialize_output nf90_def_var k2_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%q_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Q_varid), &
                                "netcdf_io_initialize_output nf90_def_var Q_varid"  )
         end if

         if (param%lenergy) then
            call netcdf_io_check( nf90_def_var(nc%id, nc%ke_orbit_varname, nc%out_type, nc%time_dimid, nc%KE_orb_varid), &
                                  "netcdf_io_initialize_output nf90_def_var KE_orb_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%ke_rot_varname, nc%out_type, nc%time_dimid, nc%KE_rot_varid), &
                                  "netcdf_io_initialize_output nf90_def_var KE_rot_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type, nc%time_dimid, nc%PE_varid), &
                                  "netcdf_io_initialize_output nf90_def_var PE_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%be_varname, nc%out_type, nc%time_dimid, nc%BE_varid), &
                                  "netcdf_io_initialize_output nf90_def_var BE_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%te_varname, nc%out_type, nc%time_dimid, nc%TE_varid), &
                                  "netcdf_io_initialize_output nf90_def_var TE_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%L_orbit_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], &
                                               nc%L_orbit_varid), &
                                  "netcdf_io_initialize_output nf90_def_var L_orbit_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%L_rot_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], &
                                               nc%L_rot_varid), &
                                  "netcdf_io_initialize_output nf90_def_var L_rot_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%L_escape_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], &
                                               nc%L_escape_varid), &
                                  "netcdf_io_initialize_output nf90_def_var L_escape_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%E_collisions_varname, nc%out_type, nc%time_dimid, nc%E_collisions_varid), &
                                  "netcdf_io_initialize_output nf90_def_var E_collisions_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%E_untracked_varname, nc%out_type, nc%time_dimid, nc%E_untracked_varid), &
                                  "netcdf_io_initialize_output nf90_def_var E_untracked_varid"  )
            call netcdf_io_check( nf90_def_var(nc%id, nc%GMescape_varname, nc%out_type, nc%time_dimid, nc%GMescape_varid), &
                                  "netcdf_io_initialize_output nf90_def_var GMescape_varid"  )
         end if

         ! Set fill mode to NaN for all variables
         call netcdf_io_check( nf90_inquire(nc%id, nVariables=nvar), "netcdf_io_initialize_output nf90_inquire nVariables" )
         do varid = 1, nvar
            call netcdf_io_check( nf90_inquire_variable(nc%id, varid, xtype=vartype, ndims=ndims), &
                                  "netcdf_io_initialize_output nf90_inquire_variable"  )
            select case(vartype)
            case(NF90_INT)
               if (varid == nc%status_varid) then
                  ! Be sure the status variable fill value is the INACTIVE symbolic value
                  call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, INACTIVE), &
                                  "netcdf_io_netcdf_initialize_output nf90_def_var_fill status variable"  )
               else
                  call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, NF90_FILL_INT), &
                                  "netcdf_io_initialize_output nf90_def_var_fill NF90_INT"  )
               end if
            case(NF90_FLOAT)
               call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, sfill), &
                                  "netcdf_io_initialize_output nf90_def_var_fill NF90_FLOAT"  )
            case(NF90_DOUBLE)
               call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, dfill), &
                                  "netcdf_io_initialize_output nf90_def_var_fill NF90_DOUBLE"  )
            case(NF90_CHAR)
               call netcdf_io_check( nf90_def_var_fill(nc%id, varid, NO_FILL, 0), &
                                  "netcdf_io_initialize_output nf90_def_var_fill NF90_CHAR"  )
            end select
         end do

         ! Set special fill mode for discard time so that we can make use of it for non-discarded bodies.
         if (param%lclose) then
            select case (vartype)
            case(NF90_FLOAT)
               call netcdf_io_check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_SP)), &
                                  "netcdf_io_initialize_output nf90_def_var_fill discard_time NF90_FLOAT"  )
            case(NF90_DOUBLE)
               call netcdf_io_check( nf90_def_var_fill(nc%id, nc%discard_time_varid, NO_FILL, huge(1.0_DP)), &
                                  "netcdf_io_initialize_output nf90_def_var_fill discard_time NF90_DOUBLE"  )
            end select
         end if

         ! Take the file out of define mode
         call netcdf_io_check( nf90_enddef(nc%id), "netcdf_io_initialize_output nf90_enddef"  )

         ! Add in the space dimension coordinates
         call netcdf_io_check( nf90_put_var(nc%id, nc%space_varid, nc%space_coords, start=[1], count=[NDIM]), &
                                  "netcdf_io_initialize_output nf90_put_var space"  )

         if (param%lrotation .and. nc%lc_lm_exists) then
            ! Populate coordinate values for l, m, and sign
            call netcdf_io_check( nf90_put_var(nc%id, nc%l_varid, [(i, i=0, nc%l_dim_max-1)], start=[1], count=[nc%l_dim_max]), &
                                 "netcdf_io_netcdf_initialize_output nf90_put_var l_varid")
            call netcdf_io_check( nf90_put_var(nc%id, nc%m_varid, [(i, i=0, nc%m_dim_max-1)], start=[1], count=[nc%m_dim_max]), &
                                 "netcdf_io_netcdf_initialize_output nf90_put_var m_varid")
            call netcdf_io_check( nf90_put_var(nc%id, nc%sign_varid, nc%sign_coords, start=[1], count=[2] ), &
                                  "netcdf_io_netcdf_initialize_output nf90_put_var sign_varid")
         end if

      end associate
      return

      667 continue
      write(*,*) "Error creating NetCDF output file. " // trim(adjustl(errmsg))
      call base_util_exit(FAILURE,param%display_unit)
   end subroutine swiftest_io_netcdf_initialize_output


   module subroutine swiftest_io_netcdf_open(self, param, readonly)
      !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton
      !!
      !! Opens a NetCDF file and does the variable inquiries to activate variable ids
      implicit none
      ! Arguments
      class(swiftest_netcdf_parameters), intent(inout) :: self     !! Parameters used to identify a particular NetCDF dataset
      class(swiftest_parameters),        intent(inout) :: param    !! Current run configuration parameters
      logical, optional,                 intent(in)    :: readonly !! Logical flag indicating that this should be open read only
      ! Internals
      integer(I4B) :: mode, status
      character(len=STRMAX) :: errmsg

      mode = NF90_WRITE
      if (present(readonly)) then
         if (readonly) mode = NF90_NOWRITE
      end if

      associate(nc => self)
         write(errmsg,*) "swiftest_io_netcdf_open nf90_open ",trim(adjustl(nc%file_name))
         call netcdf_io_check( nf90_open(nc%file_name, mode, nc%id), errmsg)
         self%lfile_is_open = .true.

         ! Dimensions
         call netcdf_io_check( nf90_inq_dimid(nc%id, nc%time_dimname, nc%time_dimid), &
                                  "swiftest_io_netcdf_open nf90_inq_dimid time_dimid"  )
         call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%time_dimid, nc%time_dimname, len=nc%max_tslot), &
                                  "swiftest_io_netcdf_open nf90_inquire_dimension max_tslot"  )
         call netcdf_io_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), &
                                  "swiftest_io_netcdf_open nf90_inq_dimid space_dimid"  )
         call netcdf_io_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), &
                                  "swiftest_io_netcdf_open nf90_inq_dimid name_dimid"  )
         call netcdf_io_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), &
                                  "swiftest_io_netcdf_open nf90_inq_dimid str_dimid"  )

         status = nf90_inq_varid(nc%id, nc%c_lm_varname, nc%c_lm_varid)
         if (status == NF90_NOERR) then
            call netcdf_io_check( nf90_inq_dimid(nc%id, nc%sign_dimname, nc%sign_dimid), &
                                    "swiftest_io_netcdf_open nf90_inq_dimid sign_dimid")
            call netcdf_io_check( nf90_inq_dimid(nc%id, nc%l_dimname, nc%l_dimid), &
                                    "swiftest_io_netcdf_open nf90_inq_dimid l_dimid")
            call netcdf_io_check( nf90_inq_dimid(nc%id, nc%m_dimname, nc%m_dimid), &
                                    "swiftest_io_netcdf_open nf90_inq_dimid m_dimid")
         end if

         ! Dimension coordinates
         call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid time_varid" )
         call netcdf_io_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid space_varid" )
         call netcdf_io_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid name_varid" )
         status = nf90_inq_varid(nc%id, nc%c_lm_varname, nc%c_lm_varid)
         if (status == NF90_NOERR) then
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%sign_dimname, nc%sign_varid), &
                                    "swiftest_io_netcdf_open nf90_inq_varid sign_varid")
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%l_dimname, nc%l_varid), &
                                    "swiftest_io_netcdf_open nf90_inq_varid l_varid")
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%m_dimname, nc%m_varid), &
                                    "swiftest_io_netcdf_open nf90_inq_varid m_varid")
         end if

         ! Required Variables
         call netcdf_io_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid name_varid" )
         call netcdf_io_check( nf90_inq_varid(nc%id, nc%Gmass_varname, nc%Gmass_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid Gmass_varid" )

         if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid rh_varid" )
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid vh_varid" )

            if (param%lgr) then
               !! check if pseudovelocity vectors exist in this file. If they are, set the correct flag so we know whe should not 
               !! do the conversion.
               status = nf90_inq_varid(nc%id, nc%gr_pseudo_vh_varname, nc%gr_pseudo_vh_varid)
               nc%lpseudo_vel_exists = (status == NF90_NOERR)
               if (param%lrestart .and. .not.nc%lpseudo_vel_exists) then
                  write(*,*) "Warning! Pseudovelocity not found in input file for GR enabled run. " &
                          // "If this is a restarted run, bit-identical trajectories are not guarunteed!"
               end if

            end if
         end if

         if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%a_varname, nc%a_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid a_varid" )
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%e_varname, nc%e_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid e_varid" )
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%inc_varname, nc%inc_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid inc_varid" )
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%capom_varname, nc%capom_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid capom_varid" )
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%omega_varname, nc%omega_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid omega_varid" )
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%capm_varname, nc%capm_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid capm_varid" )
         end if

         if (param%lclose) then
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid radius_varid" )
         end if 

         if (param%lrotation) then
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid Ip_varid" )
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid rot_varid" )

            ! rotphase may not be input by the user                      
            status = nf90_inq_varid(nc%id, nc%rotphase_varname, nc%rotphase_varid)

         end if

         if (param%lyarkovsky) then
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%albedo_varname, nc%albedo_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid albedo_varid" )
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%emissivity_varname, nc%emissivity_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid emissivity_varid" )
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_k_varname, nc%rot_k_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid rot_k_varid" )
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%gamma_varname, nc%gamma_varid), &
                                  "swiftest_io_netcdf_open nf90_inq_varid gamma_varid" )
         end if

         if (param%ltides) then
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%k2_varname, nc%k2_varid), &
                                 "swiftest_io_netcdf_open nf90_inq_varid k2_varid" )
            call netcdf_io_check( nf90_inq_varid(nc%id, nc%q_varname, nc%Q_varid), &
                                 "swiftest_io_netcdf_open nf90_inq_varid Q_varid" )
         end if

         ! Optional variables The User Doesn't Need to Know About
         status = nf90_inq_varid(nc%id, nc%mass_varname, nc%mass_varid)
         status = nf90_inq_varid(nc%id, nc%rhill_varname, nc%rhill_varid)
         status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid)
         status = nf90_inq_varid(nc%id, nc%status_varname, nc%status_varid)
         status = nf90_inq_varid(nc%id, nc%ntp_varname, nc%ntp_varid)
         status = nf90_inq_varid(nc%id, nc%j2rp2_varname, nc%j2rp2_varid)
         status = nf90_inq_varid(nc%id, nc%j4rp4_varname, nc%j4rp4_varid)
         status = nf90_inq_varid(nc%id, nc%ptype_varname, nc%ptype_varid)
         status = nf90_inq_varid(nc%id, nc%varpi_varname, nc%varpi_varid)
         status = nf90_inq_varid(nc%id, nc%lam_varname, nc%lam_varid)
         status = nf90_inq_varid(nc%id, nc%f_varname, nc%f_varid)
         status = nf90_inq_varid(nc%id, nc%cape_varname, nc%cape_varid)

         if (param%lmtiny_pl) status = nf90_inq_varid(nc%id, nc%nplm_varname, nc%nplm_varid)

         if (param%lclose) then
            status = nf90_inq_varid(nc%id, nc%origin_type_varname, nc%origin_type_varid)
            status = nf90_inq_varid(nc%id, nc%origin_time_varname, nc%origin_time_varid)
            status = nf90_inq_varid(nc%id, nc%origin_rh_varname, nc%origin_rh_varid)
            status = nf90_inq_varid(nc%id, nc%origin_vh_varname, nc%origin_vh_varid)
            status = nf90_inq_varid(nc%id, nc%collision_id_dimname, nc%collision_id_varid)
            status = nf90_inq_varid(nc%id, nc%discard_time_varname, nc%discard_time_varid)
            status = nf90_inq_varid(nc%id, nc%discard_rh_varname, nc%discard_rh_varid)
            status = nf90_inq_varid(nc%id, nc%discard_vh_varname, nc%discard_vh_varid)
            status = nf90_inq_varid(nc%id, nc%discard_body_id_varname, nc%discard_body_id_varid)
         end if

         if (param%lenergy) then
            status = nf90_inq_varid(nc%id, nc%ke_orbit_varname, nc%KE_orb_varid)
            status = nf90_inq_varid(nc%id, nc%ke_rot_varname, nc%KE_rot_varid)
            status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid)
            status = nf90_inq_varid(nc%id, nc%be_varname, nc%BE_varid)
            status = nf90_inq_varid(nc%id, nc%te_varname, nc%TE_varid)
            status = nf90_inq_varid(nc%id, nc%L_orbit_varname, nc%L_orbit_varid)
            status = nf90_inq_varid(nc%id, nc%L_rot_varname, nc%L_rot_varid)
            status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid)
            status = nf90_inq_varid(nc%id, nc%E_collisions_varname, nc%E_collisions_varid)
            status = nf90_inq_varid(nc%id, nc%E_untracked_varname, nc%E_untracked_varid)
            status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid)
         end if

         status = nf90_inq_varid(nc%id, nc%c_lm_varname, nc%c_lm_varid)

      end associate

      return
   end subroutine swiftest_io_netcdf_open


   module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask, plmmask, Gmtiny)
      !! author: David A. Minton
      !!
      !! Given an open NetCDF, returns logical masks indicating which bodies in the body arrays are active pl and tp type at the 
      !! current time. Uses the value of tslot stored in the NetCDF parameter object as the definition of current time
      use, intrinsic :: ieee_exceptions
      use, intrinsic :: ieee_arithmetic
      implicit none
      ! Arguments
      class(swiftest_netcdf_parameters),  intent(inout)          :: self    !! Parameters used to identify a particular NetCDF 
                                                                            !!   dataset
      logical, dimension(:), allocatable, intent(out)            :: plmask  !! Logical mask indicating which bodies are massive 
                                                                            !!   bodies
      logical, dimension(:), allocatable, intent(out)            :: tpmask  !! Logical mask indicating which bodies are test 
                                                                            !!   particles
      logical, dimension(:), allocatable, intent(out), optional  :: plmmask !! Logical mask indicating which bodies are fully 
                                                                            !!   interacting massive bodies
      real(DP),                           intent(in),  optional  :: Gmtiny  !! The cutoff G*mass between semi-interacting and fully
                                                                            !!   interacting massive bodies
      ! Internals
      real(DP), dimension(:), allocatable :: Gmass, a
      real(DP), dimension(:,:), allocatable :: rh
      integer(I4B), dimension(:), allocatable :: body_status 
      logical, dimension(:), allocatable :: lvalid
      integer(I4B) :: idmax, status
      logical, dimension(size(IEEE_ALL))      :: fpe_halting_modes
      character(len=NAMELEN), dimension(:), allocatable :: ctemp

      call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes)  ! Save the current halting modes so we can turn them off temporarily
      call ieee_set_halting_mode(IEEE_ALL,.false.)

      call netcdf_io_check( nf90_inquire_dimension(self%id, self%name_dimid, len=idmax), &
                                  "swiftest_io_netcdf_get_valid_masks nf90_inquire_dimension name_dimid"  )

      allocate(Gmass(idmax))
      allocate(tpmask(idmax))
      allocate(plmask(idmax))
      allocate(lvalid(idmax))
      tpmask(:) = .false.
      plmask(:) = .false.
      lvalid(:) = .false.
      associate(tslot => self%tslot)

         call netcdf_io_check( nf90_get_var(self%id, self%Gmass_varid, Gmass, start=[1,tslot], count=[idmax,1]), &
                                  "swiftest_io_netcdf_get_valid_masks nf90_getvar Gmass_varid"  )

         status = nf90_inq_varid(self%id, self%status_varname, self%status_varid) 
         if (status == NF90_NOERR) then
            allocate(body_status(idmax))
            call netcdf_io_check( nf90_get_var(self%id, self%status_varid, body_status, start=[1, tslot], count=[idmax,1]), &
                                  "swiftest_io_netcdf_get_valid_masks nf90_getvar status_varid"  )
            lvalid(:) = body_status(:) /= INACTIVE
         else
            status = nf90_inq_varid(self%id, self%rh_varname, self%rh_varid) 
            if (status == NF90_NOERR) then
               allocate(rh(NDIM,idmax))
               call netcdf_io_check( nf90_get_var(self%id, self%rh_varid, rh, start=[1, 1, tslot], count=[NDIM,idmax,1]), &
                                  "swiftest_io_netcdf_get_valid_masks nf90_getvar rh_varid"  )
               lvalid(:) = ieee_is_normal(rh(1,:)) 
            else
               status = nf90_inq_varid(self%id, self%a_varname, self%a_varid) 
               if (status == NF90_NOERR) then
                  allocate(a(idmax))
                  call netcdf_io_check( nf90_get_var(self%id, self%a_varid, a, start=[1, tslot], count=[idmax,1]), &
                                  "swiftest_io_netcdf_get_valid_masks nf90_getvar a_varid"  )
                  lvalid(:) = ieee_is_normal(a(:))
               else
                  lvalid(:) = .false.
               end if
            end if
         end if

         allocate(ctemp(idmax))
         status = nf90_get_var(self%id, self%ptype_varid, ctemp, count=[NAMELEN, idmax])
         if (status == NF90_NOERR) then
            where(ctemp(:) == PL_TYPE_NAME) plmask(:) = .true.
            where(ctemp(:) == PL_TINY_TYPE_NAME) plmask(:) = .true.
            where(ctemp(:) == TP_TYPE_NAME) tpmask(:) = .true.
         else
            plmask(:) = ieee_is_normal(Gmass(:))
            where(plmask(:)) plmask(:) = Gmass(:) > 0.0_DP
            tpmask(:) = .not. plmask(:)
            plmask(1) = .false. ! This is the central body
         end if 
         ! Select only active bodies
         plmask(:) = plmask(:) .and. lvalid(:)
         tpmask(:) = tpmask(:) .and. lvalid(:)

         if (present(plmmask) .and. present(Gmtiny)) then
            allocate(plmmask, source=plmask)
            if (status == NF90_NOERR) then
               where(ctemp(:) == PL_TINY_TYPE_NAME) plmmask(:) = .FALSE.
            else
               where(plmask(:)) plmmask = Gmass(:) > Gmtiny
            end if
         end if
         deallocate(ctemp)

         call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes)

      end associate

      return
   end subroutine swiftest_io_netcdf_get_valid_masks


   module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ierr)
      !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
      !!
      !! Read a frame (header plus records for each massive body and active test particle) from an output binary file
      implicit none
      ! Arguments
      class(swiftest_nbody_system),      intent(inout) :: self  !! Swiftest nbody_system object
      class(swiftest_netcdf_parameters), intent(inout) :: nc    !! Parameters used to identify a particular NetCDF dataset
      class(swiftest_parameters),        intent(inout) :: param !! Current run configuration parameters 
      ! Return
      integer(I4B)                                :: ierr  !! Error code: returns 0 if the read is successful
      ! Internals
      integer(I4B)                              :: i, idmax, npl_check, ntp_check, str_max, status, npl, ntp
      real(DP), dimension(:), allocatable       :: rtemp
      real(DP), dimension(:,:), allocatable     :: vectemp
      integer(I4B), dimension(:), allocatable   :: itemp
      logical, dimension(:), allocatable        :: tpmask, plmask

      associate(cb => self%cb, pl => self%pl, tp => self%tp, tslot => nc%tslot)
         call nc%open(param, readonly=.true.)
         call nc%find_tslot(self%t, tslot)
         call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%time_dimid, len=nc%max_tslot), &
                                  "netcdf_io_read_frame_system nf90_inquire_dimension time_dimid"  )
         tslot = min(tslot, nc%max_tslot)

         call self%read_hdr(nc, param)

         ! Save these values as variables as they get reset by the setup method
         npl = pl%nbody
         ntp = tp%nbody
         call pl%setup(npl, param)
         call tp%setup(ntp, param)
         call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%name_dimid, len=idmax), &
                                  "netcdf_io_read_frame_system nf90_inquire_dimension name_dimid"  )
         allocate(rtemp(idmax))
         allocate(vectemp(NDIM,idmax))
         allocate(itemp(idmax))

         call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%str_dimid, len=str_max), &
                                  "netcdf_io_read_frame_system nf90_inquire_dimension str_dimid"  )

         call nc%get_valid_masks(plmask, tpmask)

         ! Check to make sure the number of bodies is correct
         npl_check = count(plmask(:))
         ntp_check = count(tpmask(:))

         if (npl_check /= npl) then
            write(*,*) "Error reading in NetCDF file: The recorded value of npl does not match the number of active massive bodies"
            write(*,*) "Recorded: ",npl
            write(*,*) "Active  : ",npl_check
         end if

         if (ntp_check /= ntp) then
            write(*,*) "Error reading in NetCDF file: The recorded value of ntp does not match the number of active test particles"
            write(*,*) "Recorded: ",ntp
            write(*,*) "Active  : ",ntp_check
            call base_util_exit(FAILURE,param%display_unit)
         end if

         ! Now read in each variable and split the outputs by body type
         if ((param%in_form == "XV") .or. (param%in_form == "XVEL")) then
            call netcdf_io_check( nf90_get_var(nc%id, nc%rh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar rh_varid"  )
            do i = 1, NDIM
               if (npl > 0) pl%rh(i,:) = pack(vectemp(i,:), plmask(:))
               if (ntp > 0) tp%rh(i,:) = pack(vectemp(i,:), tpmask(:))
            end do

            if (param%lgr .and. nc%lpseudo_vel_exists) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%gr_pseudo_vh_varid, vectemp, start=[1, 1, tslot], &
                                                  count=[NDIM,idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar gr_pseudo_vh_varid"  )
               do i = 1, NDIM
                  if (npl > 0) pl%vh(i,:) = pack(vectemp(i,:), plmask(:))
                  if (ntp > 0) tp%vh(i,:) = pack(vectemp(i,:), tpmask(:))
               end do
            else
               call netcdf_io_check( nf90_get_var(nc%id, nc%vh_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar vh_varid"  )
               do i = 1, NDIM
                  if (npl > 0) pl%vh(i,:) = pack(vectemp(i,:), plmask(:))
                  if (ntp > 0) tp%vh(i,:) = pack(vectemp(i,:), tpmask(:))
               end do
            end if
         end if

         if ((param%in_form == "EL")  .or. (param%in_form == "XVEL")) then
            call netcdf_io_check( nf90_get_var(nc%id, nc%a_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar a_varid"  )
            if (.not.allocated(pl%a)) allocate(pl%a(npl))
            if (.not.allocated(tp%a)) allocate(tp%a(ntp))
            if (npl > 0) pl%a(:) = pack(rtemp, plmask)
            if (ntp > 0) tp%a(:) = pack(rtemp, tpmask)

            call netcdf_io_check( nf90_get_var(nc%id, nc%e_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar e_varid"  )
            if (.not.allocated(pl%e)) allocate(pl%e(npl))
            if (.not.allocated(tp%e)) allocate(tp%e(ntp))
            if (npl > 0) pl%e(:) = pack(rtemp, plmask)
            if (ntp > 0) tp%e(:) = pack(rtemp, tpmask)

            call netcdf_io_check( nf90_get_var(nc%id, nc%inc_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar inc_varid"  )
            if (.not.allocated(pl%inc)) allocate(pl%inc(npl))
            if (.not.allocated(tp%inc)) allocate(tp%inc(ntp))
            if (npl > 0) pl%inc(:) = pack(rtemp, plmask)
            if (ntp > 0) tp%inc(:) = pack(rtemp, tpmask)

            call netcdf_io_check( nf90_get_var(nc%id, nc%capom_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar capom_varid"  )
            if (.not.allocated(pl%capom)) allocate(pl%capom(npl))
            if (.not.allocated(tp%capom)) allocate(tp%capom(ntp))
            if (npl > 0) pl%capom(:) = pack(rtemp, plmask)
            if (ntp > 0) tp%capom(:) = pack(rtemp, tpmask)

            call netcdf_io_check( nf90_get_var(nc%id, nc%omega_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar omega_varid"  )
            if (.not.allocated(pl%omega)) allocate(pl%omega(npl))
            if (.not.allocated(tp%omega)) allocate(tp%omega(ntp))
            if (npl > 0) pl%omega(:) = pack(rtemp, plmask)
            if (ntp > 0) tp%omega(:) = pack(rtemp, tpmask)

            call netcdf_io_check( nf90_get_var(nc%id, nc%capm_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar capm_varid"  )
            if (.not.allocated(pl%capm)) allocate(pl%capm(npl))
            if (.not.allocated(tp%capm)) allocate(tp%capm(ntp))
            if (npl > 0) pl%capm(:) = pack(rtemp, plmask)
            if (ntp > 0) tp%capm(:) = pack(rtemp, tpmask)

         end if
      
         call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar Gmass_varid"  )

         cb%Gmass = rtemp(1)           

         ! Set initial central body mass for Helio bookkeeping
         cb%GM0 = cb%Gmass
                                  
         if (npl > 0) then
            pl%Gmass(:) = pack(rtemp, plmask)
            if (param%lmtiny_pl) pl%nplm = count(pack(rtemp,plmask) > param%GMTINY )

            status = nf90_get_var(nc%id, nc%rhill_varid, rtemp, start=[1, tslot], count=[idmax,1])
            if (status == NF90_NOERR) then
               pl%rhill(:) = pack(rtemp, plmask)
            end if
         end if

         status = nf90_inq_varid(nc%id, nc%mass_varname, nc%mass_varid)
         if (status == NF90_NOERR) then
            call netcdf_io_check( nf90_get_var(nc%id, nc%mass_varid, rtemp, start=[1,tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar mass_varid"  )
            cb%mass = rtemp(1)
            if (npl > 0) pl%mass(:) = pack(rtemp, plmask)
         else
            cb%mass = cb%Gmass / param%GU
            if (npl > 0) pl%mass(:) = pl%Gmass(:) / param%GU
         end if

         if (param%lclose) then
            call netcdf_io_check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar radius_varid"  )
            cb%radius = rtemp(1)

            if (npl > 0) pl%radius(:) = pack(rtemp, plmask)
         else
            cb%radius = param%rmin
            if (npl > 0) pl%radius(:) = 0.0_DP
         end if
         cb%R0 = cb%radius
         cb%density = cb%mass / (4.0_DP / 3.0_DP * PI * cb%radius**3)

         if (param%lrotation) then
            call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid,  vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar Ip_varid"  )
            cb%Ip(:) = vectemp(:,1)
            do i = 1, NDIM
               if (npl > 0) pl%Ip(i,:) = pack(vectemp(i,:), plmask(:))
            end do

            call netcdf_io_check( nf90_get_var(nc%id, nc%rot_varid, vectemp, start=[1, 1, tslot], count=[NDIM,idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar rot_varid"  )
            cb%rot(:) = vectemp(:,1) 
            do i = 1, NDIM
               if (npl > 0) pl%rot(i,:) = pack(vectemp(i,:), plmask(:))
            end do

            ! Set initial central body angular momentum for bookkeeping
            cb%L0(:) = cb%Ip(3) * cb%mass * cb%R0**2 * cb%rot(:) * DEG2RAD
           
            ! Check if rotphase is input by user. If not, set to 0 
            status = nf90_inq_varid(nc%id, nc%rotphase_varname, nc%rotphase_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%rotphase_varid, cb%rotphase, start=[tslot]), &
                                  "netcdf_io_read_frame_system nf90_getvar rotphase_varid"  )
            else
               cb%rotphase = 0.0_DP
            end if

            ! Check if the J2 and J4 terms are input by user. If not, set them to 0
            status = nf90_inq_varid(nc%id, nc%j2rp2_varname, nc%j2rp2_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%j2rp2_varid, cb%j2rp2, start=[tslot]), &
                                     "netcdf_io_read_frame_system nf90_getvar j2rp2_varid"  )
            else 
               cb%j2rp2 = 0.0_DP
            end if

            status = nf90_inq_varid(nc%id, nc%j4rp4_varname, nc%j4rp4_varid)   
            if (status == NF90_NOERR) then      
               call netcdf_io_check( nf90_get_var(nc%id, nc%j4rp4_varid, cb%j4rp4, start=[tslot]), &
                                     "netcdf_io_read_frame_system nf90_getvar j4rp4_varid"  )
            else 
               cb%j4rp4 = 0.0_DP
            end if
  
            ! Check if gravitational harmonics coefficients, c_lm, are set by the user. Also ensure that
            ! c_lm and j2rp2/j4rp4 are not both set.
            status = nf90_inq_varid(nc%id, nc%c_lm_varname, nc%c_lm_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%l_dimid, len = nc%l_dim_max),&
                                     "netcdf_io_read_frame_system nf90_inquire_dimension l_dimid"  )
               call netcdf_io_check( nf90_inquire_dimension(nc%id, nc%m_dimid, len = nc%m_dim_max), &
                                      "netcdf_io_read_frame_system nf90_inquire_dimension m_dimid")
              
               if (nc%l_dim_max /= nc%m_dim_max) then
                  write(*,*) "Error reading in NetCDF file: c_lm requires l_dim_max = m_dim_max"
                  call base_util_exit(FAILURE,param%display_unit) 
               end if
   
               if(.not. allocated(cb%c_lm)) allocate(cb%c_lm(nc%m_dim_max, nc%l_dim_max, 2)) 
               call netcdf_io_check( nf90_get_var(nc%id, nc%c_lm_varid, cb%c_lm, count = [nc%m_dim_max, nc%l_dim_max, 2]), &
                                     "netcdf_io_read_frame_system nf90_getvar c_lm_varid")
              
               if ((cb%j2rp2 == cb%j2rp2) .or. (cb%j4rp4 == cb%j4rp4)) then
                  if (abs(cb%j2rp2) > tiny(1.0_DP) .or. (abs(cb%j4rp4) > tiny(1.0_DP))) then
                     write(*,*) "Error reading in NetCDF file: cannot use both c_lm and j2rp2/j4rp4"
                     call base_util_exit(FAILURE,param%display_unit) 
                  end if
               endif
   
               nc%lc_lm_exists = .true.
            else
               if (allocated(cb%c_lm)) deallocate(cb%c_lm)
               nc%lc_lm_exists = .false.
            end if

         end if

         if (param%lyarkovsky) then
            call netcdf_io_check( nf90_get_var(nc%id, nc%albedo_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar albedo_varid"  )
            if (.not.allocated(pl%albedo)) allocate(pl%albedo(npl))
            if (npl > 0) pl%albedo(:) = pack(rtemp, plmask)

            call netcdf_io_check( nf90_get_var(nc%id, nc%emissivity_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar emissivity_varid"  )
            if (.not.allocated(pl%emissivity)) allocate(pl%emissivity(npl))
            if (npl > 0) pl%emissivity(:) = pack(rtemp, plmask)

            call netcdf_io_check( nf90_get_var(nc%id, nc%rot_k_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar rot_k_varid"  )
            if (.not.allocated(pl%rot_k)) allocate(pl%rot_k(npl))
            if (npl > 0) pl%rot_k(:) = pack(rtemp, plmask)

            call netcdf_io_check( nf90_get_var(nc%id, nc%gamma_varid, rtemp, start=[1, tslot], count=[idmax,1]), &
                                  "netcdf_io_read_frame_system nf90_getvar gamma_varid"  )
            if (.not.allocated(pl%gamma)) allocate(pl%gamma(npl))
            if (npl > 0) pl%gamma(:) = pack(rtemp, plmask)
         end if

         if (param%ltides) then
            call netcdf_io_check( nf90_get_var(nc%id, nc%k2_varid, rtemp, start=[1, tslot]), &
                                "netcdf_io_read_frame_system nf90_getvar k2_varid"  )
            cb%k2 = rtemp(1)
            if (npl > 0) pl%k2(:) = pack(rtemp, plmask)

            call netcdf_io_check( nf90_get_var(nc%id, nc%Q_varid,  rtemp,  start=[1, tslot]), &
                                "netcdf_io_read_frame_system nf90_getvar Q_varid"  )
            cb%Q = rtemp(1)
            if (npl > 0) pl%Q(:) = pack(rtemp, plmask)
         end if

         call self%read_particle_info(nc, param, plmask, tpmask) 

         if (param%in_form == "EL") then
            call pl%el2xv(cb)
            call tp%el2xv(cb)
         end if
         ! if this is a GR-enabled run, check to see if we got the pseudovelocities in. Otherwise, we'll need to generate them.
         if (param%lgr .and. .not.(nc%lpseudo_vel_exists)) then
            call pl%set_mu(cb)
            call tp%set_mu(cb)
            call pl%v2pv(param)
            call tp%v2pv(param)
         end if
         
      end associate

      call nc%close()

      ierr = 0
      return
   end function swiftest_io_netcdf_read_frame_system


   module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param) 
      !! author: David A. Minton
      !!
      !! Reads header information (variables that change with time, but not particle id). 
      !! This subroutine swiftest_significantly improves the output over the original binary file, allowing us to track energy, 
      !! momentum, and other quantities that previously were handled as separate output files.
      implicit none
      ! Arguments
      class(swiftest_nbody_system),      intent(inout) :: self  !! Swiftest nbody system object
      class(swiftest_netcdf_parameters), intent(inout) :: nc    !! Parameters used to for reading a NetCDF dataset to file
      class(swiftest_parameters),        intent(inout) :: param !! Current run configuration parameters
      ! Internals
      integer(I4B) :: status
      logical, dimension(:), allocatable        :: plmask, tpmask, plmmask

      associate(tslot => nc%tslot)
         call netcdf_io_check( nf90_get_var(nc%id, nc%time_varid, self%t, start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar time_varid"  )
         if (param%lmtiny_pl) then
            call nc%get_valid_masks(plmask, tpmask, plmmask, param%GMTINY)
         else
            call nc%get_valid_masks(plmask, tpmask)
         end if

         status = nf90_inq_varid(nc%id, nc%npl_varname, nc%npl_varid)
         if (status == NF90_NOERR) then
            call netcdf_io_check( nf90_get_var(nc%id, nc%npl_varid,  self%pl%nbody, start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar npl_varid"  )
         else
            self%pl%nbody = count(plmask(:))
         end if

         status = nf90_inq_varid(nc%id, nc%ntp_varname, nc%ntp_varid)
         if (status == NF90_NOERR) then
            call netcdf_io_check( nf90_get_var(nc%id, nc%ntp_varid,  self%tp%nbody, start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar ntp_varid"  )
         else
            self%tp%nbody = count(tpmask(:))
         end if

         if (param%lmtiny_pl) then
            status = nf90_inq_varid(nc%id, nc%nplm_varname, nc%nplm_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%nplm_varid,  self%pl%nplm, start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar nplm_varid"  )
            else
               self%pl%nplm = count(plmmask(:))
            end if
         end if

         if (param%lenergy) then
            status = nf90_inq_varid(nc%id, nc%ke_orbit_varname, nc%KE_orb_varid)
            if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar KE_orb_varid"  )
            status = nf90_inq_varid(nc%id, nc%ke_rot_varname, nc%KE_rot_varid)
            if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%KE_rot_varid, self%ke_rot, start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar KE_rot_varid"  )
            status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid)
            if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar PE_varid"  )
            status = nf90_inq_varid(nc%id, nc%be_varname, nc%BE_varid)
            if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, self%be, start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar BE_varid"  )
            status = nf90_inq_varid(nc%id, nc%te_varname, nc%TE_varid)
            if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%TE_varid, self%te, start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar TE_varid"  )
            status = nf90_inq_varid(nc%id, nc%L_orbit_varname, nc%L_orbit_varid)
            if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%L_orbit_varid, self%L_orbit(:), & 
                                                                         start=[1,tslot], count=[NDIM,1]), &
                                  "netcdf_io_read_hdr_system nf90_getvar L_orbit_varid"  )
            status = nf90_inq_varid(nc%id, nc%L_rot_varname, nc%L_rot_varid)
            if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%L_rot_varid, self%L_rot(:), start=[1,tslot], &
                                                            count=[NDIM,1]), &
                                  "netcdf_io_read_hdr_system nf90_getvar L_rot_varid"  )
            status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid)
            if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%L_escape(:), &
                                                                         start=[1, tslot], count=[NDIM,1]), &
                                  "netcdf_io_read_hdr_system nf90_getvar L_escape_varid"  )
            status = nf90_inq_varid(nc%id, nc%E_collisions_varname, nc%E_collisions_varid)
            if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%E_collisions_varid, self%E_collisions, &
                                                                         start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar E_collisions_varid"  )
            status = nf90_inq_varid(nc%id, nc%E_untracked_varname, nc%E_untracked_varid)
            if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%E_untracked_varid, self%E_untracked, &
                                                                          start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar E_untracked_varid"  )
            status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid)
            if (status == NF90_NOERR)  call netcdf_io_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), &
                                  "netcdf_io_read_hdr_system nf90_getvar GMescape_varid"  )
         end if

      end associate

      return
   end subroutine swiftest_io_netcdf_read_hdr_system


   module subroutine swiftest_io_netcdf_read_particle_info_system(self, nc, param, plmask, tpmask)
      !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton
      !!
      !! Reads particle information metadata from file
      implicit none
      ! Arguments
      class(swiftest_nbody_system),      intent(inout) :: self   !! Swiftest nbody system object
      class(swiftest_netcdf_parameters), intent(inout) :: nc     !! Parameters used to identify a particular NetCDF dataset
      class(swiftest_parameters),        intent(inout) :: param  !! Current run configuration parameters
      logical, dimension(:),             intent(in)    :: plmask !! Logical array indicating which index values belong to massive 
                                                                 !!     bodies
      logical, dimension(:),             intent(in)    :: tpmask !! Logical array indicating which index values belong to test 
                                                                 !!     particles

      ! Internals
      integer(I4B)                                :: i, idmax, status
      real(DP), dimension(:), allocatable         :: rtemp
      real(DP), dimension(:,:), allocatable       :: vectemp
      integer(I4B), dimension(:), allocatable     :: itemp
      character(len=NAMELEN), dimension(:), allocatable :: ctemp
      integer(I4B), dimension(:), allocatable     :: plind, tpind

      ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables
      idmax = size(plmask)
      allocate(rtemp(idmax))
      allocate(vectemp(NDIM,idmax))
      allocate(itemp(idmax))
      allocate(ctemp(idmax))

      associate(cb => self%cb, pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody)

         if (npl > 0) then
            pl%status(:) = ACTIVE
            pl%lmask(:) = .true.
            do i = 1, npl
               call pl%info(i)%set_value(status="ACTIVE")
            end do
            allocate(plind(npl))
            plind(:) = pack([(i, i = 1, idmax)], plmask(:))
         end if
         if (ntp > 0) then
            tp%status(:) = ACTIVE
            tp%lmask(:) = .true.
            do i = 1, ntp
               call tp%info(i)%set_value(status="ACTIVE")
            end do
            allocate(tpind(ntp))
            tpind(:) = pack([(i, i = 1, idmax)], tpmask(:))
         end if

         call netcdf_io_check( nf90_get_var(nc%id, nc%id_varid, itemp), "netcdf_io_read_particle_info_system nf90_getvar id_varid")
         cb%id = itemp(1)
         pl%id(:) = pack(itemp, plmask)
         tp%id(:) = pack(itemp, tpmask)
         cb%id = 0
         pl%id(:) = pack([(i,i=0,idmax-1)],plmask)
         tp%id(:) = pack([(i,i=0,idmax-1)],tpmask)

         call netcdf_io_check( nf90_get_var(nc%id, nc%name_varid, ctemp, count=[NAMELEN, idmax]), &
                                  "netcdf_io_read_particle_info_system nf90_getvar name_varid"  )
         call cb%info%set_value(name=ctemp(1))
         do i = 1, npl
            call pl%info(i)%set_value(name=ctemp(plind(i)))
         end do
         do i = 1, ntp
            call tp%info(i)%set_value(name=ctemp(tpind(i)))
         end do

         status = nf90_get_var(nc%id, nc%ptype_varid, ctemp, count=[NAMELEN, idmax])
         if (status /= NF90_NOERR) then ! Set default particle types
            call cb%info%set_value(particle_type=CB_TYPE_NAME)

            ! Handle semi-interacting bodies in SyMBA
            do i = 1, npl
               if (param%lmtiny_pl .and. (pl%Gmass(i) < param%GMTINY)) then
                  call pl%info(i)%set_value(particle_type=PL_TINY_TYPE_NAME)
               else
                  call pl%info(i)%set_value(particle_type=PL_TYPE_NAME)
               end if
            end do
            do i = 1, ntp
               call tp%info(i)%set_value(particle_type=TP_TYPE_NAME)
            end do
         else ! Use particle types defined in input file
            call cb%info%set_value(particle_type=ctemp(1))
            do i = 1, npl
               call pl%info(i)%set_value(particle_type=ctemp(plind(i)))
            end do
            do i = 1, ntp
               call tp%info(i)%set_value(particle_type=ctemp(tpind(i)))
            end do
         end if

         call cb%info%set_value(status="ACTIVE")

         if (param%lclose) then
            status = nf90_inq_varid(nc%id, nc%origin_type_varname, nc%origin_type_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%origin_type_varid, ctemp, count=[NAMELEN, idmax]), &
                                  "netcdf_io_read_particle_info_system nf90_getvar origin_type_varid"  )
            else
               ctemp = "Initial Conditions"
            end if

            call cb%info%set_value(origin_type=ctemp(1))
            do i = 1, npl
               call pl%info(i)%set_value(origin_type=ctemp(plind(i)))
            end do
            do i = 1, ntp
               call tp%info(i)%set_value(origin_type=ctemp(tpind(i)))
            end do

            status = nf90_inq_varid(nc%id, nc%origin_time_varname, nc%origin_time_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%origin_time_varid, rtemp), &
                                  "netcdf_io_read_particle_info_system nf90_getvar origin_time_varid"  )
            else
               rtemp = param%t0
            end if

            call cb%info%set_value(origin_time=rtemp(1))
            do i = 1, npl
               call pl%info(i)%set_value(origin_time=rtemp(plind(i)))
            end do
            do i = 1, ntp
               call tp%info(i)%set_value(origin_time=rtemp(tpind(i)))
            end do

            status = nf90_inq_varid(nc%id, nc%origin_rh_varname, nc%origin_rh_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%origin_rh_varid, vectemp(:,:)), &
                                  "netcdf_io_read_particle_info_system nf90_getvar origin_rh_varid"  )
            else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%rh_varid, vectemp(:,:)), &
                                  "netcdf_io_read_particle_info_system nf90_getvar rh_varid"  )
            else 
               vectemp(:,:) = 0._DP
            end if 

            do i = 1, npl
               call pl%info(i)%set_value(origin_rh=vectemp(:,plind(i)))
            end do
            do i = 1, ntp
               call tp%info(i)%set_value(origin_rh=vectemp(:,tpind(i)))
            end do

            status = nf90_inq_varid(nc%id, nc%origin_vh_varname, nc%origin_vh_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%origin_vh_varid, vectemp(:,:)), &
                                  "netcdf_io_read_particle_info_system nf90_getvar origin_vh_varid"  )
            else if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%vh_varid, vectemp(:,:)), &
                                  "netcdf_io_read_particle_info_system nf90_getvar vh_varid"  )
            else
               vectemp(:,:) = 0._DP
            end if 
            
            do i = 1, npl
               call pl%info(i)%set_value(origin_vh=vectemp(:,plind(i)))
            end do
            do i = 1, ntp
               call tp%info(i)%set_value(origin_vh=vectemp(:,tpind(i)))
            end do

            status = nf90_inq_varid(nc%id, nc%collision_id_dimname, nc%collision_id_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%collision_id_varid, itemp), &
                                  "netcdf_io_read_particle_info_system nf90_getvar collision_id_varid"  )
               do i = 1, npl
                  call pl%info(i)%set_value(collision_id=itemp(plind(i)))
               end do
               do i = 1, ntp
                  call tp%info(i)%set_value(collision_id=itemp(tpind(i)))
               end do
            end if

            status = nf90_inq_varid(nc%id, nc%discard_time_varname, nc%discard_time_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%discard_time_varid, rtemp), &
                                  "netcdf_io_read_particle_info_system nf90_getvar discard_time_varid"  )
               call cb%info%set_value(discard_time=rtemp(1))
               do i = 1, npl
                  call pl%info(i)%set_value(discard_time=rtemp(plind(i)))
               end do
               do i = 1, ntp
                  call tp%info(i)%set_value(discard_time=rtemp(tpind(i)))
               end do
            end if

            status = nf90_inq_varid(nc%id, nc%discard_rh_varname, nc%discard_rh_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%discard_rh_varid, vectemp(:,:)), &
                                  "netcdf_io_read_particle_info_system nf90_getvar discard_rh_varid"  )
               do i = 1, npl
                  call pl%info(i)%set_value(discard_rh=vectemp(:,plind(i)))
               end do
               do i = 1, ntp
                  call tp%info(i)%set_value(discard_rh=vectemp(:,tpind(i)))
               end do
            end if

            status = nf90_inq_varid(nc%id, nc%discard_vh_varname, nc%discard_vh_varid)
            if (status == NF90_NOERR) then
               call netcdf_io_check( nf90_get_var(nc%id, nc%discard_vh_varid, vectemp(:,:)), &
                                  "netcdf_io_read_particle_info_system nf90_getvar discard_vh_varid"  )
               do i = 1, npl
                  call pl%info(i)%set_value(discard_vh=vectemp(:,plind(i)))
               end do
               do i = 1, ntp
                  call tp%info(i)%set_value(discard_vh=vectemp(:,tpind(i)))
               end do
            end if
         end if

      end associate

      return
   end subroutine swiftest_io_netcdf_read_particle_info_system


   module subroutine swiftest_io_netcdf_write_frame_body(self, nc, param)
      !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton
      !!
      !! Write a frame of output of either test particle or massive body data to the binary output file
      !!    Note: If outputting to orbital elements, make sure that the conversion is done prior to calling this method
      implicit none
      ! Arguments
      class(swiftest_body),              intent(in)    :: self  !! Swiftest base object
      class(swiftest_netcdf_parameters), intent(inout) :: nc    !! Parameters used to write a NetCDF dataset to file
      class(swiftest_parameters),        intent(inout) :: param !! Current run configuration parameters 
      ! Internals
      integer(I4B)                              :: i, j, idslot, old_mode, tmp
      integer(I4B), dimension(:), allocatable   :: ind
      real(DP), dimension(NDIM)                 :: vh !! Temporary variable to store heliocentric velocity values when converting 
                                                      !! from pseudovelocity in GR-enabled runs
      real(DP)                                  :: a, e, inc, omega, capom, capm, varpi, lam, f, cape, capf
#ifdef COARRAY
      integer(I4B) :: ntp
      logical, dimension(:), allocatable        :: tpmask, plmask
#endif

      call self%write_info(nc, param)

      call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), "netcdf_io_write_frame_body nf90_set_fill"  )
      select type(self)
      class is (swiftest_body)
      select type (param)
      class is (swiftest_parameters)
         associate(n => self%nbody, tslot => nc%tslot)
            if (n == 0) return

            call util_sort(self%id(1:n), ind)

            do i = 1, n
               j = ind(i)
               call nc%find_idslot(self%id(j), idslot) 
               ! Convert from pseudovelocity to heliocentric without replacing the current value of pseudovelocity 
               if (param%lgr) call swiftest_gr_pseudovel2vel(param, self%mu(j), self%rh(:, j), self%vh(:, j), vh(:))

               if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then
                  call netcdf_io_check( nf90_put_var(nc%id, nc%rh_varid, self%rh(:, j), start=[1,idslot, tslot], count=[NDIM,1,1]),&
                                  "netcdf_io_write_frame_body nf90_put_var rh_varid"  )
                  if (param%lgr) then ! Convert from pseudovelocity to heliocentric without replacing the current value of 
                                      !  pseudovelocity
                     call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, vh(:), start=[1,idslot, tslot], count=[NDIM,1,1]), &
                                  "netcdf_io_write_frame_body nf90_put_var vh_varid"  )
                     call netcdf_io_check( nf90_put_var(nc%id, nc%gr_pseudo_vh_varid, self%vh(:, j), start=[1,idslot, tslot], &
                                                        count=[NDIM,1,1]), &
                                  "netcdf_io_write_frame_body nf90_put_var gr_pseudo_vhx_varid"  )

                  else
                     call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, self%vh(:, j), start=[1,idslot, tslot], &
                                                        count=[NDIM,1,1]), &
                                  "netcdf_io_write_frame_body nf90_put_var vh_varid"  )
                  end if
               end if

               if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then
                  if (param%lgr) then ! For GR-enabled runs, use the true value of velocity computed above
                     call swiftest_orbel_xv2el(self%mu(j), self%rh(1,j), self%rh(2,j), self%rh(3,j), &
                                       vh(1), vh(2), vh(3), &
                                       a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf)
                  else !! For non-GR runs just convert from the velocity we have
                     call swiftest_orbel_xv2el(self%mu(j), self%rh(1,j), self%rh(2,j), self%rh(3,j), &
                                       self%vh(1,j), self%vh(2,j), self%vh(3,j), &
                                       a, e, inc, capom, omega, capm, varpi, lam, f, cape, capf)
                  end if
                  call netcdf_io_check( nf90_put_var(nc%id, nc%a_varid, a, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body a_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%e_varid, e, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body e_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%inc_varid, inc, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body inc_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%capom_varid, capom, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body capom_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%omega_varid, omega, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body omega_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%capm_varid, capm, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body capm_varid"  ) 
                  call netcdf_io_check( nf90_put_var(nc%id, nc%varpi_varid, varpi, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body varpi_varid"  ) 
                  call netcdf_io_check( nf90_put_var(nc%id, nc%lam_varid, lam, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body lam_varid"  ) 
                  call netcdf_io_check( nf90_put_var(nc%id, nc%f_varid, f, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body f_varid"  ) 
                  if (e < 1.0_DP) then
                     call netcdf_io_check( nf90_put_var(nc%id, nc%cape_varid, cape, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body cape_varid"  ) 
                  else if (e > 1.0_DP) then
                     call netcdf_io_check( nf90_put_var(nc%id, nc%cape_varid, capf , start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body (capf) cape_varid"  ) 
                  end if
               end if

               select type(self)  
               class is (swiftest_pl)  ! Additional output if the passed polymorphic object is a massive body
                  call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass(j), start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body Gmass_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%mass_varid, self%mass(j), start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body mass_varid"  )
                  if (param%lrhill_present) then
                     call netcdf_io_check( nf90_put_var(nc%id, nc%rhill_varid, self%rhill(j), start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body rhill_varid"  )
                  end if
                  if (param%lclose) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, self%radius(j), &
                                                                       start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body radius_varid"  )
                  if (param%lrotation) then
                     call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:, j), start=[1,idslot, tslot], &
                                                        count=[NDIM,1,1]), &
                                  "netcdf_io_write_frame_body nf90_put_var body Ip_varid"  )
                     call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:, j), start=[1,idslot, tslot], &
                                                        count=[NDIM,1,1]), &
                                  "netcdf_io_write_frame_body nf90_put_var body rotx_varid"  )
                  end if
                  
                  if (param%lyarkovsky) then
                     call netcdf_io_check( nf90_put_var(nc%id, nc%albedo_varid, self%albedo(j), start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body albedo_varid"  )
                     call netcdf_io_check( nf90_put_var(nc%id, nc%emissivity_varid, self%emissivity(j), start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body emissivity_varid"  )
                     call netcdf_io_check( nf90_put_var(nc%id, nc%rot_k_varid, self%rot_k(j), start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body rot_k_varid"  )
                     call netcdf_io_check( nf90_put_var(nc%id, nc%gamma_varid, self%gamma(j), start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body gamma_varid"  )
                  end if

                  if (param%ltides) then
                     call netcdf_io_check( nf90_put_var(nc%id, nc%k2_varid, self%k2(j), start=[idslot, tslot]), &
                                 "netcdf_io_write_frame_body nf90_put_var body k2_varid"  )
                     call netcdf_io_check( nf90_put_var(nc%id, nc%Q_varid, self%Q(j), start=[idslot, tslot]), &
                                 "netcdf_io_write_frame_body nf90_put_var body Q_varid"  )
                  end if
               class is (swiftest_tp)
                  call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, 0.0_DP, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body Gmass_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%mass_varid, 0.0_DP, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body mass_varid"  )
                  if (param%lrhill_present) then
                     call netcdf_io_check( nf90_put_var(nc%id, nc%rhill_varid, 0.0_DP, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body rhill_varid"  )
                  end if
                  if (param%lclose) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, 0.0_DP, start=[idslot, tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var body radius_varid"  )
                  if (param%lrotation) then
                     call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, [0.0_DP,0.0_DP,0.0_DP], start=[1,idslot, tslot], &
                                                        count=[NDIM,1,1]), &
                                  "netcdf_io_write_frame_body nf90_put_var body Ip_varid"  )
                     call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, [0.0_DP,0.0_DP,0.0_DP], start=[1,idslot, tslot], &
                                                        count=[NDIM,1,1]), &
                                  "netcdf_io_write_frame_body nf90_put_var body rotx_varid"  )
                  end if
               end select
            end do
         end associate
      end select
      end select
#ifdef COARRAY
      select type(self)
      class is (swiftest_tp)
         call nc%get_valid_masks(plmask, tpmask)
         ntp = count(tpmask(:))  
         call netcdf_io_check( nf90_put_var(nc%id, nc%ntp_varid, ntp, start=[nc%tslot]), &
                                  "netcdf_io_write_frame_body nf90_put_var ntp_varid"  )
      end select
#endif   

      call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp), &
                                  "netcdf_io_write_frame_body nf90_set_fill old_mode"  )

      return
   end subroutine swiftest_io_netcdf_write_frame_body


   module subroutine swiftest_io_netcdf_write_frame_cb(self, nc, param)
      !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton
      !!
      !! Write a frame of output of the central body
      implicit none
      ! Arguments
      class(swiftest_cb),                intent(inout)    :: self  !! Swiftest base object
      class(swiftest_netcdf_parameters), intent(inout) :: nc    !! Parameters used to for writing a NetCDF dataset to file
      class(swiftest_parameters),        intent(inout) :: param !! Current run configuration parameters 
      ! Internals
      integer(I4B)                              :: idslot, old_mode, tmp

      associate(tslot => nc%tslot)
         call self%write_info(nc, param)

         call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), &
                                  "swiftest_io_netcdf_write_frame_cb nf90_set_fill"  )

         call nc%find_idslot(self%id, idslot) 
         call netcdf_io_check( nf90_put_var(nc%id, nc%status_varid, ACTIVE, start=[idslot, tslot]), &
                                  "swiftest_io_netcdf_write_frame_cb nf90_put_var cb id_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, self%Gmass, start=[idslot, tslot]), &
                                  "swiftest_io_netcdf_write_frame_cb nf90_put_var cb Gmass_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%mass_varid, self%mass, start=[idslot, tslot]), &
                                  "swiftest_io_netcdf_write_frame_cb nf90_put_var cb mass_varid"  )
         if (param%lclose) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, self%radius, start=[idslot, tslot]), &
                                  "swiftest_io_netcdf_write_frame_cb nf90_put_var cb radius_varid"  )

         if (param%lrotation) then
            call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, self%Ip(:), start=[1, idslot, tslot], count=[NDIM,1,1]), &
                                  "swiftest_io_netcdf_write_frame_cb nf90_put_var cb Ip_varid"  )
            call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, self%rot(:), start=[1, idslot, tslot], &
                                               count=[NDIM,1,1]), &
                                  "swiftest_io_netcdf_write_frame_cb nf90_put_var cb rot_varid"  )
            
            call netcdf_io_check( nf90_put_var(nc%id, nc%rotphase_varid, self%rotphase, start = [tslot]), & 
                                  "swiftest_io_netcdf_write_frame_cb nf90_put_var cb rotphase")
            
            if (nc%lc_lm_exists) then
               call netcdf_io_check( nf90_put_var(nc%id, nc%c_lm_varid, self%c_lm, count = [nc%m_dim_max, nc%l_dim_max, 2]), &
                                       "netcdf_io_write_frame_cb nf90_put_var c_lm_varid")
            else
               call netcdf_io_check( nf90_put_var(nc%id, nc%j2rp2_varid, self%j2rp2, start=[tslot]), &
                                    "swiftest_io_netcdf_write_frame_cb nf90_put_var cb j2rp2_varid" )
               call netcdf_io_check( nf90_put_var(nc%id, nc%j4rp4_varid, self%j4rp4, start=[tslot]), &
                                    "swiftest_io_netcdf_write_frame_cb nf90_put_var cb j4rp4_varid" )
            end if
         end if

         call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp), &
                                  "swiftest_io_netcdf_write_frame_cb nf90_set_fill old_mode"  )
      end associate

      return
   end subroutine swiftest_io_netcdf_write_frame_cb


   module subroutine swiftest_io_netcdf_write_frame_system(self, nc, param)
      !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
      !!
      !! Write a frame (header plus records for each massive body and active test particle) to a output binary file
      implicit none
      ! Arguments
      class(swiftest_nbody_system),      intent(inout) :: self  !! Swiftest nbody_system object
      class(swiftest_netcdf_parameters), intent(inout) :: nc    !! Parameters used to for writing a NetCDF dataset to file
      class(swiftest_parameters),        intent(inout) :: param !! Current run configuration parameters 

      call self%write_hdr(nc, param)
#ifdef COARRAY
      if (this_image() == 1) then
#endif
         call self%cb%write_frame(nc, param)
         call self%pl%write_frame(nc, param)
         select type(self)
         class is(ringmoons_nbody_system)
            call self%ring%write_frame(param)
            self%seed%maxid = self%maxid
            call self%seed%write_frame(nc, param)
            self%maxid = self%seed%maxid
         end select
#ifdef COARRAY
      end if ! this_image() == 1
#endif
      call self%tp%write_frame(nc, param)

      return
   end subroutine swiftest_io_netcdf_write_frame_system


   module subroutine swiftest_io_netcdf_write_hdr_system(self, nc, param) 
      !! author: David A. Minton
      !!
      !! Writes header information (variables that change with time, but not particle id). 
      !! This subroutine swiftest_significantly improves the output over the original binary file, allowing us to track energy,
      !! momentum, and other quantities that previously were handled as separate output files.
      implicit none
      ! Arguments
      class(swiftest_nbody_system),      intent(in)    :: self  !! Swiftest nbody system object
      class(swiftest_netcdf_parameters), intent(inout) :: nc    !! Parameters used to for writing a NetCDF dataset to file
      class(swiftest_parameters),        intent(inout) :: param !! Current run configuration parameters
      ! Internals
      integer(I4B) :: tslot

      call nc%find_tslot(self%t, tslot)
      call netcdf_io_check( nf90_put_var(nc%id, nc%time_varid, self%t, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var time_varid"  )
      call netcdf_io_check( nf90_put_var(nc%id, nc%npl_varid, self%pl%nbody, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var npl_varid"  )
#ifndef COARRAY
      call netcdf_io_check( nf90_put_var(nc%id, nc%ntp_varid, self%tp%nbody, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var ntp_varid"  )
#endif
      if (param%lmtiny_pl) call netcdf_io_check( nf90_put_var(nc%id, nc%nplm_varid, self%pl%nplm, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var nplm_varid"  )

      if (param%lenergy) then
         call netcdf_io_check( nf90_put_var(nc%id, nc%KE_orb_varid, self%ke_orbit, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var KE_orb_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%KE_rot_varid, self%ke_rot, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var KE_rot_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var PE_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%BE_varid, self%be, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var BE_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%TE_varid, self%te, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var TE_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%L_orbit_varid, self%L_orbit(:), start=[1,tslot], count=[NDIM,1]), &
                                  "netcdf_io_write_hdr_system nf90_put_var L_orbit_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%L_rot_varid, self%L_rot(:), start=[1,tslot], count=[NDIM,1]), &
                                  "netcdf_io_write_hdr_system nf90_put_var L_rot_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%L_escape_varid, self%L_escape(:), start=[1,tslot], count=[NDIM,1]), &
                                  "netcdf_io_write_hdr_system nf90_put_var L_escape_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%E_collisions_varid, self%E_collisions, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var E_collisions_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%E_untracked_varid, self%E_untracked, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var E_untracked_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), &
                                  "netcdf_io_write_hdr_system nf90_put_var GMescape_varid"  )
      end if

      return
   end subroutine swiftest_io_netcdf_write_hdr_system


   module subroutine swiftest_io_netcdf_write_info_body(self, nc, param)
      !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton
      !!
      !! Write all current particle information metadata to file
      implicit none
      ! Arguments
      class(swiftest_body),              intent(in)    :: self  !! Swiftest particle object
      class(swiftest_netcdf_parameters), intent(inout) :: nc      !! Parameters used to identify a particular NetCDF dataset
      class(swiftest_parameters),        intent(inout) :: param !! Current run configuration parameters
      ! Internals
      integer(I4B)                              :: i, j, idslot, old_mode, tmp
      integer(I4B), dimension(:), allocatable   :: ind
      character(len=NAMELEN) :: charstring
      character(len=NAMELEN), dimension(self%nbody) :: origin_type

      call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), &
                                  "netcdf_io_write_info_body nf90_set_fill NF90_NOFILL"  )

      select type(self)
         class is (swiftest_body)
         associate(n => self%nbody, tslot => nc%tslot)
            if (n == 0) return
            call util_sort(self%id(1:n), ind)
            call nc%get_idvals()

            do i = 1, n
               j = ind(i)
               call nc%find_idslot(self%id(j), idslot) 
               call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, self%id(j), start=[idslot]), &
                                  "netcdf_io_write_info_body nf90_put_var id_varid"  )
               call netcdf_io_check( nf90_put_var(nc%id, nc%status_varid, self%status(j), start=[idslot,tslot]), &
                                  "netcdf_io_write_info_body nf90_put_var status_varid"  )

               charstring = trim(adjustl(self%info(j)%name))
               call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), &
                                  "netcdf_io_write_info_body nf90_put_var name_varid"  )

               charstring = trim(adjustl(self%info(j)%particle_type))
               call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), &
                                  "netcdf_io_write_info_body nf90_put_var particle_type_varid"  )

               if (param%lclose) then
                  charstring = trim(adjustl(self%info(j)%origin_type))
                  origin_type(i) = charstring
                  call netcdf_io_check( nf90_put_var(nc%id, nc%origin_type_varid, charstring, start=[1, idslot], &
                                                     count=[NAMELEN, 1]), &
                                  "netcdf_io_write_info_body nf90_put_var origin_type_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%origin_time_varid,  self%info(j)%origin_time,  start=[idslot]), &
                                  "netcdf_io_write_info_body nf90_put_var origin_time_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%origin_rh_varid,    self%info(j)%origin_rh(:), start=[1,idslot], &
                                                     count=[NDIM,1]), &
                                  "netcdf_io_write_info_body nf90_put_var origin_rh_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%origin_vh_varid,    self%info(j)%origin_vh(:), start=[1,idslot], &
                                                     count=[NDIM,1]), &
                                  "netcdf_io_write_info_body nf90_put_var origin_vh_varid"  )
   
                  call netcdf_io_check( nf90_put_var(nc%id, nc%collision_id_varid, self%info(j)%collision_id, start=[idslot]), &
                                  "netcdf_io_write_info_body nf90_put_var collision_id_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%discard_time_varid, self%info(j)%discard_time, start=[idslot]), &
                                  "netcdf_io_write_info_body nf90_put_var discard_time_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%discard_rh_varid,   self%info(j)%discard_rh(:), start=[1,idslot], &
                                                     count=[NDIM,1]), &
                                  "netcdf_io_write_info_body nf90_put_var discard_rh_varid"  )
                  call netcdf_io_check( nf90_put_var(nc%id, nc%discard_vh_varid,   self%info(j)%discard_vh(:), start=[1,idslot], &
                                                     count=[NDIM,1]), &
                                  "netcdf_io_write_info_body nf90_put_var discard_vh_varid"  )
               end if

            end do
         end associate
      end select

      call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp) )
      return
   end subroutine swiftest_io_netcdf_write_info_body


   module subroutine swiftest_io_netcdf_write_info_cb(self, nc, param)
      !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton
      !!
      !! Write the central body info to file
      implicit none
      class(swiftest_cb),               intent(in)    :: self  !! Swiftest particle object
      class(swiftest_netcdf_parameters), intent(inout) :: nc      !! Parameters used to identify a particular NetCDF dataset
      class(swiftest_parameters),           intent(inout) :: param !! Current run configuration parameters
      ! Internals
      integer(I4B)                              :: idslot, old_mode, tmp
      character(len=NAMELEN) :: charstring

      ! This string of spaces of length NAMELEN is used to clear out any old data left behind inside the string variables
      call netcdf_io_check( nf90_set_fill(nc%id, NF90_NOFILL, old_mode), &
                                  "netcdf_io_write_info_cb nf90_set_fill NF90_NOFILL"  )

      call nc%find_idslot(self%id, idslot) 

      call netcdf_io_check( nf90_put_var(nc%id, nc%id_varid, self%id, start=[idslot]), &
                                  "netcdf_io_write_info_cb nf90_put_var id_varid"  )

      charstring = trim(adjustl(self%info%name))
      call netcdf_io_check( nf90_put_var(nc%id, nc%name_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), &
                                  "netcdf_io_write_info_cb nf90_put_var name_varid"  )

      charstring = trim(adjustl(self%info%particle_type))
      call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), &
                                  "netcdf_io_write_info_cb nf90_put_var ptype_varid"  )

      if (param%lclose) then
         charstring = trim(adjustl(self%info%origin_type))
         call netcdf_io_check( nf90_put_var(nc%id, nc%origin_type_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), &
                                  "netcdf_io_write_info_body nf90_put_var cb origin_type_varid"  )

         call netcdf_io_check( nf90_put_var(nc%id, nc%origin_time_varid, self%info%origin_time, start=[idslot]), &
                                  "netcdf_io_write_info_body nf90_put_var cb origin_time_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%origin_rh_varid, self%info%origin_rh(:), start=[1, idslot], count=[NDIM,1]), &
                                  "netcdf_io_write_info_body nf90_put_var cb origin_rh_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%origin_vh_varid, self%info%origin_vh(:), start=[1, idslot], count=[NDIM,1]), &
                                  "netcdf_io_write_info_body nf90_put_var cb origin_vh_varid"  )

         call netcdf_io_check( nf90_put_var(nc%id, nc%collision_id_varid, self%info%collision_id, start=[idslot]), &
                                  "netcdf_io_write_info_body nf90_put_var cb collision_id_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%discard_time_varid, self%info%discard_time, start=[idslot]), &
                                  "netcdf_io_write_info_body nf90_put_var cb discard_time_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%discard_rh_varid, self%info%discard_rh(:), start=[1, idslot], &
                                            count=[NDIM,1]), &
                                  "netcdf_io_write_info_body nf90_put_var cb discard_rh_varid"  )
         call netcdf_io_check( nf90_put_var(nc%id, nc%discard_vh_varid, self%info%discard_vh(:), start=[1, idslot], &
                                            count=[NDIM,1]), &
                                  "netcdf_io_write_info_body nf90_put_var cb discard_vh_varid"  )
      end if
      call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp) )

      return
   end subroutine swiftest_io_netcdf_write_info_cb


   module subroutine swiftest_io_remove_nul_char(string)
      !! author: David A. Minton
      !!
      !! Remove spaces and trailing NUL characters that are introduced by NetCDF strings
      implicit none
      ! Arguments
      character(len=*), intent(inout) :: string !! String to make upper case
      ! Internals
      integer(I4B) :: pos
  
      pos = index(string, achar(0)) - 1
      if (pos > 0) then
         string = trim(adjustl(string(1:pos)))
      else
         string = trim(adjustl(string))
      end if

      return
   end subroutine swiftest_io_remove_nul_char


   module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) 
      !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
      !!
      !! Read in parameters for the integration
      !! as the newline characters are ignored in the input file when compiled in ifort.
      !!
      !! Adapted from David E. Kaufmann's Swifter routine io_init_param.f90
      !! Adapted from Martin Duncan's Swift routine io_init_param.f
      implicit none
      ! Arguments
      class(swiftest_parameters), intent(inout) :: self       !! Collection of parameters
      integer(I4B),               intent(in)    :: unit       !! File unit number
      character(len=*),           intent(in)    :: iotype     !! Dummy argument passed to the  input/output procedure contains the 
                                                              !!    text from the char-literal-constant, prefixed with DT. 
                                                              !!    If you do not include a char-literal-constant, the iotype 
                                                              !!    argument contains only DT.
      character(len=*),           intent(in)    :: v_list(:)  !! The first element passes the integrator code to the reader
      integer(I4B),               intent(out)   :: iostat     !! IO status code
      character(len=*),           intent(inout) :: iomsg      !! Message to pass if iostat /= 0
      ! Internals
      logical                        :: tstart_set = .false.              !! Is the final time set in the input file?
      logical                        :: tstop_set = .false.               !! Is the final time set in the input file?
      logical                        :: dt_set = .false.                  !! Is the step size set in the input file?
      integer(I4B)                   :: ilength, ifirst, ilast, i         !! Variables used to parse input file
      character(STRMAX)              :: line                              !! Line of the input file
      character(len=:), allocatable  :: line_trim,param_name, param_value !! Strings used to parse the param file
      character(*),parameter         :: linefmt = '(A)'                   !! Format code for simple text string
      integer(I4B)                   :: nseeds, nseeds_from_file
      logical                        :: seed_set = .false.      !! Is the random seed set in the input file?
      real(DP)                       :: tratio, y
#ifdef COARRAY
      type(swiftest_parameters), codimension[:], allocatable :: param

      select type(self)
      type is (swiftest_parameters)
         allocate(param[*], source = self)
      end select

#else
      associate(param => self) 
#endif
         ! Parse the file line by line, extracting tokens then matching them up with known parameters if possible
         call random_seed(size = nseeds)
         if (allocated(param%seed)) deallocate(param%seed)
         allocate(param%seed(nseeds))
#ifdef COARRAY
      if (this_image() == 1) then
#endif
         open(unit = unit, file = param%param_file_name, status = 'old', err = 667, iomsg = iomsg)
         do
            read(unit = unit, fmt = linefmt, end = 1, err = 667, iomsg = iomsg) line
            line_trim = trim(adjustl(line))
            ilength = len(line_trim)
            if ((ilength /= 0)) then 
               ifirst = 1
               ! Read the pair of tokens. The first one is the parameter name, the second is the value.
               param_name = swiftest_io_get_token(line_trim, ifirst, ilast, iostat)
               if (param_name == '') cycle ! No parameter name (usually because this line is commented out)
               call swiftest_io_toupper(param_name)
               ifirst = ilast + 1
               param_value = swiftest_io_get_token(line_trim, ifirst, ilast, iostat)
               select case (param_name)
               case ("T0")
                  read(param_value, *, err = 667, iomsg = iomsg) param%t0
               case ("TSTART")
                  read(param_value, *, err = 667, iomsg = iomsg) param%tstart
                  tstart_set = .true.                  
               case ("TSTOP")
                  read(param_value, *, err = 667, iomsg = iomsg) param%tstop
                  tstop_set = .true.
               case ("DT")
                  read(param_value, *, err = 667, iomsg = iomsg) param%dt
                  dt_set = .true.
               case ("CB_IN")
                  param%incbfile = param_value
               case ("PL_IN")
                  param%inplfile = param_value
               case ("TP_IN")
                  param%intpfile = param_value
               case ("NC_IN")
                  param%nc_in = param_value
               case ("IN_TYPE")
                  call swiftest_io_toupper(param_value)
                  param%in_type = param_value
               case ("IN_FORM")
                  call swiftest_io_toupper(param_value)
                  param%in_form = param_value
               case ("ISTEP_OUT")
                  read(param_value, *) param%istep_out
               case ("NSTEP_OUT")
                  read(param_value, *) param%nstep_out
               case ("BIN_OUT")
                  param%outfile = param_value
               case ("OUT_TYPE")
                  call swiftest_io_toupper(param_value)
                  param%out_type = param_value
               case ("OUT_FORM")
                  call swiftest_io_toupper(param_value)
                  param%out_form = param_value
               case ("OUT_STAT")
                  call swiftest_io_toupper(param_value)
                  param%out_stat = param_value
               case ("DUMP_CADENCE")
                  read(param_value, *, err = 667, iomsg = iomsg) param%dump_cadence
               case ("CHK_CLOSE")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T') param%lclose = .true.
               case ("CHK_RMIN")
                  read(param_value, *, err = 667, iomsg = iomsg) param%rmin
               case ("CHK_RMAX")
                  read(param_value, *, err = 667, iomsg = iomsg) param%rmax
               case ("CHK_EJECT")
                  read(param_value, *, err = 667, iomsg = iomsg) param%rmaxu
               case ("CHK_QMIN")
                  read(param_value, *, err = 667, iomsg = iomsg) param%qmin
               case ("CHK_QMIN_COORD")
                  call swiftest_io_toupper(param_value)
                  param%qmin_coord = param_value
               case ("CHK_QMIN_RANGE")
                  read(param_value, *, err = 667, iomsg = iomsg) param%qmin_alo
                  ifirst = ilast + 2
                  param_value = swiftest_io_get_token(line, ifirst, ilast, iostat)
                  read(param_value, *, err = 667, iomsg = iomsg) param%qmin_ahi
               case ("EXTRA_FORCE")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T') param%lextra_force = .true.
               case ("BIG_DISCARD")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T' ) param%lbig_discard = .true.
               case ("RHILL_PRESENT")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T' ) param%lrhill_present = .true.
               case ("MU2KG")
                  read(param_value, *, err = 667, iomsg = iomsg) param%MU2KG
               case ("TU2S")
                  read(param_value, *, err = 667, iomsg = iomsg) param%TU2S
               case ("DU2M")
                  read(param_value, *, err = 667, iomsg = iomsg) param%DU2M
               case ("ENERGY")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T') param%lenergy = .true.
               case ("GR")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T') param%lgr = .true. 
               case ("ROTATION")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T') param%lrotation = .true. 
               case ("TIDES")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T') param%ltides = .true. 
               case("RADIATION")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T') param%lradiation = .true.
               case("YARKOVSKY")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T') param%lyarkovsky = .true.
               case("YARKOVSKY_SCHACH")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T') param%lyarkovsky_schach = .true.
               case ("INTERACTION_LOOPS")
                  call swiftest_io_toupper(param_value)
                  param%interaction_loops = param_value
               case ("ENCOUNTER_CHECK_PLPL")
                  call swiftest_io_toupper(param_value)
                  param%encounter_check_plpl = param_value
               case ("ENCOUNTER_CHECK_PLTP")
                  call swiftest_io_toupper(param_value)
                  param%encounter_check_pltp = param_value
               case ("ENCOUNTER_CHECK")
                  call swiftest_io_toupper(param_value)
                  param%encounter_check_plpl = param_value
                  param%encounter_check_pltp = param_value
               case ("FIRSTKICK")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "NO" .or. param_value == 'F') param%lfirstkick = .false. 
               case ("FIRSTENERGY")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "NO" .or. param_value == 'F') param%lfirstenergy = .false. 
               case("EORBIT_ORIG")
                  read(param_value, *, err = 667, iomsg = iomsg) param%E_orbit_orig 
               case("GMTOT_ORIG")
                  read(param_value, *, err = 667, iomsg = iomsg) param%GMtot_orig 
               case("LTOT_ORIG")
                  read(param_value, *, err = 667, iomsg = iomsg) param%L_total_orig(1)
                  do i = 2, NDIM
                     ifirst = ilast + 2
                     param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) 
                     read(param_value, *, err = 667, iomsg = iomsg) param%L_total_orig(i)
                  end do
               case("LORBIT_ORIG")
                  read(param_value, *, err = 667, iomsg = iomsg) param%L_orbit_orig(1)
                  do i = 2, NDIM
                     ifirst = ilast + 2
                     param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) 
                     read(param_value, *, err = 667, iomsg = iomsg) param%L_orbit_orig(i)
                  end do
               case("LROT_ORIG")
                  read(param_value, *, err = 667, iomsg = iomsg) param%L_rot_orig(1)
                  do i = 2, NDIM
                     ifirst = ilast + 2
                     param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) 
                     read(param_value, *, err = 667, iomsg = iomsg) param%L_rot_orig(i)
                  end do
               case("LESCAPE")
                  read(param_value, *, err = 667, iomsg = iomsg) param%L_escape(1)
                  do i = 2, NDIM
                     ifirst = ilast + 2
                     param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) 
                     read(param_value, *, err = 667, iomsg = iomsg) param%L_escape(i)
                  end do
               case("GMESCAPE")
                  read(param_value, *, err = 667, iomsg = iomsg) param%GMescape 
               case("ECOLLISIONS")
                  read(param_value, *, err = 667, iomsg = iomsg) param%E_collisions
               case("EUNTRACKED")
                  read(param_value, *, err = 667, iomsg = iomsg) param%E_untracked
               case ("COLLISION_MODEL")
                  call swiftest_io_toupper(param_value)
                  read(param_value, *) param%collision_model
               case ("GMTINY")
                  read(param_value, *) param%GMTINY
               case ("MIN_GMFRAG")
                  read(param_value, *) param%min_GMfrag
               case ("NFRAG_REDUCTION")
                  read(param_value, *) param%nfrag_reduction
               case ("ENCOUNTER_SAVE")
                  call swiftest_io_toupper(param_value)
                  read(param_value, *) param%encounter_save
               case ("COARRAY")
                  call swiftest_io_toupper(param_value)
                  if (param_value == "YES" .or. param_value == 'T') param%lcoarray = .true. 
               case("RING_FILE")
                  read(param_value, *) param%ring_file
               case("SEED")
                  read(param_value, *) nseeds_from_file
                  ! Because the number of seeds can vary between compilers/systems, we need to make sure we can handle cases in 
                  ! which the input file has a different number of seeds than the current nbody_system. If the number of seeds in 
                  ! the file is smaller than required, we will use them as a source to fill in the missing elements.
                  ! If the number of seeds in the file is larger than required, we will truncate the seed array.
                  if (nseeds_from_file > nseeds) then
                     nseeds = nseeds_from_file
                     deallocate(param%seed)
                     allocate(param%seed(nseeds))
                     do i = 1, nseeds
                        ifirst = ilast + 2
                        param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) 
                        read(param_value, *) param%seed(i)
                     end do
                  else ! Seed array in file is too small
                     do i = 1, nseeds_from_file
                        ifirst = ilast + 2
                        param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) 
                        read(param_value, *) param%seed(i)
                     end do
                     param%seed(nseeds_from_file+1:nseeds) = [(param%seed(1) - param%seed(nseeds_from_file) + i, &
                                                               i=nseeds_from_file+1, nseeds)]
                  end if
                  seed_set = .true.
               case ("RESTART")
                  if (param_value == "NO" .or. param_value == 'F') then
                     param%lrestart = .false. 
                  else if (param_value == "YES" .or. param_value == 'T') then
                     param%lrestart = .true.
                  end if 
               ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters
               case ("NPLMAX", "NTPMAX", "YORP")
               case default
                  write(*,*) "Ignoring unknown parameter -> ",param_name
               end select
            end if
         end do
         1 continue
         close(unit)
         iostat = 0

         ! Do basic sanity checks on the input values
         if ((.not. tstart_set) .or. (.not. tstop_set) .or. (.not. dt_set)) then
            write(iomsg,*) 'Valid simulation time not set'
            iostat = -1
            return
         end if
         if (param%dt <= 0.0_DP) then
            write(iomsg,*) 'Invalid timestep: '
            iostat = -1
            return
         end if
         if (param%inplfile == "") then
            write(iomsg,*) 'No valid massive body file in input file'
            iostat = -1
            return
         end if
         if ((param%in_type /= "ASCII") .and. (param%in_type /= "NETCDF_FLOAT") .and. (param%in_type /= "NETCDF_DOUBLE"))  then
            write(iomsg,*) 'Invalid input file type:',trim(adjustl(param%in_type))
            iostat = -1
            return
         end if
         if (param%istep_out <= 0) then
            write(iomsg,*) 'Invalid ISTEP_OUT. Must be a positive integer'
            iostat = -1
            return
         end if
         if (param%nstep_out <= 0) then
            param%nstep_out = int((param%tstop - param%t0) / (param%istep_out * param%dt))
            param%fstep_out = 1.0_DP ! Linear output time
            param%ltstretch = .false.
         else
            param%fstep_out = 1._DP
            tratio = (param%TSTOP - param%T0) / (param%istep_out * param%dt)
            if (int(tratio) == param%nstep_out) then
               param%ltstretch = .false.
            else
               param%ltstretch = .true.
               y = time_stretcher(param%fstep_out) 
               call solve_roots(time_stretcher,param%fstep_out)
            end if
         end if
         if (param%dump_cadence < 0) then
            write(iomsg,*) 'Invalid DUMP_CADENCE. Must be a positive integer or 0.'
            iostat = -1
            return
         end if
         if ((param%istep_out > 0) .and. (param%outfile == "")) then
            write(iomsg,*) 'Invalid outfile'
            iostat = -1
            return
         end if
         param%lrestart = (param%out_stat == "APPEND")
         if (param%outfile /= "") then
            if ((param%out_type /= "NETCDF_FLOAT") .and. (param%out_type /= "NETCDF_DOUBLE")) then
               write(iomsg,*) 'Invalid out_type: ',trim(adjustl(param%out_type))
               iostat = -1
               return
            end if
            if ((param%out_form /= "EL") .and. (param%out_form /= "XV") .and. (param%out_form /= "XVEL")) then
               write(iomsg,*) 'Invalid out_form: ',trim(adjustl(param%out_form))
               iostat = -1
               return
            end if
            if ((param%out_stat /= "NEW") .and. (param%out_stat /= "REPLACE") .and. (param%out_stat /= "APPEND")  &
            .and. (param%out_stat /= "UNKNOWN")) then
               write(iomsg,*) 'Invalid out_stat: ',trim(adjustl(param%out_stat))
               iostat = -1
               return
            end if
         end if
         if (param%qmin > 0.0_DP) then
            if ((param%qmin_coord /= "HELIO") .and. (param%qmin_coord /= "BARY")) then
               write(iomsg,*) 'Invalid qmin_coord: ',trim(adjustl(param%qmin_coord))
               iostat = -1
               return
            end if
            if ((param%qmin_alo <= 0.0_DP) .or. (param%qmin_ahi <= 0.0_DP)) then
               write(iomsg,*) 'Invalid qmin vals'
               iostat = -1
               return
            end if
         end if
         if (param%ltides .and. .not. param%lrotation) then
            write(iomsg,*) 'Tides require rotation to be turned on'
            iostat = -1
            return
         end if

         if (param%lyarkovsky .and. .not. param%lrotation) then
            write(iomsg,*) 'Yarkovsky forces require rotation to be turned on'
            iostat = -1
            return
         end if

         if (param%lyarkovsky_schach .and. .not. param%lrotation) then
            write(iomsg,*) 'Yarkovsky-Schach forces require rotation to be turned on'
            iostat = -1
            return
         end if

         if ((param%MU2KG < 0.0_DP) .or. (param%TU2S < 0.0_DP) .or. (param%DU2M < 0.0_DP)) then
            write(iomsg,*) 'Invalid unit conversion factor'
            iostat = -1
            return
         end if

         ! Calculate the G for the nbody_system units
         param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2))

         if ((param%encounter_save /= "NONE")       .and. &
               (param%encounter_save /= "TRAJECTORY") .and. &
               (param%encounter_save /= "CLOSEST")    .and. &
               (param%encounter_save /= "BOTH")) then
            write(iomsg,*) 'Invalid encounter_save parameter: ',trim(adjustl(param%out_type))
            write(iomsg,*) 'Valid options are NONE, TRAJECTORY, CLOSEST, or BOTH'
            iostat = -1
            return
         end if

         param%lenc_save_trajectory = (param%encounter_save == "TRAJECTORY") .or. (param%encounter_save == "BOTH")
         param%lenc_save_closest = (param%encounter_save == "CLOSEST") .or. (param%encounter_save == "BOTH")

         if ((param%integrator == INT_RMVS) .or. (param%integrator == INT_SYMBA) .or. (param%integrator == INT_RINGMOONS)) then
            if (.not.param%lclose) then
               write(iomsg,*) 'This integrator requires CHK_CLOSE to be enabled.'
               iostat = -1
               return
            end if
         end if

         param%lmtiny_pl = (param%integrator == INT_SYMBA) .or. (param%integrator == INT_RINGMOONS) 

         if (param%lmtiny_pl .and. param%GMTINY < 0.0_DP) then
            write(iomsg,*) "GMTINY invalid or not set: ", param%GMTINY
            iostat = -1
            return
         end if

         if ((param%collision_model /= "MERGE")       .and. &
               (param%collision_model /= "BOUNCE")    .and. &
               (param%collision_model /= "FRAGGLE")) then
            write(iomsg,*) 'Invalid collision_model parameter: ',trim(adjustl(param%out_type))
            write(iomsg,*) 'Valid options are MERGE, BOUNCE, or FRAGGLE'
            iostat = -1
            return
         end if

         if (seed_set) then
            call random_seed(put = param%seed)
         else
            call random_seed(get = param%seed)
         end if
         if (param%collision_model == "FRAGGLE" ) then
            if (param%min_GMfrag < 0.0_DP) param%min_GMfrag = param%GMTINY
            if (param%nfrag_reduction < 1.0_DP) then
               write(iomsg,*) "Warning: NFRAG_REDUCTION value invalid. Setting to 1.0" 
               param%nfrag_reduction = 1.0_DP
            end if
         end if

         ! Calculate Solar Luminosity, stefan-boltzmann constant, and inv_c2 in system units 
         if (param%lradiation .or. param%lyarkovsky) then
            param%L_SUN_sys = L_SUN / param%MU2KG / param%DU2M**2 * param%TU2S**3
            param%sigma_sys = SIGMA /param%MU2KG * param%TU2S**3 ! system units / K^4
            param%inv_c2 = einsteinC * param%TU2S / param%DU2M
            param%inv_c2 = (param%inv_c2)**(-2)
            param%yark_radius_threshold_sys = YARK_RADIUS_THRESHOLD / param%DU2M
            ! param%lgr = .true. ! placeholder for if GR is needed by future radiation or Yarkovsky models that we implement
         end if

         ! Calculate Solar Luminosity, stefan-boltzmann constant, and inv_c2 in system units 
         if (param%lradiation .or. param%lyarkovsky_schach) then
            param%L_SUN_sys = L_SUN / param%MU2KG / param%DU2M**2 * param%TU2S**3
            param%sigma_sys = SIGMA /param%MU2KG * param%TU2S**3 ! system units / K^4
            param%inv_c2 = einsteinC * param%TU2S / param%DU2M
            param%inv_c2 = (param%inv_c2)**(-2)
            ! param%lgr = .true. ! placeholder for if GR is needed by future radiation or Yarkovsky models that we implement
         end if

         ! Determine if the GR flag is set correctly for this integrator
         select case(param%integrator)
         case(INT_WHM, INT_RMVS, INT_HELIO, INT_SYMBA, INT_RINGMOONS)
         case default   
            if (param%lgr) write(iomsg, *) 'GR is not yet implemented for this integrator. This parameter will be ignored.'
            param%lgr = .false.
         end select

         if (param%lgr) then
            ! Calculate the inverse speed of light in the nbody_system units
            param%inv_c2 = einsteinC * param%TU2S / param%DU2M
            param%inv_c2 = (param%inv_c2)**(-2)
         end if

         select case(trim(adjustl(param%interaction_loops)))
         case("TRIANGULAR")
            param%lflatten_interactions = .false.
         case("FLAT")
            param%lflatten_interactions = .true.
         case default
            write(*,*) "Unknown value for parameter INTERACTION_LOOPS: -> ",trim(adjustl(param%interaction_loops))
            write(*,*) "Must be one of the following: TRIANGULAR or FLAT"
            write(*,*) "Using default value of TRIANGULAR"
            param%interaction_loops = "TRIANGULAR"
            param%lflatten_interactions = .false.
         end select

         select case(trim(adjustl(param%encounter_check_plpl)))
         case("TRIANGULAR")
            param%lencounter_sas_plpl = .false.
         case("SORTSWEEP")
            param%lencounter_sas_plpl = .true.
         case default
            write(*,*) "Unknown value for parameter ENCOUNTER_CHECK_PLPL: -> ",trim(adjustl(param%encounter_check_plpl))
            write(*,*) "Must be one of the following: TRIANGULAR or SORTSWEEP"
            write(*,*) "Using default value of TRIANGULAR"
            param%encounter_check_plpl = "TRIANGULAR"
            param%lencounter_sas_plpl = .false.
         end select

         select case(trim(adjustl(param%encounter_check_pltp)))
         case("TRIANGULAR")
            param%lencounter_sas_pltp = .false.
         case("SORTSWEEP")
            param%lencounter_sas_pltp = .true.
         case default
            write(*,*) "Unknown value for parameter ENCOUNTER_CHECK_PLTP: -> ",trim(adjustl(param%encounter_check_pltp))
            write(*,*) "Must be one of the following: TRIANGULAR or SORTSWEEP"
            write(*,*) "Using default value of TRIANGULAR"
            param%encounter_check_pltp = "TRIANGULAR"
            param%lencounter_sas_pltp = .false.
         end select

         if (param%lcoarray) then
#ifdef COARRAY
            if (num_images() == 1) then
               write(iomsg, *) "Only one Coarray image detected. Coarrays will not be used."
               param%lcoarray = .false.
            end if

            select case(param%integrator)
            case(INT_WHM, INT_RMVS, INT_HELIO)
            case default   
               write(iomsg, *) "Coarray-based parallelization of test particles are not compatible with this integrator. " &
                              // "This parameter will be ignored."
               param%lcoarray = .false.
            end select
#else
            write(iomsg,*) "Coarray capability not detected. Swiftest must be compiled with Coarrays enabled. to use this feature."
            param%lcoarray = .false.
#endif
         end if

         if ((param%integrator == INT_SYMBA) .or. (param%integrator == INT_RINGMOONS)) then
            if (.not.param%lenergy) then
               write(iomsg,*) 'This integrator requires ENERGY to be enabled.'
               iostat = -1
               return
            end if
            if (.not.param%lrotation) then
               write(iomsg,*) 'This integrator requires ROTATION to be enabled.'
               iostat = -1
               return
            end if
         end if

         iostat = 0
#ifdef COARRAY
      end if ! this_image() == 1
#else
      end associate
#endif
      select type(self)
      type is (swiftest_parameters)
#ifdef COARRAY
         call coclone(param)
         self = param
#endif
         call self%set_display(self%display_style)

         if (.not.self%lrestart) then
#ifdef COARRAY
            if (this_image() == 1 .or. self%log_output) then
#endif
               call self%writer(unit = self%display_unit, iotype = "none", v_list = [0], iostat = iostat, iomsg = iomsg)
               if (self%log_output) flush(self%display_unit) 
#ifdef COARRAY
            end if !(this_image() == 1)
            write(COLLISION_LOG_OUT,'("collision_coimage",I0.4,".log")') this_image()
#endif
            ! A minimal log of collision outcomes is stored in the following log file
            ! More complete data on collisions is stored in the NetCDF output files
            call swiftest_io_log_start(self, COLLISION_LOG_OUT, "Collision logfile")
         end if
         ! Print the contents of the parameter file to standard output
      end select

      return 
      667 continue
      write(*,*) "Error reading param file: ", trim(adjustl(iomsg))

      contains
         function time_stretcher(fstep_out) result(ans)
            !! author: David A. Minton
            !!
            !! Equation for the time stretchinf function. Solving the roots of this equation gives the time stretching factor for 
            !! non-linear file output cadence.
            implicit none
            ! Arguments
            real(DP), intent(in) :: fstep_out
            ! Result
            real(DP)             :: ans

            if (abs(fstep_out-1.0_DP) < epsilon(1.0_DP)) then
               ans = self%nstep_out - tratio
            else
               ans = (1.0_DP - fstep_out**(self%nstep_out))/ (1.0_DP - fstep_out) - tratio
            end if

            return
         end function time_stretcher

   end subroutine swiftest_io_param_reader


   module subroutine swiftest_io_param_writer(self, unit, iotype, v_list, iostat, iomsg) 
      !! author: David A. Minton
      !!
      !! Dump integration parameters to file
      !!
      !! Adapted from David E. Kaufmann's Swifter routine io_param_restart.f90
      !! Adapted from Martin Duncan's Swift routine io_param_restart.f
      implicit none
      ! Arguments
      class(swiftest_parameters),intent(in)     :: self       !! Collection of parameters
      integer, intent(in)                       :: unit       !! File unit number
      character(len=*), intent(in)              :: iotype     !! Dummy argument passed to the  input/output procedure contains the 
                                                              !! text from the char-literal-constant, prefixed with DT. 
                                                              !!    If you do not include a char-literal-constant, the iotype 
                                                              !!    argument contains only DT.
      integer, intent(in)                       :: v_list(:)  !! Not used in this procedure
      integer, intent(out)                      :: iostat     !! IO status code
      character(len=*), intent(inout)           :: iomsg      !! Message to pass if iostat /= 0
      ! Internals
      integer(I4B)                   :: nseeds

      associate(param => self)
         call io_param_writer_one("RESTART", param%lrestart, unit)
         call io_param_writer_one("T0", param%t0, unit)
         call io_param_writer_one("TSTART", param%tstart, unit)
         call io_param_writer_one("TSTOP", param%tstop, unit)
         call io_param_writer_one("DT", param%dt, unit)
         call io_param_writer_one("IN_TYPE", param%in_type, unit)
         if (param%in_type == "ASCII") then
            call io_param_writer_one("CB_IN", param%incbfile, unit)
            call io_param_writer_one("PL_IN", param%inplfile, unit)
            call io_param_writer_one("TP_IN", param%intpfile, unit)
         else 
            call io_param_writer_one("NC_IN", param%nc_in, unit)
         end if

         call io_param_writer_one("IN_FORM", param%in_form, unit)
         if (param%dump_cadence > 0) call io_param_writer_one("DUMP_CADENCE",param%dump_cadence, unit)
         if (param%istep_out > 0) then
            call io_param_writer_one("ISTEP_OUT", param%istep_out, unit)
            call io_param_writer_one("BIN_OUT", param%outfile, unit)
            call io_param_writer_one("OUT_TYPE", param%out_type, unit)
            call io_param_writer_one("OUT_FORM", param%out_form, unit)
            call io_param_writer_one("OUT_STAT", "APPEND", unit) 
         end if
         call io_param_writer_one("CHK_RMIN", param%rmin, unit)
         call io_param_writer_one("CHK_RMAX", param%rmax, unit)
         call io_param_writer_one("CHK_EJECT", param%rmaxu, unit)
         call io_param_writer_one("CHK_QMIN", param%qmin, unit)
         if (param%qmin >= 0.0_DP) then
            call io_param_writer_one("CHK_QMIN_COORD", param%qmin_coord, unit)
            call io_param_writer_one("CHK_QMIN_RANGE", [param%qmin_alo, param%qmin_ahi], unit)
         end if
         call io_param_writer_one("MU2KG", param%MU2KG, unit)
         call io_param_writer_one("TU2S", param%TU2S , unit)
         call io_param_writer_one("DU2M", param%DU2M, unit)
         call io_param_writer_one("RHILL_PRESENT", param%lrhill_present, unit)
         call io_param_writer_one("EXTRA_FORCE", param%lextra_force, unit)
         call io_param_writer_one("CHK_CLOSE", param%lclose, unit)
         call io_param_writer_one("ENERGY", param%lenergy, unit)
         call io_param_writer_one("GR", param%lgr, unit)
         call io_param_writer_one("ROTATION", param%lrotation, unit)
         call io_param_writer_one("TIDES", param%ltides, unit)
         call io_param_writer_one("RADIATION", param%lradiation, unit)
         call io_param_writer_one("YARKOVSKY", param%lyarkovsky, unit)
         call io_param_writer_one("YARKOVSKY_SCHACH", param%lyarkovsky_schach, unit)
         call io_param_writer_one("INTERACTION_LOOPS", param%interaction_loops, unit)
         call io_param_writer_one("ENCOUNTER_CHECK_PLPL", param%encounter_check_plpl, unit)
         call io_param_writer_one("ENCOUNTER_CHECK_PLTP", param%encounter_check_pltp, unit)
         call io_param_writer_one("ENCOUNTER_SAVE", param%encounter_save, unit)
         call io_param_writer_one("COARRAY", param%lcoarray, unit)

         if (param%lenergy) then
            call io_param_writer_one("FIRSTENERGY", param%lfirstenergy, unit)
         end if
         call io_param_writer_one("FIRSTKICK",param%lfirstkick, unit)

         if (param%GMTINY >= 0.0_DP) call io_param_writer_one("GMTINY",param%GMTINY, unit)
         if (param%min_GMfrag >= 0.0_DP) call io_param_writer_one("MIN_GMFRAG",param%min_GMfrag, unit)
         call io_param_writer_one("COLLISION_MODEL",param%collision_model, unit)
         if (param%collision_model == "FRAGGLE" ) then
            call io_param_writer_one("NFRAG_REDUCTION",param%nfrag_reduction, unit)
            nseeds = size(param%seed)
            call random_seed(get = param%seed)
            call io_param_writer_one("SEED", [nseeds, param%seed(:)], unit)
         end if
         if (param%integrator == INT_RINGMOONS) call io_param_writer_one("RING_FILE", param%ring_file, unit)

   
         iostat = 0
         iomsg = "UDIO not implemented"
      end associate

      return
   end subroutine swiftest_io_param_writer


   module subroutine swiftest_io_param_writer_one_char(param_name, param_value, unit)
      !! author: David A. Minton
      !!
      !! Writes a single parameter name/value pair to a file unit. 
      !! This version is for character param_value type
      implicit none
      ! Arguments
      character(len=*), intent(in) :: param_name  !! Name of parameter to print
      character(len=*), intent(in) :: param_value !! Value of parameter to print
      integer(I4B),     intent(in) :: unit        !! Open file unit number to print parameter to
      ! Internals
      character(len=NAMELEN) :: param_name_fixed_width !! Parameter label converted to fixed-width string
      character(len=STRMAX)  :: iomsg             !! Message to pass if iostat /= 0

      write(param_name_fixed_width, *) param_name
      write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name_fixed_width) // " " // trim(adjustl(param_value)) 

      return
      667 continue
      write(*,*) 'Error writing parameter: ',trim(adjustl(iomsg))
   end subroutine swiftest_io_param_writer_one_char


   module subroutine swiftest_io_param_writer_one_DP(param_name, param_value, unit)
      !! author: David A. Minton
      !!
      !! Writes a single parameter name/value pair to a file unit. 
      !! This version is for real(DP) param_value type
      implicit none
      ! Arguments
      character(len=*), intent(in)    :: param_name  !! Name of parameter to print
      real(DP),         intent(in)    :: param_value !! Value of parameter to print
      integer(I4B),     intent(in)    :: unit        !! Open file unit number to print parameter to
      ! Internals
      character(len=STRMAX) :: param_value_string   !! Parameter value converted to a string
      character(*),parameter :: Rfmt  = '(ES27.19)' !! Format label for real values 

      write(param_value_string,Rfmt) param_value
      call io_param_writer_one(param_name, param_value_string, unit)

      return
   end subroutine swiftest_io_param_writer_one_DP


   module subroutine swiftest_io_param_writer_one_DParr(param_name, param_value, unit)
      !! author: David A. Minton
      !!
      !! Writes a single parameter name/value pair to a file unit. 
      !! This version is for real(DP) arrays () param_value type
      implicit none
      ! Arguments
      character(len=*),       intent(in) :: param_name  !! Name of parameter to print
      real(DP), dimension(:), intent(in) :: param_value !! Value of parameter to print
      integer(I4B),           intent(in) :: unit        !! Open file unit number to print parameter to
      ! Internals
      character(len=STRMAX) :: param_value_string   !! Parameter value converted to a string
      character(*),parameter :: Rfmt  = '(ES27.19)' !! Format label for real values 
      character(len=27) :: arr_val
      integer(I4B) :: i, narr

      narr = size(param_value)
      do i = 1, narr
         write(arr_val, Rfmt) param_value(i)
         if (i == 1) then
            write(param_value_string, *) arr_val
         else
            param_value_string = trim(adjustl(param_value_string)) // " " // arr_val
         end if
      end do

      call io_param_writer_one(param_name, param_value_string, unit)

      return
   end subroutine swiftest_io_param_writer_one_DParr


   module subroutine swiftest_io_param_writer_one_I4B(param_name, param_value, unit)
      !! author: David A. Minton
      !!
      !! Writes a single parameter name/value pair to a file unit. 
      !! This version is for integer(I4B) param_value type
      implicit none
      ! Arguments
      character(len=*), intent(in) :: param_name  !! Name of parameter to print
      integer(I4B),     intent(in) :: param_value !! Value of parameter to print
      integer(I4B),     intent(in) :: unit        !! Open file unit number to print parameter to
      ! Internals
      character(len=STRMAX) :: param_value_string   !! Parameter value converted to a string
      character(*),parameter :: Ifmt  = '(I0)'      !! Format label for integer values

      write(param_value_string,Ifmt) param_value
      call io_param_writer_one(param_name, param_value_string, unit)

      return
   end subroutine swiftest_io_param_writer_one_I4B


   module subroutine swiftest_io_param_writer_one_I8B(param_name, param_value, unit)
      !! author: David A. Minton
      !!
      !! Writes a single parameter name/value pair to a file unit. 
      !! This version is for integer(I8B) param_value type
      implicit none
      ! Arguments
      character(len=*), intent(in)    :: param_name  !! Name of parameter to print
      integer(I8B),     intent(in)    :: param_value !! Value of parameter to print
      integer(I4B),     intent(in)    :: unit        !! Open file unit number to print parameter to
      ! Internals
      character(len=STRMAX) :: param_value_string   !! Parameter value converted to a string
      character(*),parameter :: Ifmt  = '(I0)'      !! Format label for integer values

      write(param_value_string,Ifmt) param_value
      call io_param_writer_one(param_name, param_value_string, unit)

      return
   end subroutine swiftest_io_param_writer_one_I8B


   module subroutine swiftest_io_param_writer_one_I4Barr(param_name, param_value, unit)
      !! author: David A. Minton
      !!
      !! Writes a single parameter name/value pair to a file unit. 
      !! This version is for integer(I4B) arrays param_value type
      implicit none
      ! Arguments
      character(len=*),           intent(in) :: param_name  !! Name of parameter to print
      integer(I4B), dimension(:), intent(in) :: param_value !! Value of parameter to print
      integer(I4B),               intent(in) :: unit        !! Open file unit number to print parameter to
      ! Internals
      character(len=STRMAX) :: param_value_string !! Parameter value converted to a string
      character(*),parameter :: Ifmt  = '(I0)'    !! Format label for integer values
      character(len=25) :: arr_val
      integer(I4B) :: i, narr

      narr = size(param_value)
      do i = 1, narr
         write(arr_val, Ifmt) param_value(i)
         if (i == 1) then
            write(param_value_string, *) trim(adjustl(arr_val))
         else
            param_value_string = trim(adjustl(param_value_string)) // " " // trim(adjustl(arr_val))
         end if
      end do

      call io_param_writer_one(param_name, param_value_string, unit)

      return
   end subroutine swiftest_io_param_writer_one_I4Barr


   module subroutine swiftest_io_param_writer_one_logical(param_name, param_value, unit)
      !! author: David A. Minton
      !!
      !! Writes a single parameter name/value pair to a file unit. 
      !! This version is for logical param_value type
      implicit none
      ! Arguments
      character(len=*), intent(in) :: param_name  !! Name of parameter to print
      logical,          intent(in) :: param_value !! Value of parameter to print
      integer(I4B),     intent(in) :: unit        !! Open file unit number to print parameter to
      ! Internals
      character(len=STRMAX) :: param_value_string   !! Parameter value converted to a string
      character(*),parameter :: Lfmt  = '(L1)'         !! Format label for logical values 

      write(param_value_string,Lfmt) param_value
      call io_param_writer_one(param_name, param_value_string, unit)

      return
   end subroutine swiftest_io_param_writer_one_logical

#ifdef QUADPREC
   module subroutine swiftest_io_param_writer_one_QP(param_name, param_value, unit)
      !! author: David A. Minton
      !!
      !! Writes a single parameter name/value pair to a file unit. 
      !! This version is for real(QP) param_value type
      implicit none
      ! Arguments
      character(len=*), intent(in) :: param_name  !! Name of parameter to print
      real(QP),         intent(in) :: param_value !! Value of parameter to print
      integer(I4B),     intent(in) :: unit        !! Open file unit number to print parameter to
      ! Internals
      character(len=STRMAX) :: param_value_string   !! Parameter value converted to a string
      character(*),parameter :: Rfmt  = '(ES27.19)' !! Format label for real values 

      write(param_value_string,Rfmt) param_value
      call io_param_writer_one(param_name, param_value_string, unit)

      return
   end subroutine swiftest_io_param_writer_one_QP
#endif

   module subroutine swiftest_io_read_in_body(self, param) 
      !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
      !!
      !! Read in either test particle or massive body data 
      !!
      !! Adapted from David E. Kaufmann's Swifter routine swiftest_init_pl.f90 and swiftest_init_tp.f90
      !! Adapted from Martin Duncan's Swift routine swiftest_init_pl.f and swiftest_init_tp.f
      implicit none
      ! Arguments
      class(swiftest_body),       intent(inout) :: self   !! Swiftest particle object
      class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
      ! Internals
      integer(I4B)                  :: iu = LUN
      integer(I4B)                  :: i, nbody
      character(len=:), allocatable :: infile
      character(STRMAX)             :: errmsg
      ! Internals
      integer(I4B)                                :: ierr  !! Error code: returns 0 if the read is successful

      ! Select the appropriate polymorphic class (test particle or massive body)
      if (param%in_type /= "ASCII") return ! Not for NetCDF

      select type(self)
      class is (swiftest_pl)
         infile = param%inplfile
      class is (swiftest_tp)
         infile = param%intpfile
      end select

      open(unit = iu, file = infile, status = 'old', form = 'FORMATTED', err = 667, iomsg = errmsg)
      read(iu, *, err = 667, iomsg = errmsg) nbody

      call self%setup(nbody, param)
      ierr = 0
      if (nbody > 0) then
         ierr = self%read_frame(iu, param)
         self%status(:) = ACTIVE
         self%lmask(:) = .true.
         do i = 1, nbody
            call self%info(i)%set_value(status="ACTIVE")
         end do
      end if
      close(iu, err = 667, iomsg = errmsg)

      if (ierr == 0) return

      667 continue
      write(*,*) 'Error reading in initial conditions file: ',trim(adjustl(errmsg))
      return
   end subroutine swiftest_io_read_in_body

   module subroutine swiftest_io_read_in_cb(self, param) 
      !! author: David A. Minton
      !!
      !! Reads in central body data 
      !!
      !! Adapted from David E. Kaufmann's Swifter routine swiftest_init_pl.f90
      !! Adapted from Martin Duncan's Swift routine swiftest_init_pl.f
      implicit none
      ! Arguments
      class(swiftest_cb),         intent(inout) :: self
      class(swiftest_parameters), intent(inout) :: param
      ! Internals
      integer(I4B)            :: iu = LUN
      character(len=STRMAX)   :: errmsg
      integer(I4B)            :: ierr
      character(len=NAMELEN)  :: name

      if (param%in_type /= "ASCII") return ! Not for NetCDF

      self%id = 0
      open(unit = iu, file = param%incbfile, status = 'old', form = 'FORMATTED', err = 667, iomsg = errmsg)
      read(iu, *, err = 667, iomsg = errmsg) name
      call self%info%set_value(name=name)
      read(iu, *, err = 667, iomsg = errmsg) self%Gmass
      self%mass = real(self%Gmass / param%GU, kind=DP)
      read(iu, *, err = 667, iomsg = errmsg) self%radius
      read(iu, *, err = 667, iomsg = errmsg) self%j2rp2
      read(iu, *, err = 667, iomsg = errmsg) self%j4rp4
      if (param%lrotation) then
         read(iu, *, err = 667, iomsg = errmsg) self%Ip(1), self%Ip(2), self%Ip(3)
         read(iu, *, err = 667, iomsg = errmsg) self%rot(1), self%rot(2), self%rot(3)
      end if
      ierr = 0
      close(iu, err = 667, iomsg = errmsg)

      if (ierr == 0) then
   
         if (param%rmin < 0.0) param%rmin = self%radius
         
         self%GM0 = self%Gmass
         self%dGM = 0.0_DP
         self%R0 = self%radius
         if (param%lrotation) then
            self%L0(:) = self%Ip(3) * self%mass * self%radius**2 * self%rot(:) * DEG2RAD
            self%dL(:) = 0.0_DP
         end if
      end if
      return

      667 continue
      write(*,*) "Error reading central body file: " // trim(adjustl(errmsg))
      call base_util_exit(FAILURE,param%display_unit)
   end subroutine swiftest_io_read_in_cb

   module subroutine swiftest_io_read_in_system(self, nc, param)
      !! author: David A. Minton and Carlisle A. Wishard
      !!
      !! Reads in the nbody_system from input files
      implicit none
      ! Arguments
      class(swiftest_nbody_system),      intent(inout) :: self
      class(swiftest_netcdf_parameters), intent(inout) :: nc     !! Parameters used to identify a particular NetCDF dataset
      class(swiftest_parameters),        intent(inout) :: param
      ! Internals
      integer(I4B) :: ierr, i
      class(swiftest_parameters), allocatable :: tmp_param

      if (param%in_type == "ASCII") then
         call self%cb%read_in(param)
         call self%pl%read_in(param)
         if (self%pl%nbody > 0) self%pl%id(:) = [(i, i = 1, self%pl%nbody)]
         call self%tp%read_in(param)
         if (self%tp%nbody > 0) self%tp%id(:) = [(i, i = self%pl%nbody + 1, self%pl%nbody + 1 + self%tp%nbody)]
         self%maxid = self%pl%nbody + self%tp%nbody
         ! Copy over param file variable inputs
         self%E_orbit_orig = param%E_orbit_orig
         self%GMtot_orig = param%GMtot_orig
         self%L_total_orig(:) = param%L_total_orig(:)
         self%L_orbit_orig(:) = param%L_orbit_orig(:)
         self%L_rot_orig(:) = param%L_rot_orig(:)
         self%L_escape(:) = param%L_escape(:)
         self%E_collisions = param%E_collisions
         self%E_untracked = param%E_untracked
      else
         allocate(tmp_param, source=param)
         nc%file_name = param%nc_in
         tmp_param%out_form = param%in_form
         if (.not. param%lrestart) then
            ! Turn off energy computation so we don't have to feed it into the initial conditions
            tmp_param%lenergy = .false.
         end if
         ierr = self%read_frame(nc, tmp_param)
         deallocate(tmp_param)
         if (ierr /=0) call base_util_exit(FAILURE,param%display_unit)
      end if

      param%lnon_spherical_cb = (self%cb%j2rp2 /= 0.0_DP) .or. (self%cb%j4rp4 /= 0.0_DP) .or. allocated(self%cb%c_lm)
      if (.not.param%lnon_spherical_cb) then
         if (allocated(self%pl%aobl)) deallocate(self%pl%aobl)
         if (allocated(self%tp%aobl)) deallocate(self%tp%aobl)
      else
         if (self%pl%nbody > 0) then
            if (.not. allocated(self%pl%aobl)) allocate(self%pl%aobl(NDIM,self%pl%nbody))
            self%pl%aobl(:,:) = 0.0_DP
         end if
         if (self%tp%nbody > 0) then
            if (.not. allocated(self%tp%aobl)) allocate(self%tp%aobl(NDIM,self%tp%nbody))
            self%tp%aobl(:,:) = 0.0_DP
         end if
      
      
      end if

      return
   end subroutine swiftest_io_read_in_system

   module function swiftest_io_read_frame_body(self, iu, param) result(ierr)
      !! author: David A. Minton
      !!
      !! Reads a frame of output of either test particle or massive body data from a binary output file
      !!
      !! Adapted from David E. Kaufmann's Swifter routine  io_read_frame.f90
      !! Adapted from Hal Levison's Swift routine io_read_frame.f
      implicit none
      ! Arguments
      class(swiftest_body),       intent(inout) :: self   !! Swiftest particle object
      integer(I4B),               intent(inout) :: iu     !! Unit number for the output file to write frame to
      class(swiftest_parameters), intent(inout) :: param  !! Current run configuration parameters 
      ! Result
      integer(I4B)                              :: ierr  !! Error code: returns 0 if the read is successful
      ! Internals
      character(len=STRMAX)   :: errmsg
      character(len=NAMELEN), dimension(self%nbody)  :: name
      integer(I4B) :: i
      real(QP)                      :: val

      if (self%nbody == 0) return

      if ((param%in_form /= "EL") .and. (param%in_form /= "XV")) then
         write(errmsg, *) trim(adjustl(param%in_form)) // " is not a recognized format code for input files."
         goto 667
      end if

      associate(n => self%nbody)

         if (param%in_form == "EL") then
            if (.not.allocated(self%a))     allocate(self%a(n))
            if (.not.allocated(self%e))     allocate(self%e(n))
            if (.not.allocated(self%inc))   allocate(self%inc(n))
            if (.not.allocated(self%capom)) allocate(self%capom(n))
            if (.not.allocated(self%omega)) allocate(self%omega(n))
            if (.not.allocated(self%capm))  allocate(self%capm(n))
         end if

         select case(param%in_type)
         case ("ASCII")
            do i = 1, n
               select type(self)
               class is (swiftest_pl)
                  if (param%lrhill_present) then
                     read(iu, *, err = 667, iomsg = errmsg) name(i), val, self%rhill(i)
                  else
                     read(iu, *, err = 667, iomsg = errmsg) name(i), val
                  end if
                  self%Gmass(i) = real(val, kind=DP)
                  self%mass(i) = real(val / param%GU, kind=DP)
                  if (param%lclose) read(iu, *, err = 667, iomsg = errmsg) self%radius(i)
               class is (swiftest_tp)
                  read(iu, *, err = 667, iomsg = errmsg) name(i)
               end select
               call self%info(i)%set_value(name=name(i))

               select case(param%in_form)
               case ("XV")
                  read(iu, *, err = 667, iomsg = errmsg) self%rh(1, i), self%rh(2, i), self%rh(3, i)
                  read(iu, *, err = 667, iomsg = errmsg) self%vh(1, i), self%vh(2, i), self%vh(3, i)
               case ("EL")
                  read(iu, *, err = 667, iomsg = errmsg) self%a(i), self%e(i), self%inc(i)
                  read(iu, *, err = 667, iomsg = errmsg) self%capom(i), self%omega(i), self%capm(i)
               end select

               select type (self)
               class is (swiftest_pl)
                  if (param%lrotation) then
                     read(iu, *, err = 667, iomsg = errmsg) self%Ip(1, i), self%Ip(2, i), self%Ip(3, i)
                     read(iu, *, err = 667, iomsg = errmsg) self%rot(1, i), self%rot(2, i), self%rot(3, i)
                  end if
                  if (param%ltides) then
                     read(iu, *, err = 667, iomsg = errmsg) self%k2(i)
                     read(iu, *, err = 667, iomsg = errmsg) self%Q(i)
                  end if
               end select
            end do
         end select
      end associate

      ierr = 0
      return

      667 continue
      select type (self)
      class is (swiftest_pl)
         write(*,*) "Error reading massive body file: " // trim(adjustl(errmsg))
      class is (swiftest_tp)
         write(*,*) "Error reading test particle file: " // trim(adjustl(errmsg))
      class default
         write(*,*) "Error reading body file: " // trim(adjustl(errmsg))
      end select
      call base_util_exit(FAILURE,param%display_unit)
   end function swiftest_io_read_frame_body


   module subroutine swiftest_io_read_in_param(self, param_file_name) 
      !! author: David A. Minton
      !!
      !! Read in parameters for the integration
      !!
      !! Adapted from David E. Kaufmann's Swifter routine io_init_param.f90
      !! Adapted from Martin Duncan's Swift routine io_init_param.f
      implicit none
      ! Arguments
      class(swiftest_parameters),intent(inout) :: self             !! Current run configuration parameters
      character(len=*), intent(in)             :: param_file_name !! Parameter input file name (i.e. param.in)
      ! Internals
      integer(I4B)      :: ierr = 0 !! Input error code
      character(STRMAX) :: errmsg   !! Error message in UDIO procedure

      ! Read in name of parameter file
      self%param_file_name = trim(adjustl(param_file_name))

      !! todo: Currently this procedure does not work in user-defined derived-type input mode 
      !!    as the newline characters are ignored in the input file when compiled in ifort.

      !read(LUN,'(DT)', iostat= ierr, iomsg = errmsg) self
      call self%reader(LUN, iotype= "none", v_list=[""], iostat = ierr, iomsg = errmsg)
      if (ierr == 0) return

      write(self%display_unit,*) "Error reading parameter file: " // trim(adjustl(errmsg))
      call base_util_exit(FAILURE)
   end subroutine swiftest_io_read_in_param


   module subroutine swiftest_io_set_display_param(self, display_style)
      !! author: David A. Minton
      !!
      !! Sets the display style parameters, which cause the following output behavior:
      !! "PROGRESS": A progress bar is displayed on stdout. Standard output is redirected to a log file. 
      !! "CLASSIC" : All output goes to stdout, similar to Swifter. No log file is written.
      !! "COMPACT" : A compact machine-readable output is displayed on stdout. Standard output is redirected to a log file. 
      !! "QUIET"   : No output is displayed on stdout. Standard output is redirected to a log file.
      implicit none
      ! Arguments
      class(swiftest_parameters), intent(inout) :: self            !! Current run configuration parameters
      character(*),               intent(in)    :: display_style   !! Style of the output display 
      ! Internals
      character(STRMAX) :: errmsg
      logical           :: fileExists   

      select case(display_style)
      case ('CLASSIC')
         self%display_unit = OUTPUT_UNIT !! stdout from iso_fortran_env
         self%log_output = .false.
      case ('COMPACT', 'PROGRESS', 'QUIET')
#ifdef COARRAY
         if (self%lcoarray) then
            write(SWIFTEST_LOG_FILE,'("swiftest_coimage",I0.4,".log")') this_image()
         else
            write(SWIFTEST_LOG_FILE,'("swiftest.log")')
         end if 
#endif
         inquire(file=SWIFTEST_LOG_FILE, exist=fileExists)
         if (self%lrestart.and.fileExists) then
            open(unit=SWIFTEST_LOG_OUT, file=SWIFTEST_LOG_FILE, status="OLD", position="APPEND", err = 667, iomsg = errmsg)
         else
            open(unit=SWIFTEST_LOG_OUT, file=SWIFTEST_LOG_FILE, status="REPLACE", err = 667, iomsg = errmsg)
         end if
         self%display_unit = SWIFTEST_LOG_OUT 
         self%log_output = .true.
      case default
         write(*,*) display_style, " is an unknown display style."
         call base_util_exit(USAGE)
      end select

      self%display_style = display_style

      return

      667 continue
      write(*,*) "Error opening swiftest log file: " // trim(adjustl(errmsg))
      call base_util_exit(FAILURE,self%display_unit)
   end subroutine swiftest_io_set_display_param


   module subroutine swiftest_io_toupper(string)
      !! author: David A. Minton
      !!
      !! Convert string to uppercase
      !!
      !! Adapted from David E. Kaufmann's Swifter routine: util_toupper.f90
      implicit none
      ! Arguments
      character(*), intent(inout) :: string !! String to make upper case
      ! Internals
      integer(I4B) :: i, length, idx
   
      length = len(string)
      do i = 1, length
         idx = iachar(string(i:i))
         if ((idx >= lowercase_begin) .and. (idx <= lowercase_end)) then
            idx = idx + uppercase_offset
            string(i:i) = achar(idx)
         end if
      end do
   
      return
   end subroutine swiftest_io_toupper


   module subroutine swiftest_io_initialize_output_file_system(self, nc, param)
      !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
      !!
      !! Write a frame (header plus records for each massive body and active test particle) to output binary file
      !! There is no direct file output from this subroutine
      !!
      !! Adapted from David E. Kaufmann's Swifter routine  io_write_frame.f90
      !! Adapted from Hal Levison's Swift routine io_write_frame.f
      implicit none
      ! Arguments
      class(swiftest_nbody_system),      intent(inout) :: self   !! Swiftest nbody_system object
      class(swiftest_netcdf_parameters), intent(inout) :: nc     !! Parameters used to identify a particular NetCDF dataset
      class(swiftest_parameters),        intent(inout) :: param !! Current run configuration parameters 
      ! Internals

      character(len=2*STRMAX)          :: errmsg
      logical                          :: fileExists

      associate (pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody, lfirst => self%lfirst_io)
         nc%file_name = param%outfile
         if (lfirst) then
            inquire(file=param%outfile, exist=fileExists)
#ifdef COARRAY
            if (this_image() /= 1) param%out_stat = 'APPEND'
#endif
            
            select case(param%out_stat)
            case('APPEND')
               if (.not.fileExists) then
                  errmsg = trim(adjustl(param%outfile)) // " not found! You must specify OUT_STAT = NEW, REPLACE, or UNKNOWN"
                  goto 667
               end if
               call nc%open(param)
            case('NEW')
               if (fileExists) then
                  errmsg = trim(adjustl(param%outfile))// " already exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN"
                  goto 667
               end if
               call nc%initialize(param)
            case('REPLACE', 'UNKNOWN')
               call nc%initialize(param)
            end select

            lfirst = .false.
         end if

      end associate

      return

      667 continue
      write(*,*) "Error writing nbody_system frame: " // trim(adjustl(errmsg))
      call base_util_exit(FAILURE,param%display_unit)
   end subroutine swiftest_io_initialize_output_file_system

end submodule s_swiftest_io
