output Subroutine

   impure subroutine output(t, c, ns, mx, my, lout)
   !! This routine prints the values of the individual species densities at the current time
   !! `t`, to logical unit `lout`.
   
      real(rk), intent(in) :: t
      real(rk), intent(in) :: c(ns, mx, my)
      integer, intent(in) :: ns
      integer, intent(in) :: mx
      integer, intent(in) :: my
      integer, intent(in) :: lout

      logical, save :: first = .true.
      integer :: i, jx, jy

      if (first) then
         write (lout, '(a16)', advance='no') 't'
         do i = 1, ns
            do jy = 1, my
               do jx = 1, mx
                  write (lout, '(1x, a16)', advance='no') &
                     "c("//to_string(i)//","//to_string(jx)//","//to_string(jy)//")"
               end do
            end do
         end do
         write (lout, *) ''
         first = .false.
      end if

      write (lout, '(es16.5e3)', advance="no") t
      do i = 1, ns
         do jy = 1, my
            do jx = 1, mx
               write (lout, '(1x, es16.5e3)', advance="no") c(i, jx, jy)
            end do
         end do
      end do
      write (lout, *) ''

   end subroutine output