program Useries

	implicit none

	! program to simulate the radioactive decay of U238 and U235 to Pb206
	! and Pb207.  Written by Kirsten Menking, Dept. of Geology and Geography,
	! Vassar College, Poughkeepsie, NY  12604, Aug. 7, 2002.  Based on the following readings:
	! Dalrymple, G.B., 1991, The Age of the Earth, Stanford, CA:  Stanford University Press, 
	! p. 79-90, 99-102, 115-119, and 122-124.
	! Faure, G., 1986, Principles of Isotope Geology, 2nd Edition, New York:  John Wiley and Sons,
	! p. 38-45 and p. 283-296.


	! The variables used in this model are as follows:
	! RESERVOIRS:
	! U238 = The number of atoms of the radioactive parent isotope U238.
	! U235 = The number of atoms of the radioactive parent isotope U235.
	! Pb206 = The number of atoms of the stable daughter isotope Pb206.
	! Pb207 = The number of atoms of the stable daughter isotope Pb207.
	
	! FLUXES:
	! decay238 = The transfer of atoms from U238 to Pb206 through radiometric decay.
	! decay235 = The transfer of atoms from U235 to Pb207 through radiometric decay.
	
	! CONVERTERS:
	! hlife238 = The half-life of U238.
	! hlife235 = The half-life of U235.
	! u238orig = The original number of atoms of U238, prior to decay.
	! u235orig = The original number of atoms of U235, prior to decay.
	! drate238 = The decay rate of U238.
	! drate235 = The decay rate of U235.
	! r206_238 = The ratio of Pb206 to U238.
	! r207_238 = The ratio of Pb207 to U235.

	
	! OTHER VARIABLES: 
	! i = Counting loop incrementer
	! imax = final value of counting loop
	! t = Time in millions of years.
	! dt = Time step in millions of years.
	! tmax = total length of simulation in years

	real U238, U235, Pb206, Pb207, hlife238, hlife235	!Variables that are real numbers.
	real decay238, decay235, u238orig, u235orig			!Variables that are real numbers.
	real drate238, drate235, r206_238, r207_235			!Variables that are real numbers.
	real dt, t, tmax									!Variables that are real numbers.
	integer i, imax										!Variables that are integers.
	dimension U238(20000), U235(20000), Pb206(20000), Pb207(20000)	!Variables that are arrays.

	open(unit=14,file='useries_results.txt',status='replace') !Creates the file name 
	!"useries_results.txt" for unit 14.

	!*******************Initial conditions***************************	
	U238(1)=1.e9	!Initial value = 1 billion U238 atoms
	U235(1)=1.e9	!Initial value = 1 billion U235 atoms
	u238orig=1.e9	!1 billion U238 atoms.
	u235orig=1.e9	!1 billion U235 atoms.
	Pb206(1)=0.		!Initial value = zero Pb206 atoms
	Pb207(1)=0.		!Initial value = zero Pb208 atoms
	hlife238=4468.	!Millions of years
	hlife235=703.8  !Millions of years
	drate238=log(2.)/hlife238	!Decay rate for U238 (in fortran 90, log is the natural log)
	drate235=log(2.)/hlife235	!Decay rate for U235
	dt=0.25			!Time step in millions of years.
	tmax=3500.		!Total length of simulation in years
	imax=int((tmax+dt)/dt)	!number of iterations required
	t=0.			!Starting time
	

	!Equations for Flows*****************************************************
	decay238=drate238*U238(1)*dt	!Initial equation for decay238 at time=1.
	decay235=drate235*U235(1)*dt	!Initial equation for decay235 at time=1.

 	do i=2,imax

	U238(i)=U238(i-1)-decay238		!Atoms
	U235(i)=U235(i-1)-decay235		!Atoms
	Pb206(i)=u238orig-U238(i)		!Atoms
	Pb207(i)=u235orig-U235(i)		!Atoms

	r206_238=Pb206(i)/U238(i)		!Ratio: unitless
	r207_235=Pb207(i)/U235(i)		!Ratio: unitless

	t=t+dt

	write(14,*)t,r207_235,r206_238	!Write output file

	decay238=drate238*U238(i)*dt	!Equation for decay238 for next iteration
	decay235=drate235*U235(i)*dt	!Equation for decay235 for next iteration 
	enddo

	end