write_variable_density_projection_audit Subroutine

   subroutine write_variable_density_projection_audit(mesh, flow, params, fields, energy, transport, step, time)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(case_params_t), intent(in) :: params
      type(flow_fields_t), intent(in) :: fields
      type(energy_fields_t), intent(in) :: energy
      type(transport_properties_t), intent(in) :: transport
      integer, intent(in) :: step
      real(rk), intent(in) :: time

      integer :: c, unit_cells, unit_faces
      integer :: local_worst_cell, boundary_face_count, physical_boundary_face_count
      logical, save :: projection_audit_initialized = .false.
      character(len=32) :: rank_suffix
      character(len=1024) :: cells_file, faces_file
      real(rk) :: local_worst_abs
      real(rk) :: s_c, div_c, div_res, mdot_div_c

      if (.not. params%enable_variable_density) return
      if (.not. params%write_diagnostics) return
      if (.not. allocated(transport%rho)) return

      local_worst_abs = -1.0_rk
      local_worst_cell = -1

      do c = flow%first_cell, flow%last_cell
         call compute_projection_audit_cell_balance(mesh, fields, c, div_c, mdot_div_c, &
                                                    boundary_face_count, physical_boundary_face_count)

         s_c = 0.0_rk
         if (allocated(fields%divergence_source)) s_c = fields%divergence_source(c)

         div_res = div_c - s_c
         if (abs(div_res) > local_worst_abs) then
            local_worst_abs = abs(div_res)
            local_worst_cell = c
         end if
      end do

      write(rank_suffix, '(i0)') flow%rank
      cells_file = trim(params%output_dir) // '/diagnostics/variable_density_projection_audit_cells_rank' // &
                   trim(adjustl(rank_suffix)) // '.csv'
      faces_file = trim(params%output_dir) // '/diagnostics/variable_density_projection_audit_faces_rank' // &
                   trim(adjustl(rank_suffix)) // '.csv'

      if (.not. projection_audit_initialized) then
         open(newunit=unit_cells, file=trim(cells_file), status='replace', action='write')
         write(unit_cells,'(a)') 'step,time,rank,audit_label,cell,volume,rho,T,h,S,divu,divu_minus_S,' // &
                                 'mass_flux_div,required_flux_sum,actual_flux_sum,flux_defect,' // &
                                 'mass_flux_sum,boundary_face_count,physical_boundary_face_count,pressure'
         close(unit_cells)

         open(newunit=unit_faces, file=trim(faces_file), status='replace', action='write')
         write(unit_faces,'(a)') 'step,time,rank,audit_label,cell,local_face,face_id,owner,neighbor,periodic_neighbor,' // &
                                 'is_boundary_face,is_physical_boundary_face,owner_pressure,cell_pressure,neighbor_pressure,' // &
                                 'dp_neighbor_minus_cell,face_flux,outward_face_flux,mass_flux,outward_mass_flux'
         close(unit_faces)
         projection_audit_initialized = .true.
      end if

      open(newunit=unit_cells, file=trim(cells_file), status='unknown', position='append', action='write')
      open(newunit=unit_faces, file=trim(faces_file), status='unknown', position='append', action='write')

      if (local_worst_cell > 0) then
         call write_projection_audit_one_cell(mesh, flow, fields, energy, transport, step, time, &
                                              local_worst_cell, 'local_worst', unit_cells, unit_faces)
      end if

      if (flow%first_cell <= 1 .and. 1 <= flow%last_cell) then
         if (local_worst_cell /= 1) then
            call write_projection_audit_one_cell(mesh, flow, fields, energy, transport, step, time, &
                                                 1, 'cell_1', unit_cells, unit_faces)
         end if
      end if

      close(unit_cells)
      close(unit_faces)

   end subroutine write_variable_density_projection_audit