program swid_flow

	!Governing Flow Equation: d^2h/dx^2+d^2h/dy^2=S/T*(dh/dt)-R/T
	!Initial Condition: h(x,y,t=0)=-i*x+ho
	!Boundary Conditions: 
	!Stress Period 1: d^2h/dx^2+d^2h/dy^2=0
	!h(x=0,y,0<t<=t1)=ho
	!h(x=inf,y, 0<t<=t1)=-i*inf+ho
	!h(x,y=0, 0<t<=t1)=-i*x+ho
	!h(x,y=inf, 0<t<=t1)=-i*x+ho
	!Stress Period 2: d^2h/dx^2+d^2h/dy^2=S/T*(dh/dt)-R/T
	!Same as Stress Period 1 plus: R(x=well,y=well,t1<t<=t2)=Q/(dx*dy)
	!Stress Period 3: d^2h/dx^2+d^2h/dy^2=S/T*(dh/dt)
	!Same as Stress Period 2
	!Finite Difference Equation: h[i,j,n+1]=1/A*(A*h[i,j,n]+(hbar-h[i,j,n])+(delx*dely*R)/(4*T))
	!Where A=(delx*dely*S)/(4*T*delt) & hbar=(h[i+1,j,n]+h[i-1,j,n]+h[i,j+1,n]+h[i,j-1,n])/4
	!Note 1: i is x index, j is y index, n is time (t) index
	!Note 2: (T*delt)/(S*delx*dely)<0.25 stable fully explicit scheme
	!Note 3: Stress period 1 is over relaxed at 1.5 for more efficient Gauss-Seidel iteration
	!Variables from Hoss thesis: i=-0.03, Q=42.8ft^3/day, ti=0.3125days, K=28.3ft/day, b=5ft, S=0.25, ho[well]=100ft

	implicit none
	
	! Declare, Define, Compute, Output
	
	! Declare
	
	! Stress Period 1
	
	integer :: i, j, n
	integer,parameter :: row=51, col=51, pi=3.14159
	real,dimension(row,col) :: ihead, lhead, phead
	real :: delx, dely, grad, ho 
	integer :: numit, maxit
	real :: error, maxerr, converg, omega

	! Stress Period 2 and 3
	
	real,dimension(row,col) :: hold, hnew, R
	integer :: nend, kount1, kount2
	real :: time, A, hbar, stab
	real :: delt, K, b, T, S, Q, t1, t2, t3, Vi, Vo, Mb

	! Define
	
	delx=1.
	dely=delx
	delt=0.0001
	ho=100.+0.75
	grad=-0.03
	K=28.3/3.5
	b=5.
	T=K*b
	S=0.25
	Q=42.8
	t1=1.
	t2=t1+0.3125
	t3=t2+t1
	
	omega=1.5
	numit=0
	maxit=1000

	converg=0.00005
	error=0.
	maxerr=2.*converg
	
	! Stress period 1: no pumping, steady-state flow, Gauss-Seidel, over relaxation, d^2h/dx^2+d^2h/dy^2=0
	
	! Fill heads with ho as placeholders
	
	do i=1,row
		do j=1,col
			ihead(i,j)=ho
			lhead(i,j)=ihead(i,j)
			phead(i,j)=ihead(i,j)
		end do
	end do
	
	! Boundary conditions
	
	! Left is already set at ho
	
	! Right
	
	do i=1,row
		ihead(i,col)=grad*(delx*(col-1))+ho
		lhead(i,col)=ihead(i,col)
		phead(i,col)=ihead(i,col)
	end do
	
	! Top
	
	do j=1,col
		ihead(1,j)=grad*(delx*(j-1))+ho
		lhead(1,j)=ihead(1,j)
		phead(1,j)=ihead(1,j)
	end do	
	
	! Bottom
	
	do j=1,col
		ihead(row,j)=grad*(delx*(j-1))+ho
		lhead(row,j)=ihead(row,j)
		phead(row,j)=ihead(row,j)
	end do	
	
	! Compute
	
	do
	if(maxerr<=converg .or. numit>=maxit) exit
	numit=numit+1
	maxerr=0.
		do i=2,row-1
			do j=2,col-1
				phead(i,j)=lhead(i,j)
				lhead(i,j)=(1.-omega)*lhead(i,j)+omega*(lhead(i,j+1)+lhead(i,j-1)+lhead(i+1,j)+lhead(i-1,j))/4.
				error=abs(phead(i,j)-lhead(i,j))
					if(error>maxerr) then
						maxerr=error
					end if
			end do
		end do
	end do

	! Output
	
	write(*,*) "Heads (Initial)"
	write(*,10) ((ihead(i,j),j=1,row,5),i=1,col,5)
	10 format(1x,11f8.2)
	
	write(*,*) "Heads (Stress Period 1 [0-1 days])"
	write(*,20) ((lhead(i,j),j=1,row,5),i=1,col,5)
	20 format(1x,11f8.2)
	
	write(*,30) numit, maxerr
	30 format("  numit: ",i8,/," maxerr: ",f8.5/)

	open(unit=1000, file='swid_flow_out.txt') ! Open new or overwrite old output file
	
	write(1000,100) ((lhead(i,j),j=1,row),i=1,col)
	100 format(1x,51f10.4)
	
	! Stress period 2: injection pumping, transient flow with source, explicit scheme, d^2h/dx^2+d^2h/dy^2=S/T*(dh/dt)-R/T
	
	!Finite Difference Equation: h[i,j,n+1]=1/A*(A*h[i,j,n]+(hbar-h[i,j,n])+(delx*dely*R)/(4*T))
	!Where A=(delx*dely*S)/(4*T*delt) & hbar=(h[i+1,j,n]+h[i-1,j,n]+h[i,j+1,n]+h[i,j-1,n])/4
	
	! Fill old (time=n) and new heads (time=n+1) with heads from stress period 1, new heads will be overwritten as time marches on
	
	! Define
	
	do i=1,row
		do j= 1, col
			hold(i,j)=lhead(i,j)
			hnew(i,j)=lhead(i,j)
		end do
	end do	
	
	do i=1,row
		do j=1,col
			R(i,j)=0.
		end do
	end do

	R(26,26)=Q/(delx*dely) ! Turn ON injection pump
	
	nend=3125
	time=t1
	kount1=0
	kount2=11
	Vo=0.
	Vi=0.

	! Compute & Output

		write(*,*) "Heads with Time (Stress Period 2 [1-1.3125 days])"
	
		write(*,40)"Heads","Time"
		40 format(/,44x,a5,48x,a4/)

	do n=1,nend
		
		do i=2,row-1
			do j=2,col-1
				A=(delx*dely*S)/(4.*T*delt)
				hbar=(hold(i,j+1)+hold(i,j-1)+hold(i+1,j)+hold(i-1,j))/4.
				hnew(i,j)=1./A*(A*hold(i,j)+(hbar-hold(i,j))+(delx*dely*R(i,j))/(4.*T))
			end do
		end do
		
		Vo=((((hnew(26,26)-hnew(26,26+1)))/delx)*(K))*b*delx*delt&
		+((((hnew(26,26)-hnew(26,26-1)))/delx)*(K))*b*delx*delt&
		+((((hnew(26,26)-hnew(26+1,26)))/delx)*(K))*b*dely*delt&
		+((((hnew(26,26)-hnew(26-1,26)))/delx)*(K))*b*dely*delt&
		+Vo
		
		Vi=Q*delt+Vi
		
		do i=2,row-1
			do j=2,col-1
				hold(i,j)=hnew(i,j)
			end do	
		end do
		
		time=time+delt
		kount1=1+kount1
		
		if(kount1==kount2) then
	
			write(*,50) (hnew(26,j),j=1,col,5),time
			50 format(1x,11f8.2,1f12.4)
			
			write(1000,200) (hnew(26,j),j=1,col),time
			200 format(1x,51f10.4,1f12.4)
			
			kount1=0
			
		end if
		
	end do
	
	stab=(T*delt)/(S*delx*dely)
	
	Mb=(1.+(Vo-Vi)/Vi)*100.
	
	write(*,60) stab, Vi, Vo, Mb
	60 format("   stab: ",f8.5,/,"     Vi: ",f8.5/,"     Vo: ",f8.5/,"     Mb: ",f8.5/)
	
	! Stress period 3: stop pumping, transient flow with source decay, explicit scheme, d^2h/dx^2+d^2h/dy^2=S/T*(dh/dt)-0/T
	
	!Finite Difference Equation: h[i,j,n+1]=1/A*(A*h[i,j,n]+(hbar-h[i,j,n])+(delx*dely*R)/(4*T))
	!Where A=(delx*dely*S)/(4*T*delt) & hbar=(h[i+1,j,n]+h[i-1,j,n]+h[i,j+1,n]+h[i,j-1,n])/4
	
	! Fill old (time=n) and new heads (time=n+1) with heads from stress period 2, new heads will be overwritten as time marches on
	
	! Define
	
	do i=1,row
		do j= 1, col
			hold(i,j)=hnew(i,j)
			hnew(i,j)=hnew(i,j)
		end do
	end do	
	
	do i=1,row
		do j=1,col
			R(i,j)=0.
		end do
	end do

	R(26,26)=0.0 ! Turn OFF injection pump
	
	nend=10000-3125
	time=t2
	kount1=0
	kount2=21

	! Compute & Output

		write(*,*) "Heads with Time (Stress Period 3 [1.3125-2 days])"
	
		write(*,70)"Heads","Time" !
		70 format(/,44x,a5,48x,a4/)

	do n=1,nend
		
		do i=2,row-1
			do j=2,col-1
				A=(delx*dely*S)/(4.*T*delt)
				hbar=(hold(i,j+1)+hold(i,j-1)+hold(i+1,j)+hold(i-1,j))/4.
				hnew(i,j)=1./A*(A*hold(i,j)+(hbar-hold(i,j))+(delx*dely*R(i,j))/(4.*T))
			end do
		end do
				
		do i=2,row-1
			do j=2,col-1
				hold(i,j)=hnew(i,j)
			end do	
		end do
		
		time=time+delt
		kount1=1+kount1
		
		if(kount1==kount2) then
	
			write(*,80) (hnew(26,j),j=1,col,5),time
			80 format(1x,11f8.2,1f12.4)
			
			write(1000,300) (hnew(26,j),j=1,col),time
			300 format(1x,51f10.4,1f12.4)
			
			kount1=0
			
		end if
		
	end do
	
	stab=(T*delt)/(S*delx*dely)
	
	write(*,90) stab
	90 format("   stab: ",f8.5)
	
	close(1000)
	
end program swid_flow