Program Two_D_ADI_FPDE

!Algortithm description: 
!	The 2-D fractional diffusion equation is solved here by an Implicit Euler method with 
!	Alternating Directions Implicit (ADI) method..
!	A shifted Grunwald estimate for the fractional order terms is used for stability.
!	The algorithm applies the derivation in the "problem" statement of the paper
!	titled "Finite difference methods for two-dimensional fractional dispersion equation" 
!	to appear in J. of Computational Physics.
!----------
!Start Date:		10/03/04
!Last Edit:			07/04/05
!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

!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, xterm_wt,yterm_wt, x_diffuse_wt, y_diffuse_wt, ratio, temp, u_exact, sln_err, abs_sln_err
real    :: u_int(0:400, 0:400)					!temp array of intermediate solution values prior to updates
real	:: deriv_term							!temp derivative term computation

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)

!*******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 (*,*) "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 to reduce user input errors
!BEGINNING of TEST SESSION INPUTS -->
	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, NY, NT

NT = 2 * ((NT+1)/2)		!To make sure that we have an even number of timesteps
del_t = (T_end - T_begin)/NT

!run some quick checks (not exhaustive)
if ((T_begin >= T_end) .OR. (NT <2) .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


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

!Initialize the discretized problem settings
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
yterm_wt		= del_t * G_wt_scale_y

nstep = 0					!0th time step
do while ( nstep < NT)		!main time step forward move: 2 time-steps thru this loop each time
!XXXXXXXXXXXXXXXXXXXXX-DIRECTION
	nstep = nstep+1
	tn = T_begin + nstep * del_t
	do nyfix = 1, NY-1
!		Set up matrix equations for the finite differences: Ax U_j^n+1 = U_j^n + delt * Q_j^n+1 = P_n+1 ..
!		First the NX x NX 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),tn) * xterm_wt
!			y_diffuse_wt = e(x(i), y(nyfix),tn) * yterm_wt		!not sure if best at t_n or t_{n-1}
			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
!			the high-hand side of the above equation (row i)
			rhs(i) = u(i,nyfix) + del_t * q(x(i), y(nyfix),tn)
		enddo			!do i


!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)		= xbl(y(nyfix), tn)	!this is not a consistent boundary condition, not important as we are expecting a zero boundary condition

!!
!!		deriv_term = 0.0 
!!		do i= 0,nyfix	
!!			deriv_term  = deriv_term  + Gbeta_wt(nyfix-i) * xbl(y(i+1), tn)
!!		enddo
!!		rhs(0)		= xbl(y(nyfix), tn) - yterm_wt * e(x_low, y(nyfix),tn) * deriv_term
!!
!!same as below
		rhs(0)		= 0							!intermediate sln = (1- e(x,y,t) del^\beta U/del y^\beta)U^{n+1}

		A(NX,NX-1)	= 0.0
		A(NX,NX)	= 1.0

		deriv_term = 0.0 
		do i= 0,nyfix	
			deriv_term  = deriv_term  + Gbeta_wt(nyfix-i) * xbh(y(i+1), tn)
		enddo
		rhs(NX)		= xbh(y(nyfix), tn) - yterm_wt * e(x_high, y(nyfix),tn) * deriv_term


!	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 "intermediate" solution
	u_int(0,nyfix) = rhs(0)			!u_int(0,nyfix) = low Dirichlet boundary condition
	do i = 1, NX-1
		temp = 0.0
		do j = 0,i-1	
			temp = temp + A(i,j) * u_int(j,nyfix)
		enddo	!done j
		u_int(i,nyfix) = (rhs(i) - temp) /A(i,i)
	enddo	!done i

	u_int(NX,nyfix) = rhs(NX) 		!u_int(NX,nyfix) = Dirichlet boundary condition
!
!-=-=-=-=-
!
!	Solution u(x, y ,t) is now available at this fixed y-point (nyfix). 
!	Now increment nyfix and repeat to get next y-slice
	enddo		!nyfix

!	update the solution, and then alternate the direction:
	do i = 0,NX
		do j = 1, NY-1
			u(i,j) = u_int(i,j)		!Note u(0/NX,0/NY) (u at four corners) are never used/solved for/updated in this method
		enddo
		u(i,0)	= ybl(x(i),tn)		!At the y boundaries, these values are not solved/are also not needed   weep
		u(i,NY)	= ybh(x(i),tn)
	enddo

!YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY

!	nstep = nstep + 1	!advance the timestep for the alternate direction
!	tn = T_begin + nstep * del_t
	do nxfix = 1, NX-1
!		Set up matrix equations for the finite differences: Ay U_i^n+1 =  U_i^n = rhs_i ..
!		First the NY x NY coeffs for the super-triangular matrix in Y-direction

		do i = 1 , NY-1			!compute matix coeffs for row i (y-direction implicit)
!			x_diffuse_wt = d(x(nxfix), y(i),tn) * xterm_wt		!not sure if best at t_n or t_{n-1}
			y_diffuse_wt = e(x(nxfix), y(i),tn) * yterm_wt	
			do k = 0, i+1		!column k
				A(i,k)	= - y_diffuse_wt * Gbeta_wt(i-k+1)
			enddo
			A(i,i) = A(i,i) + 1.0
!			and the right hand side of this equation
			rhs(i) = u(nxfix,i)
		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)

!	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 "intermediate" solution
		u_int(nxfix, 0) = rhs(0)			!u_int(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_int(nxfix,j)
			enddo	!done j
			u_int(nxfix,i) = (rhs(i) - temp) /A(i,i)
		enddo	!done i

		u_int(nxfix,NY) = rhs(NY) 		!Note u_int(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 the solution and then repeat the big "do while in the sky"
	do j = 0,NY
		do i = 1, NX-1
			u(i,j) = u_int(i,j)
		enddo
		u(0,j)	= xbl(y(j),tn)		!At the x boundaries, these values are not updated in this sweep
		u(NX,j)	= xbh(y(j),tn)
	enddo
enddo		!do while nstep


!Print the final results 
!*******Begin USER NOTE 2
! NOTE***_ Since for this test case 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 exact solution and the error should be deleted,
!		That is delete code for: u_exact, sln_err, abs_sln_err)
100	Format	 (5(F12.4,3x))
write(UNIT=20,FMT=*)  "x(i), y(j), u_exact, u(i,j), sln_err,T_end=", tn
abs_sln_err = 0.0 
do i = 1 ,NX-1
	do j = 1, NY-1		!do j = 0, NY
		u_exact = exp(-tn) * (x(i)**3) * (y(j)**3.6)
		sln_err = abs(u_exact - u(i,j))
		abs_sln_err = max(sln_err, abs_sln_err)
		Write(UNIT=20,FMT = 100) x(i) , y(j) , u_exact, u(i,j), sln_err
		write(*,*) "/"	
	enddo
enddo

Write(UNIT=20,FMT = *) "Absolute error = ", abs_sln_err
!*******End USER NOTE 2

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

END Program Two_D_ADI_FPDE
		
!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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		 
		 
!	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	
		
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		 
		 
!	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  =

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		 

!	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)


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		 


!	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)

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		 

!	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


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		 

!	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)

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		 

!	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


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		 


!	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)


End Function ybh