program Pbloss

	implicit none

	! Program to simulate the radioactive decay of U238 and U235 to Pb206
	! and Pb207.  This program can also take into account losses of U and Pb and a gain of U.
	! 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 = the total number of iterations
	! t = Time in millions of years
	! dt = Time step in millions of years
	! tmax = the length of the simulation in years
	! tlossgain = the time when Pb is lost or U is gained or lost



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

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


	!*******************Initial conditions***************************	
	U238(1)=1.e9	!Initial U238 value = 1 billion atoms
	U235(1)=1.e9	!Initial U235 value = 1 billion atoms
	Pb206(1)=0.		!Initial Pb206 value = zero atoms
	Pb207(1)=0.		!Initial Pb207 value = zero atoms
	hlife238=4468.	!halflife for U238 in millions of years
	hlife235=703.8  !halflife for U235 in millions of years
	dt=0.25			!timestep in millions of years
    drate238=log(2.)/hlife238	!decay rate for U238
	drate235=log(2.)/hlife235	!decay rate for U235
	t=0.			!starting time in years
	tmax=3500.		!ending time in years
	imax=int((tmax+dt)/dt)	!calculating the total number of iterations required
	tlossgain=1500.
	!***************************************************************


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

 	do i=2,imax ! go from 2 to 3500 million years

	U238(i)=U238(i-1)-decay238		!Atoms
	U235(i)=U235(i-1)-decay235		!Atoms
	Pb206(i)=Pb206(i-1)+decay238	!Atoms
	Pb207(i)=Pb207(i-1)+decay235	!Atoms
	
	 if(t.eq.tlossgain)then
	 Pb206(i)=Pb206(i-1)-(0.1*Pb206(i-1))	!Loss of Pb206 at time = 1500 million years
	 Pb207(i)=Pb207(i-1)-(0.1*Pb207(i-1))	!loss of Pb207 at time = 1500 million years
	 endif

	! if(t.eq.tlossgain)then
	! U238(i)=U238(i-1)+(0.3*U238(i-1))  !Gain of U238 at time = 1500 million years
	! U235(i)=U235(i-1)+(0.3*U235(i-1))  !Gain of U235 at time = 1500 million years
	! endif

	! if(t.eq.tlossgain)then
	! U238(i)=U238(i-1)-(0.3*U238(i-1))  !Loss of U238 at time = 1500 million years
	! U235(i)=U235(i-1)-(0.3*U235(i-1))  !Loss of U235 at time = 1500 million years
	! endif

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

	t=t+dt

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

	decay238=drate238*U238(i)*real(dt)	!Equation for decay238 at time=i
	decay235=drate235*U235(i)*real(dt)	!Equation for decay238 at time=i

	enddo

	end