
      Subroutine Artvisc(V,avab,avbc,avcd,avda,xloc,yloc,
     $ dnax,dnay,dnbx,dnby,p,lmax,mx,my)
      Implicit none
      Integer lmax,mx,my
      Real   V(lmax,mx+1,my+1)
      Real   avab(lmax,mx+1,my+1),avbc(lmax,mx+1,my+1),
     $ avcd(lmax,mx+1,my+1),avda(lmax,mx+1,my+1)
      Real  xloc(mx,my+1),yloc(mx,my+1),dnbx(mx,my+1),dnby(mx,my+1)
      Real   dnax(mx,my+1),dnay(mx,my+1)
      Real  p(mx+1,my+1)
      Real eps2(mx+1,my+1),eps4(mx+1,my+1)
      Real kp2,kp4,start,end
      Integer k, i,j,imy

!HPF$ PROCESSORS PROC(NUMBER_OF_PROCESSORS())
!HPF$ DISTRIBUTE  V *(*,*,BLOCK) onto *PROC
!HPF$ ALIGN (*,:) with *V(*,*,:) :: xloc,yloc,dnax,dnay,dnbx,dnby
!HPF$ ALIGN (*,*,:) with *V (*,*,:) ::avab,avbc,avcd,avda
!HPF$ ALIGN (*,:) with *V(*,*,:) :: p
!HPF$ ALIGN  (*,:) with V(*,*,:) :: eps2,eps4 


      kp2=4.0/4.0
      kp4=4.0/256.0

	avab=0.0
	avbc=0.0
	avcd=0.0
	avda=0.0
	
	


     
!***********Side AB & CD fluxes**********
        


	 p(2:mx+1,my+1) = p(2:mx+1,my)
      p(2:mx+1,1) = p(2:mx+1,2)
    	V(:,:,my+1)=V(:,:,my)
	  V(:,:,1)=V(:,:,2)
      Forall(k=1:lmax,i=2:mx)
	 avab(k,i,2)=0.0
       avcd(k,i,my)=0.0
	End Forall


      Forall(i=2:mx,j=3:my)

      eps2(i,j)=abs(p(i,j+1)
     $-2.0*p(i,j)+p(i,j-1))/abs(p(i,j+1)+2.0*p(i,j) +
     $ p(i,j-1))
	eps2(I,J)=kp2*MAX(eps2(i,j),abs(p(i,j)-2.0*p(i,j-1)+p(i,j-2))
     $ /abs( p(i,j)+2.0*p(i,j-1)+p(i,j-2)))
      
       eps4(i,j)=max(0.0,(kp4-eps2(i,j)))
       End Forall

       Forall(k=1:lmax,i=2:mx,j=3:my)

      	avab(k,i,j)=0.5*(V(2,i,j)/V(1,i,j)+V(2,i,j-1)/V(1,i,j-1))
     $ *dnax(i,j)

	avab(k,i,j)=avab(k,i,j)+0.5*(V(3,i,j)/V(1,i,j)+
     $V(3,i,j-1)/V(1,i,j-1))*
     $ dnay(i,j)
	avab(k,i,j)=abs(avab(k,i,j))
	avab(k,i,j)=avab(k,i,j)/sqrt(dnax(i,j)**2 + dnay(i,j)**2)
	avab(k,i,j)=avab(k,i,j)+sqrt(1.4*(p(i,j)+p(i,j-1))/
     $ (V(1,i,j)+V(1,i,j-1)))
	avab(k,i,j)=avab(k,i,j)*(eps2(i,j)*(V(k,i,j)-V(k,i,j-1))
     $ -eps4(i,j)*(V(k,i,j+1)-3.*V(k,i,j)+3.*V(k,i,j-1)-V(k,i,j-2)))
     
       End Forall 
	
    
!      Forall(i=2:mx)  eps2(i,3)=kp2*(max(abs(p(i,4)
!     $ -2.0*p(i,3)+p(i,2))/abs(p(i,4)+2.0*p(i,3) +
!     $ p(i,2)),abs(p(i,3)-1.0*p(i,2))/abs(
!     $ p(i,2)+3.0*p(i,2))))

!	eps2(:,3)=kp2*(max(abs(p(:,4)
!     $ -2.0*p(:,3)+p(:,2))/abs(p(:,4)+2.0*p(:,3) +
!     $ p(:,2)),abs(p(:,3)-1.0*p(:,2))/abs(
!     $ p(:,2)+3.0*p(:,2))))

!      Forall(i=2:mx) eps4(i,3)=max(0.0,(kp4-eps2(i,3)))

!      Forall(k=1:lmax,i=2:mx) avab(k,i,3)=(abs(0.5*(
!     $ V(2,i,3)/V(1,i,3)+ V(2,i,2)/V(1,i,2))*dnax(i,3)
!     $ + 0.5*(V(3,i,3)/V(1,i,3)+V(3,i,2)/V(1,i,2))*dnay(i,3))
!     $ /sqrt(dnax(i,3)**2+dnay(i,3)**2)
!     $ + sqrt(1.4*(p(i,3)+ p(i,2))/(V(1,i,2)+V(1,i,3))))*(eps2(i,3)
!     $ *(V(k,i,3)-V(k,i,2)) - eps4(i,3)*(V(k,i,4)-3.0*V(k,i,3)
!     $ + 2.0*V(k,i,2)))
	

!      Forall(i=2:mx) eps2(i,my)=max(2.0*eps2(i,my-1)-eps2(i,my-2),
!     $ eps2(i,my-1))
!	eps2(:,my)=max(2.0*eps2(:,my-1)-eps2(:,my-2),eps2(:,my-1))
!      Forall(i=2:mx) eps4(i,my)=max(0.0,(kp4-eps2(i,my)))
      
!       Forall(k=1:lmax,i=2:mx) avab(k,i,my)=
!     $ (abs(0.5*(V(2,i,my)/V(1,i,my)+V(2,i,my-1)/V(1,i,my-1))
!     $ *dnax(i,my) +0.5*(V(3,i,my)/V(1,i,my)+V(3,i,my-1)/
!     $ V(1,i,my-1))
!     $ *dnay(i,my))/sqrt(dnax(i,my)**2 + dnay(i,my)**2) +
!     $  sqrt(1.4*(p(i,my)+
!     $ p(i,my-1))/(V(1,i,my)+V(1,i,my-1))))*(eps2(i,my)*(V(k,i,my)
!     $ -V(k,i,my-1)) - eps4(i,my)*(-2.0*V(k,i,my)+3.0*V(k,i,my-1)
!     $ -V(k,i,my-2)))
     

      Forall(k=1:lmax,i=2:mx,j=2:my-1) avcd(k,i,j)=1.0*avab(k,i,j+1)
      
!***********now BC and DA fluxes*********

      Forall(k=1:lmax,j=2:my)
       avda(k,2,j)=0.0
       avbc(k,mx,j)=0.0
       V(k,mx+1,j)=V(k,mx,j)
	V(k,1,j)=V(k,2,j)
      End Forall

      Forall(j=1:my+1)
       p(mx+1,j)=p(mx,j)
	p(1,j)=p(2,j)
      End Forall

      	
      Forall(i=2:mx-1,j=2:my)
      eps2(i,j)=abs(p(i+1,j)
     $ -2.0*p(i,j)+p(i-1,j))/abs(p(i+1,j)+2.0*p(i,j)+p(i-1,j))

	eps2(i,j)=kp2*max(eps2(i,j),
     $ abs(p(i+2,j)-2.0*p(i+1,j)+p(i,j))/abs(p(i+2,j)+2.0*p(i+1,j)
     $ +p(i,j)))

       eps4(i,j)=max(0.0,(kp4-eps2(i,j)))
       End Forall

      Forall(k=1:lmax,i=2:mx-1,j=2:my)
      avbc(k,i,j)=
     $ 0.5*
     $ (V(2,i,j)/V(1,i,j)+V(2,i+1,j)/V(1,i+1,j))*dnbx(i,j)

	 avbc(k,i,j)=avbc(k,i,j)+
     $ 0.5*(V(3,i,j)/V(1,i,j)+V(3,i+1,j)/V(1,i+1,j))*dnby(i,j)
	avbc(k,i,j)=abs(avbc(k,i,j))

	avbc(k,i,j)=avbc(k,i,j)
     $ /sqrt(dnbx(i,j)**2+dnby(i,j)**2) 
	avbc(k,i,j)=avbc(k,i,j)+ sqrt(1.4*(p(i,j)
     $ + p(i+1,j))/(V(1,i,j)+V(1,i+1,j)))

	avbc(k,i,j)=avbc(k,i,j)*(eps2(i,j)*(V(k,i+1,j)
     $ -V(k,i,j)) -eps4(i,j)*(V(k,i+2,j)-3.0*V(k,i+1,j)
     $ + 3.0*V(k,i,j)-V(k,i-1,j)))
      End Forall

	

      Forall(k=1:lmax,i=3:mx,j=2:my) avda(k,i,j)=1.0*avbc(k,i-1,j)

	
	End subroutine

