!---------------------------------------------------------------------
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!     this function computes the norm of the difference between the
!     computed solution and the exact solution
!---------------------------------------------------------------------
      subroutine error_norm (rms)

      include 'header3d.h'
      integer  i,j,k,m,d
      double precision  xi,eta,zeta,u_exact(5),rms(5),add
      double precision  r1,r2,r3,r4,r5
      do  m = 1,5
         rms(m) = 0.0d0
      enddo  
      r1 = 0.0d0
      r2 = 0.0d0
      r3 = 0.0d0
      r4 = 0.0d0
      r5 = 0.0d0

!DVM$ region
!DVM$ PARALLEL (k,j,i) ON  u(*,i,j,k),
!DVM$& REDUCTION(SUM(r1),SUM(r2),SUM(r3),SUM(r4),SUM(r5)),
!DVM$&private(u_exact,xi,eta,zeta,m,add)
      do  k = 0,problem_size - 1
         do  j = 0,problem_size - 1
            do  i = 0,problem_size - 1
               zeta = dble (k) * dnzm1
               eta = dble (j) * dnym1
               xi = dble (i) * dnxm1

!               call exact_solution(xi, eta, zeta, u_exact)
               do  m = 1,5
                  u_exact(m) = ce(m,1) + xi * (ce(m,2) + xi * (ce(m,5) +
     & xi * (ce(m,8) + xi * ce(m,11)))) + eta * (ce(m,3) + eta * (ce(m,6
     &) + eta * (ce(m,9) + eta * ce(m,12)))) + zeta * (ce(m,4) + zeta * 
     &(ce(m,7) + zeta * (ce(m,10) + zeta * ce(m,13))))
               enddo  
               add = u(1,i,j,k) - u_exact(1)
               r1 = r1 + add * add
               add = u(2,i,j,k) - u_exact(2)
               r2 = r2 + add * add
               add = u(3,i,j,k) - u_exact(3)
               r3 = r3 + add * add
               add = u(4,i,j,k) - u_exact(4)
               r4 = r4 + add * add
               add = u(5,i,j,k) - u_exact(5)
               r5 = r5 + add * add
            enddo  
         enddo  
      enddo  

!DVM$ end region
      rms(1) = r1
      rms(2) = r2
      rms(3) = r3
      rms(4) = r4
      rms(5) = r5
      do  m = 1,5
         do  d = 1,3
            rms(m) = rms(m) / dble (grid_points(d) - 2)
         enddo  
         rms(m) = dsqrt (rms(m))
      enddo  
      return
      end


!---------------------------------------------------------------------
!---------------------------------------------------------------------
      subroutine rhs_norm (rms)

      include 'header3d.h'
      integer  i,j,k,d,m
      double precision rms(5),add,r1,r2,r3,r4,r5
      r1 = 0.0d0
      r2 = 0.0d0
      r3 = 0.0d0
      r4 = 0.0d0
      r5 = 0.0d0

!DVM$ region
!DVM$ PARALLEL (k,j,i) ON  rhs(*,i,j,k), 
!DVM$&REDUCTION(SUM(r1),SUM(r2),SUM(r3),SUM(r4),SUM(r5)),
!DVM$&private(add)
      do  k = 1,problem_size - 2
         do  j = 1,problem_size - 2
            do  i = 1,problem_size - 2
               add = rhs(1,i,j,k)
               r1 = r1 + add * add
               add = rhs(2,i,j,k)
               r2 = r2 + add * add
               add = rhs(3,i,j,k)
               r3 = r3 + add * add
               add = rhs(4,i,j,k)
               r4 = r4 + add * add
               add = rhs(5,i,j,k)
               r5 = r5 + add * add
            enddo  
         enddo  
      enddo  

!DVM$ end region
      rms(1) = r1
      rms(2) = r2
      rms(3) = r3
      rms(4) = r4
      rms(5) = r5
      do  m = 1,5
         do  d = 1,3
            rms(m) = rms(m) / dble (grid_points(d) - 2)
         enddo  
         rms(m) = dsqrt (rms(m))
      enddo  
      return
      end