PROGRAM MissouriRiver
		! This program shows how to apply the SUBROUTINE EINSTEIN to calculate the 
		!	bedload function. The problem statement can be found in the Closure of 
		!	Discussion, J. hydraul. Engrg., 132(?),???-???
		! Please fill in your information between PROGRAM and CONTAINS, and do not 
		!	change the lines after CONTAINS.
		!
		! Written by JUNKE GUO, 07/05/2005
		!
		implicit none
		real :: h, M, N, K, z, dm, E, J1, J2, qT
		
		h = 2.38 ! flow depth, m

		! Parameters of velocity and concentration distributions, see the closure, 
		!	Eqs.(23) and (24)
		M = 0.2171 ! m/s
		N = 1.8321 ! m/s
		K = 0.1775 ! kg/m^3
		z = 0.3469
	
		dm = (0.105 + 0.074)/2*1e-3		! sediment size
		E = 2*dm/h	! relative bed layer thickness
		
		CALL EINSTEIN(z,E,J1,J2)

		!unit sediment discharge, kg/s.m
		qT = K*M*h*(4.61*(1-E)**z/E**(z-1) + N/M*J1 + J2)
	
		write(*,*) 'z=',z, 'E=',E
		write(*,*) 'J1=',J1, 'J2=', J2
		write(*,*) 'qT=', qT, 'kg/s.m'
				
	CONTAINS

		SUBROUTINE EINSTEIN(z, E, J1, J2)
		! This subroutine is accoding to the references:
		!	Guo, J. and Julien, P.Y. (2004)."Efficient algorithm for computing 
		!		Einstein integrals." J. Hydraul. Engrg., 130(12), 1198-1201
		!	Guo, J. and Julien, P.Y. (2006). Closure of the discussion, JHE, 132(?)
		!
		! Written by JUNKE GUO, 07/05/2005
		!
			implicit none
			real, intent(in)  :: z, E
			real, intent(out) :: J1, J2
			real :: J10, J20, F2, S, E1
			real, parameter :: pi = 3.14159265358979, gamma = 0.57721566490153286060651
			integer :: k, n
 
			! for non-integer z
			if (abs(z-anint(z)) > 0.005) then 
				J10 = z*pi/sin(z*pi)	! J1 from zero to 1
				J1 = J10 - F1(z,E)		! J1 from E to 1

				! J2 from 0 to 1
				n = ceiling(z) + 2;
				J20 = 1-gamma - log(n-z) - 0.5/(n-z) + 1.0/12.0/(n-z)**2 - 1.0/120.0/(n-z)**4
				do k = 1, n
					J20 = J20 + 1.0/(k-z)
				end do
				J20 = -J10*J20

				! J2 from E to 1
				F2 = F1(z,E)*(log(E) + 1.0/(z-1.0))
				S = -1.0/(z-1.0)/(z-2.0)*F1(z-1.0,E)
				do k = 2, n
					S = S + (-1)**k/(z-k)/(z-k-1.0)*F1(z-k,E)
				end do
				F2 = F2 + z*S
				J2 = J20  - F2
			else
				! for z=0
				if (anint(z) == 0) then
					J1 = 1 - E
					J2 = -1 + E - E*log(E)
				else
					! for z=1				
					J1 = -1 - log(E) + E
					J2 = 1 - E + E*log(E) - 0.5*log(E)**2
				end if

				! for z=2,3,4,...
				if (anint(z) > 1) then 
					E1 = E/(1-E)
					do k = 2, anint(z)
						J1 = 1.0/(k-1.0)*E/E1**k - k/(k-1.0)*J1
						J2 = E/E1**k*log(E)/(k-1.0) - k/(k-1.0)*J2 + J1/(k-1.0)
					end do
				end if
			end if
		END SUBROUTINE	
		
		REAL FUNCTION F1(z,E)
			implicit none
			real, intent(in) :: z, E
			real :: S1, E1
			integer :: k

			E1 = E/(1.0-E)
			S1 = -1.0/(1.0-z)*E1**(1.0-z)
			do k = 2, ceiling(z) + 4
				S1 = S1 + (-1)**k/(k-z)*E1**(k-z)
			end do
				
			F1 = E/E1**z - z*S1
		END FUNCTION F1

	END PROGRAM MissouriRiver