write_boundary_mass_balance_diagnostics Subroutine

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

      integer :: p, i, f, owner, ierr, unit_id, ios
      real(rk), allocatable :: local_signed(:), global_signed(:)
      real(rk), allocatable :: local_inlet(:), global_inlet(:)
      real(rk), allocatable :: local_outlet(:), global_outlet(:)
      real(rk), allocatable :: local_area(:), global_area(:)
      real(rk), allocatable :: local_max_un(:), global_max_un(:)
      real(rk), allocatable :: local_max_umag(:), global_max_umag(:)
      real(rk) :: mdot, area, rho_face, un, umag
      character(len=1024) :: filename
      logical, save :: initialized = .false.

      if (.not. params%write_diagnostics) return
      if (.not. allocated(fields%mass_flux)) return
      if (.not. allocated(transport%rho)) return
      if (mesh%npatches <= 0) return

      allocate(local_signed(mesh%npatches), global_signed(mesh%npatches))
      allocate(local_inlet(mesh%npatches), global_inlet(mesh%npatches))
      allocate(local_outlet(mesh%npatches), global_outlet(mesh%npatches))
      allocate(local_area(mesh%npatches), global_area(mesh%npatches))
      allocate(local_max_un(mesh%npatches), global_max_un(mesh%npatches))
      allocate(local_max_umag(mesh%npatches), global_max_umag(mesh%npatches))

      local_signed = zero
      local_inlet = zero
      local_outlet = zero
      local_area = zero
      local_max_un = zero
      local_max_umag = zero

      do p = 1, mesh%npatches
         do i = 1, mesh%patches(p)%nfaces
            f = mesh%patches(p)%face_ids(i)
            if (f <= 0 .or. f > mesh%nfaces) cycle
            owner = mesh%faces(f)%owner
            if (owner < flow%first_cell .or. owner > flow%last_cell) cycle
            area = max(mesh%faces(f)%area, tiny(1.0_rk))
            mdot = fields%mass_flux(f)
            rho_face = max(transport%rho(owner), tiny(1.0_rk))
            un = mdot / (rho_face * area)
            umag = sqrt(sum(fields%u(:, owner)**2))

            local_signed(p) = local_signed(p) + mdot
            if (mdot < zero) then
               local_inlet(p) = local_inlet(p) - mdot
            else
               local_outlet(p) = local_outlet(p) + mdot
            end if
            local_area(p) = local_area(p) + area
            local_max_un(p) = max(local_max_un(p), abs(un))
            local_max_umag(p) = max(local_max_umag(p), umag)
         end do
      end do

      call MPI_Reduce(local_signed, global_signed, mesh%npatches, MPI_DOUBLE_PRECISION, MPI_SUM, 0, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI_Reduce failed for boundary signed mass flux')
      call MPI_Reduce(local_inlet, global_inlet, mesh%npatches, MPI_DOUBLE_PRECISION, MPI_SUM, 0, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI_Reduce failed for boundary inlet mass flux')
      call MPI_Reduce(local_outlet, global_outlet, mesh%npatches, MPI_DOUBLE_PRECISION, MPI_SUM, 0, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI_Reduce failed for boundary outlet mass flux')
      call MPI_Reduce(local_area, global_area, mesh%npatches, MPI_DOUBLE_PRECISION, MPI_SUM, 0, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI_Reduce failed for boundary patch area')
      call MPI_Reduce(local_max_un, global_max_un, mesh%npatches, MPI_DOUBLE_PRECISION, MPI_MAX, 0, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI_Reduce failed for boundary max normal velocity')
      call MPI_Reduce(local_max_umag, global_max_umag, mesh%npatches, MPI_DOUBLE_PRECISION, MPI_MAX, 0, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI_Reduce failed for boundary max velocity magnitude')

      if (flow%rank == 0) then
         filename = trim(params%output_dir)//'/diagnostics/boundary_mass_balance.csv'
         if (.not. initialized) then
            open(newunit=unit_id, file=trim(filename), status='replace', action='write', iostat=ios)
            if (ios == 0) then
               write(unit_id,'(a)') 'step,time,patch_id,patch_name,signed_mdot,inlet_mdot,outlet_mdot,patch_area,max_abs_normal_velocity,max_cell_velocity'
               close(unit_id)
            end if
            initialized = .true.
         end if
         open(newunit=unit_id, file=trim(filename), status='old', position='append', action='write', iostat=ios)
         if (ios == 0) then
            do p = 1, mesh%npatches
               write(unit_id,'(i0,a,es16.8,a,i0,a,a,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8,a,es16.8)') &
                  step, ',', time, ',', p, ',', trim(mesh%patches(p)%name), ',', global_signed(p), ',', &
                  global_inlet(p), ',', global_outlet(p), ',', global_area(p), ',', global_max_un(p), ',', &
                  global_max_umag(p)
            end do
            close(unit_id)
         end if
      end if

      deallocate(local_signed, global_signed, local_inlet, global_inlet, local_outlet, global_outlet, &
                 local_area, global_area, local_max_un, global_max_un, local_max_umag, global_max_umag)
   end subroutine write_boundary_mass_balance_diagnostics