program advection_2d
c
c     Written by: Steven T Zalesak
c

c     Solve the equation d/dt rho + d/dx (rho vx) + d/dy (rho vy)= 0
c     vx = +d/dy (phi), vy = -d/dx (phi); phi a streamfunction
c     on a 2-D uniform mesh with periodic boundary conditions.
c
      parameter (mx = 14, my = 14, dx = 1.0, dy = 1.0)
      real x(mx), y(my)
      real phi(mx,my), vx(mx,my), vy(mx,my)
      real rho(mx,my), flux_x(mx,my), flux_y(mx,my)
c
c     read in Courant number, nsteps
c
      write(*,*)' enter Courant number ( <= 0.25), nsteps'
      read(*,*) crn, nsteps
c
c     set initial conditions
c
      write (9,*) mx,my,nsteps,crn
      call set_ic (x, y, phi, rho, vx, vy, dx, dy, mx, my)
      call output (0,'rho',rho, mx,my)
c
c     update solution nsteps
c
      do k=1,nsteps
c
c        compute time step
c
         call dtset (vx, vy, dx, dy, crn, mx, my, dt)
c
c        compute numerical fluxes
c
         do j=3,my-2
         do i=2,mx-2
            vxh = 0.5*(vx(i,j)+vx(i+1,j))
            flux_x(i,j) = vxh*0.5*(rho(i,j)+rho(i+1,j))
     1                    -0.5*abs(vxh)*(rho(i+1,j)-rho(i,j))
            flux_x(i,j) = flux_x(i,j)*dy
         enddo
         enddo
c
         do j=2,my-2
         do i=3,mx-2
            vyh = 0.5*(vy(i,j)+vy(i,j+1))
            flux_y(i,j) = vyh*0.5*(rho(i,j)+rho(i,j+1))
     1                    -0.5*abs(vyh)*(rho(i,j+1)-rho(i,j))
            flux_y(i,j) = flux_y(i,j)*dx
         enddo
         enddo
c
c        use fluxes to update solution
c
         fact = dt/(dx*dy)
         do j=3,my-2
         do i=3,mx-2
            rho(i,j) = rho(i,j) - fact*(flux_x(i,j) - flux_x(i-1,j)
     1                               +flux_y(i,j) - flux_y(i,j-1))
         enddo
         enddo
c
c        set boundary conditions
c
         call bcset (rho, mx,my)
         write(6,*) '  step number, dt = ', k, dt
         call output (k,'rho',rho, mx,my)
c
      enddo  !   k=1,nsteps
c
      stop
      end
c
c---------------------------------------------------------------
c
      subroutine set_ic (x, y, phi, rho, vx, vy, dx, dy, mx, my)
      real x(mx), y(my)
      real phi(mx,my), vx(mx,my), vy(mx,my)
      real rho(mx,my)
c
c     set initial conditions (2D square wave)
c
      do i=1,mx
         x(i) = (float(i) - 0.5*float(mx+1))*dx
      enddo
      do j=1,my
         y(j) = (float(j) - 0.5*float(my+1))*dy
      enddo
      write(*,*)'x=',x
      write(*,*)'y=',y
      do j=1,my
      do i=1,mx

         phi(i,j) = y(j)        !This creates uniform velocity along the
                                ! +x direction
         rho(i,j) = 0.0

                 !----- Other initial conditions to try....
                 !   rsq = x(i)**2+y(j)**2  ! 1.--This would create 
                 !   phi(i,j) = rsq         !     solid body rotation
                 !
                 !   phi(i,j) = -x(i)       ! 2.--This would create uniform
                 !                          !   velocity along the +y direction
      enddo
      enddo
      do i=3*mx/4,3*mx/4+1
      do j=my/2,my/2+1
         rho(i,j) = 1.0
      enddo
      enddo
      do j=3,my-2
      do i=3,mx-2
         vx(i,j) =  (phi(i,j+1)-phi(i,j-1))/(2.0*dy)
         vy(i,j) = -(phi(i+1,j)-phi(i-1,j))/(2.0*dx)
      enddo
      enddo
      call bcset (vx, mx,my)
      call bcset (vy, mx,my)
!      write(9,*)' vx'
      call output (0,'vx ',vx,mx,my)
!      write(9,*)' vy'
      call output (0,'vy ',vy,mx,my)
c
      return
      end
c
c---------------------------------------------------------------
c
      subroutine dtset (vx, vy, dx, dy, crn, mx, my, dt)
      real vx(mx,my), vy(mx,my)
      rdt = abs(vx(3,3))/dx
      do j=3,my-2
      do i=3,mx-2
         rdt = max( abs(vx(i,j))/dx, abs(vy(i,j))/dy, rdt)
      enddo
      enddo
      dt = crn/rdt
      return
      end
c
c---------------------------------------------------------------
c
      subroutine bcset (a, mx,my)
      real a(mx,my)
      do j=3,my-2
         a(1,j)    = a(mx-3,j)
         a(2,j)    = a(mx-2,j)
         a(mx-1,j) = a(3,j)
         a(mx,j)   = a(4,j)
      enddo
      do i=1,mx
         a(i,1)    = a(i,my-3)
         a(i,2)    = a(i,my-2)
         a(i,my-1) = a(i,3)
         a(i,my)   = a(i,4)
      enddo
      return
      end
c
c---------------------------------------------------------------
c
      subroutine output (istep,name,array, mx, my)
      real array(mx,my)
      character *3, name
      write(9,9002) istep,name
      do j=1,my
        write(9,9001) (array(i,j),i=1,mx)
      enddo
      return
 9001 format(1x,14f5.2)
 9002 format(i4,1x,a3)
      end