Program Two_D_ADI_CN_EXTRAP

!Algortithm description: 
!	The 2-D fractional diffusion equation is solved here by a Crank Nicholson method using the  
!	Alternating Directions Implicit (ADI) approach.
!	A shifted Grunwald estimate for the fractional order terms is used for stability.
!	The basic method is epxected to be second order in temporal stepsize and first order in
!		the spatial meshsize. After the integration to the Tend, a Richardson extrapolation
!		is carried out to gain second order spatial accuracy.
!----------
!Start Date:		08/16/05
!Last Edit:		
!By:				Charles Tadjeran
!Acknowlegement:	This project was supported by NSF grant DMS-0139927.
!----------
!The 2-D initial boundary value FADE being solved is : 
!	del(u)/del(t) = d(x,y,t) * del^alpha(u)/del(x^alpha) + 
!					e(x,y,t) * del^alpha(u)/del(y^beta)  + q(x,y,t)
!	over the rectangular region:  x_low < x < x_high, y_low < y < y_high,
!	and time limits: T_begin < t <= T_end
!	where u = u(x,y,t)
!	
!	The initial condition is specified by: u(x,y, t=T_begin) = f(x,y)
!	Dirichlet boundary conditions are assumed on all boundaries:
!	u(x=x_low,	y,	t)	= xbl(y,t) 
!	u(x=x_high,	y,	t)	= xbh(y,t)
!	u(x, y=y_low,	t)	= ybl(x,t) 
!	u(x, y=y_high,	t)	= ybh(x,t)
!
!	Additionally:
!	d = d(x,y,t) >=0 is the	x dispersive coefficient function
!	e = e(x,y,t) >=0 is the	y dispersive coefficient function 
!	q = q(x,y,t) >=0 is the	forcing (source/sink) term
!
!	It is assumed that
!	1<= alpha <= 2 , and  1<=  beta  <= 2	fractional order of the diffusion terms
!
!	Solved using an ADI method for Crank Nicolson discretization followed by (Richardson) extrapolation.
!
!
!User to define/set variables:
!		x_low	:low limit on x
!		x_high	:high limit on x	( x_low <=  x <= x_high )
!		y_low	:low limit on y
!		y_high	:high limit on y	( y_low <=  y <= y_high )
!
!		NX		:number of subintervals in x-direction
!		NY		:number of subintervals in y-direction
!		NT		:number of subintervals (or timesteps) in time
!		alpha, beta	:the fractional orders

IMPLICIT NONE

real	:: x_low, x_high, y_low, y_high
real	:: alpha, beta	!fractional derivative orders for the diffusive terms
integer	:: NX, NY		!number of subintervals in x and y directions
integer	:: NT			!number of subintervals in time

COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha, beta 		 

real    :: T_begin		!start integration time for t
real	:: T_end		!end integration time t : on T_begin <= t <= T_end

real	:: d			!x direction diffusion coefficient d= d(x,y,t)
real	:: e			!y direction diffusion coefficient e= e(x,y,t)
real	:: q			!forcing function in the ADE where q = q(x,y,t)

real	:: f			!initial condition u(x, y,t=T_begin) = f(x,y).

real	:: xbl			!x_low (Dirichlet) boundary function	u(x=x_low,y,t) = xbl(y,t)
real	:: xbh			!x_high (Dirichlet) boundary function	u(x=x_high,y,t) = xbh(t) 

real	:: ybl			!y_low (Dirichlet) boundary function	u(x,y=y_low,t) = ybl(x,t)
real	:: ybh			!y-high (Dirichlet) boundary function	u(x,y=y_high,t) = ybh(x,t) 

real    :: del_x, del_y, del_t	! set from user inputs NX, NY, NT

integer :: npts = 401  !max number of x or y subintervals.
!	***RUNTIME CONSTRAINT: NX and NY to be <= npts. For larger values change dimensions below
real    :: x(0:400), y(0:400), u(0:400, 0:400)	!x, y and solution arrays 
real	:: A(0:400,0:400), rhs(0:400)			!coefficient matrix & high hand side vector of equations		!
real    :: Galpha_wt(0:400), Gbeta_wt(0:400)	!array of Grunwald weights- computed once here
real    :: G_wt_scale_x, G_wt_scale_y			!basically 0-th Grunwald weights

!some local working variables
integer	:: i, j, k, nxfix, nyfix, nstep
real	:: tn, tnm1, t_midpt, xterm_wt,yterm_wt, x_diffuse_wt, y_diffuse_wt, ratio, temp 
real	:: u_exact, sln_err, abs_sln_err_coarse, abs_sln_err_fine
real    :: u_temp(0:400, 0:400)					!temp array of intermediate solution values prior to updates

real    :: u_ext(0:200, 0:200),u_coarse(0:200, 0:200)			!extrapolates solution, solution on coarse grid
real	:: deriv_term							!temp derivative term computation

integer	::	NX_input, NY_input, ngrid
real	::	extrap_sln_err, abs_extrap_sln_err

integer	:: checkstatus		!file open/close status report
character(40)		:: savefilename !output file containing sim info and output at extraction well

real    :: almost_zero = 1.0E-6 ! ~machine precision (for 32-bit word)

COMMON/TEST/choose		 
	integer	:: choose

!*******Begin USER NOTE 1
!Begin SETTING UP FOR USER INPUTS below
!Read prompts for the following user inputs: x_low, x_high, y_low, y_high, T_begin, T_end, alpha, beta 

!write (*,*) "choose=?"
!read (*,*) choose

!write (*,*) "alpha, beta?"
!read (*,*) alpha, beta
!write (*,*) "Input values for  x_low, x_high, y_low, y_high ?"
!read (*,*) x_low, x_high, y_low, y_high	
!write (*,*) "Input values for  T_begin, T_end?"
!read (*,*) T_begin, T_end

!FOR GENERAL CASES, un-comment the above user read and write prompt statement, and comment out the following section 
!		which was designed to easily provide for testing and minimizes typing many inputs 
!		during each test session
!BEGINNING of TEST SESSION INPUTS -->
	choose =1
	alpha	= 1.8
	beta	= 1.6
	x_low	= 0.0
	x_high	= 1.0
	y_low	= 0.0
	y_high	= 1.0
	T_begin = 0.0
	T_end	= 1.0

!END OF TEST SESSION INPUTS  <-- which should be commented out if a new test case is constructed
!*******End  USER NOTE 1



write (*,*) "Input values for  NX, NY, NT ?"
read (*,*) NX_input, NY_input, NT


write (*,*) "Input a FILENAME to save the solution at Tend"
read(*,*) savefilename		
open(UNIT = 20, IOSTAT = checkstatus, FILE =savefilename, STATUS = "NEW")
!write(*,*) "File Open IOSTAT =", checkstatus
Write(UNIT=20,FMT=*) "alpha    =", alpha
Write(UNIT=20,FMT=*) "beta     =", beta

Write(UNIT=20,FMT=*) "x_low   =", x_low
Write(UNIT=20,FMT=*) "x_high  =", x_high
Write(UNIT=20,FMT=*) "NX       =", NX
Write(UNIT=20,FMT=*) "y_low   =", y_low
Write(UNIT=20,FMT=*) "y_high  =", y_high
Write(UNIT=20,FMT=*) "NY       =", NY
Write(UNIT=20,FMT=*) "T_begin  =", T_begin
Write(UNIT=20,FMT=*) "T_end    =", T_end
Write(UNIT=20,FMT=*) "NT       =", NT

do ngrid = 1,2		!do on two gridsizes, the second with twice the number of subintervals in x and y directions
					!ngrid=1 is the coarse grid, and ngrid=2 is the fine grid
	NX = NX_input * ngrid
	NY = NY_input * ngrid

!run some quick checks (not exhaustive)
if ((T_begin >= T_end) .OR. (NT <1) .OR. (NX <2) .OR. (NY <2) .OR. (x_low >= x_high) .OR. (y_low >= y_high) .OR. (NX>npts) .OR. (NY>npts) ) then
	Write (*,*) "*** Invalid user input. Hit <Enter> to terminate!"
	read (*,*)
	stop
endif

!Initialize the discretized problem settings
del_t = (T_end - T_begin)/NT
del_x = (x_high - x_low) / real(NX)
del_y = (y_high - y_low) / real(NY)		

!Set the gridpoints
do i = 0 , NX
    x(i)    = x_low + real(i) * del_x
enddo
do i = 0 , NY
    y(i)    = y_low + real(i) * del_y
enddo

!Set the initial condition on the gridpoints
do i = 0 , NX
	do j = 0 , NY
		u(i, j)	= f(x(i), y(j))	!set initial condition at time t = T_begin
	end do
end do

if (del_t <= almost_zero) then		!note NT is user input
    write (*,*) "Terminating integration: starting timestep is too small. Hit <Enter> to terminate!"
	read (*,*)
    stop
else
	write (*,*) " Timestep = ",del_t
end if


!We'll use normalized (ie, with 1/del_x**alpha factored out) Grunwald weights:
!The x-direction
Galpha_wt(0) = 1.0
do i = 0 , NX-1
    Galpha_wt(i+1) = - Galpha_wt(i) * (alpha - i)/real(i+1)
end do

!The y-direction
Gbeta_wt(0) = 1.0
do i = 0 , NY-1
    Gbeta_wt(i+1) = - Gbeta_wt(i) * (beta - i)/real(i+1)
end do

!  !
!For each tn, solve the resulting system: first sweep in x-direction, then in y-direction
!
!First compute some recurring weights to ease computing
G_wt_scale_x	= 1.0/(del_x**alpha)		!same as G_wt(0),factor used in Grunwald weights
G_wt_scale_y	= 1.0/(del_y**beta)	
xterm_wt		= del_t * G_wt_scale_x /2.0
yterm_wt		= del_t * G_wt_scale_y /2.0

nstep = 0
tn= T_begin

do while ( nstep < NT)		!main time step forward move: 2 time-steps thru this loop each time
!XXXXXXXXXXXXXXXXXXXXX-DIRECTION
!Note that the following is similar to Peaceman-Rachford formulation for classical CN-ADI
	nstep = nstep+1
	tnm1 = tn
	tn = T_begin + nstep * del_t
	t_midpt = 0.5 * (tnm1+tn)
	do nyfix = 1, NY-1
!		Set up matrix equations for the finite differences: (I-A) U_j^* = (I+B)U_j^n + delt * Q_j^n+1/2 = P_n+1 ..
!		First the NX x NX matrix coeffs for the super-triangular matrix in x-direction
		do i = 1 , NX-1			!compute matrix of coeffs for row i (x-direction implicit)
			x_diffuse_wt = d(x(i), y(nyfix),t_midpt) * xterm_wt			!coeffs evaluated at midpt of time interval
			do k = 0, i+1		!column k
				A(i,k)	= - x_diffuse_wt * Galpha_wt(i-k+1)
			enddo		!do k
			A(i,i) = A(i,i) + 1.0

			y_diffuse_wt = e(x(i), y(nyfix),t_midpt) * yterm_wt
			rhs(i) =  0.0	!initialize the right hand side, rhs, of the equation at x_i
			do j = 0, nyfix+1		!column k
				rhs(i) =  rhs(i) + Gbeta_wt(nyfix-j+1) * u(i,j) 	!carry out the frac deriv sum
			enddo
			rhs(i) = u(i,nyfix) + y_diffuse_wt * rhs(i) + 0.5 * del_t * q(x(i),y(nyfix),t_midpt)  !note the 0.5 weight for the source term q

		enddo			!do i

!The Dirichlet boundary condition at x=x_low for the intermediate solution
		A(0,0)		= 1.0
		A(0,1)		= 0.0
!!		see the value of rhs(NX) below for how rhs(0) is "actually" computed at x=left boundary
		rhs(0)		= 0.0			!**zero boundary ===> intermediate sln

!The Dirichlet boundary condition  at x=x_high for the  intermediate solution U^*
		A(NX,NX-1)	= 0.0
		A(NX,NX)	= 1.0

		rhs(NX)		= 0.0
		do j = 0, nyfix+1		!column k
			rhs(NX) = rhs(NX) + Gbeta_wt(nyfix-j+1) * (xbh(y(j),tnm1) - xbh(y(j),tn)) 
		enddo
		y_diffuse_wt= e(x(NX), y(nyfix),t_midpt) * yterm_wt
		rhs(NX)		= 0.5 * ( y_diffuse_wt * rhs(NX) + xbh(y(nyfix),tnm1) + xbh(y(nyfix),tn) ) !same as u_(NX,nyfix)^*

!	Solve the resulting matrix system - note A is"sum of a lower triangular and a super-diagnoal matrices
!
!	First do a backward sweep to make it a lower triangular system
!
!	First subtract a multiple of the last row with two non-zero entries	
	ratio			= A(NX-1,NX) / A(NX,NX)	
	A(NX-1,NX-1)	= A(NX-1,NX-1)- ratio * A(NX,NX-1)
	rhs(NX-1) = rhs(NX-1) - ratio * rhs(NX)
!	Remaining rows are not as sparse as the last row
	do i = NX-1, 2, -1
		if (abs(A(i,i)) <= almost_zero) then		!diagonal entry almost zero
			write(*,*) " Singular coeff matrix, at i=", i, " Press <Enter> to terminate! "
			read (*,*) 
			stop
		else
			ratio =	 A(i-1,i) / A(i,i)	
		endif

		do j = 0, i-1
			A(i-1,j) = A(i-1,j) - ratio * A(i,j)
		enddo	!do j
		rhs(i-1) = rhs(i-1) - ratio * rhs(i)
	enddo	 !do i

!	Then a forward sweep of the triangularized system to get the solution
	u_temp(0,nyfix) = rhs(0)			!u_temp(0,nyfix) = Dirichlet boundary condition at x_low
	do i = 1, NX-1
		temp = 0.0
		do j = 0,i-1	
			temp = temp + A(i,j) * u_temp(j,nyfix)
		enddo	!done j
		u_temp(i,nyfix) = (rhs(i) - temp) /A(i,i)
	enddo	!done i

	u_temp(NX,nyfix) = rhs(NX) 		!u_temp(NX,nyfix) = Dirichlet boundary condition at x_high

!	Solution u(x, y ,t) is now available at this fixed y-point (nyfix). 
!	Increment nyfix and repeat to get next y-slice
	enddo		!nyfix

!	update to the intermediate solution, and then alternate the direction:
	do i = 0,NX
		do j = 1, NY-1
			u(i,j) = u_temp(i,j)		!this u is the same as u^* (that is, the intermediate solution)
		enddo
!!		u(i,0)	= don't care		Note that u^* values are not used at y_low and at y_high (including the four corners)							
!!		u(i,NY)	= don't care		Note that u^* values are not used at y_low and at y_high (including the four corners)
	enddo

!YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY

!	Set up and solve in the alternate Y-direction
	do nxfix = 1, NX-1		!for each x slice

!		Set up the NY x NY  super-triangular coeffs matrix in Y-direction

		do j = 1 , NY-1			!compute matix coeffs for row i (y-direction implicit)
			y_diffuse_wt = e(x(nxfix), y(j),t_midpt) * yterm_wt		!for optimized code, save wt's above
			do k = 0, j+1		!column k
				A(j,k)	= - y_diffuse_wt * Gbeta_wt(j-k+1)
			enddo
			A(j,j) = A(j,j) + 1.0

!			and the right hand side of this equation
			x_diffuse_wt = d(x(nxfix), y(j),t_midpt) * xterm_wt		!for optimized code, save wt's above
			rhs(j) =  0.0		!initial the rhs at y_j for the following sum
			do i = 0, nxfix+1	!column k
				rhs(j) =  rhs(j) + Galpha_wt(nxfix-i+1) * u(i,j)	!carry out frac deriv the sum
			enddo
			rhs(j) = x_diffuse_wt * rhs(j) + u(nxfix,j) + 0.5 * del_t * q(x(nxfix), y(j),t_midpt)  !note the 0.5 weight for the source term q

		enddo

!And also the Dirichlet boundary conditions on the low/high boundaries(other entries =0)
		A(0,0)		= 1.0
		A(0,1)		= 0.0
		rhs(0)		= ybl(x(nxfix),tn)			!expecting a zero boundary condition

		A(NY,NY-1)	= 0.0
		A(NY,NY)	= 1.0
		rhs(NY)		= ybh(x(nxfix), tn)			!value of u^{n+1} at the right boundary

!	Solve the resulting matrix system - sum of a lower triangular & a super-diagnoal matrices
!	First do a backward sweep to make it a lower triangular system
!	Subtract a multiple of the last row with two non-zero entries	
	ratio			= A(NY-1,NY) / A(NY,NY)	
	A(NY-1,NY-1)	= A(NY-1,NY-1)- ratio * A(NY,NY-1)
	rhs(NY-1)		= rhs(NY-1) - ratio * rhs(NY)
!	Remaining rows are not as sparse as the last
	do i = NY-1, 2, -1
		if (abs(A(i,i)) <= almost_zero) then		!diagonal entry almost zero
			write(*,*) "Singular coeff matrix, at i=",i , "Press <Enter> to terminate"
			read(*,*) 
			stop
		else
			ratio =	 A(i-1,i) / A(i,i)	
		endif

		do j = 0, i-1
			A(i-1,j) = A(i-1,j) - ratio * A(i,j)
		enddo	!do j
		rhs(i-1) = rhs(i-1) - ratio * rhs(i)
	enddo	!do i

!		Then a forward sweep to get the solution
		u_temp(nxfix, 0) = rhs(0)			!u_temp(0,nyfix) = low Dirichlet boundary condition = 0
		do i = 1, NY-1
			temp = 0.0
			do j = 0,i-1	
				temp = temp + A(i,j) * u_temp(nxfix,j)
			enddo	!done j
			u_temp(nxfix,i) = (rhs(i) - temp) /A(i,i)
		enddo	!done i

		u_temp(nxfix,NY) = rhs(NY) 		!Note u_temp(nxfix,NY) = high Dirichlet boundary condition
!	Solution u(x, y ,t) now available at this fixed x-point (nxfix). 
	enddo		!nxfix: increment nxfix and repeat to get next x-slice

!	update to set the solution, u^{n+1},and then repeat the big "do while"
	do i = 0,NX
		do j = 1, NY-1
			u(i,j) = u_temp(i,j)	!The real McCoy! This u is the solution u^{n+1}
		enddo
		u(i,0)	= ybl(x(i),tn)		!and at the x boundaries, as these values were not updated in this sweep
		u(i,NY)	= ybh(x(i),tn)
	enddo
enddo		!do while nstep

!-----
Write(UNIT=20,FMT = *) " On ngrid =", ngrid

!Resutls from one gridsize are available now - save at Tend
if (ngrid == 1) then		!save the results for the coarse grid
	do i = 1, NX_input-1
		do j = 1, NY_input-1
			u_coarse(i,j) = u(i,j)
		enddo
	enddo
endif

enddo		!ngrid loop

!Now carry out spatial extrapolation
do i =1, NX_input-1
	do j =1, NY_input-1
		u_ext(i,j) = 2 * u(2*i, 2*j) - u_coarse(i,j)
	enddo
enddo


!Print the final results 
!*******Begin USER NOTE 2
! NOTE***_ Since for these test cases the exact solution is known, the following
!		code also computes the absolute solution error and writes the results to the output file.
!		For general problems, where the exact solution is not known, the part of the code below that 
!		computes and prints the error should be deleted/commented out,
abs_sln_err_coarse			= 0.0 
abs_sln_err_fine	= 0.0 
abs_extrap_sln_err	= 0.0


!Se the gird for the coarse  gridpoints to computee the coarse and the extrapolated solution errors 
del_x = (x_high - x_low) / real(NX_input)
del_y = (y_high - y_low) / real(NY_input)		
do i = 0 , NX_input
    x(i)    = x_low + real(i) * del_x
enddo
do i = 0 , NY_input
    y(i)    = y_low + real(i) * del_y
enddo


do i = 1 ,NX_input-1
	do j = 1, NY_input-1		!do j = 0, NY
		if (choose ==1) then		 
!			TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6;  & alpha = 1.8, beta = 1.6
			u_exact = exp(-tn) * (x(i)**3) * (y(j)**3.6)
		elseif (choose ==2) then
!			TEST CASE (2): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* x^3 * y^3.6,  & alpha=beta=2 , d=d(x,y), e=e(x,y) 
			u_exact = exp(-tn) * (x(i)**3) * (y(j)**3.6)
		elseif (choose ==3) then
!			TEST CASE (3): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2, d=const, e=const  
			u_exact=exp(-tn)* sin(x(i)) * sin(y(j))
		elseif (choose ==4) then
!			TEST CASE (4): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2 , d=const, e=const  
			u_exact=exp(-tn)* sin(x(i)) * sin(y(j))
		elseif (choose ==5) then		!time varying diffusion coeffs for classical case, d=d(t), e=e(t) 
!			TEST CASE (5): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
			u_exact=exp(-tn)* sin(x(i)) * sin(y(j))
		elseif (choose ==6) then		!!classical with time variying diffusion coeff & nonzero q
!			TEST CASE (6): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
			u_exact=exp(-tn)* sin(x(i)) * sin(y(j))
		endif

		sln_err				= abs(u_exact - u_coarse(i,j))
		abs_sln_err_coarse			= max(sln_err, abs_sln_err_coarse)
		
		sln_err				= abs(u_exact - u(2*i, 2*j))
		abs_sln_err_fine	= max(sln_err, abs_sln_err_fine)

		extrap_sln_err		= abs(u_exact - u_ext(i,j))
		abs_extrap_sln_err	= max(extrap_sln_err, abs_extrap_sln_err)
	enddo
enddo

Write(UNIT=20,FMT = *)   "Absolute error - coarse grid		= ", abs_sln_err_coarse
Write(UNIT=20,FMT = *)   "Absolute error - fine grid		= ", abs_sln_err_fine
Write(UNIT=20,FMT = *)  "Absolute Extarpolated Sln Error	= ", abs_extrap_sln_err
!*******End USER NOTE 2

Close(UNIT=20, IOSTAT= checkstatus)
write(*,*) "File Open IOSTAT =", checkstatus

END Program Two_D_ADI_CN_EXTRAP
		
!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function d(x, y, t)
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: y_low	!low limit on y
	real	:: y_high	!high limit on y
	real	:: alpha	!fractional derivative order for the x- diffusive term
	real	:: beta		!fractional derivative order for the y- diffusive term
	integer	:: NX,NY	!number of subintervals in x and y directions
	real	:: x, y, t, d

COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha, beta		 

COMMON/TEST/choose		 
	integer	:: choose

if (choose ==1) then
!	TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6,  & alpha = 1.8, beta = 1.6
	d = 0.18363375 * (x**2.8) * y		!gamma(2.2)/6 = 1.10180249088/6 =0.18363375	
elseif (choose ==2) then	!classical PDE with d=d(x,y), and e=e(x,y)
!	TEST CASE (2): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* x^3 * y^3.6,  & alpha=beta=2  
	d = 0.333333333 * x*x*y
elseif (choose ==3) then	!classical PDE with d=constant, e= constant
!	TEST CASE (3): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	d=0.5
elseif (choose == 4) then	!classical PDE with d=constant e= constant
!	TEST CASE (4): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	d=0.8
elseif (choose == 5) then		!classical with time variying diffusion coeff
!	TEST CASE (5): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	d= 0.5 * t * t + 0.25
elseif (choose == 6) then		!classical with time variying diffusion coeff & nonzero q
!	TEST CASE (6): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	d= 0.5 * t * t + 0.25
endif
		
End Function d

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function e(x, y, t)
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: y_low	!low limit on y
	real	:: y_high	!high limit on y
	real	:: alpha	!fractional derivative order for the x- diffusive term
	real	:: beta		!fractional derivative order for the y- diffusive term
	integer	:: NX,NY	!number of subintervals in x and y directions
	real	:: x, y, t, e

COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta		 

COMMON/TEST/choose		 
	integer	:: choose

if (choose ==1) then		 
!	TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6,  & alpha = 1.8, beta = 1.6
	e = 0.1494624672 * x * (y**2.6)		!2/gamma(4.6) = 2/13.3812858709  =
elseif (choose ==2) then
!	TEST CASE (2): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* x^3 * y^3.6,  & alpha=beta=2  
	e=  0.213675213675 *x*y*y		!1/(1.8*2.6)= 0.21367521367521367521367521367521
elseif (choose ==3) then
!	TEST CASE (3): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	e=0.5
elseif (choose ==4) then
!	TEST CASE (4): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	e=0.2
elseif (choose == 5) then
!	TEST CASE (5): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	e= 0.75 - 0.5 * t * t
elseif (choose == 6) then
!	TEST CASE (6): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	e= 0.75 - 0.25 * t * t
endif

End Function e

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function q(x, y, t)		!forcing function in the diff eqn, v used/added to set up test cases

	real	:: q
	real	:: x_low,  y_low	!low limit on x and y
	real	:: x_high, y_high	!high limit on x and y
	real	:: alpha, beta		!fractional derivative order for the x and y diffusive terms
	integer	:: NX, NY			!number of subintervals in x and y directions
	real	:: x, y, t, q

COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta		 

COMMON/TEST/choose		 
	integer	:: choose

if (choose ==1) then		 
!	TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6;  & alpha = 1.8, beta = 1.6
	q = -(1.0 + 2*x*y) * exp(-t) * (x**3) * (y**3.6)
elseif (choose ==2) then
!	TEST CASE (2): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* x^3 * y^3.6,  & alpha=beta=2  
	q= -(1.0+2.0*(x+y))* exp(-t) * (x**3) * (y**3.6)
elseif (choose ==3) then
!	TEST CASE (3): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	q=0.0
elseif (choose ==4) then
!	TEST CASE (4): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	q=0.0
elseif (choose ==5) then
!	TEST CASE (5): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	q=0.0
elseif (choose ==6) then
!	TEST CASE (6): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	q= 0.25 * t * t * exp(-t)* sin(x) * sin(y)
endif


End Function q

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function f(x, y)		!initial condition for u(x,y,t) at t=0
	real	:: f
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: y_low	!low limit on y
	real	:: y_high	!high limit on y
	real	:: alpha	!fractional derivative order for the x- diffusive term
	real	:: beta		!fractional derivative order for the y- diffusive term
	integer	:: NX,NY	!number of subintervals in x and y directions
	real	:: x, y, f

COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta		 

COMMON/TEST/choose		 
	integer	:: choose

if (choose ==1) then		 
!	TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6;  & alpha = 1.8, beta = 1.6
	f = (x**3) * (y**3.6)
elseif (choose ==2) then
!	TEST CASE (2): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* x^3 * y^3.6,  & alpha=beta=2  
	f = (x**3) * (y**3.6)
elseif (choose ==3) then
!	TEST CASE (3): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	f = sin(x) * sin(y)
elseif (choose ==4) then
!	TEST CASE (4): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	f = sin(x) * sin(y)
elseif (choose ==5) then
!	TEST CASE (5): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	f = sin(x) * sin(y)
elseif (choose ==6) then
!	TEST CASE (6): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	f = sin(x) * sin(y)
endif

End Function f

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function xbl(y, t)		!Dirichlet boundary condition at x_low ; u(x=x_low,y, t)
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: y_low	!low limit on y
	real	:: y_high	!high limit on y
	real	:: alpha	!fractional derivative order for the x- diffusive term
	real	:: beta		!fractional derivative order for the y- diffusive term
	integer	:: NX,NY	!number of subintervals in x and y directions
	real	:: y, t, xbl

COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta		 

COMMON/TEST/choose		 
	integer	:: choose

if (choose ==1) then		 
!	TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6;  & alpha = 1.8, beta = 1.6
	xbl = 0.0		! zero boundary function at x = x_low = 0
elseif (choose ==2) then
!	TEST CASE (2): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* x^3 * y^3.6,  & alpha=beta=2  
	xbl = 0.0		! zero boundary function at x = x_low = 0
elseif (choose ==3) then
!	TEST CASE (3): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	xbl = 0.0
elseif (choose ==4) then
!	TEST CASE (4): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	xbl = 0.0
elseif (choose ==5) then
!	TEST CASE (5): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	xbl = 0.0
elseif (choose ==6) then
!	TEST CASE (6): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	xbl = 0.0
endif


End Function xbl

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function xbh(y, t)		!Dirichlet boundary condition at x_high: u(x=x_high,y, t)
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: y_low	!low limit on y
	real	:: y_high	!high limit on y
	real	:: alpha	!fractional derivative order for the x- diffusive term
	real	:: beta		!fractional derivative order for the y- diffusive term
	integer	:: NX,NY	!number of subintervals in x and y directions
	real	:: y, t, xbh

COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta		 

COMMON/TEST/choose		 
	integer	:: choose

if (choose ==1) then
!	TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6;  & alpha = 1.8, beta = 1.6
	xbh = exp(-t) * (x_high**3) * (y**3.6)
elseif (choose ==2) then
!	TEST CASE (2): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* x^3 * y^3.6,  & alpha=beta=2  
	xbh = exp(-t) * (x_high**3) * (y**3.6)
elseif (choose ==3) then
!	TEST CASE (3): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	xbh = exp(-t) * sin(x_high) * sin(y)
elseif (choose ==4) then
!	TEST CASE (4): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	xbh = exp(-t) * sin(x_high) * sin(y)
elseif (choose ==5) then
!	TEST CASE (5): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	xbh = exp(-t) * sin(x_high) * sin(y)
elseif (choose ==6) then
!	TEST CASE (6): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	xbh = exp(-t) * sin(x_high) * sin(y)
endif

End Function xbh

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function ybl(x, t)		!lDirichlet boundary condition at y_low: u(x,y=y_low, t)
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: y_low	!low limit on y
	real	:: y_high	!high limit on y
	real	:: alpha	!fractional derivative order for the x- diffusive term
	real	:: beta		!fractional derivative order for the y- diffusive term
	integer	:: NX,NY	!number of subintervals in x and y directions
	real	:: x, t, ybl

COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta		 

COMMON/TEST/choose		 
	integer	:: choose

if (choose ==1) then
!	TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6;  & alpha = 1.8, beta = 1.6
	ybl = 0.0		! zero boundary function at y= y_low = 0
elseif (choose ==2) then
!	TEST CASE (2): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* x^3 * y^3.6,  & alpha=beta=2  
	ybl = 0.0		! zero boundary function at y= y_low = 0
elseif (choose ==3) then
!	TEST CASE (3): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	ybl = 0.0
elseif (choose ==4) then
!	TEST CASE (4): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	ybl = 0.0
elseif (choose ==5) then
!	TEST CASE (5): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	ybl = 0.0
elseif (choose ==6) then
!	TEST CASE (6): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	ybl = 0.0
endif



End Function ybl

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function ybh(x, t)		!Dirichlet boundary condition at y_high: u(x,y=y_high, t)
	real	:: x_low	!low limit on x
	real	:: x_high	!high limit on x
	real	:: y_low	!low limit on y
	real	:: y_high	!high limit on y
	real	:: alpha	!fractional derivative order for the x- diffusive term
	real	:: beta		!fractional derivative order for the y- diffusive term
	integer	:: NX,NY	!number of subintervals in x and y directions
	real	:: x, t, ybh

COMMON/equation_params/x_low, x_high, y_low, y_high, NX, NY, alpha , beta		 

COMMON/TEST/choose		 
	integer	:: choose

if (choose ==1) then
!	TEST CASE with exact solution u(x,y,t) = exp(-t)* x^3 * y^3.6;  & alpha = 1.8, beta = 1.6
	ybh =  exp(-t)* (x**3) * (y_high**3.6)
elseif (choose ==2) then
!	TEST CASE (2): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* x^3 * y^3.6,  & alpha=beta=2  
	ybh =  exp(-t)* (x**3) * (y_high**3.6)
elseif (choose ==3) then
!	TEST CASE (3): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	ybh = exp(-t) * sin(x) * sin(y_high)
elseif (choose ==4) then
!	TEST CASE (4): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	ybh = exp(-t) * sin(x) * sin(y_high)
elseif (choose ==5) then
!	TEST CASE (5): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	ybh = exp(-t) * sin(x) * sin(y_high)
elseif (choose ==6) then
!	TEST CASE (6): CLASSICAL PDE:solution u(x,y,t) = exp(-t)* sin(x) * sin(y),  & alpha=beta=2  
	ybh = exp(-t) * sin(x) * sin(y_high)
endif


End Function ybh