c-----------------------------------------------------------------------
c	deriv.f 
c	Routines for differencing 2D functions with respect to one direction.
c  
c     Routine to difference r(NT,NR) in space direction or second index
c     Convergence tested using sin(x) + cos(t) x,t in [0,Pi]. Passed


		subroutine DiffRR(r,drdR,dk,NT,NR)
			implicit none

			integer	NT, NR

			real	r(NT,NR), drdR(NT,NR)

			real  dk

c		Local Variables

			integer i, j

			real divdk

			divdk = 0.5e0 / dk

c		Difference along Rstar direction for each time
			do i = 1, NT
c			End point differencing
c				drdR(i,1) = divdk * ( -4.0e0 * r(i,1) + 7.0e0*r(i,2)
c     .						-4.0e0 * r(i,3) + r(i,4))

				drdR(i,1) = divdk *(4.0e0 *r(i,2) - r(i,3) -
     .						3.0e0 * r(i,1))

c			Centered differencing

				do j = 2, NR-1 
					drdR(i,j) = divdk * (r(i,j+1) - r(i,j-1))
				end do	

c			End point differencing
c				drdR(i,NR) = divdk * (-r(i,NR-3) + 4.0e0 * r(i,NR-2)
c     .					- 7.0e0*r(i,NR-1) + 4.0e0 * r(i,NR))

				drdR(i,NR) = divdk * (-4.0e0 * r(i,NR-1) + r(i,NR-2) +
     .							3.0e0 * r(i,NR))

			end do

			return
		end
c-----------------------------------------------------------------------
c	Differencig routine. 
c Created MFH 
c		Routine to difference r(NT,NR) in time direction or first index
c		Convergence tested using sin(x) + cos(t) x,t in [0,Pi]. Passed
c
		subroutine DiffTT(F,dFdt,dt,N1,N2)
			implicit none

			integer	N1, N2

			real	F(N1,N2), dFdt(N1,N2), dt

c		Local variables

			integer	i, j

			real 	divdt

			divdt = 0.5e0 / dt

c	Loop over 2nd index
			do i = 1, N2

c		End point differencing

c				dFdt(1,i) = divdt * ( -4.0e0 * F(1,i) + 7.0e0*F(2,i)
c     .						-4.0e0 * F(3,i) + F(4,i))

				dFdt(1,i) = divdt *(4.0e0 *F(2,i) - F(3,i) -
     .						3.0e0 * F(1,i))

c		Centered differencing
				do j = 2, N1-1
					dFdt(j,i) = divdt * ( F(j+1,i) - F(j-1,i))
				end do

c		End point differencing
					
c				dFdt(N1,i) =	divdt * (-F(N1-3,i) + 4.0e0 * F(N1-2,i)
c     .               - 7.0e0*F(N1-1,i) + 4.0e0 * F(N1,i))

				dFdt(N1,i) = divdt * (-4.0e0 * F(N1-1,i) + F(N1-2,i) +
     .							3.0e0 * F(N1,i))
				
			end do

			return
		end
c-----------------------------------------------------------------------