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