Program tso_frade_imp_Euler
!-------------------------------------------
!-------------------------------------------
!Initial release date for this code: Nov. 10, 2003
!By: Charles Tadjeran
!Acknowlegements:	This work is supported by NSF Grants DMS-0139927, DMS-0139943, AR-9980489 and
!					DOE grant DE-FG03-98ER14885
!-------------------------------------------
!Code size: 782 Executable source lines of code (SLOC)
!------------------------------------------
!**Program Description: This Fortran program solves the Fractional Radial Advection Dispersion Equation(FRADE) problem.
!	The mobile-immobile FRADE initial-boundary value differential equation is of the form:
!	del(c)/del(t) + theta * del^beta (c)/del(t^beta) = - v(r) * del(c)/del(r) + d(r) * del^alpha(c)/del(r^alpha) + q(r,t)
!	over the radial range:  r_left < r < r_right , and time limits: T_begin < t <= T_end
!	where c = c(r,t) is the	solute concentration (M/L^3)
!	
!	The initial condition is specified by: c(r, t=T_begin) = f(r)
!	The boundary condition on the left boundary is a Dirichelt condition: c(r=r_left, t) = lb(t) 
!	The boundary condition on the right boundary is a Neumann condition: del(c(r=right,t))/del(r) = rb(t)
!
!	Additionally:
!	v = v(r) >=0 is advective (pore water) velocity function	(L/T)
!	d = d(r) >=0 is the	dispersive coefficient function			(L^2/T)
!	q = q(r,t)   is the non-homogenous source/sink function		(M/[L^3·T])
!
!-------------------------------------------
!**Algorithm description:
!   The first part of this program solves the mobile FRADE problem (i.e., for which
!	value of theta = 0 in the above equation). 
!	Note that a value of theta = 0 corresponds to a FRADE of the form:
!	del(c)/del(t) = - v(r) * del(c)/del(r) + d(r) * del^alpha(c)/del(r^alpha) + q(r,t)
!	which is the differential equation sovled first by this program. This problem is solved
!	by an implicit Euler method, applying a shifted Grunwald finite difference formula to estimate the fractional 
!	radial derivative (the dispersive term). An upstream finite estimate is used to estimate
!	the advective term.
!	The second part of this program, solves the mobile-immobile FRADE by applying a temporal 
!	subordination (TSO) to the mobile case solution obtained aboove,
!	to estimate the solute concentration at the extraction/observation well.
!	Also note that this method of TSO to estimate the concentration of the extracted solute 
!	at the right boundary relies on the observation/requirement that del(c(r=right,t))/del(r) =0.
!	See "Fractal Mobile/Immobile Solute Transport", by Schumer et al.[2003] for additional details 
!	on subordination.
!-------------------------------------------
!
!	The 'mobile' FRADE being solved is : 
!	del(c)/del(t) = - v * del(c)/del(r) + d * del^alpha (c)/del(r^alpha) + q(r,t), 
!	for r_left < r < r_right ,  T_begin < t <= T_end
!	and where 
!	c = c(r,t) ,			solute concentration
!	v = v(r) = Vconst/|r|	advective velocity function
!	d = d(r) = Dconst/|r|	dispersion coefficient function
!	q = q(r,t)				"non-homogenous"forcing (e.g., injection well) function 
!	1 <= alpha <= 2			fractional order derivatvie diffusion term
!	c(r, 0) = f(r)			initial condition function
!	c(r_left,t) = lb(t)		Dirichlet left boundary function (lb(t) expected to be = 0 for the fractional PDEs)
!	del(c(r_right,t))/del(r)= rb(t) Neumann boundry function on the right edge of r domain
!
!NOTE:
!	This formulation of the problem assumes a flow from left to the right, 
!	and such that the injection well is to the left of the extraction well.
!	The Tracer extraction/output well, with a radius equal to |r_right|, is on 
!	the right end of r-interval and is assumed to be centered at r =0. 
!	The left boundary should be chosen/set at far enough from the injection
!	well, located at r = r_inj_well, so that no solute, if any, is expected to reach the left 
!	boundary during the time interval of the interest.
!	For EXAMPLE, for the Shoal test, with an extraction well of radius = 0.127 m,
!	we set r_right = -0.127m, r = r_inj_well = -30.127m (i.e., the distance between 
!	injection and extraction is 30m), and we arbitrarily set 
!	r_left =-60.127m (injection well far away from left boundary).
!
!-----------------------------------
!-----------------------------------
!Following are the USER INPUT VALUES (user will be prompted at run time):
!		r_left		:left limit on r
!		r_right		:right limit on r	( r_left < r < r_right <0)
!		N			:number of sunintervals in r-direction   (defines del_r = (r_right - r_left)/N  )
!		alpha		:the fractional order
!		T_begin		:start time (defaults to zero, change here if different start time is desired)
!		T_end		:final integration time t : on T_begin <= t <= T_end
!		M			:M number of time steps    (defines del_t = (T_end - T_begin)/M  )
!		r_inj_well	:location of injection well
!		inject_time	:injection time of solute 
!		mass_injected:total mass of the injected solute into acquifer
!		Qwell		:average extraction rate  (L^3/T)
!		thick		:aquifer thickness 
!		por			:aquifer porosity 
!		disper		:aquifer dispersivity [i.e., defines d(r) = disper * v(r) ]
!
!		savefilename:a file name to save the solute concentration at the extraction well
!						this file will contain the mobile phase result, and when needed,
!						the output from mobile/immobile results using TSO method.
!
!Additional USER INPUT VALUES, if immobile phase is also present (i.e, TSO will be done):
!		beta		: fractional order of time derivative (for immobile phase)
!			**---->	: a value of beta <= 0, will skip the temporal subordination (no immobile phase)
!	IF beta >0, then
!		n_of_betas	:number of different beta values to run the sim for
!		beta_incr	:increments of beta to compute immobile phase
!		theta_start	:a starting value theta to run the sim for 
!		n_of_thetas	:number of different theta values to run the sim for
!		theta_incr	:increments of theta to compute immobile phase
!-----------------------------------
!-----------------------------------
!Following are the USER DEFINED FUNCTIONS (User to provide these functions, then compile and run)
!		Function v(r)	:transport velocity, (The example here is defined by v(r) = Vconst/r)
!		Function d(r)	:dispersion coefficient (The example here is defined by d(r) = Dconst/r)
!		Function f(r)	:initial condition - expected to VANISH when r is "close" to left boundary
!			**---->		:also for the TSO implimented here, we require c(r=r_right,t=0) = f(r_right) = 0
!		Function q(r,t)	:forcing function (Example in the function here is the Shoal injection well)
!		Function lb(t)	:left (Dirichelt) boundary function c(r=rleft,t) = lb(t) (theoretically =0)
!		Function rb(t)	:right (Neumann) boundary function - del(c(r=r_right,t))/del(r) = rb(t)  (theoretically =0)
!
!-----------------------------------
! Comment: For example, the Shoal problem was solved with the following user inputs:
!		r_left = -60.127		(meters)
!		r_right = -0.127		(meters, extraction well wall)
!		N = 60
!		alpha = 1.8
!		T_begin	=0
!		T_end = 600				(days)
!		M = 4800				(large M used to get accurate results from TSO algorithm)
!		r_inj_well = -30.128	(meters)
!		inject_time	= 4.0 days
!		mass_injected = 20.81	(kg)
!		Qwell = 12.4			(meters^3/day)
!		thickness = 35			(meters) 
!		porosity = 0.011	
!		disperersivity = 0.34	(meters)
!		beta = 0.4
!		n_of_betas = 1
!		beta_incr = 0
!		theta_start	0.179 
!		n_of_thetas = 1	
!		theta_incr = 0
!
!-----------------------------------
!
	real	:: alpha		!fractional derivative order for the diffusive term
	real	:: Vconst		!advection velocity term in v(r) = Vconst/ abs(r)
	real	:: Dconst		!dispersion coeff term for for d(r) = Dconst/ abs(r)
COMMON/eqn_params/alpha, Vconst, Dconst

	real	:: T_begin =0.0	!beginning integration time: defaults to 0.0
	real	:: T_end		!final integration time t : on T_begin <= t <= T_end
	integer	:: M			!number of time steps to take
	real	:: del_t		!time step size. Note the method is unconditionally stable.
COMMON/t_params/T_begin, Tend,    M, del_t

	real	:: r_left		!left limit on r
	real	:: r_right		!right limit on r
	integer	:: N			!number of subintervals in r direction
	real    :: del_r		!radial subinterval. del_r is set by user input N
COMMON/r_params/r_left,  r_right, N, del_r

	real	:: r_inj_well	!location of injection well
	real	:: inflow_rate	!computed/normalized injected solute flow rate
	real	:: inject_time	!injection time of solute 
COMMON/inject_well_params/ r_inj_well, inflow_rate, inject_time 		 

	real	:: mass_injected	!total mass for injecting solute into acquifer

!Following are user input values for aquifer parameters 
	real	:: qwell, thick, por, disper	!Qwell, thickness, porosity, dispersion

!The following declarations are made to avoid runtime errors in Fortran!
	real	:: f			!initial concentration c(r,t=T_begin) = f(r).
	real	:: q			!forcing(source/sink) function in the ADE where q = q(r,t)

	real	:: lb			!left boundary function c(r=r_left,t) = lb(t)   :Dirichelt 
	real	:: rb			!right boundary function - del(c(r=r_right,t))/del(r) = rb(t)  (Neumann) 

!Start of program variables

!NOTE: Fortran compiler performs no array out of bounds checking.
!Following array dimension, and N_max, have to be bumped for a larger number of r-subintervals is desired
	integer, parameter:: N_max = 1000	!max number of r-subintervals. Note that N is REQUIRED to be <= N_max
	real    :: r(0:N_max), c(0:N_max), vv(0:N_max), dd(0:N_max)
	real	:: A(0:N_max,0:N_max), P(0:N_max)	!coefficient matrix and the right hand side vector of equations		!
	real    :: G_wt(0:N_max)				!array of Grunwald weights - computed onlt once during the integration

	integer, parameter:: M_max = 10000	!M_max should be increased if larger number of timesteps is desired
	real	:: c_extract(0:M_max)		!solute time history extracted at the right bounday

	real	:: csum, tn_hours

	integer	:: ntime			!basically counts timesteps
	integer	:: checkstatus		!file open/close status report
	character(40)		:: savefilename !output file containing sim info and output at extraction well
	logical	:: done				!a logical flag to finish tasks

!-----Following declarations for the Temporal SubOrdination(TSO) section
	real	::	beta , theta			!for ADE  del(c)/del(t) + theta * del^beta(t)/del(t^beta) = L c(x,t)
	real	::	beta_start, beta_incr	!to do a number of sims for different values of beta
	real	::	theta_start, theta_incr	!to do a number of sims for different values of theta
	integer	::	n_of_betas, n_of_thetas

	integer	::	i_beta, i_theta, ntso, n_u
	real	::	beta_term, tso, h_tso, sum_term, u, theta_u_term ! some temporary working variables
	real	::	g_tso_arg, error_estimate_tso, tso_hrs
!-----End temporal subordiantion declarations


	real , parameter	:: pi = 3.141592
	real , parameter    :: almost_zero = 1.0E-6 ! ~machine precision (for 32-bit word)		
 
	done = .FALSE.
	do while (.NOT. done) 
		write (*,*) "Input value for T_begin ?"
		read (*,*) T_begin

		write (*,*) "Input value for T_end ?"
		read (*,*) T_end
 
		write (*,*) "Input the number of time steps M ?"
		read (*,*) M

		if( (T_begin >= T_end) .OR. (M <= 0 ) )then
			write(*,*) " Improper input for time limits/steps! Try again!"
		else
			del_t = (T_end - T_begin)/real(M)		!sets timestep size
			if (del_t <= almost_zero) then
				write(*,*) "Timestep size= ", del_t, "is too small."
				write(*,*) "Try again! "
			else
				write (*,*) "Timestep = ", del_t
				done= .TRUE.
			end if
		end if
	
	end do ! while not done

	if ( M > M_max) then
		write(*,*) "Doing compiler's job! Maximum array size is exceeded."
		write(*,*) "In the code, set M_max =  ", M, "and then re-compile/re-run the program."
		write(*,*) "Press <Enter> to terminate program execution."
		read(*,*)
		stop
	endif

	done = .FALSE.
	do while (.NOT. done) 
		write (*,*) "Input value for left boundary r_left? "
		read (*,*) r_left

		write (*,*) "Input value for location of injection well r_inj_well ?"
		read (*,*) r_inj_well

		write (*,*) "Input value for right boundary r_right ?"
		read (*,*) r_right

		if( ((r_left*r_right) <= 0.0) .OR. (r_inj_well <= r_left) .OR. (r_inj_well >= r_right) )then
			write(*,*) "Improper r-boundary limits or injection well location! Try again! "
		else	
			done= .TRUE.
		endif
	end do ! while not done

	done = .FALSE.
	do while (.NOT. done) 
		write (*,*) "Input the number of subintervals in r-direction ?"
		read (*,*) N
		if (N <= 0) then
			write(*,*) "Improper number of subintervals! Try again! "	
		else	
			done= .TRUE.
		endif
	end do ! while not done

	if ( N > N_max) then
		write(*,*) "Doing compiler's job! Maximum array size is exceeded."
		write(*,*) "In the code, set N_max =  ", N, "and then re-compile/re-run the program."
		write (*,*) "Press <Enter> to terminate program execution."
		read (*,*)
		stop
	endif

	done = .FALSE.
	do while (.NOT. done) 
		write (*,*) "Input value for alpha ?"
		read (*,*) alpha
		if ((alpha < 1.0) .OR. (alpha > 2.0)) then		!Note for other values of alpha, modify/delete this statement	
			write(*,*) "Expected 1<= alpha <= 2! Try again! "	
		else	
			done= .TRUE.
		endif
	end do ! while not done


	done = .FALSE.
	do while (.NOT. done) 
		write (*,*) "Input value for average fluid extraction rate Qwell ?"
		read(*,*) qwell
		if(qwell <= 0.0) then
			write(*,*) " Improper input for fluid extraction! Try again!"
		else	
			done= .TRUE.
		endif
	end do ! while not done

	done = .FALSE.
		do while (.NOT. done) 
		write (*,*) "Input value for aquifer thickness ?"
		read(*,*) thick
		if(thick <= 0.0) then
			write(*,*) " Improper thickness input! Try again!"
		else	
			done= .TRUE.
		endif
	end do ! while not done

	done = .FALSE.
	do while (.NOT. done) 
		write (*,*) "Input value for aquifer porosity ?"
		read(*,*) por
		if ( (por <= 0.0) .OR. (por>=1.0) )then
			write(*,*) " Improper input for porposity! Try again!"
		else	
			done= .TRUE.
		endif
	end do ! while not done

	done = .FALSE.
	do while (.NOT. done) 
		write (*,*) "Input value hydraulic dispersivity"
		read(*,*) disper
		if(disper <= 0.0) then
			write(*,*) " Improper input for dispersivity! Try again!"
		else	
			done= .TRUE.
		endif
	end do ! while not done

	done = .FALSE.
	do while (.NOT. done) 
		write (*,*) "Input total amount of solute injected ?"
		read(*,*) mass_injected
		if(mass_injected <= almost_zero) then
			write (*,*) " Improper value for solute injected (none)! Try again!"
		else	
			done= .TRUE.
		endif
	end do ! while not done

	done = .FALSE.
	do while (.NOT. done) 
		write (*,*) "Input number of days solute is injected ?"
		read(*,*) inject_time
		if(inject_time <= del_t) then
			write(*,*) " Injection time smaller than timestep --> Numerically no solute is injected!"
			write(*,*) " Improper injection duration or timestep size. Please reconsider the input parameters!"
			write(*,*) " Press <Enter> to terminate program execution."
			read(*,*)
			stop 
		else	
			done= .TRUE.
		endif
	end do ! while not done

	write (*,*) "Input a FILENAME to save the extracted solute data at the extraction well ?"
	read(*,*) savefilename		

!Initialize the discretized problem settings
	del_r = (r_right - r_left)/real(N)		!sets meshsize in r-direction
	Vconst	= qwell/(2.0*pi*thick*por)		!sets advective velocity term  for v = Vconst/r
	Dconst	= Vconst*disper					!sets dispersion constant term for d = Dconst/r
	inflow_rate = mass_injected/(inject_time* pi*thick*por*((r_inj_well - 0.5*del_r)**2.-(r_inj_well + 0.5*del_r)**2.))
										!above sets solute injection mass rate

	open(UNIT = 20, IOSTAT = checkstatus, FILE =savefilename, STATUS = "NEW")
	write(*,*) "File Open IOSTAT =", checkstatus
	Write(UNIT=20,FMT=*) "r_left             =" , r_left
	Write(UNIT=20,FMT=*) "r_inj_well         =" , r_inj_well
	Write(UNIT=20,FMT=*) "r_right            =" , r_right
	Write(UNIT=20,FMT=*) "N(r-sunintervals)  =" , N
	Write(UNIT=20,FMT=*) "del_r              =" , del_r
	Write(UNIT=20,FMT=*) "T_begin            =" , T_begin
	Write(UNIT=20,FMT=*) "T_end              =" , T_end
	Write(UNIT=20,FMT=*) "M(timesteps)       =" , M
	Write(UNIT=20,FMT=*) "del_t              =" , del_t
	Write(UNIT=20,FMT=*) "alpha              =" , alpha
	Write(UNIT=20,FMT=*) "Qwell              =" , qwell
	Write(UNIT=20,FMT=*) "Thickness          =" , thick
	Write(UNIT=20,FMT=*) "Porosity           =" , por
	Write(UNIT=20,FMT=*) "Dispersivity       =" , disper
	Write(UNIT=20,FMT=*) "Inject_Time        =" , inject_time
	Write(UNIT=20,FMT=*) "Mass_Injected      =" , mass_injected
	Write(UNIT=20,FMT=*) "Vconst             =" , Vconst
	Write(UNIT=20,FMT=*) "Dconst             =" , Dconst
	Write(UNIT=20,FMT=*) "Solute_Inflow_Rate =" , inflow_rate

!Below it is assumed that V and D are indepedent of time. If not, then this block should be repositioned
	do i = 0 , N
		r(i)    = r_left + real(i) * del_r 
		vv(i)   = v(r(i))				!set convection constant at grid points
		dd(i)	= d(r(i))				!set diffusion constant at grid points
		c(i)	= f(r(i))				!sets initial condition at time t = T_begin
	end do

!Compute some weights for computational aid
	del_r_inv	= real(N)/(r_right - r_left)	!value of 1/del_r
	G_wt_scale	= del_r_inv ** alpha			!the 1/del_r**alpha,same as G_wt(0),factor used in Grunwald weights

!The Grunwald weights, defined below, are computed only once for the entire r-domain
!	G_wt(0) = 1.0 / (del_r**alpha)
!	do i = 0 , N-1
!		G_wt(i+1) = - G_wt(i) * (alpha - i)/real(i+1)
!	end do
!We'll use normalized (with 1/del_r**alpha factored out) Grunwald weights instead:
	G_wt(0) = 1.0
	do i = 0 , N-1
		G_wt(i+1) = - G_wt(i) * (alpha - i)/real(i+1)
	end do
	

!	Set up matrix equations for the finite differences: in the interior A C^{m+1} = C^m + del_t * Q^{m+1} ==> P^{m+1} ..
!	First the the N x N coeffs matrix (for interior C^{m}= transpose[c(1,m),c(2,m)....c(N,m)]) for shifted(=1) Grunwald formulation

    do i = 1 , N-1			!compute matix coeffs for row i -- ** Compute only once
!		convect_wght = vv(i) * del_t /(abs(r(i)) * del_r)
		convect_wght = vv(i) * del_t / del_r
		diffuse_wght = dd(i) * del_t / del_r**alpha

		do j = 0, i+1		!column j
			A(i,j)	= - G_wt(i-j+1)	* diffuse_wght
		enddo

		A(i,i-1)	=     - convect_wght + A(i,i-1)
		A(i,i)		= 1.0 + convect_wght + A(i,i)
	enddo

!And the Neumann boundary conditons with the one-sided finitie difference equation gives
	A(N,N-1)= -1.0
	A(N,N)	= 1.0

!	Solve the resulting matrix system - note A is"sum of a lower tri-angular and a super-diagnoal matrices
!
!   Since the coeff matrix is independent of time, invert/solve only once(save coeffs)
!	Start backward sweep to make it a lower triangular system
!
!	First subtract a multiple of the last row with two non-zero entries	
	ratio		= A(N-1,N) / A(N,N)				!With A(N,N)=1, ratio is same as A(N-1,N)
	A(N-1,N-1)	= A(N-1,N-1)- ratio * A(N,N-1)	!Note A(N,N)=1 for 1-sided finite diff for Neumann boundary
	A(N-1,N)	= ratio							!Just to make the formula universal across all N's
!	Remaining rows are not so sparse as the last
	do i = N-1, 2, -1
		if (abs(A(i,i)) <= almost_zero) then		!diagnal entry almost zero
			write(*,*) "Singular coeff matrix A, at i = ",i
			write (*,*) " 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	!done j

		A(i-1,i) = ratio		!store in an "un-used" location for backward sweeps of rhs
	enddo		!done i
!
!	now matrix A is a lower triangular system
!	**Computational note/enhance later: Array elements A(i,j) for j=i+2, all i's, are not used at all!
!
!Initialize step one of integration
	tn			= T_begin	!starting integration time t
	c_extract(0)	= c(N)		!solute extracted at right boundray

	Write(UNIT=20,FMT=*) "**Results** "
	Write(UNIT=20,FMT=*) "Days,  Hours, C(r_right,tn)"

!Then apply Implicit Euler method at each timestep.
!Use a right-shifted Grunwald formula to get stability.

	do ntime = 1, M
		tnp1		= T_begin + del_t *  ntime       !next integration time       
		do i = 1 , N-1			!P is the right hand-side of the matrix equation A C = P
			P(i) = c(i) + del_t * q(r(i),tnp1)  !-A(i,0) * c(r=r_left,t=tnp1) but with c(.)=0,.. 
		enddo

		P(0) = lb(tnp1)				!from Dirichlet condition on left boundary ('expected' to be =0)
		P(N) = del_r * rb(tnp1)		!from Neumann boundary condition ('expected' to be =0)

!Solving the resulting matrix system
!Start backward sweep to "adjust" the rhs

		do i = N-1, 1, -1
			P(i) = P(i)	- A(i,i+1) * P(i+1)
		enddo

!	Start forward sweep to get the solution
		c(0) = P(0)	!Note c(0) is the left Dirichlet boundary condition
		do i = 1, N-1
			temp = 0.0
			do j = 0,i-1	
				temp = temp + A(i,j) * c(j)
			enddo	!done j
			c(i) = (P(i) - temp) /A(i,i)
		enddo	!done i

		c(N) = P(N) +  c(N-1)		!and the Neumann boundary cond on the right

!Solution c(r,t) is now available at time tnp1

		tn = tnp1
		c_extract(ntime) = c(N)					!solute extracted at right boundary-store for TSO
		tn_hours = tn * 24.0					!time in hours for plotting purposes
		Write(UNIT=20,FMT=*) tn, tn_hours, c(N)	!Output right boundary results to file

	enddo !ntime
	write (*, *) "*** Conservation Integration Completed ***"


!=============== Start Temporal Subordination (tso) ===============
	write (*, *) "*** Starting Temporal Subordination ***"
	write(UNIT=20,FMT=*) "**TSO Data**"

!
!Note: the equation being solved is del(c)/del(t) + theta * del^beta(c)/del(t^beta) = L c(t) 
!Note: The following section of the code computes the temporal subordination accroding to
! h_tso(r=r_right, t= tso) = 
!	 intg_{u=0}^{u=tso} {c(r_right,u)*g_beta[(tso-u)/(theta*u)**(1/beta)]/ (theta*u)**(1/beta)}du
!Note: computations in this section are not optimized for performance.	
!-----------
!Note: will automatically cycle through a fixed set of values for beta and theta for computaional speed		
!Following for ease of use/avoid to type at terminal		
	write(*,*) "Input starting value for fractional order of time derivative beta (beta=0 means no TSO)"
	read(*,*)	beta_start
	if(beta_start > 0) then
		write(*,*) "Input number of beta values (>= 1) ?"
		read(*,*) n_of_betas

		write(*,*) "Input value for beta increment (>= 0.0) ?"
		read(*,*) beta_incr

		write(*,*) "Input starting value for theta (>= 0.0)?"
		read(*,*) theta_start

		write(*,*) "Input number of theta values (>= 1) ?"
		read(*,*) n_of_thetas

		write(*,*) "Input value for theta increment (>= 0.0) ?"
		read(*,*) theta_incr

		do i_beta =	 1, n_of_betas		!**--> gives a pre-defined set of values for beta 
			do i_theta = 1, n_of_thetas

				Write(UNIT=20,FMT=*) "Vconst     = " , Vconst
				Write(UNIT=20,FMT=*) "Dconst     = " , Dconst
				Write(UNIT=20,FMT=*) "alpha      = " , alpha

				Write(UNIT=20,FMT=*) "r_left     =", r_left
				Write(UNIT=20,FMT=*) "r_inj_well =", r_inj_well
				Write(UNIT=20,FMT=*) "r_right    =", r_right
				Write(UNIT=20,FMT=*) "N          =", N
				Write(UNIT=20,FMT=*) "del_r      =", del_r

				Write(UNIT=20,FMT=*) "T_begin    =", T_begin
				Write(UNIT=20,FMT=*) "T_end      =", T_end
				Write(UNIT=20,FMT=*) "M          =", M
				Write(UNIT=20,FMT=*) "del_t      =", del_t
	
				Write(UNIT=20,FMT=*) "beta,  theta,  hours,  concentration"

				beta	= (i_beta-1) * beta_incr + beta_start	! for multiple beta values
				theta	= (i_theta-1)* theta_incr + theta_start	! for multiple theta valuestheta 
				write(*,*) "Now computing TSO for beta = ", beta ," theta = ", theta
!				write(*,*) " Press <Enter> to continue"
!				read (*,*)

				beta_term = -1.0/beta
	
				do ntso = 1, M								!compute tso'ed values for each t
					tso = ntso * del_t						!time at which to compute tso'ed solute concentration
!					Note: h_tso = h(r=r_right,t=tso)		!this is what we want to compute below
!					The values of gamma(1-beta) do not matter here, since we have assumed f(r_right) =0
!					h_tso = 0.5*f(r_right)*gamma(1-beta)..	!initially zero concentration at time u=0
					h_tso = 0.0								!just what the above is at time u=0

					do n_u = 1, ntso  
						u				= n_u * del_t
						theta_u_term	= (theta * u) ** beta_term
						g_tso_arg		= (tso - u) * theta_u_term
						sum_term		= g_tso(beta, g_tso_arg) * theta_u_term * c_extract(n_u) 
						h_tso			= h_tso + sum_term
					end do ! n_u
					h_tso = (h_tso - 0.5 * sum_term) * del_t  !and adjust for the last point in trapezoid integration

!				Note : the error estimate below just for informational purpose	
!				This is "an" error estimate based on how fast the integrand dies off - ie, effect of an additional term for u-->infinity
					error_estimate_tso = del_t * sum_term

					tso_hrs = 24.0 * tso
!					Write(UNIT=20,FMT=*) beta, theta, tso_hrs , h_tso,  error_estimate_tso
					Write(UNIT=20,FMT=*) beta, theta, tso_hrs , h_tso
				end do	!ntso for each time t
			end do ! i_theta 	
		end do ! i_beta 	
		
!===========   End Temporal Subordination
	else
		write(*,*) "No temporal subordination was perforemd (beta=0). "
	endif		!beta_start =0 means no TSO

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

	end Program  tso_frade_imp_Euler


!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function v(r)
	real	:: rb
	real	:: r		!radial distance variable

	real	:: alpha	!fractional derivative order for the diffusive term
	real	:: Vconst	!advection velocity term V(r) = Vconst/ abs(r)
	real	:: Dconst	!dispersion coeff term for for D(r) = Dconst/ abs(r)
COMMON/eqn_params/alpha, Vconst, Dconst

v = Vconst / abs(r)
End Function v			

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function d(r)
	real	:: d
	real	:: r		!radial distance variable

	real	:: alpha	!fractional derivative order for the diffusive term
	real	:: Vconst	!advection velocity term V(r) = Vconst/ abs(r)
	real	:: Dconst	!dispersion coeff term for for D(r) = Dconst/ abs(r)
COMMON/eqn_params/alpha, Vconst, Dconst

d = Dconst/ abs(r)
End Function d

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function q(r, t)		!forcing function in the diff eqn, v used/added to set up test cases
	real	:: q
	real	:: r		!radial distance variable
	real	:: t		!time varible

	real	:: r_left	!left limit on r
	real	:: r_right	!right limit on r
	integer	:: N		!number of subintervals in r direction
	real    :: del_r	!radial subinterval. del_r is set by user input N
COMMON/r_params/r_left,  r_right, N, del_r

	real	:: r_inj_well	!location of injection well
	real	:: inflow_rate	!computed/normalized injected solute flow rate
	real	:: inject_time	!injection time of solute 
COMMON/inject_well_params/ r_inj_well, inflow_rate, inject_time 		 

q = 0.0
if (t <= inject_time + 1e-6) then		!This is an example from Shoal injection well
	if (abs(r-r_inj_well) < (0.5*del_r)) then
		q = inflow_rate		! units of q are in concentration (kg/m^3 water)
	endif
endif

End Function q


!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function f(r)			!initial condition for c(r,t) at t=0
	real	:: f
	real	:: r		!radial distance variable

	real	:: r_left	!left limit on r
	real	:: r_right	!right limit on r
	integer	:: N		!number of subintervals in r direction
	real    :: del_r	!radial subinterval. del_r is set by user input N
COMMON/r_params/r_left,  r_right, N, del_r
 
f = 0.0              

End Function f

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function lb(t)			!left Dirichlet boundary condition c(r=r_left,t) = lb(t)
	real	:: lb
	real	:: t		!time variable

	real	:: T_begin		!beginning integration time
	real	:: T_end		!final integration time t : on T_begin <= t <= T_end
	integer	:: M			!M number of time steps i user input
	real	:: del_t		!time step size. Note the method is unconditionally stable.
COMMON/t_params/T_begin, Tend,    M, del_t

lb = 0.0		! zero boundary function
End Function lb

!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function rb(t)			!right Neumann boundary condition del (c(r=r_right,t))/del(r) = rb(t)

	real	:: rb
	real	:: t		!time variable


	real	:: T_begin		!beginning integration time
	real	:: T_end		!final integration time t : on T_begin <= t <= T_end
	integer	:: M			!M number of time steps i user input
	real	:: del_t		!time step size. Note the method is unconditionally stable.
COMMON/t_params/T_begin, Tend,    M, del_t

rb = 0.0		! zero boundary function
End Function rb


!-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
Function g_tso(beta,t)			!Inverse Lapalace transform by Kanter
!
!Description:
!This function compute the function g_tso(beta,t) for user specified values of beta and t
!	and error tolerance. All computations are single precision in this code.
!	

	IMPLICIT NONE  
 
	real	:: g_tso
	real	:: beta, t , err_tol = 0.00001	!use err_tol for a satisfactory criterion for inverse Lapalace transform

	real , parameter	:: pi = 3.141592
	real , parameter    :: too_small = 1.0E-6 ! ~machine precision (for 32-bit word)		

	real	:: g_intgl_wt, del_u, beta_term, t_term, t_term2, u ,au, est_intgl_fine, est_intgl_coarse
	integer :: i, n_subints , n_subints_max=2000, tolerance_not_ok

	if( t .ge. too_small) then	 
		beta_term = 1.0 / (1.0-beta)
		t_term = (1.0 / t)**beta_term
		t_term2 = (1.0 / t) ** (beta * beta_term)
		g_intgl_wt = beta * beta_term * t_term /pi		!factor outside the integral for g_beta(t)

		n_subints = 10			!start the inegral estimate with this many subintervals
		est_intgl_coarse = 0.0	!some starting wag on the integral value
		tolerance_not_ok = 1

		do while (tolerance_not_ok)
			del_u = pi/ n_subints

!Now apply the trapezoidal rule to estimate the integral
			au = beta **(beta * beta_term) * (1.0-beta)			!limit as u-->0 for a(u=0)
			est_intgl_fine = 0.5 * au * exp(- t_term2 * au)

			do i = 1, (n_subints-1)
				u = i * del_u
				au = (sin(beta*u))**(beta * beta_term) * sin((1.0-beta)*u)/(sin(u))** beta_term
				est_intgl_fine = est_intgl_fine + au * exp(- t_term2 * au)
       		end do

!Note: As u -->pi, a(u) -->+infty, but integrand -->0, so we ignore the last term in the sum

			est_intgl_fine = g_intgl_wt * est_intgl_fine * del_u

!test to see if estimate is good enough
			if(abs(est_intgl_fine - est_intgl_coarse) < err_tol) then
				tolerance_not_ok = 0
			else					!get ready to refine the grid and re-compute
				est_intgl_coarse = est_intgl_fine
				n_subints = 2 * n_subints
			endif

			if(n_subints > n_subints_max) then
				write(*,*) "Number of subintervals exceeds max allowed. n_subints =", n_subints
				write(*,*) "Estimate for g_beta(t) =", est_intgl_fine
				write(*,*) "Press <Enter> to terminate"
				read(*,*)
				stop
			endif
		end do !while(tolerance_not_ok)
!Report the results of successful integration 
!		write(*,*) " Successful completion. g_beta(t)= ", est_intgl_fine
!		write(*,*) " Number of subintervals used, n_subints=", n_subints

		g_tso = est_intgl_fine

	else		! t is too small
		g_tso = 0.0
	endif

End Function g_tso