c---------------------------------------------------------------------
c---------------------------------------------------------------------

       subroutine error_norm(rms)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c this function computes the norm of the difference between the
c computed solution and the exact solution
c---------------------------------------------------------------------

       include 'header.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),private(zeta,eta,xi,add,u_exact,m)
!DVM$& ,reduction(SUM(r1),SUM(r2),SUM(r3),SUM(r4),SUM(r5))
! DVM$& ,shadow_renew(u, rhs)
       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
        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))))
        end do
		
		    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			
             end do
          end do
       end do
!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)
          end do
          rms(m) = dsqrt(rms(m))
       end do

       return
       end



       subroutine rhs_norm(rms)

       include 'header.h'

       integer i, j, k, d, m
       double precision rms(5), add

       do m = 1, 5
          rms(m) = 0.0d0
       enddo	   	   
	   
!DVM$ region
!DVM$ parallel (k,j,i) on u(*,i,j,k),private(add)
!DVM$& ,reduction(SUM(rms))
       do k = 1, nz2
          do j = 1, ny2
             do i = 1, nx2
               add = rhs(1,i,j,k)
               rms(1) = rms(1) + add*add
               add = rhs(2,i,j,k)
               rms(2) = rms(2) + add*add
               add = rhs(3,i,j,k)
               rms(3) = rms(3) + add*add
               add = rhs(4,i,j,k)
               rms(4) = rms(4) + add*add
               add = rhs(5,i,j,k)
               rms(5) = rms(5) + add*add			   
             end do 
          end do 
       end do 
!DVM$ end region   
       do   m = 1, 5
          do   d = 1, 3
             rms(m) = rms(m) / dble(grid_points(d)-2)
          end do
          rms(m) = dsqrt(rms(m))
       end do

       return
       end