program opgave1
c
      implicit none
c
c solution of instationary diffusion equation with several methods
c
      integer N,i
c
      parameter (N=10)
      double precision Cold(0:N)
      double precision Cnew(0:N)
      double precision x(0:N)
      double precision dx
      double precision RHS(1:N)
      double precision time, dt
      double precision Kdif
      double precision length
      double precision U
      integer imethod
      integer istep, nstep
c
      open(unit=10,file='param')
      read(10,*)nstep
      read(10,*)dt
      read(10,*)imethod
      close(10)
      Kdif = 0.0001d0
      U    = 25.d0

      call grid(x,dx,N)
      call init(Cnew,x,N)
      call copy(Cold,Cnew,N)
      time = 0.d0
c
      
c
c check stability
c
         call checkst(dx,x,dt,U,Kdif,N)
      do istep=1,nstep
         call explicit(Cnew,Cold,RHS,x,dx,dt,U,Kdif,imethod,N)
         call bound(Cnew,N)
         time = time  + dt
      enddo
c
      call output(Cnew,x,U,Kdif, time,  N)
c
      stop
      end
 
      subroutine grid(x,dx,N)
c
      implicit none
      integer N
      double precision strech
      double precision x(0:N)
      double precision dx
      double precision length
c
      double precision xend
       integer i
c 
       length = 1.d0
       dx = length/N
       do i=0,N
          x(i) = i*dx
       enddo
c
       return
       end
c
       subroutine init(C,x,N)
c
       implicit none
c
       integer N, i
       double precision C(0:N)
       double precision x(0:N)
       double precision pi
       pi = 4.d0*atan(1.d0)
c
       do i=0,N
          C(i) = 0.d0
       enddo
c
       return
       end
       subroutine copy(Cold,Cnew,N)
c
       implicit none
c
       integer N
c
       double precision Cold(0:N)
       double precision Cnew(0:N)
c
       integer i
c
       do i=0,N
          Cold(i) = Cnew(i)
       enddo
c
       return
       end
       subroutine checkst(dx,x,dt,U,Kdif,N)
c
c implement your stability criterion in this routine
c
       implicit none
c
       integer N
       double precision dx
       double precision x (0:N)
       double precision dt, U,Kdif, Pe
       integer i
       double precision difnum, sounum
c
          difnum = 0.00001
          if(difnum.gt.0.5d0)then
             write(6,*)'diffusion number gt 0.5 ' , 
     &       'difnum = ', difnum
          endif
       return
       end

       subroutine explicit(Cnew,Cold,RHS,x,dx,dt,U,Kdif,imethod,N)

       implicit none
c
       integer N
c
      double precision Cold(0:N)
      double precision Cnew(0:N)
      double precision RHS (1:N)
      double precision x(0:N)
      double precision dx
      double precision dt, U,Kdif
      integer  imethod
c
      integer i
      if(imethod.eq.1)then
c
c PUT YOUR OWN STUFF HERE
      do i=0,N
         Cold(i) = 0.d0
      enddo
c
c PUT YOUR OWN STUFF HERE
      do i=1,N-1
         Cnew(i) = 0.d0
      enddo
      elseif(imethod.eq.2)then
c
c PUT YOUR OWN STUFF HERE
      do i=0,N
         Cold(i) = 0.d0
      enddo
c
c PUT YOUR OWN STUFF HERE
      do i=1,N-1
         Cnew(i) = 0.d0
      enddo
      else
         write(6,*)'ERROR: NO valid method chosen'
         write(6,*)'ERROR: imethod must be 1 or 2'
         stop
      endif
c
       return
       end
c
       subroutine bound(C,N)
c
       implicit none
c
       integer N
c
       double precision C(0:N)
c
c PUT YOUR OWN STUFF HERE
       C(0) = 0.d0
c PUT YOUR OWN STUFF HERE
       C(N) = 0.d0
c
       return
       end


       subroutine output(C,x,U,Kdif,time,N)
c
       implicit none
c
       integer N
c
       double precision C(0:N)
       double precision x(0:N)
       double precision U,Kdif, time
       double precision exact, pi, ermax
c
       integer i
       pi = 4.d0*atan(1.d0)
c
       open(unit=10,file='output')
c
       ermax = 0.d0
c
       do i=0,N
C PUT YOUR OWN STUFF HERE (calculate exact, calculate maximum absolute error)
          exact = 0.d0
          write(10,*)x(i),C(i),exact
       enddo
       write(6,*)'time, ermax = ', time, ermax
c
       close(10)
c
       return
       end