program test_cn
	use crank_nicholson
	implicit none
	real :: D, h, t
	logical :: ld, rd
	real, dimension(401) :: u, unew
	real, dimension(2) :: ub
	integer :: ierr,i,j,counter,frame
	character(*),parameter :: filename_base = 'cn_data/step_'
	character(len=30) :: filename
	
	D = 1.0
	h = 0.05
	t = 0.000625
	ld = .true.
	rd = .true.
	call initialize_cn(D,h,ld,rd)
	
	ub = (/0.0,0.0/)
	call set_boundary_condition(ub)
	u = 0.0
	forall(i = 2:81) u(i) = 1.0
	unew = 0.0
	write(filename,'(a,i0,a)') filename_base,0,'.data'
	open(unit=10,file=trim(filename))
	write(10,'(f9.6)') 0.0
	write(10,'(f9.6,tr2,f9.6)') ((i-1)*h,u(i),i=1,401)
	close(10)
	open (unit=20,file='cn_data/lightcurve.data')
	counter = 1
	frame = 1
 	do j = 1,8000
		call step_cn(u,t,unew,ierr)
		if (counter == 32) then
			write(filename,'(a,i0,a)') filename_base,frame,'.data'
			open(unit=10,file=trim(filename))
			write(10,'(f9.6)') j*t
			write(10,'(f9.6,tr2,f9.6)') ((i-1)*h,unew(i),i=1,401)
			close(10)
			write(20,'(i4,tr2,f9.6,tr2,f9.6)') frame,j*t,D*(unew(2)-unew(1))/h
			counter = 0
			frame = frame + 1
		end if
		counter = counter+1
		u = unew
	end do
	close(20)
end program test_cn