C
C   1D CIP  test program 
C		2002/09/11	M. YAMADA
C
C----------------------------------------
C  parameter input
C----------------------------------------
      program main
      include "parameter.f"

C     IM: GRID
C     f:  Value		fn:  next step of f
C     gf: df/dx		gfn: next step of gf
C     u:  velocity
      common /VALUE/  f(0:IM),  gf(0:IM)
      common /NEXT/   fn(0:IM), gfn(0:IM)
      common /VELO/   u(0:IM)
      common /param/  dx, dt 

      real  stime, etime, time
      integer  estep

C----------------------------------------
C  set initial value
C----------------------------------------
      time = 0.0
      etime = 5.0
      estep = int(etime/dt)

      dx = 2./(real(IN))
      dt = 0.2*dx
 
      call initial
C----------------------------------------
C  main loop
C----------------------------------------

c     do 100 i=0,  estep
      do 100 i=0,  200
        call boundary
        call advection
        call shift

        time = time + dt
  100 end do

        do 1000 j=0, in
          write(6,*) j, f(j),gf(j),u(j)
 1000   end do

C----------------------------------------
C
C----------------------------------------

      stop
      end

C============================================================
C   SUBROUTINES
C============================================================
C----------------------------------------
C  shift
C----------------------------------------
      subroutine shift
      include "parameter.f"
      common /VALUE/  f(0:IM),  gf(0:IM)
      common /NEXT/   fn(0:IM), gfn(0:IM)
      do 100 i=0, IM
         f(i)=fn(i)
         gf(i)=gfn(i)
  100 end do
      return
      end

C----------------------------------------
C  advection
C----------------------------------------
      subroutine advection
      include "parameter.f"
      common /VALUE/  f(0:IM),  gf(0:IM)
      common /NEXT/   fn(0:IM), gfn(0:IM) 
      common /VELO/   u(0:IM) 
      common /param/  dx, dt 

      do 100 i=1, IN

        xi = - u(i)*dt
        if(u(i).le.0.)then
          isgn=1
        else
          isgn=-1
        end if

C--------------------
C  eq. (6)
C--------------------
        A =  (gf(i)+gf(i+isgn))/(dx*dx)
     &      + 2.*(f(i)-f(i+isgn))/isgn/(dx*dx*dx)
C--------------------
C  eq. (7)
C--------------------
        B =  3.*(f(i+isgn)-f(i))/(dx*dx)
     &      - (2.*gf(i)+gf(i+isgn))/isgn/dx
C--------------------
C  eq. (10)
C--------------------
        fn(i)  = ((A*xi+B)*xi+gf(i))*xi+f(i)
C--------------------
C  eq. (11)
C--------------------
        gfn(i) = (3.*A*xi+2.*B)*xi+gf(i)
  100 end do
      return
      end

C----------------------------------------
C  subloutine initial 
C----------------------------------------
      subroutine initial
      include "parameter.f"
      common /VALUE/  f(0:IM), gf(0:IM)
      common /NEXT/   fn(0:IM), gfn(0:IM) 
      common /VELO/   u(0:IM)
      common /param/  dx, dt 

      do i=0, IN	
        f(i) = 0.0
        if(i.ge.0.2*IN .and. i.le.0.4*IN)then
          f(i) = 2.0
        endif

        u(i)=1.
      end do

      do i=1, IN	
        gf(i) = (f(i+1)-f(i-1))/(2.0*dx)
      end do

      return
      end

C----------------------------------------
C  subloutine boundary
C----------------------------------------
      subroutine boundary
      include "parameter.f"
      common /VALUE/  f(0:IM), gf(0:IM)
      common /NEXT/   fn(0:IM), gfn(0:IM) 
      common /VELO/   u(0:IM)
      common /param/  dx, dt 

      f(0)=f(IN)
      f(IM)=f(1)

      gf(0)=gf(IN)
      gf(IM)=gf(1)

      return
      end
