      Subroutine James(V,xloc,yloc,dnax,dnay,dnbx,dnby,min,
     $ lmax,mx,my)
      Implicit none
      Integer lmax,mx,my
      Real V(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 Fluxab(lmax,mx+1,my+1),Fluxbc(lmax,mx+1,my+1)
      Real Fluxcd(lmax,mx+1,my+1),Fluxda(lmax,mx+1,my+1)
      Real Res(lmax,mx,my+1),h(lmax,mx+1,my+1)
      Real Q(lmax,mx+1,my+1)
      Real p(mx+1,my+1)
      Real avab(lmax,mx+1,my+1),avbc(lmax,mx+1,my+1),
     $ avcd(lmax,mx+1,my+1)
      Real avda(lmax,mx+1,my+1)
      Real dt(mx,my+1)
      Real area(mx,my+1)
      Real pback,cfl,min,start,end
      Integer l,i,j,z
      Integer exit,iter,itermax

!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 (*,*,:):: Q,Fluxab,Fluxbc,Fluxcd,Fluxda,h,
     $ avab,avbc,avcd,avda,res
!HPF$ ALIGN (*,:) with V (*,*,:) :: p,dt,area


      write(6,*) 'Please enter pback'
      Read(5,*) pback
      write(6,*) 'Please enter CFl number'
      read(5,*) cfl
      write(6,*) 'enter exit condition: 1=subsonic, 2=supersonic'
      read(5,*) exit
      iter=1

	if (exit .eq. 2) then
	  itermax=120.0
	end if

	if (exit .eq. 1) then
		itermax=10000 
	end if

	call timer(start,0.0)
		res(3,4,5)=0.4
	avab=0.0
	avbc=0.0
	avcd=0.0
	avda=0.0
	
	  p=pback
	open(unit=8,file='err.dat',status='unknown')
	open(unit=9,file='in.dat',status='unknown')
	
	
	do i=1,mx
	 write(9,*) xloc(i,2),V(2,i,2)/V(1,i,2)/sqrt(p(i,2)*1.4
     $ /V(1,i,2)),i
	end do

      Forall(i=2:mx,j=2:my) area(i,j)=0.5*abs((xloc(i,j)
     $ -xloc(i-1,j-1))*(yloc(i-1,j)-yloc(i,j-1)) - 
     $ (xloc(i-1,j)-xloc(i,j-1))*(yloc(i,j)-yloc(i-1,j-1)))

!-----------Start looping-------------
!General form of the loop:
!        -Compute time step at each axial location
!        -Compute flux vectors for each "volume"
!        -Compute dissipation vectors
!        -Compute residual
!        -Compute stage of Runge Kutta
!        -Repeat for all 4 stages
!        -Update solution 


      Do While(iter .le. itermax)
         Call time_step(V,dt,xloc,yloc,lmax,mx,my)
	
         Call Flux(V,Fluxab,Fluxbc,Fluxcd,Fluxda,xloc,yloc,p,lmax,mx,my)
	
	 Call Artvisc(V,avab,avbc,avcd,avda,xloc,yloc,
     $ dnax,dnay,dnbx,dnby,p,lmax,mx,my)

!         Call source(V,h,p,yloc,lmax,mx,my)
	
         Forall(l=1:lmax,i=2:mx,j=2:my)
	 Res(l,i,j)=-1.0/area(i,j)
     $ * (Fluxab(l,i,j)+Fluxbc(l,i,j)+Fluxcd(l,i,j)+Fluxda(l,i,j))
	Res(l,i,j)=Res(l,i,j) + 
     $ + 1.0/area(i,j)*(avbc(l,i,j)*sqrt(dnbx(i,j)**2
     $ + dnby(i,j)**2)-avda(l,i,j)*sqrt(dnbx(i-1,j)**2
     $ + dnby(i-1,j)**2))
	Res(l,i,j)=Res(l,i,j) +1.0/area(i,j)*( avcd(l,i,j)
     $ *sqrt(dnax(i,j+1)**2
     $ + dnay(i,j+1)**2)-avab(l,i,j)*sqrt(dnax(i,j)**2+dnay(i,j)
     $ **2)) + h(l,i,j)
	 Q(l,i,j)=V(l,i,j)
     $         + 0.25 *Res(l,i,j)*dt(i,j)
	End forall
         Forall(i=2:mx,j=2:my)
	 p(i,j)=0.4*(Q(4,i,j)-0.5/Q(1,i,j)
     $ *(Q(2,i,j)**2 + Q(3,i,j)**2))
	 End Forall

	  Call bound(Q,p,min,pback,lmax,mx,my,exit)
	 	           
         
          Call time_step(Q,dt,xloc,yloc,lmax,mx,my)
         Call Flux(Q,Fluxab,Fluxbc,Fluxcd,Fluxda,xloc,yloc,p,lmax,mx,my)
         Call Artvisc(Q,avab,avbc,avcd,avda,xloc,yloc,
     $ dnax,dnay,dnbx,dnby,p,lmax,mx,my)

!         Call source(Q,h,p,yloc,lmax,mx,my)
          Forall(l=1:lmax,i=2:mx,j=2:my)

         Res(l,i,j)=-1.0/area(i,j)
     $ * (Fluxab(l,i,j)+Fluxbc(l,i,j)+Fluxcd(l,i,j)+Fluxda(l,i,j))
        Res(l,i,j)=Res(l,i,j) + 
     $ + 1.0/area(i,j)*(avbc(l,i,j)*sqrt(dnbx(i,j)**2
     $ + dnby(i,j)**2)-avda(l,i,j)*sqrt(dnbx(i-1,j)**2
     $ + dnby(i-1,j)**2))
        Res(l,i,j)=Res(l,i,j) +1.0/area(i,j)*( avcd(l,i,j)
     $ *sqrt(dnax(i,j+1)**2
     $ + dnay(i,j+1)**2)-avab(l,i,j)*sqrt(dnax(i,j)**2+dnay(i,j)
     $ **2)) + h(l,i,j)


         Q(l,i,j)=V(l,i,j)
     $   + 1.0/3.0 *Res(l,i,j)*dt(i,j)
	End Forall
         Forall(i=2:mx,j=2:my)
	 p(i,j)=0.4*(Q(4,i,j)-0.5/Q(1,i,j)
     $ *(Q(2,i,j)**2 + Q(3,i,j)**2))
	 End Forall
         Call bound(Q,p,min,pback,lmax,mx,my,exit)

         Call time_step(Q,dt,xloc,yloc,lmax,mx,my)
         Call Flux(Q,Fluxab,Fluxbc,Fluxcd,Fluxda,xloc,yloc,p,lmax,mx,my)
         Call Artvisc(Q,avab,avbc,avcd,avda,xloc,yloc,
     $ dnax,dnay,dnbx,dnby,p,lmax,mx,my)
!         Call source(Q,h,p,yloc,lmax,mx,my)
         Forall(l=1:lmax,i=2:mx,j=2:my)

 	Res(l,i,j)=-1.0/area(i,j)
     $ * (Fluxab(l,i,j)+Fluxbc(l,i,j)+Fluxcd(l,i,j)+Fluxda(l,i,j))
        Res(l,i,j)=Res(l,i,j) + 
     $ + 1.0/area(i,j)*(avbc(l,i,j)*sqrt(dnbx(i,j)**2
     $ + dnby(i,j)**2)-avda(l,i,j)*sqrt(dnbx(i-1,j)**2
     $ + dnby(i-1,j)**2))
        Res(l,i,j)=Res(l,i,j) +1.0/area(i,j)*( avcd(l,i,j)
     $ *sqrt(dnax(i,j+1)**2
     $ + dnay(i,j+1)**2)-avab(l,i,j)*sqrt(dnax(i,j)**2+dnay(i,j)
     $ **2)) + h(l,i,j)



         Q(l,i,j)=V(l,i,j)
     $         + 1.0/2.0 *Res(l,i,j)*dt(i,j)
	End Forall
         Forall(i=2:mx,j=2:my)
	 p(i,j)=0.4*(Q(4,i,j)-0.5/Q(1,i,j)
     $ *(Q(2,i,j)**2 + Q(3,i,j)**2))
	End Forall
         Call bound(Q,p,min,pback,lmax,mx,my,exit)


         Call time_step(Q,dt,xloc,yloc,lmax,mx,my)
         Call Flux(Q,Fluxab,Fluxbc,Fluxcd,Fluxda,xloc,yloc,p,lmax,mx,my)
         Call Artvisc(Q,avab,avbc,avcd,avda,xloc,yloc,
     $ dnax,dnay,dnbx,dnby,p,lmax,mx,my)
!         Call source(Q,h,p,yloc,lmax,mx,my)
          Forall(l=1:lmax,i=2:mx,j=2:my)
	  Res(l,i,j)=-1.0/area(i,j)

     $ * (Fluxab(l,i,j)+Fluxbc(l,i,j)+Fluxcd(l,i,j)+Fluxda(l,i,j))
        Res(l,i,j)=Res(l,i,j) + 
     $ + 1.0/area(i,j)*(avbc(l,i,j)*sqrt(dnbx(i,j)**2
     $ + dnby(i,j)**2)-avda(l,i,j)*sqrt(dnbx(i-1,j)**2
     $ + dnby(i-1,j)**2))
        Res(l,i,j)=Res(l,i,j) +1.0/area(i,j)*( avcd(l,i,j)
     $ *sqrt(dnax(i,j+1)**2
     $ + dnay(i,j+1)**2)-avab(l,i,j)*sqrt(dnax(i,j)**2+dnay(i,j)
     $ **2)) + h(l,i,j)


         V(l,i,j)=V(l,i,j)
     $    + 1.0 *Res(l,i,j)*dt(i,j)
	End Forall
         Forall(i=2:mx,j=2:my)
	 p(i,j)=0.4*(V(4,i,j)-0.5/V(1,i,j)
     $ *(V(2,i,j)**2 + V(3,i,j)**2))
	End Forall

         Call bound(V,p,min,pback,lmax,mx,my,exit)
        
	iter=iter+1
	
	End do
	call timer(end,start)
	write(6,*) 'Execution took:',end,'seconds'
	call Contour(V,xloc,yloc,lmax,mx,my )
	
	
	do i=1,mx
	 write(8,*) xloc(i,8),sqrt((V(2,i,17)/V(1,i,17))**2
     $ + (V(3,i,17)/V(1,i,17))**2)/
     $ sqrt(1.4*p(i,17)/V(1,i,17)), sqrt((V(2,i,2)/V(1,i,2))**2
     $ + (V(3,i,2)/V(1,i,2))**2)/
     $ sqrt(1.4*p(i,2)/V(1,i,2)),v(3,i,2)/V(1,i,2)
	end do
	write(8,*) 'and this was the last residual',maxval(abs(res))
         End subroutine
