write_variable_density_diagnostics Subroutine

   subroutine write_variable_density_diagnostics(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, f, lf, fid, owner, nb, ierr, unit_id
      integer :: local_worst_cell, local_winner_rank, global_winner_rank, global_worst_cell
      integer :: mpi_int_send(1), mpi_int_recv(1)
      logical :: file_exists
      logical, save :: vd_diag_initialized = .false.
      logical, save :: vd_worst_initialized = .false.
      character(len=1024) :: filename
      real(rk) :: rho_c, rho_old_c, s_c, div_c, div_res, vol
      real(rk) :: flux_out, mdot_out, mdot_div_c, face_area_value, mass_flux_value
      real(rk) :: local_rho_min, local_rho_max
      real(rk) :: local_rho_old_min, local_rho_old_max
      real(rk) :: local_S_min, local_S_max
      real(rk) :: local_div_min, local_div_max
      real(rk) :: local_res_max, local_res_l2_num
      real(rk) :: local_mdot_div_min, local_mdot_div_max, local_mdot_div_l2_num
      real(rk) :: local_mass, local_net_mdot, local_volume
      real(rk) :: local_h_min, local_h_max
      real(rk) :: local_T_min, local_T_max
      real(rk) :: local_worst_abs, local_worst_divu, local_worst_S
      real(rk) :: local_worst_mdot_div, local_worst_rho, local_worst_T, local_worst_h
      real(rk) :: global_worst_abs, global_worst_divu, global_worst_S
      real(rk) :: global_worst_mdot_div, global_worst_rho, global_worst_T, global_worst_h
      real(rk) :: rho_min, rho_max
      real(rk) :: rho_old_min, rho_old_max
      real(rk) :: S_min, S_max
      real(rk) :: div_min, div_max
      real(rk) :: res_max, res_l2_num, res_l2
      real(rk) :: mdot_div_min, mdot_div_max, mdot_div_l2_num, mdot_div_l2
      real(rk) :: mass_integral, net_boundary_mdot, volume_total
      real(rk) :: h_min, h_max
      real(rk) :: T_min, T_max
      real(rk) :: mpi_reduce_send(1), mpi_reduce_recv(1)

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

      local_rho_min = huge(1.0_rk)
      local_rho_max = -huge(1.0_rk)
      local_rho_old_min = huge(1.0_rk)
      local_rho_old_max = -huge(1.0_rk)
      local_S_min = huge(1.0_rk)
      local_S_max = -huge(1.0_rk)
      local_div_min = huge(1.0_rk)
      local_div_max = -huge(1.0_rk)
      local_res_max = 0.0_rk
      local_res_l2_num = 0.0_rk
      local_mdot_div_min = huge(1.0_rk)
      local_mdot_div_max = -huge(1.0_rk)
      local_mdot_div_l2_num = 0.0_rk
      local_mass = 0.0_rk
      local_net_mdot = 0.0_rk
      local_volume = 0.0_rk
      local_h_min = huge(1.0_rk)
      local_h_max = -huge(1.0_rk)
      local_T_min = huge(1.0_rk)
      local_T_max = -huge(1.0_rk)
      local_worst_abs = -1.0_rk
      local_worst_cell = -1
      local_worst_divu = 0.0_rk
      local_worst_S = 0.0_rk
      local_worst_mdot_div = 0.0_rk
      local_worst_rho = 0.0_rk
      local_worst_T = 0.0_rk
      local_worst_h = 0.0_rk

      do c = flow%first_cell, flow%last_cell
         vol = mesh%cells(c)%volume
         rho_c = transport%rho(c)
         rho_old_c = rho_c
         if (allocated(transport%rho_old)) rho_old_c = transport%rho_old(c)

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

         div_c = 0.0_rk
         mdot_div_c = 0.0_rk

         if (allocated(fields%face_flux)) then
            do lf = 1, mesh%ncell_faces(c)
               fid = mesh%cell_faces(lf, c)
               if (mesh%faces(fid)%owner == c) then
                  flux_out = fields%face_flux(fid)
               else
                  flux_out = -fields%face_flux(fid)
               end if
               div_c = div_c + flux_out / vol
            end do
         end if

         if (allocated(fields%mass_flux)) then
            do lf = 1, mesh%ncell_faces(c)
               fid = mesh%cell_faces(lf, c)
               if (mesh%faces(fid)%owner == c) then
                  mdot_out = fields%mass_flux(fid)
               else
                  mdot_out = -fields%mass_flux(fid)
               end if
               mdot_div_c = mdot_div_c + mdot_out / vol
            end do
         end if

         div_res = div_c - s_c

         local_rho_min = min(local_rho_min, rho_c)
         local_rho_max = max(local_rho_max, rho_c)
         local_rho_old_min = min(local_rho_old_min, rho_old_c)
         local_rho_old_max = max(local_rho_old_max, rho_old_c)
         local_S_min = min(local_S_min, s_c)
         local_S_max = max(local_S_max, s_c)
         local_div_min = min(local_div_min, div_c)
         local_div_max = max(local_div_max, div_c)
         local_res_max = max(local_res_max, abs(div_res))
         local_res_l2_num = local_res_l2_num + div_res * div_res * vol
         local_mdot_div_min = min(local_mdot_div_min, mdot_div_c)
         local_mdot_div_max = max(local_mdot_div_max, mdot_div_c)
         local_mdot_div_l2_num = local_mdot_div_l2_num + mdot_div_c * mdot_div_c * vol
         local_mass = local_mass + rho_c * vol
         local_volume = local_volume + vol

         if (allocated(energy%h)) then
            local_h_min = min(local_h_min, energy%h(c))
            local_h_max = max(local_h_max, energy%h(c))
         end if
         if (allocated(energy%T)) then
            local_T_min = min(local_T_min, energy%T(c))
            local_T_max = max(local_T_max, energy%T(c))
         end if

         if (abs(div_res) > local_worst_abs) then
            local_worst_abs = abs(div_res)
            local_worst_cell = c
            local_worst_divu = div_c
            local_worst_S = s_c
            local_worst_mdot_div = mdot_div_c
            local_worst_rho = rho_c
            if (allocated(energy%T)) local_worst_T = energy%T(c)
            if (allocated(energy%h)) local_worst_h = energy%h(c)
         end if
      end do

      if (.not. allocated(fields%mass_flux)) then
         local_mdot_div_min = 0.0_rk
         local_mdot_div_max = 0.0_rk
         local_mdot_div_l2_num = 0.0_rk
      end if
      if (.not. allocated(energy%h)) then
         local_h_min = 0.0_rk
         local_h_max = 0.0_rk
      end if
      if (.not. allocated(energy%T)) then
         local_T_min = 0.0_rk
         local_T_max = 0.0_rk
      end if

      if (allocated(fields%mass_flux)) then
         do f = 1, mesh%nfaces
            owner = mesh%faces(f)%owner
            nb = mesh%faces(f)%neighbor
            if (nb <= 0) nb = mesh%faces(f)%periodic_neighbor

            if (nb <= 0) then
               if (owner >= flow%first_cell .and. owner <= flow%last_cell) then
                  local_net_mdot = local_net_mdot + fields%mass_flux(f)
               end if
            end if
         end do
      end if

      mpi_reduce_send(1) = local_rho_min
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      rho_min = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd rho_min')

      mpi_reduce_send(1) = local_rho_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      rho_max = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd rho_max')

      mpi_reduce_send(1) = local_rho_old_min
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      rho_old_min = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd rho_old_min')

      mpi_reduce_send(1) = local_rho_old_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      rho_old_max = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd rho_old_max')

      mpi_reduce_send(1) = local_S_min
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      S_min = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd S_min')

      mpi_reduce_send(1) = local_S_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      S_max = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd S_max')

      mpi_reduce_send(1) = local_div_min
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      div_min = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd div_min')

      mpi_reduce_send(1) = local_div_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      div_max = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd div_max')

      mpi_reduce_send(1) = local_res_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      res_max = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd residual max')

      mpi_reduce_send(1) = local_res_l2_num
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      res_l2_num = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd residual l2')

      mpi_reduce_send(1) = local_mdot_div_min
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      mdot_div_min = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd mdot div min')

      mpi_reduce_send(1) = local_mdot_div_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      mdot_div_max = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd mdot div max')

      mpi_reduce_send(1) = local_mdot_div_l2_num
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      mdot_div_l2_num = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd mdot div l2')

      mpi_reduce_send(1) = local_mass
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      mass_integral = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd mass')

      mpi_reduce_send(1) = local_net_mdot
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      net_boundary_mdot = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd net mdot')

      mpi_reduce_send(1) = local_volume
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      volume_total = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd volume')

      mpi_reduce_send(1) = local_h_min
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      h_min = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd h_min')

      mpi_reduce_send(1) = local_h_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      h_max = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd h_max')

      mpi_reduce_send(1) = local_T_min
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      T_min = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd T_min')

      mpi_reduce_send(1) = local_T_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      T_max = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd T_max')

      mpi_reduce_send(1) = local_worst_abs
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      global_worst_abs = mpi_reduce_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd worst residual')

      if (abs(local_worst_abs - global_worst_abs) <= max(1.0e-12_rk, 1.0e-10_rk * max(1.0_rk, global_worst_abs))) then
         local_winner_rank = flow%rank
      else
         local_winner_rank = huge(1)
      end if
      mpi_int_send(1) = local_winner_rank
      call MPI_Allreduce(mpi_int_send, mpi_int_recv, 1, MPI_INTEGER, MPI_MIN, flow%comm, ierr)
      global_winner_rank = mpi_int_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing vd winner rank')

      if (flow%rank == global_winner_rank) then
         mpi_int_send(1) = local_worst_cell
      else
         mpi_int_send(1) = 0
      end if
      call MPI_Allreduce(mpi_int_send, mpi_int_recv, 1, MPI_INTEGER, MPI_SUM, flow%comm, ierr)
      global_worst_cell = mpi_int_recv(1)

      mpi_reduce_send(1) = merge(local_worst_divu, 0.0_rk, flow%rank == global_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_worst_divu = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_worst_S, 0.0_rk, flow%rank == global_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_worst_S = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_worst_mdot_div, 0.0_rk, flow%rank == global_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_worst_mdot_div = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_worst_rho, 0.0_rk, flow%rank == global_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_worst_rho = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_worst_T, 0.0_rk, flow%rank == global_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_worst_T = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_worst_h, 0.0_rk, flow%rank == global_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_worst_h = mpi_reduce_recv(1)

      if (volume_total > 0.0_rk) then
         res_l2 = sqrt(res_l2_num / volume_total)
         mdot_div_l2 = sqrt(mdot_div_l2_num / volume_total)
      else
         res_l2 = 0.0_rk
         mdot_div_l2 = 0.0_rk
      end if

      if (flow%rank == 0) then
         filename = trim(params%output_dir) // '/diagnostics/variable_density_diagnostics.csv'

         if (.not. vd_diag_initialized) then
            file_exists = .false.
            open(newunit=unit_id, file=trim(filename), status='replace', action='write')
            vd_diag_initialized = .true.
         else
            file_exists = .true.
            open(newunit=unit_id, file=trim(filename), status='unknown', position='append', action='write')
         end if

         if (.not. file_exists) then
            write(unit_id,'(a)') 'step,time,rho_min,rho_max,rho_old_min,rho_old_max,S_min,S_max,' // &
                                'divu_min,divu_max,divu_minus_S_current_max,divu_minus_S_current_l2,' // &
                                'mass_flux_div_min,mass_flux_div_max,mass_flux_div_l2,' // &
                                'mass_integral,net_boundary_mass_flux,h_min,h_max,T_min,T_max,' // &
                                'worst_residual_abs,worst_residual_rank,worst_residual_cell,' // &
                                'worst_divu,worst_S,worst_mass_flux_div,worst_rho,worst_T,worst_h'
         end if

         write(unit_id,'(i0,21(",",ES26.16E4),2(",",i0),6(",",ES26.16E4))') step, time, rho_min, rho_max, &
              rho_old_min, rho_old_max, S_min, S_max, div_min, div_max, res_max, res_l2, &
              mdot_div_min, mdot_div_max, mdot_div_l2, mass_integral, net_boundary_mdot, &
              h_min, h_max, T_min, T_max, global_worst_abs, global_winner_rank, global_worst_cell, &
              global_worst_divu, global_worst_S, global_worst_mdot_div, global_worst_rho, &
              global_worst_T, global_worst_h
         close(unit_id)
      end if

      if (.not. vd_worst_initialized) then
         if (flow%rank == 0) then
            filename = trim(params%output_dir) // '/diagnostics/variable_density_worst_cell.csv'
            open(newunit=unit_id, file=trim(filename), status='replace', action='write')
            write(unit_id,'(a)') 'step,time,rank,cell,x,y,z,volume,rho,T,h,S,divu,divu_minus_S,mass_flux_div,abs_residual'
            close(unit_id)

            filename = trim(params%output_dir) // '/diagnostics/variable_density_worst_cell_faces.csv'
            open(newunit=unit_id, file=trim(filename), status='replace', action='write')
            write(unit_id,'(a)') 'step,time,rank,cell,local_face,face_id,owner,neighbor,face_area,face_flux,mass_flux,outward_face_flux,outward_mass_flux'
            close(unit_id)
         end if
         vd_worst_initialized = .true.
      end if

      call MPI_Barrier(flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure in vd worst-cell IO barrier')

      if (flow%rank == global_winner_rank .and. local_worst_cell > 0) then
         filename = trim(params%output_dir) // '/diagnostics/variable_density_worst_cell.csv'
         open(newunit=unit_id, file=trim(filename), status='unknown', position='append', action='write')
         write(unit_id,'(i0,",",ES26.16E4,",",i0,",",i0,12(",",ES26.16E4))') step, time, flow%rank, local_worst_cell, &
              0.0_rk, 0.0_rk, 0.0_rk, mesh%cells(local_worst_cell)%volume, local_worst_rho, local_worst_T, &
              local_worst_h, local_worst_S, local_worst_divu, local_worst_divu - local_worst_S, &
              local_worst_mdot_div, local_worst_abs
         close(unit_id)

         filename = trim(params%output_dir) // '/diagnostics/variable_density_worst_cell_faces.csv'
         open(newunit=unit_id, file=trim(filename), status='unknown', position='append', action='write')
         do lf = 1, mesh%ncell_faces(local_worst_cell)
            fid = mesh%cell_faces(lf, local_worst_cell)
            face_area_value = 0.0_rk

            if (allocated(fields%mass_flux)) then
               mass_flux_value = fields%mass_flux(fid)
            else
               mass_flux_value = 0.0_rk
            end if

            if (mesh%faces(fid)%owner == local_worst_cell) then
               flux_out = fields%face_flux(fid)
               mdot_out = mass_flux_value
            else
               flux_out = -fields%face_flux(fid)
               mdot_out = -mass_flux_value
            end if

            write(unit_id,'(i0,",",ES26.16E4,6(",",i0),5(",",ES26.16E4))') &
                 step, time, flow%rank, local_worst_cell, lf, fid, mesh%faces(fid)%owner, &
                 mesh%faces(fid)%neighbor, face_area_value, fields%face_flux(fid), &
                 mass_flux_value, flux_out, mdot_out
         end do
         close(unit_id)
      end if
      call write_variable_density_boundary_residual_scan(mesh, flow, params, fields, energy, transport, step, time)
      call write_variable_density_projection_audit(mesh, flow, params, fields, energy, transport, step, time)
      call write_variable_density_transport_conservation_diagnostics(mesh, flow, params, fields, energy, transport, step, time)
      call write_variable_density_compatibility_diagnostics(mesh, flow, params, fields, step, time)
      call write_variable_density_continuity_residual_diagnostics(mesh, flow, params, fields, transport, step, time)

   end subroutine write_variable_density_diagnostics