write_vtu_ascii Subroutine

   subroutine write_vtu_ascii(params, flow, mesh, fields, species, energy, transport, step, time)
      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
      real(rk), intent(in) :: time

      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)//'/VTK/'//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(ES26.16E4,1x))') fields%u(1,c), fields%u(2,c), fields%u(3,c)
      end do
      write(unit_id,'(a)') '        </DataArray>'

      write(unit_id,'(a)') '        <DataArray type="Float64" Name="U_mag" format="ascii">'
      do n = 1, n_owned
         c = owned_indices(n)
         write(unit_id,'(ES26.16E4)') sqrt(dot_product(fields%u(:,c), fields%u(:,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,'(ES26.16E4)') 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,'(ES26.16E4)') 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,'(ES26.16E4)') fields%div(c)
      end do
      write(unit_id,'(a)') '        </DataArray>'

      if (.not. allocated(transport%rho) .or. .not. allocated(transport%nu)) then
         call fatal_error('output', 'transport rho/nu output requested but transport arrays are not allocated')
      end if

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

      write(unit_id,'(a)') '        <DataArray type="Float64" Name="nu" format="ascii">'
      do n = 1, n_owned
         c = owned_indices(n)
         write(unit_id,'(ES26.16E4)') transport%nu(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,'(ES26.16E4)') 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,'(ES26.16E4)') 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,'(ES26.16E4)') energy%qrad(c)
         end do
         write(unit_id,'(a)') '        </DataArray>'


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

         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,'(ES26.16E4)') 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,'(ES26.16E4)') 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,'(ES26.16E4)') 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,'(ES26.16E4)') energy%rho_thermo(c)
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end if

         if (allocated(energy%chem_heat_release_rate)) then
            write(unit_id,'(a)') '        <DataArray type="Float64" Name="chem_heat_release_rate" format="ascii">'
            do n = 1, n_owned
               c = owned_indices(n)
               write(unit_id,'(ES26.16E4)') energy%chem_heat_release_rate(c)
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end if
         if (allocated(energy%chem_dT_dt)) then
            write(unit_id,'(a)') '        <DataArray type="Float64" Name="chem_dT_dt" format="ascii">'
            do n = 1, n_owned
               c = owned_indices(n)
               write(unit_id,'(ES26.16E4)') energy%chem_dT_dt(c)
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end if
         if (allocated(energy%chem_dY_max_dt)) then
            write(unit_id,'(a)') '        <DataArray type="Float64" Name="chem_dY_max_dt" format="ascii">'
            do n = 1, n_owned
               c = owned_indices(n)
               write(unit_id,'(ES26.16E4)') energy%chem_dY_max_dt(c)
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end if
         if (allocated(energy%chem_rel_rho_change_dt)) then
            write(unit_id,'(a)') '        <DataArray type="Float64" Name="chem_rel_rho_change_dt" format="ascii">'
            do n = 1, n_owned
               c = owned_indices(n)
               write(unit_id,'(ES26.16E4)') energy%chem_rel_rho_change_dt(c)
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end if
         if (allocated(energy%chem_source_alpha)) then
            write(unit_id,'(a)') '        <DataArray type="Float64" Name="chem_source_alpha" format="ascii">'
            do n = 1, n_owned
               c = owned_indices(n)
               write(unit_id,'(ES26.16E4)') energy%chem_source_alpha(c)
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end if
         if (allocated(energy%chem_active)) then
            write(unit_id,'(a)') '        <DataArray type="Float64" Name="chem_active" format="ascii">'
            do n = 1, n_owned
               c = owned_indices(n)
               write(unit_id,'(ES26.16E4)') energy%chem_active(c)
            end do
            write(unit_id,'(a)') '        </DataArray>'
         end if

         if (params%enable_variable_density) then
            call write_energy_reconciliation_vtu_arrays(unit_id, mesh, flow, energy, transport)
         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,'(ES26.16E4)') 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,'(ES26.16E4)') 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,'(ES26.16E4)') 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>'


      ! `fields%mass_flux` is face-centered and owner-to-neighbor oriented.
      ! The volume VTU exports cell-centered rho*u and div(rho*u) diagnostics
      ! so mass-flow information appears in ParaView/PVTU.

      block

         integer :: mf_c, mf_f, mf_owner, mf_nb

         real(rk), allocatable :: mf_div(:)

         allocate(mf_div(mesh%ncells))

         mf_div = 0.0_rk

      

         if (allocated(fields%mass_flux)) then

            do mf_f = 1, mesh%nfaces

               mf_owner = mesh%faces(mf_f)%owner

               mf_nb = mesh%faces(mf_f)%neighbor

               if (mf_nb <= 0) mf_nb = mesh%faces(mf_f)%periodic_neighbor

      

               if (mf_owner > 0) then

                  mf_div(mf_owner) = mf_div(mf_owner) + fields%mass_flux(mf_f) / mesh%cells(mf_owner)%volume

               end if

               if (mf_nb > 0) then

                  mf_div(mf_nb) = mf_div(mf_nb) - fields%mass_flux(mf_f) / mesh%cells(mf_nb)%volume

               end if

            end do

         end if

      

         write(unit_id,'(a)') '        <DataArray type="Float64" Name="mass_flux_vector" NumberOfComponents="3" format="ascii">'

         do mf_c = 1, mesh%ncells

            if (.not. flow%owned(mf_c)) cycle

            write(unit_id,'(3(1x,ES26.16E4))') &

               transport%rho(mf_c) * fields%u(1,mf_c), &

               transport%rho(mf_c) * fields%u(2,mf_c), &

               transport%rho(mf_c) * fields%u(3,mf_c)

         end do

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

      

         write(unit_id,'(a)') '        <DataArray type="Float64" Name="mass_flux_divergence" format="ascii">'

         do mf_c = 1, mesh%ncells

            if (.not. flow%owned(mf_c)) cycle

            write(unit_id,'(1x,ES26.16E4)') mf_div(mf_c)

         end do

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

      

         deallocate(mf_div)

      end block


      ! Current target low-Mach divergence source for div(u)=S projection.

      write(unit_id,'(a)') '        <DataArray type="Float64" Name="lowmach_divergence_source" format="ascii">'

      do c = 1, mesh%ncells

         if (.not. flow%owned(c)) cycle

         if (allocated(fields%divergence_source)) then

            write(unit_id,'(1x,ES26.16E4)') fields%divergence_source(c)

         else

            write(unit_id,'(1x,ES26.16E4)') 0.0_rk

         end if

      end do

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

      if (params%enable_variable_density) then

         call write_lowmach_debug_vtu_arrays(unit_id, mesh, flow, params, fields, transport)

      end if

      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(ES26.16E4,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)

    end subroutine write_vtu_ascii