write_vtu_unstructured Subroutine

   subroutine write_vtu_unstructured(params, flow, mesh, fields, species, energy, transport, step)
      type(case_params_t), intent(in) :: params
      type(flow_mpi_t), intent(inout) :: flow
      type(mesh_t), intent(in) :: mesh
      type(flow_fields_t), intent(in) :: fields
      type(species_fields_t), intent(in) :: species
      type(energy_fields_t), intent(in) :: energy
      type(transport_properties_t), intent(in) :: transport
      integer, intent(in) :: step

      integer :: unit_id
      integer :: p, c, n, m, k, n_owned
      character(len=path_len + 64) :: filename, local_filename
      integer, allocatable :: owned_indices(:)

      if (.not. params%write_vtu) return

      ! Count owned cells and collect their indices
      n_owned = 0
      do c = 1, mesh%ncells
         if (flow%owned(c)) n_owned = n_owned + 1
      end do

      allocate(owned_indices(n_owned))
      n_owned = 0
      do c = 1, mesh%ncells
         if (flow%owned(c)) then
            n_owned = n_owned + 1
            owned_indices(n_owned) = c
         end if
      end do

      ! Every rank writes its own piece file
      write(local_filename,'("flow_",i6.6,"_P",i4.4,".vtu")') step, flow%rank
      filename = trim(params%output_dir)//'/'//trim(local_filename)
      open(newunit=unit_id, file=trim(filename), status='replace', action='write')

      write(unit_id,'(a)') '<?xml version="1.0"?>'
      write(unit_id,'(a)') '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">'
      write(unit_id,'(a)') '  <UnstructuredGrid>'
      write(unit_id,'(a,i0,a,i0,a)') '    <Piece NumberOfPoints="', mesh%npoints, &
                                      '" NumberOfCells="', n_owned, '">'

      ! ------------------------------------------------------------
      ! Correct XML VTK Piece order:
      !
      ! PointData
      ! CellData
      ! Points
      ! Cells
      ! ------------------------------------------------------------

      write(unit_id,'(a)') '      <PointData>'
      write(unit_id,'(a)') '      </PointData>'

      write(unit_id,'(a)') '      <CellData Scalars="pressure" Vectors="velocity">'

      write(unit_id,'(a)') '        <DataArray type="Float64" Name="velocity" NumberOfComponents="3" format="ascii">'
      do n = 1, n_owned
         c = owned_indices(n)
         write(unit_id,'(3(es24.16,1x))') fields%u(1,c), fields%u(2,c), fields%u(3,c)
      end do
      write(unit_id,'(a)') '        </DataArray>'

      write(unit_id,'(a,a,a)') '        <DataArray type="Float64" Name="pressure" format="ascii">'
      do n = 1, n_owned
         c = owned_indices(n)
         write(unit_id,'(es24.16)') fields%p(c)
      end do
      write(unit_id,'(a)') '        </DataArray>'


      write(unit_id,'(a)') '        <DataArray type="Float64" Name="thermo_pressure" format="ascii">'
      do n = 1, n_owned
         write(unit_id,'(es24.16)') params%background_press
      end do
      write(unit_id,'(a)') '        </DataArray>'

      write(unit_id,'(a,a,a)') '        <DataArray type="Float64" Name="divergence" format="ascii">'
      do n = 1, n_owned
         c = owned_indices(n)
         write(unit_id,'(es24.16)') fields%div(c)
      end do
      write(unit_id,'(a)') '        </DataArray>'

      if (params%enable_energy) then
         if (.not. allocated(energy%T) .or. .not. allocated(energy%h) .or. &
             .not. allocated(energy%qrad)) then
            call fatal_error('output', 'energy output requested but energy arrays are not allocated')
         end if

         write(unit_id,'(a)') '        <DataArray type="Float64" Name="temperature" format="ascii">'
         do n = 1, n_owned
            c = owned_indices(n)
            write(unit_id,'(es24.16)') energy%T(c)
         end do
         write(unit_id,'(a)') '        </DataArray>'

         write(unit_id,'(a)') '        <DataArray type="Float64" Name="enthalpy" format="ascii">'
         do n = 1, n_owned
            c = owned_indices(n)
            write(unit_id,'(es24.16)') energy%h(c)
         end do
         write(unit_id,'(a)') '        </DataArray>'

         write(unit_id,'(a)') '        <DataArray type="Float64" Name="qrad" format="ascii">'
         do n = 1, n_owned
            c = owned_indices(n)
            write(unit_id,'(es24.16)') energy%qrad(c)
         end do
         write(unit_id,'(a)') '        </DataArray>'

         if (allocated(energy%cp)) then
            write(unit_id,'(a)') '        <DataArray type="Float64" Name="cp" format="ascii">'
            do n = 1, n_owned
               c = owned_indices(n)
               write(unit_id,'(es24.16)') energy%cp(c)
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end if

         if (allocated(energy%lambda)) then
            write(unit_id,'(a)') '        <DataArray type="Float64" Name="thermal_conductivity" format="ascii">'
            do n = 1, n_owned
               c = owned_indices(n)
               write(unit_id,'(es24.16)') energy%lambda(c)
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end if


         if (allocated(energy%lambda) .and. allocated(energy%cp)) then
            write(unit_id,'(a)') '        <DataArray type="Float64" Name="thermal_diffusivity" format="ascii">'
            do n = 1, n_owned
               c = owned_indices(n)
               write(unit_id,'(es24.16)') energy%lambda(c) / max(params%rho * energy%cp(c), tiny(1.0_rk))
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end if

         if (allocated(energy%rho_thermo)) then
            write(unit_id,'(a)') '        <DataArray type="Float64" Name="rho_thermo" format="ascii">'
            do n = 1, n_owned
               c = owned_indices(n)
               write(unit_id,'(es24.16)') energy%rho_thermo(c)
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end if

      end if

      if (params%enable_species .and. params%nspecies > 0) then

         do k = 1, params%nspecies
            write(unit_id,'(a,a,a)') '        <DataArray type="Float64" Name="Y_', &
               trim(params%species_name(k)), '" format="ascii">'
            do n = 1, n_owned
               c = owned_indices(n)
               write(unit_id,'(es24.16)') species%Y(k,c)
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end do

         write(unit_id,'(a)') '        <DataArray type="Float64" Name="sum_Y" format="ascii">'
         do n = 1, n_owned
            c = owned_indices(n)
            write(unit_id,'(es24.16)') sum(species%Y(:,c))
         end do
         write(unit_id,'(a)') '        </DataArray>'

         if (allocated(transport%diffusivity)) then
            do k = 1, params%nspecies
               write(unit_id,'(a,a,a)') '        <DataArray type="Float64" Name="D_', &
                  trim(params%species_name(k)), '" format="ascii">'
               do n = 1, n_owned
                  c = owned_indices(n)
                  write(unit_id,'(es24.16)') transport%diffusivity(k,c)
               end do
               write(unit_id,'(a)') '        </DataArray>'
            end do
         end if

      end if

      ! Helpful debug scalar so you can color by cell_id in ParaView.
      write(unit_id,'(a)') '        <DataArray type="Int32" Name="cell_id" format="ascii">'
      do n = 1, n_owned
         c = owned_indices(n)
         write(unit_id,'(i0)') c
      end do
      write(unit_id,'(a)') '        </DataArray>'

      write(unit_id,'(a)') '      </CellData>'

      write(unit_id,'(a)') '      <Points>'
      write(unit_id,'(a)') '        <DataArray type="Float64" NumberOfComponents="3" format="ascii">'
      do p = 1, mesh%npoints
         write(unit_id,'(3(es24.16,1x))') mesh%points(1,p), mesh%points(2,p), mesh%points(3,p)
      end do
      write(unit_id,'(a)') '        </DataArray>'
      write(unit_id,'(a)') '      </Points>'

      write(unit_id,'(a)') '      <Cells>'

      write(unit_id,'(a)') '        <DataArray type="Int32" Name="connectivity" format="ascii">'
      do n = 1, n_owned
         c = owned_indices(n)
         write(unit_id,'(8(i0,1x))') (mesh%cells(c)%nodes(m) - 1, m = 1, 8)
      end do
      write(unit_id,'(a)') '        </DataArray>'

      write(unit_id,'(a)') '        <DataArray type="Int32" Name="offsets" format="ascii">'
      do n = 1, n_owned
         write(unit_id,'(i0)') 8 * n
      end do
      write(unit_id,'(a)') '        </DataArray>'

      write(unit_id,'(a)') '        <DataArray type="UInt8" Name="types" format="ascii">'
      do n = 1, n_owned
         write(unit_id,'(i0)') 12
      end do
      write(unit_id,'(a)') '        </DataArray>'

      write(unit_id,'(a)') '      </Cells>'

      write(unit_id,'(a)') '    </Piece>'
      write(unit_id,'(a)') '  </UnstructuredGrid>'
      write(unit_id,'(a)') '</VTKFile>'

      close(unit_id)

      deallocate(owned_indices)

      ! Rank 0 writes the master PVTU file
      if (flow%rank == 0) then
         call write_pvtu_master(params, flow, step)
      end if
   end subroutine write_vtu_unstructured