c     $Header$

#include "cctk.h"

      subroutine apply_dissipation (var, rhs, ni, nj, nk, dx, dt, order, epsdis)
      implicit none
      
      integer   ni, nj, nk
      CCTK_REAL var(ni,nj,nk), rhs(ni,nj,nk)
      CCTK_REAL dx(3), dt
      CCTK_INT  order
      CCTK_REAL epsdis
      
      integer   i, j, k
      
      if (order .eq. 1) then
         
         do k = 2, nk-1
            do j = 2, nj-1
               do i = 2, ni-1
                  
                  rhs(i,j,k) = rhs(i,j,k) + epsdis
     $                 * (+ (var(i-1,j,k) - 2*var(i,j,k) - var(i+1,j,k)) / dx(1)
     $                    + (var(i,j-1,k) - 2*var(i,j,k) - var(i,j+1,k)) / dx(2)
     $                    + (var(i,j,k-1) - 2*var(i,j,k) - var(i,j,k+1)) / dx(3))
                  
               end do
            end do
         end do
         
      else if (order .eq. 3) then
         
         do k = 3, nk-2
            do j = 3, nj-2
               do i = 3, ni-2
                  
                  rhs(i,j,k) = rhs(i,j,k) - epsdis / 16
     $                 * (+ (var(i-2,j,k) - 4*var(i-1,j,k) + 6*var(i,j,k) - 4*var(i+1,j,k) + var(i+2,j,k)) / dx(1)
     $                    + (var(i,j-2,k) - 4*var(i,j-1,k) + 6*var(i,j,k) - 4*var(i,j+1,k) + var(i,j+2,k)) / dx(2)
     $                    + (var(i,j,k-2) - 4*var(i,j,k-1) + 6*var(i,j,k) - 4*var(i,j,k+1) + var(i,j,k+2)) / dx(3))
                  
               end do
            end do
         end do
         
      else
         
         call CCTK_WARN (0, "internal error")
         
      end if
      
      end