output Subroutine

   subroutine output(action)
   !! Auxiliary routine to save results to file.
      integer, intent(in) :: action
         !! parameter to select output action

      character(*), parameter :: folder = ".\output\example2\"
      real(rk) :: time_elapsed
      integer, save :: funit_x(size(nc)) = 0, funit_u = 0
      integer :: i, j

      select case (action)

         ! Open files and write headers and grid
      case (1)

         write (stdout, '(1x, a)') "Running example2 ..."
         write (stdout, '(1x, a, i3)') "Max # threads: ", omp_get_max_threads()
         
         ! Write grid
         do concurrent(j=1:size(nc))
            open (newunit=funit_x(j), file=folder//"x"//to_string(j)//".txt", &
                  status="replace", action="write", position="rewind")

            write (funit_x(j), '(a5, 2(1x, a15))') "i", "x(i)", "dx(i)"
            do i = 1, nc(j)
               write (funit_x(j), '(i5, 2(1x, es15.5))') i, gx(j)%center(i), gx(j)%width(i)
            end do
         end do
         
         call timer()

         ! Write header u
         open (newunit=funit_u, file=folder//"u.txt", status="replace", &
               action="write", position="rewind")

         write (funit_u, '(a16)', advance="no") "t"
         do i = 1, size(u)
            write (funit_u, '(1x, a16)', advance="no") "u("//to_string(i)//")"
         end do
         write (funit_u, *) ""

         ! Write values u(x,t)
      case (2)
         write (funit_u, '(es16.5e3)', advance="no") time
         do i = 1, size(u)
            write (funit_u, '(1x, es16.5e3)', advance="no") u(i)
         end do
         write (funit_u, *) ""

         ! Close files
      case (3)
         do concurrent(j=1:size(nc))
            close (funit_x(j))
         end do
         close (funit_u)
         call timer(time_elapsed)
         write (stdout, '(1x, a, 1x, f5.2)') "Elapsed time (s) :", time_elapsed

      end select

   end subroutine output