      Subroutine Flux(V,Fluxab,Fluxbc,Fluxcd,Fluxda,xloc,yloc,
     $ p,lmax,mx,my)
      Implicit none
      Integer lmax,mx,my
      Real  V(lmax,mx+1,my+1)
      Real Fluxab(lmax,mx+1,my+1),Fluxbc(lmax,mx+1,my+1),
     $ Fluxcd(lmax,mx+1,my+1),Fluxda(lmax,mx+1,my+1)
      Real  xloc(mx,my+1),yloc(mx,my+1)
      Real  p(mx+1,my+1)
      Integer l,i,j

!HPF$ PROCESSORS PROC(NUMBER_OF_PROCESSORS())
!HPF$ DISTRIBUTE  V *(*,*,BLOCK) onto *PROC
!HPF$ ALIGN (*,:) with *V(*,*,:) :: xloc,yloc
!HPF$ ALIGN (*,*,:) with *V (*,*,:) ::Fluxab,Fluxbc,Fluxcd,Fluxda
!HPF$ ALIGN (*,:) with *V (*,*,:) :: p
!************Compute Boundary Fluxes***************
! --------------Top (solid) Boundary---------------
      Fluxcd=0.0
      Forall(i=2:mx)
	 Fluxcd(2,i,my)=p(i,my)*(yloc(i-1,my)-yloc(i,my))
	
         Fluxcd(3,i,my)=-1.0*p(i,my)*(xloc(i-1,my)
     $ -xloc(i,my))
         End Forall
! -------------Bottom (symmetry) Boundary-----------------
      Fluxab=0.0

      Forall(i=2:mx) 
      Fluxab(1,i,2)=V(2,i,2)*(yloc(i,1)-yloc(i-1,1))
       Fluxab(2,i,2)=(V(2,i,2)**2/V(1,i,2)
     $  + p(i,2))*(yloc(i,1)-yloc(i-1,1))
       Fluxab(3,i,2)=-1.0*p(i,2)
     $  *(xloc(i,1)-xloc(i-1,1))
       Fluxab(4,i,2)=(V(4,i,2)+p(i,2))
     $ *V(2,i,2)/V(1,i,2)*(yloc(i,1)-yloc(i-1,1))
       End Forall

!****************Start Interior Points*************

	Forall(i=2:mx,j=3:my)
        Fluxab(1,i,j)=0.5*(V(2,i,j)+
     $ V(2,i,j-1))*(yloc(i,j-1)-yloc(i-1,j-1)) 
     $ -0.5*(V(3,i,j)+V(3,i,j-1))*(xloc(i,j-1)-xloc(i-1,j-1))
	End forall
	
	

	Forall(i=2:mx,j=3:my)
         Fluxab(2,i,j)=0.5*(V(2,i,j)**2
     $ / V(1,i,j) +p(i,j) + V(2,i,j-1)**2/V(1,i,j-1) +p(i,j-1))
     $ *(yloc(i,j-1)-yloc(i-1,j-1))
	Fluxab(2,i,j)=Fluxab(2,i,j) -0.5*(V(2,i,j)*V(3,i,j)/V(1,i,j)
     $ + V(2,i,j-1)*V(3,i,j-1)/V(1,i,j-1))
     $ *(xloc(i,j-1)-xloc(i-1,j-1)) 
        End Forall

	Forall(i=2:mx,j=3:my)
         Fluxab(3,i,j)=0.5*(V(2,i,j)
     $  *v(3,i,j)/V(1,i,j) + V(2,i,j-1)*V(3,i,j-1)/V(1,i,j-1)
     $ )*(yloc(i,j-1)-yloc(i-1,j-1))
	Fluxab(3,i,j)=Fluxab(3,i,j) -0.5*(V(3,i,j)**2/V(1,i,j)
     $ +p(i,j) +p(i,j-1)+ V(3,i,j-1)**2/V(1,i,j-1))
     $ *(xloc(i,j-1)-xloc(i-1,j-1))
	 end forall

	Forall(i=2:mx,j=3:my)
         Fluxab(4,i,j)=(V(4,i,j)
     $ + p(i,j))*V(2,i,j)/V(1,i,j)
	 Fluxab(4,i,j)=Fluxab(4,i,j) + (V(4,i,j-1) +
     $ p(i,j-1))*V(2,i,j-1)/V(1,i,j-1)
	Fluxab(4,i,j)=0.5*Fluxab(4,i,j)*(yloc(i,j-1)-yloc(i-1,j-1))
	Fluxab(4,i,j)=Fluxab(4,i,j)
     $ -0.5*(V(4,i,j)+p(i,j))*V(3,i,j)/V(1,i,j)
     $ *(xloc(i,j-1)-xloc(i-1,j-1))
	 Fluxab(4,i,j)=Fluxab(4,i,j)-
     $ 0.5*(V(4,i,j-1)+p(i,j-1))*V(3,i,j-1)/V(1,i,j-1)
     $ * (xloc(i,j-1)-xloc(i-1,j-1))

	End Forall 
		
! Now bc fluxes

	Forall(i=1:mx,j=2:my)
        Fluxbc(1,i,j)=0.5*(V(2,i,j)
     $ + V(2,i+1,j))*(yloc(i,j)-yloc(i,j-1)) -0.5*(V(3,i,j)
     $ + V(3,i+1,j))*(xloc(i,j)-xloc(i,j-1))
	End Forall

	Forall(i=1:mx,j=2:my)
         Fluxbc(2,i,j)=0.5*(V(2,i,j)**2
     $ /V(1,i,j) + p(i,j)+V(2,i+1,j)**2/V(1,i+1,j)+p(i+1,j))
     $ *(yloc(i,j)-yloc(i,j-1))
	Fluxbc(2,i,j)=Fluxbc(2,i,j)
     $	 -0.5*(V(3,i,j)*V(2,i,j)/V(1,i,j)
     $ + V(3,i+1,j)*V(2,i+1,j)/V(1,i+1,j))
     $ * (xloc(i,j)-xloc(i,j-1))
	end Forall
	
	Forall(i=1:mx,j=2:my)
         Fluxbc(3,i,j)=0.5*(V(2,i,j)
     $  *V(3,i,j)/V(1,i,j) + V(2,i+1,j)*V(3,i+1,j)/V(1,i+1,j))
     $ *(yloc(i,j)-yloc(i,j-1)) 
	Fluxbc(3,i,j)=Fluxbc(3,i,j)-0.5*(V(3,i,j)**2/V(1,i,j)
     $ + V(3,i+1,j)**2/V(1,i+1,j) +p(i,j) +p(i+1,j))
     $ *(xloc(i,j)-xloc(i,j-1))
	End Forall

	Forall(i=1:mx,j=2:my)
         Fluxbc(4,i,j)=(V(4,i,j)
     $ + P(i,j))*V(2,i,j)/V(1,i,j)
	Fluxbc(4,i,j)=Fluxbc(4,i,j) + (V(4,i+1,j)+p(i+1,j))
     $ *V(2,i+1,j)/V(1,i+1,j)
	Fluxbc(4,i,j)=0.5*Fluxbc(4,i,j)*(yloc(i,j)-yloc(i,j-1))
     $  -0.5*(V(4,i,j)
     $ + P(i,j))*V(3,i,j)/V(1,i,j)*(xloc(i,j)-xloc(i,j-1))
        Fluxbc(4,i,j)=Fluxbc(4,i,j)-0.5* (V(4,i+1,j)+p(i+1,j))
     $ *V(3,i+1,j)/V(1,i+1,j)*(xloc(i,j)-xloc(i,j-1))
	End Forall


        Forall(l=1:lmax,i=2:mx,j=2:my-1)
        fluxcd(l,i,j)=
     $ -1.0*fluxab(l,i,j+1)
        End FORALL
	

        Forall(l=1:lmax,i=2:mx,j=2:my)
        fluxda(l,i,j)=
     $ -1.0*fluxbc(l,i-1,j)
        End Forall
	
    
	End Subroutine
 

