write_pvtu_master Subroutine

   subroutine write_pvtu_master(params, flow, step)
      type(case_params_t), intent(in) :: params
      type(flow_mpi_t), intent(in) :: flow
      integer, intent(in) :: step

      integer :: unit_id, r, k
      character(len=path_len + 64) :: filename
      character(len=64) :: piece_name

      write(filename,'(a,"/VTK/flow_",i6.6,".pvtu")') trim(params%output_dir), step
      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="PUnstructuredGrid" version="0.1" byte_order="LittleEndian">'
      write(unit_id,'(a)') '  <PUnstructuredGrid GhostLevel="0">'

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

      write(unit_id,'(a)') '    <PCellData Scalars="pressure" Vectors="velocity">'
      write(unit_id,'(a)') '      <PDataArray type="Float64" Name="velocity" NumberOfComponents="3" format="ascii"/>'
      write(unit_id,'(a)') '      <PDataArray type="Float64" Name="U_mag" format="ascii"/>'
      write(unit_id,'(a)') '      <PDataArray type="Float64" Name="pressure" format="ascii"/>'
      write(unit_id,'(a)') '      <PDataArray type="Float64" Name="thermo_pressure" format="ascii"/>'
      write(unit_id,'(a)') '      <PDataArray type="Float64" Name="divergence" format="ascii"/>'
      write(unit_id,'(a)') '      <PDataArray type="Float64" Name="rho" format="ascii"/>'
      write(unit_id,'(a)') '      <PDataArray type="Float64" Name="nu" format="ascii"/>'
      if (params%enable_energy) then
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="temperature" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="enthalpy" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="qrad" format="ascii"/>'
         if (params%enable_species_enthalpy_diffusion) then
            write(unit_id,'(a)') '      <PDataArray type="Float64" Name="species_enthalpy_diffusion" format="ascii"/>'
         end if
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="cp" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="thermal_conductivity" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="thermal_diffusivity" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="rho_thermo" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="chem_heat_release_rate" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="chem_dT_dt" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="chem_dY_max_dt" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="chem_rel_rho_change_dt" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="chem_source_alpha" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="chem_active" format="ascii"/>'
         if (params%enable_variable_density) then
            write(unit_id,'(a)') '      <PDataArray type="Float64" Name="rho_h_output_state" format="ascii"/>'
            write(unit_id,'(a)') '      <PDataArray type="Float64" Name="rho_h_operator_consistent" format="ascii"/>'
            write(unit_id,'(a)') '      <PDataArray type="Float64" Name="rho_h_density_reconciliation" format="ascii"/>'
            write(unit_id,'(a)') '      <PDataArray type="Float64" Name="relative_rho_h_density_reconciliation" format="ascii"/>'
         end if
      end if
      if (params%enable_species .and. params%nspecies > 0) then
         do k = 1, params%nspecies
            write(unit_id,'(a,a,a)') '      <PDataArray type="Float64" Name="Y_', &
               trim(params%species_name(k)), '" format="ascii"/>'
         end do
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="sum_Y" format="ascii"/>'
         do k = 1, params%nspecies
            write(unit_id,'(a,a,a)') '      <PDataArray type="Float64" Name="D_', &
               trim(params%species_name(k)), '" format="ascii"/>'
         end do
      end if
      write(unit_id,'(a)') '      <PDataArray type="Int32" Name="cell_id" format="ascii"/>'
      write(unit_id,'(a)') '        <PDataArray type="Float64" Name="mass_flux_vector" NumberOfComponents="3"/>'
      write(unit_id,'(a)') '        <PDataArray type="Float64" Name="mass_flux_divergence"/>'
      write(unit_id,'(a)') '        <PDataArray type="Float64" Name="lowmach_divergence_source"/>'

      if (params%enable_variable_density) then
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="lowmach_source_current" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="lowmach_source_projection" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="lowmach_source_difference" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="divu_recomputed" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="divu_minus_S_projection" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="divu_minus_S_current" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="rho_current" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="rho_projection" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="rho_current_minus_projection" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="mass_flux_divergence_recomputed" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="lowmach_source_history_estimate" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="lowmach_source_advective_density" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="u_dot_grad_rho" format="ascii"/>'
         write(unit_id,'(a)') '      <PDataArray type="Float64" Name="continuity_residual_estimate" format="ascii"/>'
      end if
      write(unit_id,'(a)') '    </PCellData>'

      write(unit_id,'(a)') '    <PPoints>'
      write(unit_id,'(a)') '      <PDataArray type="Float64" NumberOfComponents="3" format="ascii"/>'
      write(unit_id,'(a)') '    </PPoints>'

      do r = 0, flow%nprocs - 1
         write(piece_name,'("flow_",i6.6,"_P",i4.4,".vtu")') step, r
         write(unit_id,'(a,a,a)') '    <Piece Source="', trim(piece_name), '"/>'
      end do

      write(unit_id,'(a)') '  </PUnstructuredGrid>'
      write(unit_id,'(a)') '</VTKFile>'
      close(unit_id)
   end subroutine write_pvtu_master