implicit double precision(A-H,O-Z) real*8 t(83),f(83),s(83), & rm(86),sm(86),pdis(1590,86),psto(1590),alpha(5),beta(5) integer*4 itemp(86) character*80 infile,ofile integer ix, iy, iz, ndate(5) common /randc/ ix, iy, iz ix=1234 iy=2345 iz=3456 open(1,file='cal2.txt') open(2,file='test.txt') open(3,file='wjeff.txt') read(3,*) ndate(1) do 5 i=1,ndate(1) read(3,200) rm(i),sm(i) 200 format(T25,F4.0,T32,F3.0) 5 continue close(3) open(3,file='mndv1.txt') read(3,*) ndate(2) do 6 i=ndate(1)+1,ndate(1)+ndate(2) read(3,200) rm(i),sm(i) 6 continue close(3) infile='mndv2.txt' open(3,file=infile) read(3,*) ndate(3) do 7 i=ndate(1)+ndate(2)+1,ndate(1)+ndate(2)+ndate(3) read(3,200) rm(i),sm(i) 7 continue close(3) infile='mndv3.txt' open(3,file=infile) read(3,*) ndate(4) do 8 i=ndate(1)+ndate(2)+ndate(3)+1, & ndate(1)+ndate(2)+ndate(3)+ndate(4) read(3,200) rm(i),sm(i) 8 continue close(3) infile='mndv4.txt' open(3,file=infile) read(3,*) ndate(5) do 9 i=ndate(1)+ndate(2)+ndate(3)+ndate(4)+1, & ndate(1)+ndate(2)+ndate(3)+ndate(4)+ndate(5) read(3,200) rm(i),sm(i) 9 continue close(3) C Read calibration curve data do 10 i=1,83 read(1,*) t(i),f(i),s(i) t(i)=1950.-t(i) 10 continue do 11 i=1,5 beta(i)=1600 alpha(i)=251 11 continue do 20 i=1,86 call calibrate(t,f,s,rm(i),sm(i),psto) do 30 j=1,1590 pdis(j,i)=psto(j) 30 continue 20 continue C Burn-in for 5000 cycles do 90 k=1,5000 call mcmc(ndate,beta,alpha,pdis) 90 continue C Now do 50,000 iterations, writing every 50th one do 100 k=1,50000 do 110 l=1,50 call mcmc(ndate,beta,alpha,pdis) 110 continue write(2,99) idnint(beta(5)),idnint(alpha(5)), & idnint(beta(4)),idnint(alpha(4)), & idnint(beta(3)),idnint(alpha(3)), & idnint(beta(2)),idnint(alpha(2)), & idnint(beta(1)),idnint(alpha(1)) print *,k 99 format(10(i4,1x)) c99 format(4(1x,F7.1)) 100 continue stop end ************************************************************************ subroutine MCMC(ndate,beta,alpha,pdis) ************************************************************************ implicit double precision(A-H,O-Z) real*8 beta(5),alpha(5),pdis(1590,86),a,b,zlo,zhi,a1,b2,td(1348) real random integer itemp(86),ipoint,ndate(5) integer ix, iy, iz common /randc/ ix, iy, iz do 2 i=1,1348 2 td(i)=0.d0 1 a1=250.5d0 b2=1600.d0 C Draw dates within boundaries for first (earlier) phase do 910 j=1,ndate(1) zlo=pdis(ipoint(beta(5)),j) zhi=pdis(ipoint(alpha(5)),j) u=dble(random()) z=u*(zhi-zlo)+zlo call bsearch(z,pdis,j,idate,beta(5),alpha(5)) itemp(j)=1601-idate 910 continue C Draw dates within boundaries for second (later) phase do 920 j=ndate(1)+1,ndate(1)+ndate(2) zlo=pdis(ipoint(beta(4)),j) zhi=pdis(ipoint(alpha(4)),j) u=dble(random()) z=u*(zhi-zlo)+zlo call bsearch(z,pdis,j,idate,beta(4),alpha(4)) itemp(j)=1601-idate 920 continue C Draw dates within boundaries for third (later) phase do 930 j=ndate(1)+ndate(2)+1,ndate(1)+ndate(2)+ndate(3) zlo=pdis(ipoint(beta(3)),j) zhi=pdis(ipoint(alpha(3)),j) u=dble(random()) z=u*(zhi-zlo)+zlo call bsearch(z,pdis,j,idate,beta(3),alpha(3)) itemp(j)=1601-idate 930 continue C Draw dates within boundaries for fourth (later) phase do 940 j=ndate(1)+ndate(2)+ndate(3)+1,ndate(1)+ndate(2)+ndate(3) & +ndate(4) zlo=pdis(ipoint(beta(2)),j) zhi=pdis(ipoint(alpha(2)),j) u=dble(random()) z=u*(zhi-zlo)+zlo call bsearch(z,pdis,j,idate,beta(2),alpha(2)) itemp(j)=1601-idate 940 continue C Draw dates within boundaries for fifth (later) phase do 950 j=ndate(1)+ndate(2)+ndate(3)+ndate(4)+1, & ndate(1)+ndate(2)+ndate(3)+ndate(4)+ndate(5) zlo=pdis(ipoint(beta(1)),j) zhi=pdis(ipoint(alpha(1)),j) u=dble(random()) z=u*(zhi-zlo)+zlo call bsearch(z,pdis,j,idate,beta(1),alpha(1)) itemp(j)=1601-idate 950 continue C=============This Section for Getting New West Jefferson Boundaries m=ndate(1) power=1.d0/(1.d0-dble(m)) imin=2000 do 93 j=1,ndate(1) if(itemp(j).lt.imin) imin=itemp(j) 93 continue imax=0 do 94 j=1,ndate(1) if(itemp(j).gt.imax) imax=itemp(j) 94 continue b=dble(imin)-0.5d0 a=dble(imax)+0.5d0 if(beta(1).gt.a) a=beta(1) C Sample beta for West Jefferson c u=dble(random()) c diff1=b2-alpha(5) c diff2=a-alpha(5) cumsum=0.d0 do 951 i=imax,1348+alpha(1)+1 t=dble(i) d=(t-alpha(5))**ndate(1)*(t-alpha(1))**8*(1349.d0-t+alpha(1)) if(d.le.0.d0) then print *,i,alpha(1) pause endif d=1/d td(i-imax+1)=d cumsum=cumsum+d 951 continue do 952 i=imax,1348+alpha(1)+1 td(i-imax+1)=td(i-imax+1)/cumsum 952 continue print *,td pause beta(5)=alpha(5)+((1-u)*diff1**(1-m)+u*diff2**(1-m))**power C Sample alpha for West Jefferson u=dble(random()) diff1=beta(5)-b diff2=beta(5)-alpha(4) alpha(5)=beta(5)-((1-u)*diff1**(1-m)+u*diff2**(1-m))**power C=============This Section for Getting New Moundville I Boundaries m=ndate(2) power=1.d0/(1.d0-dble(m)) imin=2000 do 95 j=ndate(1)+1,ndate(1)+ndate(2) if(itemp(j).lt.imin) imin=itemp(j) 95 continue imax=0 do 96 j=ndate(1)+1,ndate(1)+ndate(2) if(itemp(j).gt.imax) imax=itemp(j) 96 continue b=dble(imin)-0.5d0 if(alpha(5).lt.b) b=alpha(5) a=dble(imax)+0.5d0 C Sample beta for Moundville I u=dble(random()) diff1=beta(5)-alpha(4) diff2=a-alpha(4) beta(4)=alpha(4)+((1-u)*diff1**(1-m)+u*diff2**(1-m))**power C Sample alpha for Moundville I u=dble(random()) diff1=beta(4)-b diff2=beta(4)-a1 alpha(4)=beta(4)-((1-u)*diff1**(1-m)+u*diff2**(1-m))**power C=============This Section for Getting New Moundville II Boundaries m=ndate(3) power=1.d0/(1.d0-dble(m)) imin=2000 do 97 j=ndate(1)+ndate(2)+1,ndate(1)+ndate(2)+ndate(3) if(itemp(j).lt.imin) imin=itemp(j) 97 continue imax=0 do 98 j=ndate(1)+ndate(2)+1,ndate(1)+ndate(2)+ndate(3) if(itemp(j).gt.imax) imax=itemp(j) 98 continue b=dble(imin)-0.5d0 if(alpha(4).lt.b) b=alpha(4) a=dble(imax)+0.5d0 C Sample beta for Moundville II u=dble(random()) diff1=beta(4)-alpha(3) diff2=a-alpha(3) beta(3)=alpha(3)+((1-u)*diff1**(1-m)+u*diff2**(1-m))**power C Sample alpha for Moundville II u=dble(random()) diff1=beta(3)-b diff2=beta(3)-a1 alpha(3)=beta(3)-((1-u)*diff1**(1-m)+u*diff2**(1-m))**power C=============This Section for Getting New Moundville III Boundaries m=ndate(4) power=1.d0/(1.d0-dble(m)) imin=2000 do 99 j=ndate(1)+ndate(2)+ndate(3)+1,ndate(1)+ndate(2)+ndate(3) & +ndate(4) if(itemp(j).lt.imin) imin=itemp(j) 99 continue imax=0 do 100 j=ndate(1)+ndate(2)+ndate(3)+1,ndate(1)+ndate(2)+ndate(3) & +ndate(4) if(itemp(j).gt.imax) imax=itemp(j) 100 continue b=dble(imin)-0.5d0 if(alpha(3).lt.b) b=alpha(3) a=dble(imax)+0.5d0 C Sample beta for Moundville III u=dble(random()) diff1=beta(3)-alpha(2) diff2=a-alpha(2) beta(2)=alpha(2)+((1-u)*diff1**(1-m)+u*diff2**(1-m))**power C Sample alpha for Moundville III u=dble(random()) diff1=beta(2)-b diff2=beta(2)-a1 alpha(2)=beta(2)-((1-u)*diff1**(1-m)+u*diff2**(1-m))**power C=============This Section for Getting New Moundville IV Boundaries m=ndate(5) power=1.d0/(1.d0-dble(m)) imin=2000 do 101 j=ndate(1)+ndate(2)+ndate(3)+ndate(4)+1,ndate(1)+ndate(2) & + ndate(3)+ndate(4)+ndate(5) if(itemp(j).lt.imin) imin=itemp(j) 101 continue imax=0 do 102 j=ndate(1)+ndate(2)+ndate(3)+ndate(4)+1,ndate(1)+ndate(2) & + ndate(3)+ndate(4)+ndate(5) if(itemp(j).gt.imax) imax=itemp(j) 102 continue b=dble(imin)-0.5d0 if(alpha(2).lt.b) b=alpha(2) a=dble(imax)+0.5d0 C Sample beta for Moundville IV u=dble(random()) diff1=beta(2)-alpha(1) diff2=a-alpha(1) beta(1)=alpha(1)+((1-u)*diff1**(1-m)+u*diff2**(1-m))**power C Sample alpha for Moundville IV u=dble(random()) diff1=beta(1)-b diff2=beta(1)-a1 alpha(1)=beta(1)-((1-u)*diff1**(1-m)+u*diff2**(1-m))**power return end ************************************************************************ function ipoint(zdate) ************************************************************************ real*8 zdate ipoint=1601-idnint(zdate) return end ************************************************************************ subroutine calibrate(t,f,s,rm,sm,pdis) ************************************************************************ implicit double precision(A-H,O-Z) real*8 t(83),f(83),s(83),df,dt,dj,djp1, & a,b,df2,dj2,djp12,a2,b2,sto,cal,cals,p,rm,sm,pmax, & pdis(1590),ldate,udate,z integer*4 ip ip=0 pmax=0.d0 do 20 j=2,83 dt=t(j-1)-t(j) b=(f(j)-f(j-1))/(t(j)-t(j-1)) a=f(j)-b*t(j) do 35 k=0,idnint(dt)-1 sto=t(j-1)-dble(k) cal=a+b*sto p=dexp(-0.5d0*(rm-cal)**2/(sm**2)) ip=ip+1 if(ip.eq.1) then pdis(ip)=p else pdis(ip)=p+pdis(ip-1) endif 35 continue 20 continue do 40 j=1,1590 40 pdis(j)=pdis(j)/pdis(1590) return end ************************************************************************ subroutine calibrate2(rm,sm,pdis) ************************************************************************ * * WARNING: This subroutine never called * ************************************************************************ implicit double precision(A-H,O-Z) real*8 pdis(1590) integer*4 ip cal=1600.d0 p=dexp(-0.5d0*(rm-cal)**2/(sm**2)) pdis(1)=p do 20 jcal=2,1590 cal=1600.d0-dble(jcal) p=dexp(-0.5d0*(rm-cal)**2/(sm**2)) pdis(jcal)=p+pdis(jcal-1) 20 continue do 40 j=1,1590 40 pdis(j)=pdis(j)/pdis(1590) return end ************************************************************************ subroutine bsearch(z,pdis,ndate,idate,beta,alpha) ************************************************************************ implicit double precision(A-H,O-Z) real*8 pdis(1590,86),z,beta,alpha integer a,b,c,ndate ic=0 a=ipoint(beta) c=ipoint(alpha) b=(a+c)/2 1 ic=ic+1 if(z.le.pdis(b,ndate)) then c=b b=(a+b)/2 else a=b b=(b+c)/2 endif if(a.eq.b) then idate=a return endif if(b.eq.c) then idate=b return endif goto 1 return end ******************************************************************* real function random() ******************************************************************* c c Algorithm AS 183 Appl. Statist. (1982) vol.31, no.2 c c Returns a pseudo-random number rectangularly distributed c between 0 and 1. The cycle length is 6.95E+12 (See page 123 c of Applied Statistics (1984) vol.33), not as claimed in the c original article. c c IX, IY and IZ should be set to integer values between 1 and c 30000 before the first entry. c c Integer arithmetic up to 30323 is required. c integer ix, iy, iz common /randc/ ix, iy, iz c ix = 171 * mod(ix, 177) - 2 * (ix / 177) iy = 172 * mod(iy, 176) - 35 * (iy / 176) iz = 170 * mod(iz, 178) - 63 * (iz / 178) c if (ix .lt. 0) ix = ix + 30269 if (iy .lt. 0) iy = iy + 30307 if (iz .lt. 0) iz = iz + 30323 c c If integer arithmetic up to 5212632 is available, the preceding c 6 statements may be replaced by: c c ix = mod(171 * ix, 30269) c iy = mod(172 * iy, 30307) c iz = mod(170 * iz, 30323) c random = mod(float(ix) / 30269. + float(iy) / 30307. + + float(iz) / 30323., 1.0) return end