c----------------------------------------------------------------------
c	This is a general purpose time marching 2-D Navier-Stokes solver
c
c	input data must be in a file 'chan.dat' in the run directory
c
c	program allows changes to occur during run time to
c	Reynolds number, the program can be instructed to
c	save everything and then be restarted later.
c
c	Uses multigrid to solve for stream function
c	from vorticity and different methods can be used
c	for time marching the vorticity
c
c-------------------------------------------------------------channel--
      program channel
      include 'step.h'

      call start
1     continue
      call runtimesequence
2     call whatnext(key)
      if(key.gt.0)then
         if(key.eq.5)then
            goto 2
         else
            goto 1
         endif
      endif
      call shutdown
      stop
      end
c-------------------------------------------------------------whatnext-
      subroutine whatnext(key)
      include 'step.h'
      write(logfile,22)
22    format(' 0=end,1=reflect,2=change Re,3=set dt,4=method,5=archive')
      read(inputfile,*)key
      if     (key.eq.swapflow)then
           write(logfile,23)
23         format('swapping:')
           call swap
      else if(key.eq.changereynolds)then
           read(inputfile,*)reynolds
           write(logfile,24)reynolds
24         format('reynolds changed to:',f10.2)
	    pex=dt/(reynolds*strouhal*dx*dx)
	    pey=dt/(reynolds*strouhal*dy*dy)

      else if(key.eq.changedt)then
           read(inputfile,*)newdt
           write(logfile,25)newdt
25         format('dt changed to:',f10.4)
	    dtratio=dt/newdt
	    printsteps=dtratio*printsteps
	    dtratio=newdt/dt
	    do 30 i=0,ix1
	    do 30 j=-iy1,channeltype*iy1
	     oldvort(i,j)=dtratio*oldvort(i,j)+(1.d0-dtratio)*vort(i,j)
30	    continue

	    dt=newdt
	    cx=dt/(strouhal*dx)
	    cy=dt/(strouhal*dy)
	    pex=dt/(reynolds*strouhal*dx*dx)
	    pey=dt/(reynolds*strouhal*dy*dy)


      else if(key.eq.changemethod)then
	    read(inputfile,*)method
           write(logfile,35)method
35         format('method changed to:',i4)
      else if(key.eq.archiveflow)then
	    call archive
           write(logfile,45)reynolds
45         format('writing to archive re=',f10.2)
      endif
      return
      end
c----------------------------------------------------------------------
      subroutine swap
      include 'step.h'
c
c	reflect y values
c
      do 10 i=0,ix1
         do 10 j=-iy1,iy1/2-1
            temp=psi(i,j)
            psi(i,j)=-psi(i,iy1-j)
            psi(i,iy1-j)=-temp

            temp=newvort(i,j)
            newvort(i,j)=-newvort(i,iy1-j)
            newvort(i,iy1-j)=-temp

            temp=vort(i,j)
            vort(i,j)=-vort(i,iy1-j)
            vort(i,iy1-j)=-temp

            temp=oldvort(i,j)
            oldvort(i,j)=-oldvort(i,iy1-j)
            oldvort(i,iy1-j)=-temp

            temp=u(i,j)
            u(i,j)=u(i,iy1-j)
            u(i,iy1-j)=temp

            temp=v(i,j)
            v(i,j)=-v(i,iy1-j)
            v(i,iy1-j)=-temp
10    continue
      return
      end
c-----------------------------------------------------------------------
      subroutine shutdown
      include 'step.h'
      write(psifile,30)ix1/2+1,(1+channeltype)*iy1/2+1
      write(vortfile,30)ix1/2+1,iy1+1
30    format('aspect:.2',/,'height:0. -1. 1.',/,'contour ',2i4)
      write(psifile,31)((psi(i,j),i=0,ix1,2),j=-iy1,channeltype*iy1,2)
      write(vortfile,31)((vort(i,j),i=0,ix1,2),j=-iy1,channeltype*iy1,2)
31	format(8e14.5)
       close(unit=psifile)
       close(unit=vortfile)
       close(unit=sumfile)
       close(unit=logfile)
       close(unit=inputfile)
       return
       end

c-----------------------------------------------------------------------
      subroutine runtimesequence
      include 'step.h'
      do 10 k=1,ns1
       call onetimestep
       call export(1)
	if((k.gt.10).and.(vincmax.lt.1.e-6))then
	  write(logfile,5)k,reynolds
5	  format('converged to steady flow after ',i5,'steps : Re =',f8.2)
	goto 15
	endif
10    continue
15     do 20 k=1,ns2
        call onetimestep
        call export(2)
20	continue
      return
      end

c-----------------------------------------------------------------------
      subroutine onetimestep
      include 'step.h'
           currenttime=currenttime+dt
           call setflux
	    if     (method.eq.central)then
                call stepcentral
	    else if(method.eq.upwind)then
		call stepupwind
	    else if(method.eq.ncupwind)then
		call stepnonc
	    else if(method.eq.laxwendroff)then
		call steplax
	    else if(method.eq.quickest)then
		call stepquick
	    endif
        call calculatepsi
        call calculatevelocities
        call calculateboundaryvorticity
        do 50 i=0,ix1
           do 50 j=-iy1,channeltype*iy1
              oldvort(i,j)=vort(i,j)
              vort(i,j)=newvort(i,j)
50      continue

      return
      end
c--------------------------------------------------------------------
      subroutine export(sequence)
      include 'step.h'
      integer sequence
      s1=0.
      s2=0.
      do 10 i=ixc1,ix1-1
         if(vort(i,-iy1)*vort(i+1,-iy1).lt.0.)then
           fac=(0.-vort(i,-iy1))/(vort(i+1,-iy1)-vort(i,-iy1))
           s1=(i+fac-ixc1)*dx
         endif
         if(vort(i,channeltype*iy1)*vort(i+1,channeltype*iy1).lt.0.)then
           fac=(0.-vort(i,channeltype*iy1))/
     +            (vort(i+1,channeltype*iy1)-vort(i,channeltype*iy1))
           s2=(i+fac-ixc1)*dx
         endif
10    continue
         pmax=0.
         pmin=0.
         ipmax=0
         ipmin=0
         jpmin=0
         jpmax=0
         do 12 i=0,ix1
         do 12 j=-iy1,channeltype*iy1
            if(psi(i,j).gt.pmax)then
               pmax=psi(i,j)
               ipmax=i
               jpmax=j
            endif
            if(psi(i,j).lt.pmin)then
               pmin=psi(i,j)
               ipmin=i
               jpmin=j
            endif
12          continue
	write(screen,15)currenttime,psi(0,iy1),reynolds,
     +   vort(ixc1,0),s1,ipmin,jpmin,pmin,
     +   ipmax,jpmax,pmax,vincmax
15	format(f8.4,f8.3,f8.2,f8.1,f8.3,2(2i4,f10.5),e14.4)
       if(sequence.eq.2)then
	  call calculateasymmetry(asym,vasym)
	  write(logfile,17)currenttime,reynolds,
     +   s1,ipmin,jpmin,pmin,
     +   s2,ipmax,jpmax,pmax
         write(sumfile,16)reynolds,s1,s2,asym,vasym
16       format(f8.2,2f10.4,2e14.4)
17       format(f8.3,f8.2,2(f8.4,2i4,f10.5))
       endif
      return
      end

c------------------------------------------------------------------
      subroutine calculateasymmetry(asym,vasym)
	include 'step.h'

         area=ixc1*dx+3*(ix1-ixc1)*dx
         asym=0.
         do 20 i=1,ix1
            do 10 j=iy1/2,channeltype*iy1
               f=psi(i,j)+psi(i,iy1-j)
               asym=asym+f*f
10          continue
20       continue
         asym=sqrt(asym*dx*dy/area)
	 vasym=(1.-vfac)*v(jv,iy1/2)+vfac*v(jv+1,iy1/2)
         asym=sign(asym,vasym)
      return
      end
c---------------------------------------------------------------------
      subroutine initialiseall()
      include 'step.h'
         do 20 i=0,ix1
            do 10 j=-iy1,channeltype*iy1
               psi(i,j)=0.
               vort(i,j)=0.
               oldvort(i,j)=0.
               newvort(i,j)=0.
               u(i,j)=0.
               v(i,j)=0.
10          continue
20       continue
      return
      end
c----------------------------------------------------------------------
      subroutine steplax
      include 'step.h'
	real *8 peu,pev
	
	vincmax=0.
        do 10 i=1,ixc1
        do 10 j=1,iy1-1
              au=-0.5*cx*u(i-1,j)
              bu=0.
              cu=0.5*cx*u(i+1,j)
              av=-0.5*cy*v(i,j-1)
              bv=0.
              cv=0.5*cy*v(i,j+1)
	      peu=pex + cx*cx*u(i,j)*u(i,j)/2
	      pev=pey + cy*cy*v(i,j)*v(i,j)/2
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-cu*vort(i+1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +peu*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pev*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
10         continue
        do 20 i=ixc1+1,ix1-1
        do 20 j=-iy1+1,channeltype*iy1-1
              au=-0.5*cx*u(i-1,j)
              bu=0.
              cu=0.5*cx*u(i+1,j)
              av=-0.5*cy*v(i,j-1)
              bv=0.
              cv=0.5*cy*v(i,j+1)
	      peu=pex + cx*cx*u(i,j)*u(i,j)/2
	      pev=pey + cy*cy*v(i,j)*v(i,j)/2
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-cu*vort(i+1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +peu*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pev*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
20      continue
	i=ix1
        do 30 j=-iy1+1,channeltype*iy1-1
              au=-0.5*cx*u(i-1,j)
              bu=0.
              cu=0.5*cx*u(i-1,j)
              av=-0.5*cy*v(i,j-1)
              bv=0.
              cv=0.5*cy*v(i,j+1)
	      peu=pex + cx*cx*u(i,j)*u(i,j)/2
	      pev=pey + cy*cy*v(i,j)*v(i,j)/2
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-cu*vort(i-1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +peu*(vort(i-1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pev*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
30      continue

      return
      end
c----------------------------------------------------------------------
      subroutine stepcentral
      include 'step.h'
	vincmax=0.
        do 10 i=1,ixc1
        do 10 j=1,iy1-1
              au=-0.5*cx*u(i-1,j)
              bu=0.
              cu=0.5*cx*u(i+1,j)
              av=-0.5*cy*v(i,j-1)
              bv=0.
              cv=0.5*cy*v(i,j+1)
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-cu*vort(i+1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
10         continue
        do 20 i=ixc1+1,ix1-1
        do 20 j=-iy1+1,channeltype*iy1-1
              au=-0.5*cx*u(i-1,j)
              bu=0.
              cu=0.5*cx*u(i+1,j)
              av=-0.5*cy*v(i,j-1)
              bv=0.
              cv=0.5*cy*v(i,j+1)
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-cu*vort(i+1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
20      continue
	i=ix1
        do 30 j=-iy1+1,channeltype*iy1-1
              au=-0.5*cx*u(i-1,j)
              bu=0.
              cu=0.5*cx*u(i-1,j)
              av=-0.5*cy*v(i,j-1)
              bv=0.
              cv=0.5*cy*v(i,j+1)
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-cu*vort(i-1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +pex*(vort(i-1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
30      continue

      return
      end
c----------------------------------------------------------------------
      subroutine stepnonc
      include 'step.h'
	vincmax=0.
        do 10 i=1,ixc1
        do 10 j=1,iy1-1
	   if(u(i,j).ge.0.)then
              au=-cx*u(i,j)
              bu= cx*u(i,j)
              cu= 0.
	   else
              au=0 .
              bu=-cx*u(i,j)
              cu= cx*u(i,j)
	   endif
	   if(v(i,j).ge.0.)then
              av=-cy*v(i,j)
              bv= cy*v(i,j)
              cv= 0.
	   else
              av= 0.
              bv=-cy*v(i,j)
              cv= cy*v(i,j)
	   endif
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i+1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
10         continue
        do 20 i=ixc1+1,ix1-1
        do 20 j=-iy1+1,channeltype*iy1-1
	   if(u(i,j).ge.0.)then
              au=-cx*u(i,j)
              bu= cx*u(i,j)
              cu= 0.
	   else
              au=0 .
              bu=-cx*u(i,j)
              cu= cx*u(i,j)
	   endif
	   if(v(i,j).ge.0.)then
              av=-cy*v(i,j)
              bv= cy*v(i,j)
              cv= 0.
	   else
              av= 0.
              bv=-cy*v(i,j)
              cv= cy*v(i,j)
	   endif
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i+1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
20      continue
	i=ix1
        do 30 j=-iy1+1,channeltype*iy1-1
	   if(u(i,j).ge.0.)then
              au=-cx*u(i,j)
              bu= cx*u(i,j)
              cu= 0.
	   else
              au=0 .
              bu=-cx*u(i,j)
              cu= cx*u(i,j)
	   endif
	   if(v(i,j).ge.0.)then
              av=-cy*v(i,j)
              bv= cy*v(i,j)
              cv= 0.
	   else
              av= 0.
              bv=-cy*v(i,j)
              cv= cy*v(i,j)
	   endif
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i-1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +pex*(vort(i-1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
30      continue

      return
      end
c----------------------------------------------------------------------
      subroutine stepupwind
      include 'step.h'
	vincmax=0.
        do 10 i=1,ixc1
        do 10 j=1,iy1-1
	   if(u(i,j).ge.0.)then
              au=-cx*u(i-1,j)
              bu= cx*u(i,j)
              cu= 0.
	   else
              au=0 .
              bu=-cx*u(i,j)
              cu= cx*u(i+1,j)
	   endif
	   if(v(i,j).ge.0.)then
              av=-cy*v(i,j-1)
              bv= cy*v(i,j)
              cv= 0.
	   else
              av= 0.
              bv=-cy*v(i,j)
              cv= cy*v(i,j+1)
	   endif
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i+1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
10         continue
        do 20 i=ixc1+1,ix1-1
        do 20 j=-iy1+1,channeltype*iy1-1
	   if(u(i,j).ge.0.)then
              au=-cx*u(i-1,j)
              bu= cx*u(i,j)
              cu= 0.
	   else
              au=0 .
              bu=-cx*u(i,j)
              cu= cx*u(i+1,j)
	   endif
	   if(v(i,j).ge.0.)then
              av=-cy*v(i,j-1)
              bv= cy*v(i,j)
              cv= 0.
	   else
              av= 0.
              bv=-cy*v(i,j)
              cv= cy*v(i,j+1)
	   endif
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i+1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +pex*(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
20      continue
	i=ix1
        do 30 j=-iy1+1,channeltype*iy1-1
	   if(u(i,j).ge.0.)then
              au=-cx*u(i-1,j)
              bu= cx*u(i,j)
              cu= 0.
	   else
              au=0 .
              bu=-cx*u(i,j)
              cu= cx*u(i-1,j)
	   endif
	   if(v(i,j).ge.0.)then
              av=-cy*v(i,j-1)
              bv= cy*v(i,j)
              cv= 0.
	   else
              av= 0.
              bv=-cy*v(i,j)
              cv= cy*v(i,j+1)
	   endif
	      vincrement(i,j)=
     +         -au*vort(i-1,j)-(bu+bv)*vort(i,j)-cu*vort(i-1,j)
     +         -av*vort(i,j-1)-cv*vort(i,j+1)
     +         +pex*(vort(i-1,j)-2*vort(i,j)+vort(i-1,j))
     +         +pey*(vort(i,j-1)-2*vort(i,j)+vort(i,j+1))
              newvort(i,j)=vort(i,j)+vincrement(i,j)
	temp=abs(vincrement(i,j))
	if(temp.gt.vincmax)vincmax=temp
30      continue

      return
      end
c------------------------------------------------------------------
      subroutine stepquick
      include 'step.h'
        vincmax=0.
        do 10 i=1,ixc1
        do 10 j=1,iy1-1
              call quick1(i,j,vinc)
              newvort(i,j)=vort(i,j)+vinc
              temp=abs(vinc)
              if(temp.gt.vincmax)vincmax=temp
10      continue
        do 20 i=ixc1+1,ix1-1
        do 20 j=-iy1+1,channeltype*iy1-1
              call quick2(i,j,vinc)
              newvort(i,j)=vort(i,j)+vinc
              temp=abs(vinc)
              if(temp.gt.vincmax)vincmax=temp
20      continue
        i=ix1
        do 30 j=-iy1+1,channeltype*iy1-1
              call quick3(i,j,vinc)
              newvort(i,j)=vort(i,j)+vinc
              temp=abs(vinc)
              if(temp.gt.vincmax)vincmax=temp
30      continue

      return
      end
c----------------------------------------------------------------------
        subroutine quick1(i,j,vortinc)
        include "step.h"
c in inlet section, i<= ixc1
        if(i.eq.1)then
                wleft=vort(0,j)
        else
                wleft=vort(i-2,j)
        endif
        if(j.eq.1)then
                wsouth=2*vort(i,j-1)-vort(i,j)
        else
                wsouth=vort(i,j-2)
        endif
        if(j.eq.(iy1-1))then
                wnorth=2*vort(i,j+1)-vort(i,j)
        else
                wnorth=vort(i,j+2)
        endif
        u0=u(i,j)
        v0=v(i,j)
        cox=u0*cx
        coy=v0*cy
        cox2=cox*cox
        coy2=coy*coy
        if(u0.ge.0.)then
                wxxx=(vort(i+1,j)-3*vort(i,j)+3*vort(i-1,j)-wleft)/6
        else
                wxxx=(vort(i+2,j)-3*vort(i+1,j)+3*vort(i,j)-vort(i-1,j))/6
        endif
        if(v0.ge.0.)then
                wyyy=(vort(i,j+1)-3*vort(i,j)+3*vort(i,j-1)-wsouth)/6
        else
                wyyy=(wnorth-3*vort(i,j+1)+3*vort(i,j)-vort(i,j-1))/6
        endif
        wxxy=((vort(i+1,j+1)-2*vort(i,j+1)+vort(i-1,j+1))-
     +        (vort(i+1,j-1)-2*vort(i,j-1)+vort(i-1,j-1)))/4

        wxyy=((vort(i+1,j+1)-2*vort(i+1,j)+vort(i+1,j-1))-
     +        (vort(i-1,j+1)-2*vort(i-1,j)+vort(i-1,j-1)))/4
        wxx=(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))/2
        wyy=(vort(i,j+1)-2*vort(i,j)+vort(i,j-1))/2
        wxy=(vort(i+1,j+1)-vort(i+1,j-1)-vort(i-1,j+1)+vort(i-1,j-1))/4
        wx=(vort(i+1,j)-vort(i-1,j))/2-wxxx
        wy=(vort(i,j+1)-vort(i,j-1))/2-wyyy
        vortinc=-cox*wx-coy*wy+(2*pex+cox2)*wxx+(2*pey+coy2)*wyy+
     +       cox*coy*wxy-(6*pex+cox2)*cox*wxxx-coy*(2*pex+cox2)*wxxy-
     +       cox*(2*pey+coy2)*wxyy-(6*pey+coy2)*coy*wyyy
        return
        end
        subroutine quick2(i,j,vortinc)
        include "step.h"
c in channel section, i> ixc1
        wleft=vort(i-2,j)
        if(i.eq.(ixc1+1).and.((j.lt.0).or.(j.gt.iy1)))then
                wleft=2*vort(i-1,j)-vort(i,j)
        endif
        if(i.lt.(ix1-1))then
                wright=vort(i+2,j)
        else
                wright=vort(i,j)
        endif

        if(j.eq.(1-iy1))then
                wsouth=2*vort(i,j-1)-vort(i,j)
        else
                wsouth=vort(i,j-2)
        endif
        if(j.eq.(channeltype*iy1-1))then
                wnorth=2*vort(i,j+1)-vort(i,j)
        else
                wnorth=vort(i,j+2)
        endif
        u0=u(i,j)
        v0=v(i,j)
        cox=u0*cx
        coy=v0*cy
        cox2=cox*cox
        coy2=coy*coy
        if(u0.ge.0.)then
                wxxx=(vort(i+1,j)-3*vort(i,j)+3*vort(i-1,j)-wleft)/6
        else
                wxxx=(wright-3*vort(i+1,j)+3*vort(i,j)-vort(i-1,j))/6
        endif
        if(v0.ge.0.)then
                wyyy=(vort(i,j+1)-3*vort(i,j)+3*vort(i,j-1)-wsouth)/6
        else
                wyyy=(wnorth-3*vort(i,j+1)+3*vort(i,j)-vort(i,j-1))/6
        endif
        wxxy=((vort(i+1,j+1)-2*vort(i,j+1)+vort(i-1,j+1))-
     +        (vort(i+1,j-1)-2*vort(i,j-1)+vort(i-1,j-1)))/4

        wxyy=((vort(i+1,j+1)-2*vort(i+1,j)+vort(i+1,j-1))-
     +        (vort(i-1,j+1)-2*vort(i-1,j)+vort(i-1,j-1)))/4
        wxx=(vort(i+1,j)-2*vort(i,j)+vort(i-1,j))/2
        wyy=(vort(i,j+1)-2*vort(i,j)+vort(i,j-1))/2
        wxy=(vort(i+1,j+1)-vort(i+1,j-1)-vort(i-1,j+1)+vort(i-1,j-1))/4
        wx=(vort(i+1,j)-vort(i-1,j))/2-wxxx
        wy=(vort(i,j+1)-vort(i,j-1))/2-wyyy
        vortinc=-cox*wx-coy*wy+(2*pex+cox2)*wxx+(2*pey+coy2)*wyy+
     +       cox*coy*wxy-(6*pex+cox2)*cox*wxxx-coy*(2*pex+cox2)*wxxy-
     +       cox*(2*pey+coy2)*wxyy-(6*pey+coy2)*coy*wyyy
        return
        end
        subroutine quick3(i,j,vortinc)
        include "step.h"
c in channel section, i=ix1
        wleft=vort(i-2,j)
        wright=vort(i-2,j)

        if(j.eq.(1-iy1))then
                wsouth=2*vort(i,j-1)-vort(i,j)
        else
                wsouth=vort(i,j-2)
        endif
        if(j.eq.(channeltype*iy1-1))then
                wnorth=2*vort(i,j+1)-vort(i,j)
        else
                wnorth=vort(i,j+2)
        endif
        u0=u(i,j)
        v0=v(i,j)
        cox=u0*cx
        coy=v0*cy
        cox2=cox*cox
        coy2=coy*coy
        if(u0.ge.0.)then
                wxxx=(vort(i-1,j)-3*vort(i,j)+3*vort(i-1,j)-wleft)/6
        else
                wxxx=(wright-3*vort(i-1,j)+3*vort(i,j)-vort(i-1,j))/6
        endif
        if(v0.ge.0.)then
                wyyy=(vort(i,j+1)-3*vort(i,j)+3*vort(i,j-1)-wsouth)/6
        else
                wyyy=(wnorth-3*vort(i,j+1)+3*vort(i,j)-vort(i,j-1))/6
        endif
        wxxy=((vort(i-1,j+1)-2*vort(i,j+1)+vort(i-1,j+1))-
     +        (vort(i-1,j-1)-2*vort(i,j-1)+vort(i-1,j-1)))/4

        wxyy=((vort(i-1,j+1)-2*vort(i-1,j)+vort(i-1,j-1))-
     +        (vort(i-1,j+1)-2*vort(i-1,j)+vort(i-1,j-1)))/4
        wxx=(vort(i-1,j)-2*vort(i,j)+vort(i-1,j))/2
        wyy=(vort(i,j+1)-2*vort(i,j)+vort(i,j-1))/2
        wxy=(vort(i-1,j+1)-vort(i-1,j-1)-vort(i-1,j+1)+vort(i-1,j-1))/4
        wx=(vort(i-1,j)-vort(i-1,j))/2-wxxx
        wy=(vort(i,j+1)-vort(i,j-1))/2-wyyy
        vortinc=-cox*wx-coy*wy+(2*pex+cox2)*wxx+(2*pey+coy2)*wyy+
     +       cox*coy*wxy-(6*pex+cox2)*cox*wxxx-coy*(2*pex+cox2)*wxxy-
     +       cox*(2*pey+coy2)*wxyy-(6*pey+coy2)*coy*wyyy
        return
        end

c------------------------------------------------------------------
      subroutine calculatepsi
      include 'step.h'
         do 10 i=0,ix1
            do 10 j=-iy1,channeltype*iy1

               r1(i,j)=0.-newvort(i,j)
10       continue
        call mgrid(psi,r1,residual)
      return
      end
c----------------------------------------------------------------------
      subroutine calculateboundaryvorticity()
      include 'step.h'
	real *8 gxc,gyc
	 gy=2./(dy*dy)
	 gx=2./(dx*dx)
         gxc=2.*gx/3.
	 gyc=2.*gy/3.
         do 20 i=1,ixc1-1
            	newvort(i,0)=gy*(psi(i,0)-psi(i,1))
         	newvort(i,iy1)=gy*(psi(i,iy1)-psi(i,iy1-1))
20       continue
	 i=ixc1
         newvort(i,0)=gyc*(psi(i,0)-psi(i,1))+gxc*(psi(i,0)-psi(i+1,0))
	if(channeltype.eq.1)then
		gyc=gy
		gxc=0.
	endif
         newvort(i,iy1)=
     +      gyc*(psi(i,iy1)-psi(i,iy1-1))+gxc*(psi(i,iy1)-psi(i+1,iy1))
	 do 30 j=-iy1+1,-1
         	newvort(i,j)=gx*(psi(i,j)-psi(i+1,j))
30	continue
	 do 35 j=iy1+1,channeltype*iy1-1
         	newvort(i,j)=gx*(psi(i,j)-psi(i+1,j))
35	continue
         do 40 i=ixc1+1,ix1
            	newvort(i,-iy1)=gy*(psi(i,-iy1)-psi(i,-iy1+1))
         	newvort(i,channeltype*iy1)=
     +            gy*(psi(i,channeltype*iy1)-psi(i,channeltype*iy1-1))
40       continue
      return
      end
c------------------------------------------------------------------
      subroutine calculatevelocities()
      include 'step.h'
      gx=1./dx
      gy=1./dy
	do 10 i=1,ixc1
         do 10 j=1,iy1-1
               u(i,j)=gy*(psi(i,j+1)-psi(i,j-1))/2
               v(i,j)=gx*(psi(i-1,j)-psi(i+1,j))/2
10          continue
	do 20 i=ixc1+1,ix1-1
         do 20 j=-iy1+1,channeltype*iy1-1
               u(i,j)=gy*(psi(i,j+1)-psi(i,j-1))/2
               v(i,j)=gx*(psi(i-1,j)-psi(i+1,j))/2
20          continue
	i=ix1
         do 30 j=-iy1+1,channeltype*iy1-1
               u(i,j)=gy*(psi(i,j+1)-psi(i,j-1))/2
               v(i,j)=gx*(psi(i-1,j)-psi(i-1,j))/2
30          continue
      return
      end
c--------------------------------------------------------------------
      subroutine setflux()
      include 'step.h'
         if(currenttime.lt.0.25)then
           f=sin(2*3.14159*currenttime)
         else
           f=mean+(1.-mean)*sin(2*3.14159*currenttime)
         endif
         do 20 i=0,ixc1
            psi(i,0)=0.-f
            psi(i,iy1)=f
20       continue
	do 30 j=-iy1,0
		psi(ixc1,j)=0.-f
30	continue
	do 35 j=iy1,channeltype*iy1
		psi(ixc1,j)=f
35	continue
	do 40 i=ixc1+1,ix1
		psi(i,-iy1)=0.-f
                psi(i,channeltype*iy1)=f
40	continue
       do 50 j=0,iy1
          y=j*dy
          u(0,j)=12*f*y*(1.-y)
	  v(0,j)=0.
          newvort(0,j)=-12*f*(1-2*y)
          psi(0,j)=f*(6*y*y-4*y*y*y-1)
50        continue

      return
      end
c--------------------------------------------------------------------

      subroutine mgrid(z1,r1,difference)
      include 'step.h'
         difference=1.
         k=0
         do while((difference.gt.tolerance).and.(k.lt.20))
            call gauss   (z1,r1,ixm1,iym1,iyp1,ix1,iy1,1,1)
            call residual(z1,r1,ixm1,iym1,iyp1,ix1,iy1,
     +                       r2,ixm2,iym2,iyp2,ix2,iy2,1  )
            call zero    (z2,   ixm2,iym2,iyp2,ix2,iy2    )
            call gauss   (z2,r2,ixm2,iym2,iyp2,ix2,iy2,1,2)
            call residual(z2,r2,ixm2,iym2,iyp2,ix2,iy2,
     +                       r3,ixm3,iym3,iyp3,ix3,iy3,2  )
            call zero    (z3,   ixm3,iym3,iyp3,ix3,iy3    )
	 if(iy1.gt.8)then
            call gauss   (z3,r3,ixm3,iym3,iyp3,ix3,iy3,1,4)
            call residual(z3,r3,ixm3,iym3,iyp3,ix3,iy3,
     +                       r4,ixm4,iym4,iyp4,ix4,iy4,4  )
            call zero    (z4,   ixm4,iym4,iyp4,ix4,iy4    )
            call gauss   (z4,r4,ixm4,iym4,iyp4,ix4,iy4,4,8)
            call interp  (z3,   ixm3,iym3,iyp3,ix3,iy3,
     +                    z4,   ixm4,iym4,iyp4,ix4,iy4    )
	 endif
            call gauss   (z3,r3,ixm3,iym3,iyp3,ix3,iy3,2,4)
            call interp  (z2,   ixm2,iym2,iyp2,ix2,iy2,
     +                    z3,   ixm3,iym3,iyp3,ix3,iy3    )
            call gauss   (z2,r2,ixm2,iym2,iyp2,ix2,iy2,2,2)
            call interp  (z1,   ixm1,iym1,iyp1,ix1,iy1,
     +                    z2,   ixm2,iym2,iyp2,ix2,iy2    )
            call gauss   (z1,r1,ixm1,iym1,iyp1,ix1,iy1,2,1)
            call diff    (z1,r1,ixm1,iym1,iyp1,ix1,iy1,difference)
            k=k+1
         end do
         if(k.eq.20)then
            write(screen,30)
30          format('warning: multigrid not converged')
         endif
      return
      end
c---------------------------------------------------------------------
      subroutine gauss(z,r,ixm,iym,iyp,ix,iy,repeats,level)
      implicit real*8(a-h,o-z)
	implicit integer (i-n)
      common/geometry/dx,dy,vfac,jv,ixc1,channeltype
      integer repeats,channeltype
      real *8 z(0:ixm,iym:iyp),r(0:ixm,iym:iyp)
	ixc=ixc1/level
	gx=1./(level*dx*level*dx)
	gy=1./(level*dy*level*dy)
	gxy=1./(2*(gx+gy))
         do 30 k=1,repeats
            do 10 i=1,ixc
               do 10 j=1,iy-1
                  z(i,j)=gxy*( gx*(z(i+1,j)+z(i-1,j))+
     +                         gy*(z(i,j+1)+z(i,j-1))-r(i,j) )
10             continue
	    do 20 i=ixc+1,ix-1
		do 20 j=-iy+1,channeltype*iy-1
                  z(i,j)=gxy*( gx*(z(i+1,j)+z(i-1,j))+
     +                         gy*(z(i,j+1)+z(i,j-1))-r(i,j) )
20          continue
	    i=ix
		do 25 j=-iy+1,channeltype*iy-1
                  z(i,j)=gxy*( gx*(z(i-1,j)+z(i-1,j))+
     +                         gy*(z(i,j+1)+z(i,j-1))-r(i,j) )
25          continue
30       continue
      return
      end
c----------------------------------------------------------------------
      subroutine zero(z,ixm,iym,iyp,ix,iy)
         implicit integer (i-n)
         real *8 z(0:ixm,iym:iyp)
	
         do 10 i=0,ix
            do 10 j=-iy,2*iy
               z(i,j)=0.
10       continue
      return
      end
c----------------------------------------------------------------------
      subroutine interp(z1,ixm1,iym1,iyp1,ix1,iy1,
     +                  z2,ixm2,iym2,iyp2,ix2,iy2)
         implicit integer (i-n)
         real *8 z1(0:ixm1,iym1:iyp1),z2(0:ixm2,iym2:iyp2)
	
         do 10 i=0,ix2
         do 10 j=-iy2+1,2*iy2-1
            z1(2*i,2*j)=z1(2*i,2*j)+z2(i,j)
10       continue
         do 20 i=1,ix1-1,2
         do 20 j=-iy1,2*iy1,2
            z1(i,j)=z1(i,j)+0.5*(z2((i+1)/2,j/2)+z2((i-1)/2,j/2))
20       continue
         do 30 i=0,ix1,2
         do 30 j=-iy1+1,2*iy1-1,2
            z1(i,j)=z1(i,j)+0.5*(z2(i/2,(j+1)/2)+z2(i/2,(j-1)/2))
30       continue
         do 40 i=1,ix1-1,2
         do 40 j=-iy1+1,2*iy1-1,2
            z1(i,j)=z1(i,j)+0.25*(z2((i+1)/2,(j+1)/2)
     +              +z2((i+1)/2,(j-1)/2)+z2((i-1)/2,(j+1)/2)
     +              +z2((i-1)/2,(j-1)/2))
40       continue
      return
      end
c-----------------------------------------------------------------------
      subroutine residual(z1,r1,ixm1,iym1,iyp1,ix1,iy1,
     +                       r2,ixm2,iym2,iyp2,ix2,iy2,level)
	implicit real *8(a-h,o-z)
	implicit integer (i-n)
       common/geometry/dx,dy,vfac,jv,ixc1,channeltype
      integer channeltype
	real *8 z1(0:ixm1,iym1:iyp1),r1(0:ixm1,iym1:iyp1),
     +                              r2(0:ixm2,iym2:iyp2)
	gx=1./(level*dx*level*dx)
	gy=1./(level*dy*level*dy)
         do 20 j=-iy1+2,2*iy1-2,2
           do 10 i=2,ix1-2,2
            r2(i/2,j/2)=r1(i,j)-gx*(z1(i+1,j)-2*z1(i,j)+z1(i-1,j))
     +                         -gy*(z1(i,j+1)-2*z1(i,j)+z1(i,j-1))
10         continue
            i=ix1
            r2(i/2,j/2)=r1(i,j)-gx*(z1(i-1,j)-2*z1(i,j)+z1(i-1,j))
     +                         -gy*(z1(i,j+1)-2*z1(i,j)+z1(i,j-1))
20       continue
      return
      end
c-----------------------------------------------------------------------
      subroutine diff(z1,r1,ixm1,iym1,iyp1,ix1,iy1,w)
      implicit real *8(a-h,o-z)
	implicit integer (i-n)
      common/geometry/dx,dy,vfac,jv,ixc1,channeltype
      integer channeltype

      real *8 z1(0:ixm1,iym1:iyp1),r1(0:ixm1,iym1:iyp1)

	gx=1./(dx*dx)
	gy=1./(dy*dy)
         w=0.
         d=0.
         do 10 i=1,ixc1
         do 10 j=1,iy1-1
             d= r1(i,j)-
     +          gx*( z1(i+1,j)-2*z1(i,j)+z1(i-1,j) )
     +         -gy*( z1(i,j+1)-2*z1(i,j)+z1(i,j-1) )
              w=w+abs(d)
10          continue
         do 20 i=ixc1+1,ix1-1
         do 20 j=-iy1+1,channeltype*iy1-1
             d= r1(i,j)-
     +          gx*( z1(i+1,j)-2*z1(i,j)+z1(i-1,j) )
     +         -gy*( z1(i,j+1)-2*z1(i,j)+z1(i,j-1) )
              w=w+abs(d)
20      continue
	i=ix1
         do 30 j=-iy1+1,channeltype*iy1-1
             d= r1(i,j)-
     +          gx*( z1(i-1,j)-2*z1(i,j)+z1(i-1,j) )
     +         -gy*( z1(i,j+1)-2*z1(i,j)+z1(i,j-1) )
              w=w+abs(d)
30      continue
              d=abs(z1(0,iy1)-z1(0,0))
            if(d.gt.0.)then
               w=w/(d*ix1*iy1)
            else
               w=w/(ix1*iy1)
            endif
      return
      end
c----------------------------------------------------------------------
      subroutine start()
      include 'step.h'

      keyboard =5
      screen   =6
      arcfile  =8
      logfile  =9
      inputfile=10
      sumfile  =11
      psifile  =12
      vortfile =13
      open(unit=inputfile,file='chan.dat',form='formatted')
      open(unit=logfile  ,file='chan.log',form='formatted')
      open(unit=psifile  ,file='chan.psi',form='formatted')
      open(unit=vortfile ,file='chan.vor',form='formatted')
      open(unit=sumfile  ,file='chan.sum',form='formatted')
c	if method>4 then start from archive with method-5 as method
	central		=0
	upwind		=1
	ncupwind	=2
	laxwendroff	=3
	quickest	=4
	finish		=0
	swapflow	=1
	changereynolds	=2
	changedt	=3
	changemethod	=4
	archiveflow	=5

      read(inputfile,1)channeltype
      read(inputfile,1)method
      read(inputfile,1)ix1
      read(inputfile,1)iy1
      read(inputfile,1)ixc1
      read(inputfile,1)ns1
      read(inputfile,1)ns2
      read(inputfile,2)reynolds
      read(inputfile,2)strouhal
      read(inputfile,2)tolerance
      read(inputfile,2)mean
      read(inputfile,2)dt
1     format(i5)
2     format(f10.3)


	ix2=ix1/2
	ix3=ix2/2
	ix4=ix3/2
	iy2=iy1/2
	iy3=iy2/2
	iy4=iy3/2

	dy=1./iy1
	dx=dy
	jv=6.38/dx
	vfac=(6.38-jv*dx)/dx
	jv=jv+ixc1
	cx=dt/(strouhal*dx)
	cy=dt/(strouhal*dy)
	pex=dt/(reynolds*strouhal*dx*dx)
	pey=dt/(reynolds*strouhal*dy*dy)
	if(channeltype.eq.1)then
		write(screen, 3)
		write(logfile,3)
	elseif(channeltype.eq.2)then
		write(screen,4)
		write(logfile,4)
	else
		write(screen,5)
		write(logfile,5)
	endif
3	format('Backstep')
4	format('Symmetric Expansion')
5	format('Panic: unknown channel type')
	if((method.eq.0).or.((method-5).eq.0))then
		write(screen,6)
		write(logfile,6)
	else if((method.eq.1).or.((method-5).eq.1))then
		write(screen,7)
		write(logfile,7)
	else if((method.eq.2).or.((method-5).eq.2))then
		write(screen,8)
		write(logfile,8)
	else if((method.eq.3).or.((method-5).eq.3))then
		write(screen,9)
		write(logfile,9)
	else if((method.eq.4).or.((method-5).eq.4))then
		write(screen,10)
		write(logfile,10)
	else 
		write(screen,11)
		write(logfile,11)
	endif
6	format('Method: central difference')
7	format('Method: upwind difference')
8	format('Method: nonconservative upwind difference')
9	format('Method: Lax Wendroff difference')
10	format('Method: Quickest difference')
11	format('Method: Panic: unkown method')



      write(screen,20)ns1,ns2,reynolds,strouhal,tolerance,mean,dt
      write(logfile,20)ns1,ns2,reynolds,strouhal,tolerance,mean,dt
20    format(i6,'steps then',i5,' steps per print',/,
     +     'reynolds number = ',f8.1,/,
     +     'strouhal number = ',f8.4,/,
     +     'tolerance       = ',f8.6,/,
     +     'mean flow       = ',f8.3,/,
     +     'time step       = ',f8.6)
      write(screen,30)dx*ix1,dx*(ix1-ixc1),cx,cy,pex,pey
      write(logfile,30)dx*ix1,dx*(ix1-ixc1),cx,cy,pex,pey
30    format('channel length',f10.2,/,
     +      'lee length    ',f10.2,/,
     +      'cx,cy         ',2f10.4,/,
     +      'pex,pey       ',2f10.4)
	if(method.ge.5)then
		method=method-5
		call unarchive()
	else
	        currenttime=0.
      		call initialiseall
	endif
      return
      end
c---------------------------------------------------------------------
c----------------------------------------------------------------------
      subroutine archive()
      include 'step.h'
      open(unit=arcfile  ,file='chan.arc',form='unformatted')
      write(arcfile)channeltype
      write(arcfile)method
      write(arcfile)ix1
      write(arcfile)iy1
      write(arcfile)ixc1
      write(arcfile)ns1
      write(arcfile)ns2
      write(arcfile)reynolds
      write(arcfile)strouhal
      write(arcfile)tolerance
      write(arcfile)mean
      write(arcfile)dt
      write(arcfile)currenttime
      write(arcfile)psi
      write(arcfile)oldvort
      write(arcfile)vort
       close(unit=arcfile)
      return
      end
c----------------------------------------------------------------------
      subroutine unarchive()
      include 'step.h'
	integer channeltypearc,methodarc,ix1arc,iy1arc,ixc1arc,
     +    ns1arc,ns2arc
	real*8 reynoldsarc,strouhalarc,tolerancearc,dtarc
     
c
c	want to be able to unarchive from coarser resolutions
c	but assuming global array sizes are same. Use oldvort
c	as dummy array in this case
c
      open(unit=arcfile  ,file='chan.arc',form='unformatted')
      read(arcfile)channeltypearc
      read(arcfile)methodarc
      read(arcfile)ix1arc
      read(arcfile)iy1arc
      read(arcfile)ixc1arc
      read(arcfile)ns1arc
      read(arcfile)ns2arc
      read(arcfile)reynoldsarc
      read(arcfile)strouhalarc
      read(arcfile)tolerancearc
      read(arcfile)meanarc
      read(arcfile)dtarc
      read(arcfile)currenttime
      	if(iy1arc.eq.iy1)then
		write(screen,5)
5		format('unarchiving normal data')
      		read(arcfile)psi
      		read(arcfile)oldvort
     		read(arcfile)vort
	else if(2*iy1arc.eq.iy1)then
		write(screen,8)
8		format('unarchiving from coarse data')

		read(arcfile)oldvort
c		1)set values from coarse array
		do 10 i=0,ix1arc
		do 10 j=-iy1arc,2*iy1arc
			psi(2*i,2*j)=oldvort(i,j)
10		continue
c		2)interpolate cell by cell
		do 20 i=0,ix1-2,2
		do 20 j=-iy1,2*iy1-2,2
			psi(i+1,j)=0.5*(psi(i,j)+psi(i+2,j))
			psi(i,j+1)=0.5*(psi(i,j)+psi(i,j+2))
			psi(i+1,j+1)=0.25*(psi(i,j)+psi(i+2,j)+
     +                                   psi(i,j+2)+psi(i+2,j+2))
20		continue
		do 30 i=0,ix1-2,2
			psi(i+1,2*iy1)=0.5*(psi(i,2*iy1)+psi(i+2,2*iy1))
30		continue
		do 40 j=-iy1,2*iy1-2,2
			psi(ix1,j+1)=0.5*(psi(ix1,j)+psi(ix1,j+2))
40		continue
		read(arcfile)oldvort
		read(arcfile)oldvort
c		1)set values from coarse array
		do 50 i=0,ix1arc
		do 50 j=-iy1arc,2*iy1arc
			vort(2*i,2*j)=oldvort(i,j)
50		continue
c		2)interpolate cell by cell
		do 60 i=0,ix1-2,2
		do 60 j=-iy1,2*iy1-2,2
			vort(i+1,j)=0.5*(vort(i,j)+vort(i+2,j))
			vort(i,j+1)=0.5*(vort(i,j)+vort(i,j+2))
			vort(i+1,j+1)=0.25*(vort(i,j)+vort(i+2,j)+
     +                                   vort(i,j+2)+vort(i+2,j+2))
60		continue
		do 70 i=0,ix1-2,2
			vort(i+1,2*iy1)=0.5*(vort(i,2*iy1)+vort(i+2,2*iy1))
70		continue
		do 80 j=-iy1,2*iy1-2,2
			vort(ix1,j+1)=0.5*(vort(ix1,j)+vort(ix1,j+2))
80		continue
		do 90 i=0,ix1
		do 90 j=-iy1,2*iy1
			oldvort(i,j)=vort(i,j)
90		continue
	else if(2*iy1.eq.iy1arc)then
		write(screen,92)
92	format('unarchiving fine data')

		read(arcfile)oldvort
		do 100 i=0,ix1
		do 100 j=-iy1,2*iy1
			psi(i,j)=oldvort(2*i,2*j)
100		continue
		read(arcfile)newvort
		do 110 i=0,ix1
		do 110 j=-iy1,2*iy1
			oldvort(i,j)=newvort(2*i,2*j)
110		continue
		read(arcfile)newvort
		do 120 i=0,ix1
		do 120 j=-iy1,2*iy1
			vort(i,j)=newvort(2*i,2*j)
120		continue
	else
		write(screen,95)ix1,iy1,ix1arc,iy1arc
95		format('PANIC in restorearc:',4i6)
		
	endif
	call calculatevelocities
       close(unit=arcfile)
      return
      end



