C PROGRAM TO AUTOMATICALLY CALCULATE A THERMAL HISTORY C BOTH TEMPERATURE ENDS ARE FIXED. INITIAL TEMP = 600 C ************ COMMENTS ************ C MC = # OF COEFF. THAT DEFINE TH C MFIT = # OF ADJUSTED COEFFICIENTS C MPS = # OF TIME SUBINTERVALS C NTEMPS = # OF TEMP INTERVALS C ************ END COMMENTS *********** implicit double precision (a-h,o-z) parameter(ns=200,nc=100,ntst=1001,mxi=350,nn=319,cht0=1.d-2, $nd=10,xlambd=0.0005543d00,a=0.d00,perc=0.01d00,mc=15, $nwcy=10,mps=40,ntemps=22) dimension c(nc),c0(mc),tt(2,0:ntst),d0(nd),r39(0:ns), $xmat(nd,0:ns,mxi),lc(mxi),vc(nd),tilab(ns),telab(ns), $ys(0:ns),wt(ns),lista(nc),covar(nc,nc),alpha(nc,nc), $chq(nwcy),xi(0:2*nc,nd,mxi),matx1(ns,ns),matx2(ns,ns),matx(ns) common /var/ tti(nwcy,2,ntst),tto(nwcy,2,ntst), $agei(nwcy,2,ns),ageo(nwcy,2,ns) common /cont/ ncyc,nmaxi(nwcy),nmaxo(nwcy) dimension dyda(nc) external ftemp character tab1*9,yes*1 tab1=char(09) print *, 'Insert number of runnings' read *, nrun chqmin = 1.d20 ncyc=1 print *, "Insert max. plateau age" read *,tt(1,1) print *, "Insert actual sample temperature" read *,temp0 idem = -tt(1,1)*nrun mcyc = 0 mfit=mc-2 C SET OF CONTOUR ROUTINE VARIABLES print *,"Do you want to create contour matrices? (y or n)" read *,yes if(yes.eq."Y".or.yes.eq."y")then print *, "Insert minimum age" read *,agein agend = tt(1,1) C NUMBER OF TIME INTERVALS nsteps = int(tt(1,1)) if(nsteps.gt.100)nsteps = 100 dage = (tt(1,1) - agein)/nsteps C OUTPUT FILES open(unit=25,file='matx1.dat') open(unit=26,file='matx2.dat') open(unit=27,file='matx.dat') endif open(unit=20,file='param-free.out') open(unit=32,file='chisq.dat') open(unit=50,file='wtage.dat') open(unit=21,file='agenew-sd.samp') open(unit=40,file='chistall-inp.dat') open(unit=42,file='chistall-out.dat') open(unit=44,file='agesall-inp.dat') open(unit=46,file='agesall-out.dat') C INPUT FILES open(unit=17,file='temstep.in',status='old') open(unit=22,file='age.in',status='old') open(unit=24,file='sig.in',status='old') open(unit=30,file='arr-me.in',status='old') C INITIALIZATION OF OUTPUT FILES FOR XVGR SETS kc0=48 kc1=48 call wgraph(40,kc1) call wgraph(42,kc1) call wgraph(44,kc1) call wgraph(46,kc1) read(17,*)ni do 12 j=1,ni r39(j)=0.d0 read(17,*)telab(j),tilab(j) 12 tilab(j)=tilab(j)/5.256d11 close(17) np = 3 tt(2,1)= 873.d00 tt(2,2)=573.d00 tt(1,2)=tt(1,1)/2.d0 tt(2,3)=temp0 + 273.d0 tt(1,3)= 0.d0 C READ DOMAIN DISTRIBUTION PARAMETERS read(30,*)nst do 14 js=1,nst read(30,*)e,d0(js),vc(js) d0(js)= 10.d0**d0(js) 14 d0(js)=d0(js)/4.d0 *(24.d0 *3600.d0 *365d06) close(30) C CALCULATION OF LC nsq=2 nsum=0 m=-1 do 32 mi=1,nn if(nsum.gt.20)then nsum=1 nsq=2*nsq endif m=m+nsq if(m.gt.200)nsum=nsum+1 lc(mi)=m 32 continue call lab(r39,d0,e,vc,telab,tilab,lc,xmat,ni,nst,nn) call agesamp(r39,ys,wt,ni) call chebft(a,tt(1,1),c,nc,ftemp,tt,np) tt(2,1)=tt(2,1)-273.d0 tt(1,0)=tt(1,1)+1.d0 tt(2,0)=tt(2,1) sa = c(1)/2.d0 sb = c(1)/2.d0 c0(1)=c(1) do 13 i=1,mfit lista(i)=i+2 13 continue do 15 i=2,mfit+2 c0(i)=c(i) sa = sa + c(i) sb = sb + c(i)*(-1)**(i+1) 15 continue 34 call ranchist(c,c0,mc,sa,sb,idem) C TEST OF RANCHIST CONSTRAINT SOLUTIONS sa1 = c(1)/2.d0 sb1 = c(1)/2.d0 do 16 i=2,mfit+2 sa1 = sa1 + c(i) sb1 = sb1 + c(i)*(-1)**(i+1) 16 continue print *,sa-sa1,sb-sb1 C END TEST dch=0.d00 alamda=-1 call mrqmin(r39,ys,wt,ni,c,mc,lista,mfit,covar,alpha, $ nc,chisq,alamda,tt,xmat,lc,vc,d0,E,nst,nn,sa,sb) if (chisq.eq.0) then nfail=nfail + 1 goto 34 endif k=1 itst=0 nck2=0 write(20,*)' chi-squared: ',ni*chisq alam=alamda 1 write(20,*)'iteration #= ',k,' alamda: ',alamda if(dch.gt.0.d0)then write(20,*)' chi-squared: ',ni*chisq endif k=k+1 ochisq=chisq call mrqmin(r39,ys,wt,ni,c,mc,lista,mfit,covar,alpha, $ nc,chisq,alamda,tt,xmat,lc,vc,d0,E,nst,nn,sa,sb) dch=abs(ochisq-chisq) cht=cht0*chisq if (dch.le.cht.and.alamda.lt.alam) then itst=itst+1 write(20,*)'itst= ',itst elseif(dch.gt.cht)then itst=0 endif C SECOND CRITERIUS TO STOP if(dch.gt.cht.and.nchini.gt.nwcy)then dch1 = abs(chisq-chqmin) if(dch1/dch.gt.50.)then nck2=nck2 + 1 write(20,*)'nck2= ',nck2 endif endif alam=alamda if (itst.lt.3.and.nck2.lt.3) goto 1 write(20,*)'difference= ',(ochisq-chisq)*ni chq(ncyc)=chisq aq = 0.5d0 * (ni-2.d0) q = gammq(aq,0.5*ni*chisq) write(20,*)'chisq= ',ni*chisq,' GOF= ',q write(32,130) nchini+1,ni*chisq,chisq,q,k alamda=0.d0 c40=1.d0-dexp(-xlambd*tt(1,0)) call zita(c,mc,lista,mfit,lc,xi,tt(1,1),d0,E,tt,nst $ ,nn,perc,alamda) do 42 i=1,ni call age(i,c,ymod,dyda,mc,lista,mfit,vc,xmat,r39,xi, $ c40,nn,nst,ni,perc,alamda) 42 continue alamda=0.d0 call mrqmin(r39,ys,wt,ni,c,mc,lista,mfit,covar,alpha, $ nc,chisq,alamda,tt,xmat,lc,vc,d0,E,nst,nn,sa,sb) if (nchini.eq.0)chqmin = chisq C if (nchini.lt.nrun) then nchini = nchini + 1 write(20,*) write(20,*) 'New cycle:',nchini write(20,*) if (ncyc.eq.nwcy.or.nchini.eq.nrun)then open(unit=10,file='chist-inp.dat') open(unit=12,file='chist-out.dat') open(unit=18,file='ages-inp.dat') open(unit=19,file='ages-out.dat') if(nchini.gt.nwcy)then call readfile(10) call readfile(12) call readfile(18) call readfile(19) else call wgraph(10,kc0) call wgraph(12,kc0) call wgraph(18,kc0) call wgraph(19,kc0) endif if(ncyc.eq.nchini)then do 300 j=1,ncyc if(chq(j).lt.chqmin)chqmin = chq(j) 300 continue endif do 310 j=1,ncyc kcyc = nchini-ncyc+j C WRITE OF ALL SETS do 362 j1=1,ni write(44,100)agei(j,1,j1),tab1,agei(j,2,j1) write(44,100)agei(j,1,j1+1),tab1,agei(j,2,j1) write(46,100)ageo(j,1,j1),tab1,ageo(j,2,j1) write(46,100)ageo(j,1,j1+1),tab1,ageo(j,2,j1) 362 continue do 322 j1=1,nmaxi(j) write(40,100)tti(j,1,j1),tab1,tti(j,2,j1) 322 continue do 342 j1=1,nmaxo(j) write(42,100)tto(j,1,j1),tab1,tto(j,2,j1) 342 continue if(kcyc.ne.nrun)then write(40,*)'& cycle #',kcyc write(42,*)'& cycle #',kcyc,' chisq= ',chisq*ni write(44,*)'& cycle #',kcyc write(46,*)'& cycle #',kcyc endif if(kcyc/30*30-kcyc.eq.0.and.kcyc.ne.nrun)then kc1=kc1+1 call wgraph(40,kc1) call wgraph(42,kc1) call wgraph(44,kc1) call wgraph(46,kc1) endif C WRITE OF SELECTED SETS if(chq(j).lt.chqmin*2.d0)then mcyc=mcyc+1 do 320 j1=1,nmaxi(j) write(10,100)tti(j,1,j1),tab1,tti(j,2,j1) 320 continue do 340 j1=1,nmaxo(j) write(12,100)tto(j,1,j1),tab1,tto(j,2,j1) 340 continue C CONTOUR MATRIX CALCULATION if(yes.eq."y".or.yes.eq."Y")then mi = nmaxi(j) mo = nmaxo(j) do 345 jc = 1,nsteps do 345 kc = 1,mps xk = kc xage = agein + (jc-1)*dage + xk/mps*dage call contour(xage,jc,tti,j,mi,matx1) call contour(xage,jc,tto,j,mo,matx2) 345 continue endif do 360 j1=1,ni write(18,100)agei(j,1,j1),tab1,agei(j,2,j1) write(18,100)agei(j,1,j1+1),tab1,agei(j,2,j1) write(19,100)ageo(j,1,j1),tab1,ageo(j,2,j1) write(19,100)ageo(j,1,j1+1),tab1,ageo(j,2,j1) 360 continue write(10,*)'& cycle #',kcyc write(12,*)'& cycle #',kcyc,' chisq= ',chisq*ni write(18,*)'& cycle #',kcyc write(19,*)'& cycle #',kcyc if(mcyc/30*30-mcyc.eq.0.and.kcyc.ne.nrun)then kc0=kc0+1 call wgraph(10,kc0) call wgraph(12,kc0) call wgraph(18,kc0) call wgraph(19,kc0) endif endif 310 continue ncyc=1 close(10) close(12) close(18) close(19) else ncyc=ncyc + 1 endif if(nchini.ne.nrun)goto 34 C WRITE CONTOUR MATRICES if(yes.eq."y".or.yes.eq."Y")then do 60 k2=1,ntemps write(25,110)(matx1(k1,k2),k1=1,nsteps) write(26,110) (matx2(k1,k2),k1=1,nsteps) do 65 k1=1,nsteps matx(k1) = matx2(k1,k2)-matx1(k1,k2) if(matx(k1).lt.0)matx(k1)=0.d0 65 continue write(27,110)(matx(k1),k1=1,nsteps) 60 continue endif C REMOVE LAST '&' FROM OUTPUT FILES ??? write(20,*)'End of CH cycles, total = ',nrun, 'completed =',mcyc write(20,*)'initial failed = ',nfail stop 'End of CH cycles' 100 format(1x,6(g20.8,a1)) 140 format(1x,4(f12.4,a1)) 130 format(1x,i2,' chisq= ',f12.2,f12.4, $' GOF= ',f9.5,' Iter#= ',i4) 120 format(1x,f5.2,6(f12.4)) 110 format(100I5) end C SUBROUTINE WGRAPH subroutine wgraph(n,k) character *2,graph if(k.lt.48.or.k.gt.57)return graph = "g"//char(k) write(n,120)"@WITH "//graph write(n,130)"@"//graph//" ON" 120 format(a8) 130 format(a6) return end C SUBROUTINE CONTOUR subroutine contour(xage,jc,t,j,n,mat) implicit double precision (a-h,o-z) C MIN TEMP 50C, MAX TEMP 600C parameter(ns=200,nt=1001,ny=10,tin=50.d0,tend=600.d0,ntemps=22) dimension t(ny,2,nt),mat(ns,ns) external temfc dtp = (tend-tin)/ntemps dtemp = temfc(t,j,n,xage) if(dtemp.gt.tend.or.dtemp.lt.tin)return ntem = (dtemp-tmpin)/dtp + 1 mat(jc,ntem)=mat(jc,ntem) + 1 return end C FUNCTION TEMFC double precision function temfc(t,jc,mi,x) implicit double precision (a-h,o-z) parameter(ntst=1001,nwcy=10) dimension t(nwcy,2,ntst) do 10 j=mi,1,-1 if(t(jc,1,j).ge.x.and.t(jc,1,j+1).lt.x)then slope = (t(jc,2,j+1)-t(jc,2,j))/(t(jc,1,j+1)-t(jc,1,j)) temfc=t(jc,2,j)+ slope * (x-t(jc,1,j)) mi = j return endif 10 continue end C SUBROUTINE RANCHIST subroutine ranchist(c,c0,mc,sa,sb,idem) implicit double precision (a-h,o-z) dimension c(mc),c0(mc) sma = 0.d0 smb = 0.d0 do 30 l=3,mc perc = 5.d0 * ran1(idem)*(0.6 - ran1(idem)) c(l)=c0(l) + perc sma = sma + c(l)*(1+(-1)**(l+1)) smb = smb + c(l)*(1+(-1)**l)/2.d0 30 continue c(1)=sa+sb-sma c(2)=(sa-sb)/2.d0-smb 20 continue return end C FUNCTION RAN1 double precision function ran1(idum) implicit double precision (a-h,o-z) dimension r(97) parameter (m1=259200,ia1=7141,ic1=54773,rm1=3.8580247e-6) parameter (m2=134456,ia2=8121,ic2=28411,rm2=7.4373773e-6) parameter (m3=243000,ia3=4561,ic3=51349) data iff /0/ if (idum.lt.0.or.iff.eq.0) then iff=1 ix1=mod(ic1-idum,m1) ix1=mod(ia1*ix1+ic1,m1) ix2=mod(ix1,m2) ix1=mod(ia1*ix1+ic1,m1) ix3=mod(ix1,m3) do 11 j=1,97 ix1=mod(ia1*ix1+ic1,m1) ix2=mod(ia2*ix2+ic2,m2) r(j)=(float(ix1)+float(ix2)*rm2)*rm1 11 continue idum=1 endif ix1=mod(ia1*ix1+ic1,m1) ix2=mod(ia2*ix2+ic2,m2) ix3=mod(ia3*ix3+ic3,m3) j=1+(97*ix3)/m3 if(j.gt.97.or.j.lt.1)then write(20,*)'ERROR(RAN1): (j.gt.97.or.j.lt.1)' stop 'ERROR(RAN1): (j.gt.97.or.j.lt.1)' endif ran1=r(j) r(j)=(float(ix1)+float(ix2)*rm2)*rm1 return end C SUBROUTINE AGESAMP subroutine agesamp(r39,ys,wt,ni) implicit double precision (a-h,o-z) parameter (ns=200) dimension xs(0:ns),ys(0:ni),r39(0:ni+1),ya(0:ns) dimension wt(ni),sig(ns) character*9 tab1 tab1=char(09) do 10 j=1,ni read(24,*)sig(j) if(sig(j).le.0.)then write(20,*)'ERROR(AGESAMP): (sig(j).le.0)' stop 'ERROR(AGESAMP): (sig(j).le.0)' endif read(22,*)xs(j),ya(j) write(21,100)xs(j-1)*100.d0,tab1,ya(j)+sig(j) write(21,100)xs(j),tab1,ya(j)+sig(j) xs(j)=xs(j)/100.d00 wt(j)=1.d00/dsqrt(r39(j)-r39(j-1))*sig(j) 10 continue do 90 k=1,ni write(21,100)xs(ni-k+1)*100.,tab1,ya(ni-k+1)-sig(ni-k+1) write(21,100)xs(ni-k)*100.,tab1,ya(ni-k+1)-sig(ni-k+1) 90 continue write(21,100)xs(0),tab1,ya(1)+sig(1) ya(0)=ya(1) xs(ni+1)=1.d00 ys(ni+1)=ys(ni) jold = 1 do 30 k=1,ni ys(k)=0.d0 do 20 j=jold,ni+1 if(xs(j-1).lt.r39(k).and.xs(j).gt.r39(k))then if(jold.eq.j)then ys(k)=ya(j) else do 34 jm=jold+1,j-1 ys(k)=ya(jm)*(xs(jm)-xs(jm-1))+ys(k) 34 continue ys(k)=ys(k)+ya(jold)*(xs(jold)-r39(k-1))+ya(j)*(r39(k)-xs(j-1)) ys(k)=ys(k)/(r39(k)-r39(k-1)) jold=j endif goto 32 endif 20 continue 32 write(50,100)r39(k-1)*100.d0,tab1,ys(k) write(50,100)r39(k)*100.d0,tab1,ys(k),tab1,wt(k) 30 continue return 100 format(1x,4(f14.8,a1)) end C SUBROUTINE LAB subroutine lab(r39,d0,e,vc,telab,tilab,lc,xmat,ni,nst,nn) implicit double precision (a-h,o-z) parameter(pi=3.14159265359d00,imp=2, $b=8.d0,ns=200,R= 1.987d-03,nd=10) dimension xmat(1:nd,0:ns,1:nn),lc(nn),zita(0:ns),r39(0:ni) dimension tilab(ni),telab(ni),d0(nst),vc(nst) zita(0)=0.d00 nsq=2 nsq1=1 do 60 j=1,nn an2=(lc(j)*pi)**2 xmat(1,ni+1,j)=b/an2 if(j.ge.121)then if((j-121)/20*20-(j-121).eq.0)then nsq=nsq*2 nsq2=nsq/imp endif do 64 mi=1,nsq1-1 an=((lc(j)-mi*imp)*pi)**2 xmat(1,ni+1,j)=xmat(1,ni+1,j)+b/an*(nsq1-mi)/nsq1 64 continue if(j.eq.nn)goto 60 do 68 mi=1,nsq2-1 an=((lc(j)+mi*imp)*pi)**2 xmat(1,ni+1,j)=xmat(1,ni+1,j)+b/an*(nsq2-mi)/nsq2 68 continue nsq1=nsq2 endif 60 continue do 10 k=1,nst do 10 nt=1,ni zita(nt)=zita(nt-1)+d0(k)*tilab(nt)*dexp(-e/(r*telab(nt))) nsq=2 nsq1=1 do 48 j=1,nn an2=(lc(j)*pi)**2 if(an2*zita(nt).gt.60)then do 70 jx=j,nn xmat(k,nt,jx)=xmat(1,ni+1,jx) 70 continue goto 72 else xmat(k,nt,j)=xmat(1,ni+1,j)-b*dexp(-an2*zita(nt))/an2 endif if(j.ge.121)then if((j-121)/20*20-(j-121).eq.0)then nsq=nsq*2 nsq2=nsq/imp endif do 42 mi=1,nsq1-1 an=((lc(j)-mi*imp)*pi)**2 xmat(k,nt,j)=xmat(k,nt,j)-b*dexp(-an*zita(nt)) $ /an* (nsq1-mi)/nsq1 42 continue if(j.eq.nn)goto 50 do 44 mi=1,nsq2-1 an=((lc(j)+mi*imp)*pi)**2 xmat(k,nt,j)=xmat(k,nt,j)-b*dexp(-an*zita(nt)) $ /an*(nsq2-mi)/nsq2 44 continue nsq1=nsq2 endif 50 continue 48 continue 72 do 20 k1=1,80000,imp an2=(k1*pi)**2 r39(nt)=r39(nt)+b/an2*vc(k) if (zita(nt)*(k1*pi)**2.gt.100.) then ab = 0.d00 else ab= dexp(-zita(nt)*(k1*pi)**2)/(k1*pi)**2 r39(nt)=r39(nt)-b*ab*vc(k) endif 20 continue 10 continue return end C SUBROUTINE CHEBFT subroutine chebft(a,b,c,n,func,tt,np) implicit double precision (a-h,o-z) parameter (nmax=100, pi=3.141592653589793d0) dimension c(n),f(nmax),tt(2,0:np) bma=0.5d00*(b-a) bpa=0.5d00*(b+a) do 11 k=1,n y=dcos(pi*(k-0.5d0)/n) f(k)=func(tt,np,y*bma+bpa) 11 continue fac=2.d0/n do 13 j=1,n sum=0.d0 do 12 k=1,n sum=sum+f(k)*dcos((pi*(j-1))*((k-0.5d0)/n)) 12 continue c(j)=fac*sum 13 continue return end C DOUBLE PRECISION FUNCTION TEMP double precision function ftemp(a,n,x) implicit double precision (a-h,o-z) dimension a(2,0:n) do 10 j=2,n if(a(1,j).lt.x.and.a(1,j-1).gt.x)then ftemp=a(2,j-1)-(a(1,j-1)-x)*(a(2,j-1)-a(2,j))/(a(1,j-1)-a(1,j)) ftemp= dsqrt(ftemp - 273.d0) return endif 10 continue end C SUBROUTINE MRQMIN subroutine mrqmin(x,y,sig,ndata,a,ma,lista,mfit, * covar,alpha,nca,chisq,alamda,tt,xmat, * lc,vc,d0,E,nst,nn,sa,sb) implicit double precision (a-h,o-z) parameter (mmax=100,ntst=1001,nd=10,ns=200) dimension x(0:ndata),y(0:ndata),sig(ndata),a(ma),lista(ma), * covar(nca,nca),alpha(nca,nca),atry(mmax),beta(mmax),da(mmax), * xmat(1:nd,0:ns,1:nn),d0(nst),vc(nst),tt(2,0:ntst),lc(nn) * ,covaux(mmax,mmax),daux(mmax) t0 = tt(1,1) if(alamda.lt.0.)then chisq=0.d0 kk=mfit+1 do 12 j=1,ma ihit=0 do 11 k=1,mfit if(lista(k).eq.j)ihit=ihit+1 11 continue if (ihit.eq.0) then lista(kk)=j kk=kk+1 else if (ihit.gt.1) then write(20,*)'ERROR(MRQMIN): improper permutation in lista-1' stop 'ERROR(MRQMIN): improper permutation in lista-1' endif 12 continue if (kk.ne.(ma+1))then write(20,*)'ERROR(MRQMIN): improper perm. in lista-2' stop 'ERROR(MRQMIN): improper perm. in lista-2' endif if(chebev(0.,t0,a,ma,0.)**2.gt.tt(2,1))then write(20,*)'BAD INITIAL: temp > Maxtemp at age=0' chisq=0.d0 return endif call mrqcof(x,y,sig,ndata,a,ma,lista,mfit,alpha,beta,nca,chisq, *tt,xmat,lc,vc,d0,e,nst,nn,alamda) if(chisq.eq.0)return alamda=0.001d0 ochisq=chisq do 13 j=1,ma atry(j)=a(j) 13 continue endif do 15 j=1,mfit do 14 k=1,mfit covar(j,k)=alpha(j,k) 14 continue covar(j,j)=alpha(j,j)*(1.d0+alamda) da(j)=beta(j) 15 continue ng = 1 call gaussj(covar,mfit,nca,da,1,1) if(ng.eq.-1)stop 'Singular Matrix' if(alamda.eq.0.)then call covsrt(covar,nca,ma,lista,mfit) return endif sma=0.d0 smb=0.d0 C CHANGE THE "ATRY(1) AND ATRY(2)" INITIAL AND FINAL TEMP CONSTRAINT do 16 j=1,mfit if(alamda.ge.1.00d20)then da(j)=0.d0 endif atry(lista(j))=a(lista(j))+da(j) sma = sma + atry(lista(j))*(1+(-1)**(lista(j)+1)) smb = smb + atry(lista(j))*(1+(-1)**lista(j))/2.d0 16 continue atry(1)=sa+sb-sma atry(2)=(sa-sb)/2.d0-smb call mrqcof(x,y,sig,ndata,atry,ma,lista,mfit,covar,da,nca,chisq, *tt,xmat,lc,vc,d0,e,nst,nn,alamda) if(chisq.le.ochisq)then C CHECK IF COVAR IS NOT SINGULAR do 19 j=1,mfit do 20 k=1,mfit covaux(j,k)=covar(j,k) 20 continue daux(j)=da(j) 19 continue call gaussj(covaux,mfit,nca,daux,ng,1) if(ng.eq.-1)then alamda=10.d0*alamda chisq=ochisq ng = 1 return endif C END CHECKING FOR SINGULAR MATRIX alamda=0.1d0*alamda ochisq=chisq do 18 j=1,mfit do 17 k=1,mfit alpha(j,k)=covar(j,k) 17 continue beta(j)=da(j) a(lista(j))=atry(lista(j)) 18 continue a(1)=atry(1) a(2)=atry(2) else alamda=10.d0*alamda chisq=ochisq endif return end C SUBROUTINE MRQCOF subroutine mrqcof(x,y,sig,ndata,a,ma,lista,mfit,alpha,beta,nalp, *chisq,tt,xmat,lc,vc,d0,e,nst,nn,al) implicit double precision (a-h,o-z) parameter (mmax=100,ntst=1001,perc=0.01d0,xlambd=0.0005543d0, *nc=100,mxi=500,nd=10,ns=200) dimension x(0:ndata),y(0:ndata),sig(ndata),alpha(nalp,nalp), * dyda(mmax),lista(mfit),a(ma),xmat(1:nd,0:ns,1:nn),d0(nst), * vc(nst),tt(2,0:ntst),lc(nn),beta(ma),xi(0:2*nc,nd,mxi) t0=tt(1,1) c40=1.d0-dexp(-xlambd*tt(1,0)) do 12 j=1,mfit do 11 k=1,j alpha(j,k)=0.d0 11 continue beta(j)=0.d0 12 continue call zita(a,ma,lista,mfit,lc,xi,t0,d0,E,tt,nst,nn,perc,al) if(tt(2,ntst).eq.-1.)then tt(2,ntst)=0.d0 chisq=chisq*2.d0 return endif chisq=0.d0 do 15 i=1,ndata call age(i,a,ymod,dyda,ma,lista,mfit,vc,xmat,x,xi,c40,nn, * nst,ndata,perc,al) sig2i=1.d0/(sig(i)*sig(i)) dy=y(i)-ymod do 14 j=1,mfit wt=dyda(lista(j))*sig2i do 13 k=1,j alpha(j,k)=alpha(j,k)+wt*dyda(lista(k)) 13 continue beta(j)=beta(j)+dy*wt 14 continue chisq=chisq+dy*dy*sig2i 15 continue do 17 j=2,mfit do 16 k=1,j-1 alpha(k,j)=alpha(j,k) 16 continue 17 continue return end C SUBROUTINE COVSRT subroutine covsrt(covar,ncvm,ma,lista,mfit) implicit double precision (a-h,o-z) dimension covar(ncvm,ncvm),lista(mfit) do 12 j=1,ma-1 do 11 i=j+1,ma covar(i,j)=0.d0 11 continue 12 continue do 14 i=1,mfit-1 do 13 j=i+1,mfit if(lista(j).gt.lista(i)) then covar(lista(j),lista(i))=covar(i,j) else covar(lista(i),lista(j))=covar(i,j) endif 13 continue 14 continue swap=covar(1,1) do 15 j=1,ma covar(1,j)=covar(j,j) covar(j,j)=0.d0 15 continue covar(lista(1),lista(1))=swap do 16 j=2,mfit covar(lista(j),lista(j))=covar(1,j) 16 continue do 18 j=2,ma do 17 i=1,j-1 covar(i,j)=covar(j,i) 17 continue 18 continue return end C SUBROUTINE GAUSSJ subroutine gaussj(a,n,np,b,m,mp) implicit double precision (a-h,o-z) parameter (nmax=100) dimension a(np,np),b(np,mp),ipiv(nmax),indxr(nmax),indxc(nmax) do 11 j=1,n ipiv(j)=0 11 continue do 22 i=1,n big=0.d0 do 13 j=1,n if(ipiv(j).ne.1)then do 12 k=1,n if (ipiv(k).eq.0) then if (dabs(a(j,k)).ge.big)then big=dabs(a(j,k)) irow=j icol=k endif else if (ipiv(k).gt.1) then write(20,*)'ERROR(GAUSS): singular matrix-1' m= -1 return endif 12 continue endif 13 continue ipiv(icol)=ipiv(icol)+1 if (irow.ne.icol) then do 14 l=1,n dum=a(irow,l) a(irow,l)=a(icol,l) a(icol,l)=dum 14 continue do 15 l=1,m dum=b(irow,l) b(irow,l)=b(icol,l) b(icol,l)=dum 15 continue endif indxr(i)=irow indxc(i)=icol if (a(icol,icol).eq.0.) then write(20,*)'ERROR(GAUSS): singular matrix-2' m= -1 return endif pivinv=1.d0/a(icol,icol) a(icol,icol)=1.d0 do 16 l=1,n a(icol,l)=a(icol,l)*pivinv 16 continue do 17 l=1,m b(icol,l)=b(icol,l)*pivinv 17 continue do 21 ll=1,n if(ll.ne.icol)then dum=a(ll,icol) a(ll,icol)=0.d0 do 18 l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum 18 continue do 19 l=1,m b(ll,l)=b(ll,l)-b(icol,l)*dum 19 continue endif 21 continue 22 continue do 24 l=n,1,-1 if(indxr(l).ne.indxc(l))then do 23 k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum 23 continue endif 24 continue return end C SUBROUTINE ZITA subroutine zita(c,mc,lista,mfit,lc,xi,t0,d0,E,tt,nst, $ nn,perc,al) implicit double precision (a-h,o-z) parameter(ntst=1001,a=0,nc=100,ns=200) dimension lc(nn),d0(nst),lista(mfit) dimension c(mc),tt(2,0:ntst),cder(nc),xi(0:2*nc,1:nst,1:nn) Parameter (nwcy=10) common /var/ tti(nwcy,2,ntst),tto(nwcy,2,ntst), $agei(nwcy,2,ns),ageo(nwcy,2,ns) common /cont/ ncyc,nmaxi(nwcy),nmaxo(nwcy) n=0 call chder(a,t0,c,cder,mc) call sched(c,cder,mc,tt,imax) if(imax.eq.0)then tt(2,ntst)=-1.d0 return endif call geosec(imax,d0,E,tt,xi,nn,lc,nst,n) write(20,*)'# of CH points= ',imax if(al.le.0.)then if (al.lt.0) then nmaxi(ncyc)=imax+1 else nmaxo(ncyc)=imax+1 endif do 40 j=1,imax+1 if (al.lt.0) then tti(ncyc,1,j)=tt(1,j-1) tti(ncyc,2,j)=tt(2,j-1) else tto(ncyc,1,j)=tt(1,j-1) tto(ncyc,2,j)=tt(2,j-1) endif 40 continue endif do 20 k=1,mfit csk=c(lista(k)) jm=2+lista(k)/2*2-lista(k) csj=c(jm) dc=perc * dabs(c(lista(k))) do 30 js=1,2 n=n+1 c(lista(k))=c(lista(k))+(-1)**js*dc c(jm)=c(jm)-dc*(-1)**js*2.d0/jm do 42 j=2,imax tt(2,j)=chebev(a,t0,c,mc,tt(1,j))**2 42 continue call geosec(imax,d0,E,tt,xi,nn,lc,nst,n) c(lista(k))=csk c(jm)=csj 30 continue 20 continue return end C SUBROUTINE AGE subroutine age(nt,c,y,dyda,mc,lista,mfit,vc,xmat,r39,xi,c40, *nn,nst,ndata,perc,al) implicit double precision (a-h,o-z) parameter(xlambd=0.0005543d0,nc=100,nd=10,ns=200,ntst=1001) dimension xmat(1:nd,0:ns,1:nn),r39(0:ndata),lista(mfit) dimension c(mc),dyda(mc),vc(nst),xi(0:2*nc,1:nst,1:nn) parameter(nwcy=10) common /var/ tti(nwcy,2,ntst),tto(nwcy,2,ntst), $agei(nwcy,2,ns),ageo(nwcy,2,ns) common /cont/ ncyc,nmaxi(nwcy),nmaxo(nwcy) c39wd=1.d0 - c40 s1=0.d0 do 52 k=1,nst do 52 m=1,nn s1=s1 + vc(k)*xi(0,k,m)*(xmat(k,nt,m)-xmat(k,nt-1,m)) 52 continue dr39=r39(nt)-r39(nt-1) dr40=dr39*(c40-1)+s1 y=0.d0 if(dr40/dr39.gt.0.)then y=dlog(1.d0+dr40/dr39/c39wd)/xlambd dya=1.d0/(dr39*c39wd+dr40)/xlambd endif if(al.lt.0)then agei(ncyc,1,nt+1)=r39(nt)*100.d0 agei(ncyc,2,nt)=y else if (al.eq.0.) then ageo(ncyc,1,nt+1)=r39(nt)*100.d0 ageo(ncyc,2,nt)=y endif 42 continue do 56 ls=1,mfit l=(lista(ls)-2)*2-1 dc=perc*dabs(c(lista(ls))) s1=0.d0 do 54 k=1,nst do 54 m=1,nn dxi= (xi(l+1,k,m)-xi(l,k,m))/(2*dc) s1=s1 + vc(k)*(xmat(k,nt,m)-xmat(k,nt-1,m))*dxi 54 continue dyda(lista(ls))=s1*dya 56 continue return end C SUBROUTINE GEOSEC subroutine geosec(imax,d0,E,tt,xi,nn,lc,nst,j1) implicit double precision (a-h,o-z) parameter(xlambd=5.543d-04,R=1.987d-03,max=1002, $pi=3.14159265359d0,nc=100) dimension xi(0:2*nc,1:nst,1:nn),tt(2,0:imax),dzita(0:max) dimension lc(nn),d0(nst),d(0:max) do 80 k=1,nst nend=imax-1 do 10 j=0,nend avtemp=(tt(2,j+1)+tt(2,j))/2.d0 + 273.d0 if(e/r/avtemp.gt.80)then d(j)=0.d0 goto 10 endif d(j)=d0(k)*dexp(-e/r/avtemp) 10 continue 60 dzita(nend+1)=0 do 20 j=nend,0,-1 20 dzita(j)=dzita(j+1)+dabs(tt(1,j+1)- $tt(1,j))*d(j) C COMPUTO OF XI(M) - M<80000 --> nn=319 do 30 mi=1,nn xlogm=2*dlog(lc(mi)*pi) sum=0.d0 do 40 n=0,nend if(d(n).eq.0.)goto 40 uplus=(lc(mi)*pi)**2*dzita(n+1)+xlambd*(tt(1,0)-tt(1,n+1)) if (uplus-xlogm.gt.25.)goto 40 al=(lc(mi)*pi)**2*d(n)-xlambd xal= al *(tt(1,n)-tt(1,n+1)) if (dabs(xal).gt.30)then camal=1.d0/al else if (dabs(xal).gt.0.001)then camal=(1.d0-dexp(-xal))/al else camal=(tt(1,n)-tt(1,n+1))*(1-xal/2+xal**2/6.d0) endif endif sum=sum+d(n)*dexp(-uplus)*camal 40 continue xi(j1,k,mi)=sum*(pi*lc(mi))**2 30 continue 80 continue return end C SUBROUTINE SCHED (CH PARTITIONING) subroutine sched(c,cder,m,tt,i) implicit double precision (a-h,o-z) parameter (ntst=1001,xacc=0.01d0,ti=50.d0) dimension tt(2,0:ntst),c(m),cder(m) C DTP = TEMPERATURE STEP OF THE PARTITION C CHECK FOR TEMPERATURE GREATER THAN MAX TEMP. dti = tt(1,1)/10000.d0 a= 0.d0 t0=tt(1,1) ki=1 ntp=10 call tmax(c,cder,m,tt,tini,tpini) if(tini.lt.tt(1,1))then tt(1,2)=tini tt(2,2)=tpini I=3 else I=2 endif C CALCULATION OF THE CH PARTITION 1 continue do 10 k=1,10000 x=tt(1,i-1) - k*dti if(x.lt.0)then tt(1,i)=0.d0 tt(2,i)=chebev(a,t0,c,m,tt(1,i))**2 return endif tchk =dabs(chebev(a,t0,c,m,x)**2-tt(2,i-1)) if(tchk.gt.ntp)then sign = (chebev(a,t0,c,m,x)**2-tt(2,i-1))/tchk tchk=tt(2,i-1)+sign*ntp tt(1,i)=rtsafe(x,tt(1,i-1),xacc,c,cder,m,tt(1,1),tchk) tt(2,i)=chebev(a,t0,c,m,tt(1,i))**2 if(tt(2,i).lt.400.d0.and.tt(2,i).gt.100.d0)then ntp = 4 else ntp = 10 endif i= i + 1 goto 20 endif 10 continue 20 if(i.gt.ntst-1)then i=0 return endif goto 1 end C SUBROUTINE TMAX (Check for temp > tempmax = 600C) subroutine tmax(c,cder,m,tt,age,temp) implicit double precision (a-h,o-z) parameter (ntst=1001,xacc=0.01d0) dimension tt(2,0:ntst),c(m),cder(m) a = 0.d0 ag4 = tt(1,1) t0 = tt(1,1) do 10 j=1,1000 age = tt(1,1)/1000.d0*j temp = chebev(a,t0,c,m,age)**2 if(temp.eq.tt(2,1))return if (temp.gt.tt(2,1)) then age= rtsafe(a,age,xacc,c,cder,m,tt(1,1),tt(2,1)) temp = chebev(a,t0,c,m,age)**2 return endif 10 continue age = tt(1,1) temp = tt(2,1) return end C SUBROUTINE CHEBEV function chebev(a,b,c,m,x) implicit double precision (a-h,o-z) dimension c(m) if ((x-a)*(x-b).gt.0.)then write(20,*)'ERROR(CHEBEV): x not in range' stop 'ERROR(CHEBEV): x not in range' endif d=0.d0 dd=0.d0 y=(2.d0*x-a-b)/(b-a) y2=2.d0*y do 11 j=m,2,-1 sv=d d=y2*d-dd+c(j) dd=sv 11 continue chebev=y*d-dd+0.5d0*c(1) return end C SUBROUTINE CHDER subroutine chder(a,b,c,cder,n) implicit double precision (a-h,o-z) dimension c(n),cder(n) cder(n)=0.d0 cder(n-1)=2*(n-1)*c(n) if(n.ge.3)then do 11 j=n-2,1,-1 cder(j)=cder(j+2)+2*j*c(j+1) 11 continue endif con=2.d0/(b-a) do 12 j=1,n cder(j)=cder(j)*con 12 continue return end C DOUBLE PRECISION FUNCTION RTSAFE double precision function rtsafe(x1,x2,xacc,c,cder,mc,t0,te) implicit double precision (a-h,o-z) parameter (maxit=100) dimension c(mc),cder(mc) call funcd(x1,fl,df,c,cder,t0,mc,te) call funcd(x2,fh,df,c,cder,t0,mc,te) if(fl*fh.ge.0.) then write(20,*)'ERROR(RTSAFE): root must be bracketed' stop 'ERROR(RTSAFE): root must be bracketed' endif if(fl.lt.0.)then xl=x1 xh=x2 else xh=x1 xl=x2 swap=fl fl=fh fh=swap endif rtsafe=.5d0*(x1+x2) dxold=abs(x2-x1) dx=dxold call funcd(rtsafe,f,df,c,cder,t0,mc,te) do 11 j=1,maxit if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).ge.0. * .or. abs(2.d0*f).gt.abs(dxold*df) ) then dxold=dx dx=0.5d0*(xh-xl) rtsafe=xl+dx if(xl.eq.rtsafe)return else dxold=dx dx=f/df temp=rtsafe rtsafe=rtsafe-dx if(temp.eq.rtsafe)return endif if(abs(dx).lt.xacc) return call funcd(rtsafe,f,df,c,cder,t0,mc,te) if(f.lt.0.) then xl=rtsafe fl=f else xh=rtsafe fh=f endif 11 continue write(20,*)'ERROR(RTSAFE): exceeding maximum iterations' stop 'ERROR(RTSAFE): exceeding maximum iterations' end C SUBROUTINE FUNCD subroutine funcd(x,fn,df,c,cder,t0,m,te) implicit double precision (a-h,o-z) dimension c(m), cder(m) fn=chebev(0.d0,t0,c,m,x)**2-te + tsf df=2.* chebev(0.d0,t0,cder,m,x)*chebev(0.d0,t0,c,m,x) return end double precision function gammq(a,x) implicit double precision (a-h,o-z) if(x.lt.0..or.a.le.0.)then write(20,*)'ERROR(GAMMQ): (x.lt.0..or.a.le.0.)' stop 'ERROR(GAMMQ): (x.lt.0..or.a.le.0.)' endif if(x.lt.a+1.)then call gser(gamser,a,x,gln) gammq=1.-gamser else call gcf(gammcf,a,x,gln) gammq=gammcf endif return end subroutine gcf(gammcf,a,x,gln) implicit double precision (a-h,o-z) parameter (itmax=100,eps=3.e-7) gln=gammln(a) gold=0. a0=1. a1=x b0=0. b1=1. fac=1. do 11 n=1,itmax an=float(n) ana=an-a a0=(a1+a0*ana)*fac b0=(b1+b0*ana)*fac anf=an*fac a1=x*a0+anf*a1 b1=x*b0+anf*b1 if(a1.ne.0.)then fac=1./a1 g=b1*fac if(dabs((g-gold)/g).lt.eps)go to 1 gold=g endif 11 continue write(20,*)'ERROR(GCF): a too large, itmax too small' stop 'ERROR(GCF): a too large, itmax too small' 1 gammcf=dexp(-x+a*dlog(x)-gln)*g return end subroutine gser(gamser,a,x,gln) implicit double precision (a-h,o-z) parameter (itmax=100,eps=3.e-7) gln=gammln(a) if(x.le.0.)then if(x.lt.0.)then write(20,*)'ERROR(GSER): (x.lt.0.)' stop 'ERROR(GSER): (x.lt.0.)' endif gamser=0. return endif ap=a sum=1./a del=sum do 11 n=1,itmax ap=ap+1. del=del*x/ap sum=sum+del if(dabs(del).lt.dabs(sum)*eps)go to 1 11 continue write(20,*) 'ERROR(GSER): a too large, itmax too small' stop 'ERROR(GSER): a too large, itmax too small' 1 gamser=sum*dexp(-x+a*dlog(x)-gln) return end function gammln(xx) C real*8 cof(6),stp,half,one,fpf,x,tmp,ser implicit double precision (a-h,o-z) dimension cof(6) data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0, * -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/ data half,one,fpf/0.5d0,1.0d0,5.5d0/ x=xx-one tmp=x+fpf tmp=(x+half)*dlog(tmp)-tmp ser=one do 11 j=1,6 x=x+one ser=ser+cof(j)/x 11 continue gammln=tmp+dlog(stp*ser) return end subroutine readfile(l) 10 read(l,*,end=20) goto 10 20 return end