read_restart_file Subroutine

   subroutine read_restart_file(mesh, flow, params, fields, species, energy, transport, step, time, chemistry_accum_dt)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      type(case_params_t), intent(inout) :: params
      type(flow_fields_t), intent(inout) :: fields
      type(species_fields_t), intent(inout) :: species
      type(energy_fields_t), intent(inout) :: energy
      type(transport_properties_t), intent(inout) :: transport
      integer, intent(inout) :: step
      real(rk), intent(inout) :: time
      real(rk), intent(inout), optional :: chemistry_accum_dt

      character(len=32) :: magic
      integer :: unit_id, ios, version
      integer :: file_ncells, file_nfaces, file_nspecies, file_step, writer_nprocs
      integer :: has_species, has_energy, has_variable_density, rhs_old_valid_int
      integer :: species_rhs_valid_int, energy_rhs_valid_int
      real(rk) :: file_time, file_dt, file_chemistry_accum_dt
      real(rk), allocatable :: tmp_matrix(:,:)
      real(rk), allocatable :: energy_diag(:)
      logical :: load_species, load_energy

      if (.not. params%restart_from_file) return
      if (len_trim(params%restart_file) == 0) then
         call fatal_error('restart', 'restart_from_file=.true. requires restart_file')
      end if

      open(newunit=unit_id, file=trim(params%restart_file), access='stream', form='unformatted', &
           status='old', action='read', iostat=ios)
      if (ios /= 0) call fatal_error('restart', 'could not open restart_file: '//trim(params%restart_file))

      read(unit_id, iostat=ios) magic
      if (ios /= 0) call fatal_error('restart', 'failed reading restart magic')
      if (trim(magic) /= trim(restart_magic)) then
         call fatal_error('restart', 'restart file magic mismatch; expected '//trim(restart_magic))
      end if

      read(unit_id) version
      if (version /= restart_version) call fatal_error('restart', 'unsupported restart version')
      read(unit_id) file_ncells, file_nfaces, file_nspecies, file_step, writer_nprocs
      read(unit_id) file_time, file_dt, file_chemistry_accum_dt
      read(unit_id) has_species, has_energy, has_variable_density, rhs_old_valid_int, &
                    species_rhs_valid_int, energy_rhs_valid_int

      if (file_ncells /= mesh%ncells) call fatal_error('restart', 'restart ncells does not match current mesh')
      if (file_nfaces /= mesh%nfaces) call fatal_error('restart', 'restart nfaces does not match current mesh')
      if (params%enable_species .and. file_nspecies > params%nspecies) then
         call fatal_error('restart', 'restart nspecies exceeds current case nspecies; cannot load restart')
      end if
      if (params%enable_species .and. flow%rank == 0 .and. file_nspecies /= params%nspecies) then
         write(output_unit,'(a,i0,a,i0,a)') &
            'restart: WARNING: file nspecies=', file_nspecies, &
            ' < current nspecies=', params%nspecies, &
            '; remaining species initialized to zero'
      end if
      if (has_energy == 1 .and. .not. params%enable_energy) then
         call fatal_error('restart', 'restart contains energy state but current case has enable_energy=.false.')
      end if
      if (has_species == 1 .and. .not. params%enable_species) then
         call fatal_error('restart', 'restart contains species state but current case has enable_species=.false.')
      end if
      if (has_variable_density == 1 .and. .not. params%enable_variable_density) then
         call fatal_error('restart', 'restart was written from variable-density mode but current case is not variable-density')
      end if

      read(unit_id) fields%u
      read(unit_id) fields%u_old
      read(unit_id) fields%p
      read(unit_id) fields%phi
      read(unit_id) fields%divergence_source
      read(unit_id) fields%projection_rho
      read(unit_id) fields%projection_divergence_source
      read(unit_id) fields%rhs_old
      read(unit_id) fields%face_flux
      read(unit_id) fields%mass_flux

      read(unit_id) transport%rho
      read(unit_id) transport%rho_old
      read(unit_id) transport%mu
      read(unit_id) transport%nu
      read(unit_id) transport%lambda

      if (file_nspecies > 0) then
         allocate(tmp_matrix(file_nspecies, mesh%ncells))
         read(unit_id) tmp_matrix
         if (allocated(transport%diffusivity)) then
            transport%diffusivity(1:min(file_nspecies, params%nspecies), :) = &
               tmp_matrix(1:min(file_nspecies, params%nspecies), :)
         end if
         deallocate(tmp_matrix)
      end if

      load_species = has_species == 1 .and. file_nspecies > 0
      if (load_species) then
         allocate(tmp_matrix(file_nspecies, mesh%ncells))
         read(unit_id) tmp_matrix
         species%Y(1:min(file_nspecies, params%nspecies), :) = &
            tmp_matrix(1:min(file_nspecies, params%nspecies), :)
         if (allocated(species%Y_old)) species%Y_old = species%Y
         read(unit_id) tmp_matrix
         if (allocated(species%rhs_old)) then
            species%rhs_old(1:min(file_nspecies, params%nspecies), :) = &
               tmp_matrix(1:min(file_nspecies, params%nspecies), :)
         end if
         species%rhs_old_valid = species_rhs_valid_int /= 0
         deallocate(tmp_matrix)
      else
         species%rhs_old_valid = .false.
      end if

      load_energy = has_energy == 1
      if (load_energy) then
         read(unit_id) energy%T
         read(unit_id) energy%T_old
         read(unit_id) energy%h
         read(unit_id) energy%h_old
         read(unit_id) energy%rhs_old
         energy%rhs_old_valid = energy_rhs_valid_int /= 0
         read(unit_id) energy%qrad
         read(unit_id) energy%cp
         read(unit_id) energy%lambda
         read(unit_id) energy%rho_thermo
         allocate(energy_diag(12))
         read(unit_id) energy_diag
         call restore_energy_diagnostics(energy, energy_diag)
         deallocate(energy_diag)
      end if

      close(unit_id)

      fields%rhs_old_valid = rhs_old_valid_int /= 0
      step = file_step
      time = file_time
      if (params%restart_use_file_dt) params%dt = file_dt
      if (params%restart_nsteps_are_additional) params%nsteps = file_step + params%nsteps
      if (present(chemistry_accum_dt)) chemistry_accum_dt = file_chemistry_accum_dt

      if (flow%rank == 0) then
         write(output_unit,'(a,a,a,i0,a,es13.5,a,es13.5,a,es13.5,a,i0,a,i0,a,l1,a,l1)') &
            'read restart: ', trim(params%restart_file), ' step=', step, &
            ' time=', time, ' file_dt=', file_dt, ' run_dt=', params%dt, &
            ' nsteps=', params%nsteps, ' writer_ranks=', writer_nprocs, &
            ' use_file_dt=', params%restart_use_file_dt, &
            ' nsteps_additional=', params%restart_nsteps_are_additional
      end if
   end subroutine read_restart_file