program chem3 implicit none save include 'chem3.dek' c..chemical evolution code; version 3 c..this version supports an arbitrary number of isotopes, c..arbitrary number of zones, solar masses in solar masses out formulation, c..metallicity and mass dependent nucleosynthesis, c..different types of type1a supernova and classical novae, c..gauss-legendre quadrature, and a cash-karp stepper. c..local variables external imf,gasevo,rkqs integer nok,nbad,i,j,k,nz double precision tbeg,stpmin,mtot,mdot,sum,xnow,ggg,zzz, 1 xxx(nqmax),www(nqmax) c..format statements 01 format(5(1x,1pe10.2,' ',a,' ')) 02 format(1x,a,i6,a,i6,a,i6) 03 format(1x,'cpu time : ',i2,' hrs ',i2,' min ',i2,' sec') 04 format(1x,a,i10,a,i10) c..read the users input deck call rchem c..get todays date and time call today(adat,atim) c..start the clock call zsecond(timtot) c..read in the stellar nucleosynthesis call rnuc neq = nzone * niso write(6,02) ' ' write(6,02) 'number of equations :',neq write(6,02) ' ' if (neq .gt. gphys) then write(6,02) 'neq',neq,' gphys',gphys stop 'neq > gphys' end if c..get the imf normalization call gauleg(imfmin,mbmin,xxx,www,nnhi) call qgaus(imf,xxx,www,nnhi,ximf_norm) if (mbmax .lt. imfmax) then call gauleg(mbmin,mbmax,xxx,www,nnhi) call qgaus(imf,xxx,www,nnhi,sum) ximf_norm = ximf_norm + sum call gauleg(mbmax,imfmax,xxx,www,nnhi) call qgaus(imf,xxx,www,nnhi,sum) ximf_norm = ximf_norm + sum else call gauleg(mbmin,imfmax,xxx,www,nnhi) call qgaus(imf,xxx,www,nnhi,sum) ximf_norm = ximf_norm + sum end if c..initialize the time history table hit counters khist = 0 kmiss = 0 khyes = 0 khno = 0 kt1miss = 0 kt1yes = 0 kt1no = 0 knvmiss = 0 knvyes = 0 knvno = 0 c..initialize the solution with bigbang abundances tbeg = 0.0d0 stpmin = 1.0e-5 * stptry kount = 1 do nz=1,nzone nowzon = nz call surfden(tbeg,mtot,mdot) do i=1,niso j = i + (nz - 1)*niso init(j) = mtot * bigb(i) enddo enddo c..now integrate the odes contained in routine gasevo c..for no changes in the infall composition c..integrate from the begginning to the end if (thalo .ge. tend) then call chemint(tbeg,stpmin,tend,toler,neq,nok,nbad,gasevo,rkqs) c..allow for a change in the infall composition at time thalo else if (thalo .lt. tend) then c..integrate from tbeg to thalo call chemint(tbeg,stpmin,thalo,toler,neq,nok,nbad,gasevo,rkqs) write(6,02) 'nok =',nok,' nbad=',nbad,' kount=',kount c..change the infall composition at thalo write(6,02) ' ' write(6,02) 'changing infall composition:' nowzon = nrcomp sum = 0.0d0 do i=1,niso if (i .ne. nal26 .and. i .ne. nfe60) then nowiso = i call getgas(thalo,ggg,zzz,bigb(i)) sum = sum + bigb(i) end if enddo xnow = 1.0d0 - sum write(6,01) xnow,'mass non conservation of new infall' c..normalize the new infall composition do i=1,niso bigb(i) = bigb(i)/sum enddo c..write out the new infall composition write(6,*) 'new infall composition:' write(6,01) (bigb(i),namnuc(i), i=1,niso) c..and finish off the integration call chemint(thalo,stpmin,tend,toler,neq,nok,nbad,gasevo,rkqs) end if c..output the model call modout c..write out the table hit statistics write(6,*) write(6,04) 'maximum history storage depth =',khist write(6,04) 'direct history hits =',khyes write(6,04) 'searched history hits =',khno write(6,04) 'missed histories =',kmiss if (sn1fac .ne. 0.0) then write(6,*) write(6,04) 'direct type1a history hits =',kt1yes write(6,04) 'searched type1a history hits =',kt1no write(6,04) 'missed typ1a histories =',kt1miss end if if (novfac .ne. 0.0) then write(6,*) write(6,04) 'direct nova history hits =',knvyes write(6,04) 'searched nova history hits =',knvno write(6,04) 'missed nova histories =',knvmiss end if c..close up shop call zsecond(timtot) call timlap(timtot,hours,minuts,secs,msecs) write(6,*) write(6,03) hours,minuts,secs stop 'normal termination' end subroutine chemint(start,stpmin,stopp, 1 eps,ylogi,nok,nbad,derivs,steper) implicit none save include 'chem3.dek' c..ode integration driver c..input: c..start = c..stpmin = c..stopp = c..eps = c..ylogi = c..derivs = c..steper = c..output: c..nok = c..nbad = c..declare external derivs,steper integer nok,nbad,stpmax, 1 ylogi,i,j,nstp parameter (stpmax=50000) double precision stpmin,eps,stopp,start, 1 yscal(gphys),y(gphys),dydx(gphys), 2 x,h,hdid,hnext,zero parameter (zero=0.0d0) c..here are the format statements for printouts as we integrate 100 format(1x,i5,1p2e12.4) 101 format(1x,1p12e10.2) c..initialize if (ylogi .gt. gphys) stop 'ylogi > yphys in routine chemint' x = start h = sign(stptry,stopp-start) nok = 0 nbad = 0 c..store the first step do i=1,ylogi y(i) = init(i) enddo c..take at most stpmax steps do nstp=1,stpmax call derivs(x,y,dydx) c..just for counting purposes if (nstp .eq. 1) kount = kount - 1 c..scaling vector used to monitor accuracy do i=1,ylogi yscal(i) = max(y(i),odescal) enddo c..store intermediate results if ( kount .lt. (tphys-1) ) then c if ( (modal26 .eq. 0) .or. c 1 (modal26 .ne. 0 .and. mod(nstp,modal26) .eq. 1) .or. c 2 (kount .lt. 30) ) then kount = kount+1 time(kount) = x do i=1,ylogi grk(i,kount) = y(i) dgrk(i,kount) = dydx(i) enddo write(6,100) kount,time(kount),grk(1,kount) if (iprint .eq. 1) then write(6,101) (grk(j,kount), j=1,ylogi) end if c end if end if c..if the step can overshoot the stop point then cut it if ((x+h-stopp)*(x+h-start).gt.zero) h = stopp - x c..do an integration step call steper(y,dydx,ylogi,x,h,eps,yscal,hdid,hnext,derivs) if (hdid.eq.h) then nok = nok+1 else nbad = nbad+1 end if c..this is the normal exit point, save the final step if ( (nstp .eq. stpmax) .or. 1 ((x-stopp)*(stopp-start) .ge. zero) .or. 2 (kount+1 .eq. tphys) ) then do i=1,ylogi init(i) = y(i) enddo kount = kount+1 time(kount) = x do i=1,ylogi grk(i,kount) = y(i) dgrk(i,kount) = dydx(i) enddo write(6,100) kount,time(kount),grk(1,kount) if (iprint .eq. 1) then write(6,101) (grk(j,kount), j=1,ylogi) end if return end if c..set the step size for the next iteration; stay above stpmin h=hnext if (abs(hnext).lt.stpmin) stop 'hnext < stpmin in chemint' c..back for another iteration or death enddo write(6,*) '> than stpmax steps required in chemint' return end double precision function taum(m) implicit none save c..given the stellar mass, this routine returns the main-sequence lifetime c..declare the pass double precision m c..the talbott & arnett 1971 homology relation c.. double precision m,tsun,tinf c.. parameter (tsun = 11.7d0, tinf = 1.0d-3) c.. taum = tsun / (m*m) + tinf c..ww95 for high mass stars, schaller et al for low mass stars integer iat,mtmax parameter (mtmax=27) double precision a,b,mstar(mtmax),tstar(mtmax) data mstar / 100.0d0, 40.217d0, 35.190d0, 30.163d0, 25.136d0, 1 22.119d0, 20.109d0, 19.000d0, 18.098d0, 15.081d0, 2 13.071d0, 12.065d0, 11.065d0, 9.0d0, 7.0d0, 3 5.0d0, 4.0d0, 3.0d0, 2.5d0, 2.0d0, 4 1.7d0, 1.5d0, 1.25d0, 1.0d0, 0.9d0, 5 0.8d0, 0.08d0/ data tstar / 3.104d-3, 5.605d-3, 6.132d-3, 6.886d-3, 8.069d-3, 1 9.039d-3, 9.996d-3, 1.044d-2, 1.093d-2, 1.357d-2, 2 1.637d-2, 1.838d-2, 2.097d-2, 2.913d-2, 4.791d-2, 3 1.068d-1, 1.909d-1, 4.387d-1, 7.303d-1, 1.356d0, 4 1.827d0, 2.694d0, 3.948d0, 9.961d0, 1.550d1, 5 2.5027d1, 1.000d3/ data iat/10/ c..hard wired two point linear fit; only get the index i if need be if (m .gt. mstar(iat) .or. m .lt. mstar(iat+1)) then call locate(mstar,mtmax,m,iat) end if b = mstar(iat) - mstar(iat+1) b = (tstar(iat) - tstar(iat+1)) / b a = tstar(iat+1) - b * mstar(iat+1) taum = b*m + a c..the instantaneous recycling approximation c.. taum = tstar(2) return end double precision function mtau(t) implicit none save include 'chem3.dek' c..given a main-sequence lifetime, this routine returns the stellar mass c..declare the pass double precision t c..local variables double precision t100 c..the talbott & arnett 1972 homolgy relation c.. double precision tsun,tinf c.. parameter (tsun = 11.7d0, tinf = 1.0d-3) c.. mtau = mmax c.. if (t .ge. t100) mtau = sqrt(tsun/(t - tinf)) c..ww95 for high mass stars, schaller et al for low mass stars integer iat,mtmax,ifirst parameter (mtmax=27) double precision taum,a,b,mstar(mtmax),tstar(mtmax) data mstar / 100.0d0, 40.217d0, 35.190d0, 30.163d0, 25.136d0, 1 22.119d0, 20.109d0, 19.000d0, 18.098d0, 15.081d0, 2 13.071d0, 12.065d0, 11.065d0, 9.0d0, 7.0d0, 3 5.0d0, 4.0d0, 3.0d0, 2.5d0, 2.0d0, 4 1.7d0, 1.5d0, 1.25d0, 1.0d0, 0.9d0, 5 0.8d0, 0.08d0/ data tstar / 3.104d-3, 5.605d-3, 6.132d-3, 6.886d-3, 8.069d-3, 1 9.039d-3, 9.996d-3, 1.044d-2, 1.093d-2, 1.357d-2, 2 1.637d-2, 1.838d-2, 2.097d-2, 2.913d-2, 4.791d-2, 3 1.068d-1, 1.909d-1, 4.387d-1, 7.303d-1, 1.356d0, 4 1.827d0, 2.694d0, 3.948d0, 9.961d0, 1.550d1, 5 2.5027d1, 1.00d3/ data iat/10/, ifirst/0/ c..initialize if (ifirst .eq. 0) then t100 = taum(mmax) ifirst = 1 end if c..hard wired two point linear fit; only get the index i if need be mtau = mmax if (t .lt. t100) return if (t .lt. tstar(iat) .or. t .gt. tstar(iat+1)) then call locate(tstar,mtmax,t,iat) endif b = tstar(iat) - tstar(iat+1) b = (mstar(iat) - mstar(iat+1)) / b a = mstar(iat+1) - b * tstar(iat+1) mtau = b*t + a c..the instantaneous recycling c.. mtau = mmin + 1.0e-3 return end double precision function imf(m) implicit none save include 'chem3.dek' c..call the imf of choice c..declare the pass double precision m c..local variables double precision salp,smilin,sc86in c..for a salpeter imf if (wimf .eq. 1) then imf = salp(m) c..for a scalo-miller imf else if (wimf .eq. 2) then imf = smilin(m) c..scalo 86 imf else if (wimf .eq. 3) then imf = sc86in(m) end if return end double precision function brate(mtot,gtot) implicit none save include 'chem3.dek' c..given the total mass mtot and a gas mass gtot, this routine computes the c..schmidt stellar birth rate. c..declare the pass double precision mtot,gtot c..to try and avoid the power function, break into 3 cases if (kk .eq. 1.0) then brate = nu * gtot else if (kk .eq. 2.0) then brate = nu * gtot * gtot/mtot else brate = nu * mtot * (gtot/mtot)**kk end if return end subroutine gasevo(t,g,dgdt) implicit none save include 'chem3.dek' c..this routine sets up the right hand side of the ode's c..declare the pas double precision t,g(1),dgdt(1) c..local variables external gasbak,sn1int,novaint integer i,j,loff,nz,ii,ifirst double precision brate,birth,death,mtau,mlo,mtot, 1 mdot,sum,gastot,thalf_k40,tmean_k40, 2 brnch1,brnch2,fac, 3 thalf_al26,tmean_al26,thalf_fe60,tmean_fe60, 4 bdum c..half and mean lifetimes of k40, al26 and fe60 parameter (thalf_k40 = 1.28d0, 1 tmean_k40 = thalf_k40/0.693147181d0, 2 brnch1 = 0.893d0, 3 brnch2 = 0.107d0, 4 thalf_al26 = 7.40d-4, 5 tmean_al26 = thalf_al26/0.693147181d0, 6 thalf_fe60 = 1.50d-3, 7 tmean_fe60 = thalf_fe60/0.693147181d0) c..for the quadrature schemes double precision xa1(nqmax),wa1(nqmax),xa2(nqmax),wa2(nqmax), 1 xa3(nqmax),wa3(nqmax),xa4(nqmax),wa4(nqmax), 2 xa5(nqmax),wa5(nqmax) c..one-time initialization flag data ifirst/0/ c..initialize some of the quadrature weights if (ifirst .eq. 0) then ifirst = 1 call gauleg(mbmin,mbmax,xa1,wa1,nnmid) call gauleg(mbmax,mmax,xa2,wa2,nnhi) end if c..put this time point in the common block; get the mass this c..time point corresponds to, and then get the remaining c..the gaussian absiccas and weights timcom = t mlo = mtau(t) if (mlo .lt. mmax) then if (mlo.gt.mbmax) call gauleg(mlo,mmax,xa3,wa3,nnhi) if (mlo.gt.mbmin) call gauleg(mlo,mbmax,xa4,wa4,nnmid) if (mlo.gt.mmin) call gauleg(mlo,mbmin,xa5,wa5,nnlo) end if c..no negative gas masses do i=1,neq g(i) = max(g(i),0.0d0) enddo c..for each zone do nz=1,nzone nowzon = nz loff = (nz-1)*niso c..total mass and accretion rate call surfden(t,mtot,mdot) c..total mass of the gas gastot = 0.0d0 do i=1,niso gastot = gastot + g(i + loff) enddo c..the total birthrate bdum = brate(mtot,gastot) c..for each element do nowiso=1,niso j = nowiso + loff c..the birth rate birth = bdum * g(j)/gastot c..the death rate; if we have a valid integral range death = 0.0d0 if (mlo .lt. mmax) then c..from mbmax to mmax of the mass range if (mlo.gt.mbmax) then call qgaus(gasbak,xa3,wa3,nnhi,death) c..else from mbmin to mmax of the mass range else if (mlo.gt.mbmin) then c..with type 1a from mlo to mbmax if (sn1fac .ne. 0.0) then call qgaus(sn1int,xa4,wa4,nnmid,sum) death = death + sum * sn1fac end if c..with nova from mlo to mbmax if (novfac .ne. 0.0) then call qgaus(novaint,xa4,wa4,nnmid,sum) death = death + sum * novfac * trecur end if c..with or without type 1a or nova from mlo to mbmax call qgaus(gasbak,xa4,wa4,nnmid,sum) death = death + (1.0d0 - sn1fac - novfac) * sum c..from mbmax to mmax call qgaus(gasbak,xa2,wa2,nnhi,sum) death = death + sum c..else from mlo to mmax else if (mlo.gt.mmin) then call qgaus(gasbak,xa5,wa5,nnlo,death) c..with type 1a from mbmin to mbmax if (sn1fac .ne. 0.0) then call qgaus(sn1int,xa1,wa1,nnmid,sum) death = death + sum * sn1fac end if c..with nova from mbmin to mbmax if (novfac .ne. 0.0) then call qgaus(novaint,xa1,wa1,nnmid,sum) death = death + sum * novfac * trecur end if c..with or without type 1a or nova from mbmin to mbmax call qgaus(gasbak,xa1,wa1,nnmid,sum) death = death + (1.0d0 - sn1fac - novfac)*sum c..from mbmax to mmax; end of mlo lt mmax if block call qgaus(gasbak,xa2,wa2,nnhi,sum) death = death + sum end if end if c..the odes and the end of the isotope loop dgdt(j) = death - birth + mdot*bigb(nowiso) enddo c..decay k40 if (nk40 .ne. 0 .and. nca40 .ne. 0 .and. nar40 .ne. 0) then if (g(nk40 + loff) .gt. 0.0) then fac = g(nk40+loff)/tmean_k40 dgdt(nk40+loff) = dgdt(nk40+loff) - fac dgdt(nca40+loff) = dgdt(nca40+loff) + brnch1 * fac dgdt(nar40+loff) = dgdt(nar40+loff) + brnch2 * fac end if end if c..if modal26 is larger than zero then decay both al26 and fe60. setting this c..will induce small time steps, on the order of 0.001 gyr. c..decay al26 if (modal26 .gt. 0) then if (nal26 .ne. 0 .and. nmg26 .ne. 0) then if (g(nal26+loff) .gt. 0.0) then fac = g(nal26+loff)/tmean_al26 dgdt(nal26+loff) = dgdt(nal26+loff) - fac dgdt(nmg26+loff) = dgdt(nmg26+loff) + fac end if end if c..decay fe60 if (nfe60 .ne. 0 .and. nni60 .ne. 0) then if (g(nfe60+loff) .gt. 0.0) then fac = g(nfe60+loff)/tmean_fe60 dgdt(nfe60+loff) = dgdt(nfe60+loff) - fac dgdt(nni60+loff) = dgdt(nni60+loff) + fac end if end if end if c..end of zone loop enddo return end subroutine getgas(tdif,gastot,mettot,xiso) implicit none save include 'chem3.dek' c..interpolation/extrapolation routine for the abundance histories. c..input through the pass: c..tdif = retarded time c..input through common: c..nowiso = the isotope under consideration c..output through the pass: c..gastot = total gas mass c..mettot = metallicity mass fraction c..xiso = mass fraction of nowiso c..declare the pass double precision tdif,gastot,mettot,xiso c..local variables integer ihmax parameter (ihmax=nqmax*nqmax + 4*nqmax) integer i,j,ll,loff,ka,iat,okount,ihfind,kb,ke double precision y,thist(ihmax),ghist(ihmax),zhist(ihmax), 1 xdum(gphys),xhist(gphys,ihmax),a,b,denom,dum1 c..initialize counters data ka/0/, okount/-20/, iat/1/, ihfind/0/ c..initialize the output gastot = 0.0d0 mettot = 0.0d0 xiso = 0.0d0 c..if we have had a succesful timestep since the last visit c..to this routine, then reset ka if (okount .ne. kount) then okount = kount ka = 0 end if c..for the first isotope of every zone if (nowiso .eq. 1) then c..find the time index iat to use, if the last one is not good if (tdif .lt. time(iat) .or. tdif .gt. time(iat+1)) then call locate(time,kount,tdif,iat) if (iat .eq. 0) iat = 1 if (iat .eq. kount) iat = kount - 1 end if c..only do the work if the time point is different followed by a c..hard wired two point linear fit, y is the interpolated value ka = ka + 1 if (ka .gt. ihmax) stop 'ka > ihmax in routine getgas' if (ka.eq.1 .or. thist(max(1,ka-1)) .ne. tdif) then denom = time(iat) - time(iat + 1) if (denom .eq. 0) stop 'denominator 0 in routine getgas' denom = 1.0d0/denom loff = (nowzon-1)*niso kb = 1 + loff ke = niso + loff do j=kb,ke dum1 = grk(j,iat+1) b = (grk(j,iat) - dum1) * denom a = dum1 - b * time(iat + 1) y = b*tdif + a gastot = gastot + y xdum(j)= y enddo kb = nzbeg + loff ke = nzend + loff do j=kb,ke mettot = mettot + xdum(j) enddo xiso = xdum(nowiso + loff) c..store tdif, gastot and mettot khist = max(ka,khist) thist(ka) = tdif ghist(ka) = gastot zhist(ka) = mettot kb = 1 + loff ke = niso + loff do j=kb,ke xhist(j,ka) = xdum(j) enddo c..otherwise either ka is not one and the time point is the same else ka = ka - 1 end if c..else this is not the first isotope of a zone else c..this history point i should be a good guess; c..khyes counts the direct hits i = min(ihfind+1,ka) if (tdif .eq. thist(i)) then khyes = khyes + 1 ihfind = ihfind + 1 if (i .eq. ka) ihfind = 0 gastot = ghist(i) mettot = zhist(i) j = nowiso + (nowzon-1)*niso xiso = xhist(j,i) if (nowiso .eq. niso .and. i .eq. ka) ka = 0 goto 100 c..not a good guess, find the stored history; c..khno counts the bad table hits else khno = khno + 1 do i=1,ka if (tdif .eq. thist(i)) then ihfind = i gastot = ghist(i) mettot = zhist(i) j = nowiso + (nowzon-1)*niso xiso = xhist(j,i) if (nowiso .eq. niso .and. i .eq. ka) ka = 0 goto 100 end if enddo end if c..we should never enter this section c..kmiss counts the number of times we do enter this section c..we get here if this is not the first isotope in a zone c..and the history point was not found in storage. do it straight. kmiss = kmiss + 1 call locate(time,kount,tdif,iat) if (iat .eq. 0) iat = 1 if (iat .eq. kount) iat = kount - 1 denom = time(iat) - time(iat+1) if (denom .eq. 0) stop 'denominator 0 in routine getgas' denom = 1.0d0/denom loff = (nowzon-1)*niso kb = 1 + loff ke = niso + loff do j=kb,ke dum1 = grk(j,iat+1) b = (grk(j,iat) - dum1) * denom a = dum1 - b * time(iat+1) y = b*tdif + a gastot = gastot + y xdum(j) = y enddo kb = nzbeg + loff ke = nzend + loff do j=kb,ke mettot = mettot + xdum(j) enddo xiso = xdum(nowiso + loff) end if c..metal and isotope mass fraction 100 mettot = mettot/gastot xiso = xiso/gastot return end subroutine getgas2(tdif,gastot) implicit none save include 'chem3.dek' c..this routine is only used for the type2 supernova, c..type1 supernova, and classical nova rate calculations. c..it is much like the last portion of routine getgas, but simpler. c..input: c..tdif = the retarted time c..output c..gastot = total gas mass at tdif c..declare the pass double precision tdif,gastot c..local variables integer j,kb,ke,iat,loff double precision y,a,b,denom,dum1 data iat/1/ c..initialize gastot = 0.0d0 c..find the location in the stored time history call locate(time,kount,tdif,iat) if (iat .eq. 0) iat = 1 if (iat .eq. kount) iat = kount - 1 denom = time(iat) - time(iat+1) if (denom .eq. 0) stop 'denominator 0 in routine getgas2' denom = 1.0d0/denom loff = (nowzon-1)*niso c..sum the mass in each isotope kb = 1 + loff ke = niso + loff do j=kb,ke dum1 = grk(j,iat+1) b = (grk(j,iat) - dum1) * denom a = dum1 - b * time(iat+1) y = b*tdif + a gastot = gastot + y enddo return end double precision function gasbak(m) implicit none save include 'chem3.dek' c..the integro part of the problem. c..this routine sets the integrand and should be called c..from an integration routine. c..declare the pass double precision m c..local variables double precision taum,tdif,y,imf,mtot,mdot,brate,gastot, 1 mettot,xiso c..retarded time tdif = timcom - taum(m) c..retarded total mass and accretion rate call surfden(tdif,mtot,mdot) c..retarded total gas mass, metallicity and mass fraction call getgas(tdif,gastot,mettot,xiso) c..get the yield of this element at this metallicity call getyeld2(mettot,m,y) c..integrand gasbak = brate(mtot,gastot) * imf(m)/ximf_norm * y/m return end double precision function type2rn(m) implicit none save include 'chem3.dek' c..this routine is only used for the type2 supernova rate by number c..calculations. it is much like routine gasbak, but simpler. c..declare the pass double precision m c..local variables double precision taum,tdif,imf,mtot,mdot,brate,gastot c..the retarted time tdif = timcom - taum(m) c..total mass at the retarted time call surfden(tdif,mtot,mdot) c..total gas mass at the retarted time call getgas2(tdif,gastot) c..the integrand type2rn = brate(mtot,gastot)/m * imf(m)/ximf_norm return end double precision function sn1int(m) implicit none save include 'chem3.dek' c..this routine computes the double integral for the c..type 1a nucleosynthesis contributions. input through c..the pass is the mass of the binary system. input c..through common is the isotope under consideration nowiso. c..output is the mass of isotope nowiso ejected by type1a events. c..declare the pass double precision m c..local variables external gasbak2 integer jj,jjmax,i,ihfind parameter (jjmax=4*nqmax) double precision m2,mtau,mulo,sum,xxx2(nqmax),sum2,half, 1 www2(nqmax),tiny,muhist(jjmax),sumhist(jjmax), 2 gasbak2 parameter (half = 0.5d0, tiny = 1.0e-20) data jj/0/, ihfind/0/ c..initialize c..the time corresponds to the lifetime of the secondary c..mulo is the lower limit of integration sn1int = 0.0d0 mbin = m m2 = mtau(timcom) mulo = max(m2/mbin,(mbin-half*mbmax)/mbin) + tiny if (mulo .ge. half) return c..if this is the first isotope integrate and store the histories if (nowiso .eq. 1) then jj = jj + 1 if (jj .gt. jjmax) stop 'jj > jjmax in routine sn1int' call gauleg2(mulo,half,xxx2,www2,nnhi) call qgaus2(gasbak2,xxx2,www2,nnhi,sum) muhist(jj) = mulo sumhist(jj) = sum c..else this is not the first isotope else c..this should be a pretty good direct table hit c..kt1yes counts the number of table hits i = min(ihfind+1,jj) if (mulo .eq. muhist(i)) then kt1yes = kt1yes + 1 sum = sumhist(i) ihfind = ihfind + 1 if (i .eq. jj) ihfind = 0 if (nowiso .eq. niso .and. i .eq. jj) jj = 0 goto 30 c..not a direct table hit. do a search. c..kt1no counts the number of table misses else kt1no = kt1no + 1 do i=1,jj if (mulo .eq. muhist(i)) then sum = sumhist(i) ihfind = i if (nowiso .eq. niso .and. i .eq. jj) jj = 0 goto 30 end if enddo end if c..we should never enter this section, but if we do c..do the integration and mark how many times we eneter this section kt1miss = kt1miss + 1 call gauleg2(mulo,half,xxx2,www2,nnhi) call qgaus2(gasbak2,xxx2,www2,nnhi,sum) end if c..sum over the different types of type1a supernova 30 continue sum2 = 0.0d0 do i=1,nt1a sum2 = sum2 + fract1a(i) * mt1a(i,nowiso) enddo sn1int = sum * sum2 return end double precision function gasbak2(phi) implicit none save include 'chem3.dek' c..this routine is the integrand for the type 1a nucleosynthesis c..contributions. this routine sets the integrand and should c..be called from an integration routine. c..input is mass fraction of the binary's secondary c..declare the pass double precision phi c..local variables double precision taum,tdif,imf,mtot,mdot,brate,gastot, 1 mettot,xiso c..the retarded time tdif = timcom - taum(phi*mbin) c..the total mass at the retarded time call surfden(tdif,mtot,mdot) c..the total gas mass at the retarded time call getgas(tdif,gastot,mettot,xiso) c..integrand; final factor is the greggio and renzini distribution c..function phi**gamma with gamma = 2 and normalized over 0 to one half c.. c.. gasbak2 = brate(mtot,gastot) * imf(mbin)/ximf_norm * c.. 1 mt1a(nowiso)/(phi*mbin) * 24.0*phi*phi c..it more efficient to pull the type1a yields out of the integrand c..and cancel a phi gasbak2 = brate(mtot,gastot) 1 * imf(mbin)/ximf_norm 2 * 24.0d0*phi/mbin return end double precision function sn1rint(m) implicit none save include 'chem3.dek' c..this routine is only for the typ1a supernova rate calculations. c..input is the mass of the binary. c..declare the pass double precision m c..local variables external type1arn double precision m2,mtau,mulo,sum,xxx2(nqmax),www2(nqmax), 1 half,tiny parameter (half = 0.5d0, tiny = 1.0e-20) c..initialize c..the time corresponds to the lifetime of the secondary c..mulo is the lower limit of integration sn1rint = 0.0d0 mbin = m m2 = mtau(timcom) mulo = max(m2/mbin,(mbin-half*mbmax)/mbin) + tiny if (mulo .ge. half) return c..integrate call gauleg2(mulo,half,xxx2,www2,nnhi) call qgaus2(type1arn,xxx2,www2,nnhi,sum) sn1rint = sum return end double precision function type1arn(phi) implicit none save include 'chem3.dek' c..this routine is only used for the type1a supernova rate calculations, c..and sets the integrand for routine sn1rint. this routine is much c..like routine gasbak2, but simpler. input through the pass is phi, c..the mass fraction of the binary's secondary. c..declare the pass double precision phi c..local variables double precision taum,tdif,imf,mtot,mdot,brate,gastot c..the retarded time tdif = timcom - taum(phi*mbin) c..the total mass at the retarded time call surfden(tdif,mtot,mdot) c..the gas mass at the retarded time call getgas2(tdif,gastot) c..the integrand type1arn = brate(mtot,gastot) 1 * imf(mbin)/ximf_norm 2 * 24.0d0*phi/mbin return end double precision function novaint(m) implicit none save include 'chem3.dek' c..this routine computes the double integral for the c..nova nucleosynthesis contributions. input through c..the pass is the mass of the binary system. input c..through common is the isotope under consideration nowiso. c..output is the mass of isotope nowiso ejected by nova events. c..declare external gasbak2 integer jj,ll,jjmax,i,ihfind parameter (jjmax=50) double precision m,m2,mtau,mulo,sum,xxx2(nqmax),www2(nqmax), 1 half,tiny,muhist(jjmax),sumhist(jjmax),sum2, 2 gasbak2 parameter (half = 0.5d0, tiny = 1.0e-20) data jj/0/, ihfind/0/ c..initialize c..the time corresponds to the lifetime of the secondary c..mulo is the lower limit of integration novaint = 0.0d0 mbin = m m2 = mtau(timcom) mulo = max(m2/mbin,(mbin-half*mbmax)/mbin) + tiny if (mulo .ge. half) return c..if this is the first isotope integrate and store the histories if (nowiso .eq. 1) then jj = jj + 1 if (jj .gt. jjmax) stop 'jj > jjmax in routine novaint' call gauleg2(mulo,half,xxx2,www2,nnhi) call qgaus2(gasbak2,xxx2,www2,nnhi,sum) muhist(jj) = mulo sumhist(jj) = sum c..else this is not the first isotope else c..this should be a pretty good direct table hit c..knvyes counts the number of table hits i = min(ihfind+1,jj) if (mulo .eq. muhist(i)) then knvyes = knvyes + 1 sum = sumhist(i) ihfind = ihfind + 1 if (i .eq. jj) ihfind = 0 if (nowiso .eq. niso .and. i .eq. jj) jj = 0 goto 30 c..not a direct table hit. do a search. c..kt1no counts the number of table misses else knvno = knvno + 1 do i=1,jj if (mulo .eq. muhist(i)) then sum = sumhist(i) ihfind = i if (nowiso .eq. niso .and. i .eq. jj) jj = 0 goto 30 end if enddo end if c..we should never enter this section, but if we do c..do the integration and mark how many times we eneter this section knvmiss = knvmiss + 1 call gauleg2(mulo,half,xxx2,www2,nnhi) call qgaus2(gasbak2,xxx2,www2,nnhi,sum) end if c..sum over the different types of nova 30 continue sum2 = 0.0d0 do i=1,nnova sum2 = sum2 + fracnova(i) * mnova(i,nowiso) enddo novaint = sum * sum2 return end double precision function novrint(m) implicit none save include 'chem3.dek' c..this routine is only for the nova rate calculations. c..input is the mass of the binary. c..declare the pass double precision m c..local variables external novarn double precision m2,mtau,mulo,sum,xxx2(nqmax),www2(nqmax), 1 half,tiny parameter (half = 0.5d0, tiny = 1.0e-20) c..initialize c..the time corresponds to the lifetime of the secondary c..mulo is the lower limit of integration novrint = 0.0d0 mbin = m m2 = mtau(timcom) mulo = max(m2/mbin,(mbin-half*mbmax)/mbin) + tiny if (mulo .ge. half) return c..integrate call gauleg2(mulo,half,xxx2,www2,nnhi) call qgaus2(novarn,xxx2,www2,nnhi,sum) novrint = sum return end double precision function novarn(phi) implicit none save include 'chem3.dek' c..this routine is only used for the nova rate calculations, c..and sets the integrand for routine novrint. this routine is much c..like routine gasbak2, but simpler. input through the pass is phi, c..the mass fraction of the binary's secondary. c..declare the pass double precision phi c..local variables double precision taum,tdif,imf,mtot,mdot,brate,gastot c..the retarded time tdif = timcom - taum(phi*mbin) c..total mass at the retarded time call surfden(tdif,mtot,mdot) c..total gas at the retarded time call getgas2(tdif,gastot) c..the integrand novarn = brate(mtot,gastot) 1 * imf(mbin)/ximf_norm 2 * 24.0d0*phi/mbin return end subroutine surfden(t,mtot,mdot) implicit none save include 'chem3.dek' c..this routine computes the total mass and accretion rate, c..and sets the units for the whole calculation. c..input: c..t = time in gyr c..output: c..mtot = total mass in solar masses per square parsec c..mdot = accretion rate in solar masses per square parsec per gigayear c..profile for 0 < r < rbreak kpc is an inverse square c.. mtot = anucl/(r + rnucl)**2 c.. c..profile for r > rbreak kpc is an exponential decay c.. mtot = adisk * exp(-r/rdisk) c.. c..all modulo c.. mtot = mtot * kdisk * exp(-time/tdisk) c..declare the pass double precision t,mtot,mdot c..local variables external rgfunc integer ifirst,i,niter double precision rdisk,adisk,rnucl,anucl,xx, 1 lo,hi,rgfunc,zbrent,fact,acrete,trdisk, 2 tol,rtab1(nzmax),rtab2(nzmax),radnow parameter (tol=1.0d-6) data ifirst/0/ c..initialize radial parameters if (ifirst .eq. 0) then ifirst = 1 lo = 0.0d0 hi = 25.0e3 rnucl = zbrent(rgfunc,lo,hi,tol,niter) anucl = mcenter * rnucl * rnucl rdisk = 0.5d0 * (rbreak + rnucl) adisk = massol * exp(radsol/rdisk) c write(6,01) 'anucl=',anucl,'rnucl=',rnucl, c 1 'adisk=',adisk,'rdisk=',rdisk c01 format(1x,4(a,1pe11.3,' ')) c..precompute most of the exponentials do i=1,nzone rtab1(i) = exp(-rpoint(i)/rdisk) trdisk = atdisk * 1.0d-3 * rpoint(i) + btdisk rtab2(i) = exp(-tgal/trdisk) enddo end if c..here is the spatial dependence radnow = rpoint(nowzon) if (radnow .le. rbreak) then xx = radnow + rnucl mtot = anucl/(xx*xx) else mtot = adisk * rtab1(nowzon) end if c..here is the time dependence trdisk = atdisk * 1.0d-3 * radnow + btdisk acrete = (mtot-sigzer) / (trdisk*(1.0d0 - rtab2(nowzon))) fact = exp(-t/trdisk) c..the total mass and the accretion rate mtot = sigzer + acrete * trdisk * (1.0d0 - fact) mdot = acrete * fact return end double precision function rgfunc(rn) implicit none save include 'chem3.dek' c..root find routine used by routine surfden c.. c..declare and go double precision rn,xx xx = rbreak + rn rgfunc = massol * exp(2.0d0*(radsol - rbreak)/xx) 1 - mcenter * rn * rn/(xx*xx) return end subroutine modout implicit none save include 'chem3.dek' c..this routine outputs the results c..declare integer i,j,nz,loff double precision gas,zsum c..get the normalizing solar and ism abundances call get_solar c..for each radial point c..compute and store the gas and total metals do nz=1,nzone nowzon = nz loff = (nz - 1)*niso do i=1,kount gas = 0.0d0 do j=1,niso gas = gas + grk(j + loff,i) enddo gashist(i) = gas zsum = 0.0d0 do j=nzbeg,nzend zsum = zsum + grk(j + loff,i) enddo zsumhist(i) = zsum enddo c..write out lots of files call isotope_histories(nz) call ratio_histories(nz,1) call ratio_histories(nz,2) call supernova_histories(nz) call kinematics(nz) call gdwarf_distribution(nz) call solar_distribution(nz) c..go back for another radial zone enddo c..now do the radial integrals, if appropriate. if (nzone .lt. 3) return call radial_sums return end subroutine get_solar implicit none save include 'chem3.dek' c..local variables integer i double precision xsol,xelem,z,a,andgrev,gas,zsumf,xnow,tiny parameter (tiny = 1.0d-100) c..for easy zeroing of the solar and mean ism elements integer nelem parameter (nelem = 30) double precision elemsol(nelem) equivalence (elemsol(1),hsol) double precision elemism(nelem) equivalence (elemism(1),hism) c..zero the solar and mean ism elemental quantities do i=1,nelem elemsol(i) = tiny elemism(i) = tiny enddo c..set tagout to be at least as long as the final time integration point if (tagout .gt. time(kount)) tagout = 0.9d0 * time(kount) c..get the solar and mean ism elemental mass fractions nowzon = iznorm do i=1,niso if (i .eq. nal26) then xsol = 1.0d0 xelem = 1.0d0 z = 13.0 a = 26.0 else if (i .eq. nfe60) then xsol = 1.0d0 xelem = 1.0d0 z = 26.0 a = 60.0 else xsol = andgrev(namnuc(i),z,a,xelem) end if solnorm(i) = xsol zsolnorm(i) = z asolnorm(i) = a nowiso = i call getgas(tagout,gas,zsumf,xnow) ismnorm(i) = xnow c..sum the isotopic contributions to an element if (i .eq. 1 .or. i .eq. nh2) then hsol = xelem hism = hism + xnow else if (i .eq. nhe3 .or. i .eq. nhe4) then hesol = xelem heism = heism + xnow else if (i .eq. nli6 .or. i .eq. nli7) then lisol = xelem liism = liism + xnow else if (i .eq. nbe9) then besol = xelem beism = xnow else if (i .eq. nc12 .or. i .eq. nc13) then csol = xelem cism = cism + xnow else if (i .eq. nn14 .or. i .eq. nn15) then nsol = xelem nism = nism + xnow else if (i .eq. no16 .or. i .eq. no17 .or. i .eq. no18) then osol = xelem oism = oism + xnow else if (i .eq. nf19) then fsol = xelem fism = xnow else if (i .eq. nne20 .or. i .eq. nne21 .or. i .eq. nne22) then nesol = xelem neism = neism + xnow else if (i .eq. nna23) then nasol = xelem naism = xnow else if (i .eq. nal27) then alsol = xelem alism = xnow else if (i .eq. nmg24 .or. i .eq. nmg25 .or. i .eq. nmg26) then mgsol = xelem mgism = mgism + xnow else if (i .eq. nsi28 .or. i .eq. nsi29 .or. i .eq. nsi30) then sisol = xelem siism = siism + xnow else if (i .eq. np31) then psol = xelem pism = xnow else if (i .eq. ns32 .or. i .eq. ns33 .or. 1 i .eq. ns34 .or. i .eq. ns36) then susol = xelem suism = suism + xnow else if (i .eq. ncl35 .or. i .eq. ncl37) then clsol = xelem clism = clism + xnow else if (i .eq. nar36 .or. i .eq. nar38 .or. i .eq. nar40) then arsol = xelem arism = arism + xnow else if (i .eq. nk39 .or. i .eq. nk40 .or. i .eq. nk41) then ksol = xelem kism = kism + xnow else if (i .eq. nca40 .or. i .eq. nca42 .or. 1 i .eq. nca43 .or. i .eq. nca44 .or. 2 i .eq. nca46 .or. i .eq. nca48) then casol = xelem caism = caism + xnow else if (i .eq. nsc45) then scsol = xelem scism = xnow else if (i .eq. nti46 .or. i .eq. nti47 .or. 1 i .eq. nti48 .or. i .eq. nti49 .or. 2 i .eq. nti50 ) then tisol = xelem tiism = tiism + xnow else if (i .eq. nv50 .or. i .eq. nv51) then vsol = xelem vism = vism + xnow else if (i .eq. ncr50 .or. i .eq. ncr52 .or. 1 i .eq. ncr53 .or. i .eq. ncr54) then crsol = xelem crism = crism + xnow else if (i .eq. nmn55) then mnsol = xelem mnism = xnow else if (i .eq. nfe54 .or. i .eq. nfe56 .or. 1 i .eq. nfe57 .or. i .eq. nfe58) then fesol = xelem feism = feism + xnow else if (i .eq. nco59) then cosol = xelem coism = xnow else if (i .eq. nni58 .or. i .eq. nni60 .or. 1 i .eq. nni61 .or. i .eq. nni62 .or. 2 i .eq. nni64 ) then nisol = xelem niism = niism + xnow else if (i .eq. ncu63 .or. i .eq. ncu65) then cusol = xelem cuism = cuism + xnow else if (i .eq. nzn64 .or. i .eq. nzn66 .or. 1 i .eq. nzn67 .or. i .eq. nzn68 .or. 2 i .eq. nzn70 ) then znsol = xelem znism = znism + xnow end if c..back for another isotope enddo c..compute the metallicity of the gas zsol = 1.0d0 - hsol - hesol zism = 1.0d0 - hism - heism return end subroutine isotope_histories(nz) implicit none save include 'lunits.dek' include 'parze.dek' include 'chem3.dek' c..writes out the isotopic histories for zone nz c..declare the pass integer nz c..local variables integer i,k,lop,jrem,ilop,kb,ke,loff double precision gas,zsum,zsumf,zln,hsum,fesum,x,z,tiny parameter (tiny = 1.0d-100) c..popular format statements 01 format(' *',/, 1 ' ',t5,a,t17,a,t29,a,t41,a,t53,a,t65,a, 2 t77,a,t89,a,t101,a,t113,a,t125,a,t137,a) 02 format(1x,1p14e12.4) 09 format(a20,'_z',i4,'_',i4,a) c..set the zone offset loff = (nz - 1)*niso c..write out the solution in blocks of 8 isotopes c..lop is how many groups of 8 exist; jrem is the remainder lop = int(niso/8) jrem = niso - 8*lop do ilop = 1,lop+1 kb = 1 + 8*(ilop-1) ke = 8 + 8*(ilop-1) if (ilop .eq. lop+1 .and. jrem .eq. 0) goto 33 if (ilop .eq. lop+1) ke = niso c..raw solar masses per square parsec file write(string,09) xmasfil,nz,ilop,'.raw' call sqeeze(string) call bbb(redwrt,fil2,string,err) call header(fil2) write(fil2,01) 'time','ztot',(namnuc(k), k=kb,ke) c..mass fraction file write(string,09) xmasfil,nz,ilop,'.dat' call sqeeze(string) call bbb(redwrt,fil3,string,err) call header(fil3) write(fil3,01) 'time','z/gas','[fe/h]',(namnuc(k), k=kb,ke) c..mass fraction relative to solar write(string,09) xmasfil,nz,ilop,'.sol' call sqeeze(string) call bbb(redwrt,fil8,string,err) call header(fil8) write(fil8,01) 'time','z/zsol','[fe/h]',(namnuc(k), k=kb,ke) c..mass fraction relative to mean ism at solar birth write(string,09) xmasfil,nz,ilop,'.ism' call sqeeze(string) call bbb(redwrt,fil9,string,err) call header(fil9) write(fil9,01) 'time','z/zism','[fe/h]',(namnuc(k), k=kb,ke) c..raw solar masses per square parsec per gigayear file write(string,09) xmasfil(1:19)//'i',nz,ilop,'.raw' call sqeeze(string) call bbb(redwrt,fil11,string,err) call header(fil11) write(fil11,01) 'time','z/gas',(namnuc(k), k=kb,ke) c..for every time step do i=1,kount gas = gashist(i) zsum = zsumhist(i) zsumf = zsum/gas zln = log10(max(tiny,zsumf/zsol)) c..total hydrogen hsum = grk(1 + loff,i)/gas c..[fe/h] fesum = tiny if (nfe54 .ne. 0) fesum = fesum + grk(nfe54 + loff,i) if (nfe56 .ne. 0) fesum = fesum + grk(nfe56 + loff,i) if (nfe57 .ne. 0) fesum = fesum + grk(nfe57 + loff,i) if (nfe58 .ne. 0) fesum = fesum + grk(nfe58 + loff,i) fesum = fesum/gas c..relative to solar x = log10(max(tiny,fesum/fesol * hsol/hsum)) c..relative to mean ism z = log10(max(tiny,fesum/feism * hism/hsum)) c..solution vectors if (zln .gt. -5.0) then write(fil2,02) time(i),zsum,(grk(k+loff,i),k=kb,ke) write(fil3,02) time(i),zsumf,z,(grk(k+loff,i)/gas,k=kb,ke) write(fil8,02) time(i),zsumf/zsol,x, 1 (grk(k+loff,i)/gas/solnorm(k),k=kb,ke) write(fil9,02) time(i),zsumf/zism,z, 1 (grk(k+loff,i)/gas/ismnorm(k),k=kb,ke) write(fil11,02) time(i),zsumf,(dgrk(k+loff,i),k=kb,ke) end if c..end of time step loop enddo c..close all the files and c..go back for another block of 8 isotopes call bbb(clse,fil2,string,k) call bbb(clse,fil3,string,k) call bbb(clse,fil8,string,k) call bbb(clse,fil9,string,k) call bbb(clse,fil11,string,k) 33 continue enddo return end subroutine ratio_histories(nz,itype) implicit none save include 'lunits.dek' include 'parze.dek' include 'chem3.dek' c..writes out the ratio histories for zone nz c..declare the pass integer nz,itype c..local variables character*4 ctype integer i,j,k,loff double precision gas,zsum,zsumf,zln,hsum,fesum,aa,dd,tiny parameter (tiny = 1.0d-100) c..element ratio storage integer neratio parameter (neratio=35) character*7 ceratio(neratio) double precision eratio(neratio) equivalence (eratio(1),h2h) double precision fe_norm, h_norm, c_norm, n_norm, o_norm, f_norm, 1 ne_norm, na_norm, al_norm, mg_norm, si_norm, 2 p_norm, su_norm, cl_norm, ar_norm, k_norm, 3 ca_norm, sc_norm, ti_norm, v_norm, cr_norm, 4 mn_norm, co_norm, ni_norm, cu_norm, zn_norm c..initialize the element ratio names data ceratio/'d/h ','he3/h ','n(li) ','n(be) ','n(b) ', 1 '[fe/h] ','[c/fe] ','[n/fe] ','[o/fe] ','[f/o] ', 2 '[ne/fe]','[na/fe]','[al/fe]','[mg/fe]','[na/mg]', 3 '[al/mg]','[si/fe]','[p/fe] ','[s/fe] ','[cl/fe]', 4 '[ar/fe]','[k/fe] ','[ca/fe]','[sc/fe]','[ti/fe]', 5 '[v/fe] ','[cr/fe]','[mn/fe]','[co/fe]','[ni/fe]', 6 '[cu/fe]','[zn/fe]','delta y','delydlz','dydz '/ c..popular format statements 01 format(' *',/, 1 ' ',t5,a,t17,a,t29,a,t41,a,t53,a,t65,a, 2 t77,a,t89,a,t101,a,t113,a,t125,a,t137,a) 02 format(1x,1p14e12.4) 09 format(a20,'_z',i4,'_',i4,a) c..set the zone offset loff = (nz - 1)*niso if (itype .eq. 1) then ctype = '.sol' fe_norm = fesol h_norm = hsol c_norm = csol n_norm = nsol o_norm = osol f_norm = fsol ne_norm = nesol na_norm = nasol al_norm = alsol mg_norm = mgsol si_norm = sisol p_norm = psol su_norm = susol cl_norm = clsol ar_norm = arsol k_norm = ksol ca_norm = casol sc_norm = scsol ti_norm = tisol v_norm = vsol cr_norm = crsol mn_norm = mnsol co_norm = cosol ni_norm = nisol cu_norm = cusol zn_norm = znsol else if (itype .eq. 2) then ctype = '.ism' fe_norm = feism h_norm = hism c_norm = cism n_norm = nism o_norm = oism f_norm = fism ne_norm = neism na_norm = naism al_norm = alism mg_norm = mgism si_norm = siism p_norm = pism su_norm = suism cl_norm = clism ar_norm = arism k_norm = kism ca_norm = caism sc_norm = scism ti_norm = tiism v_norm = vism cr_norm = crism mn_norm = mnism co_norm = coism ni_norm = niism cu_norm = cuism zn_norm = znism endif c..open the element ratio relative to solar or ism files write(string,09) decfil,nz,1,ctype call sqeeze(string) call bbb(redwrt,fil2,string,i) call header(fil2) write(fil2,01) 'time','z','[fe/h]',(ceratio(k),k=1,7) write(string,09) decfil,nz,2,ctype call sqeeze(string) call bbb(redwrt,fil3,string,i) call header(fil3) write(fil3,01) 'time','z','[fe/h]',(ceratio(k),k=8,14) write(string,09) decfil,nz,3,ctype call sqeeze(string) call bbb(redwrt,fil4,string,i) call header(fil4) write(fil4,01) 'time','z','[fe/h]',(ceratio(k),k=15,21) write(string,09) decfil,nz,4,ctype call sqeeze(string) call bbb(redwrt,fil8,string,i) call header(fil8) write(fil8,01) 'time','z','[fe/h]',(ceratio(k),k=22,28) write(string,09) decfil,nz,5,ctype call sqeeze(string) call bbb(redwrt,fil9,string,i) call header(fil9) write(fil9,01) 'time','z','[fe/h]',(ceratio(k),k=29,neratio) c..for every time step do i=1,kount gas = gashist(i) zsum = zsumhist(i) zsumf = zsum/gas zln = log10(max(tiny,zsumf/zsol)) c..initial element ratios do j=1,neratio eratio(j) = tiny enddo c..hydrogen mass fraction hsum = grk(1 + loff,i)/gas c..d/h by number if (nh2 .ne. 0) h2h = h2h + grk(nh2 + loff,i) h2h = h2h/(2.0d0 * gas * hsum) c..he3/h by number if (nhe3 .ne. 0) he3h = he3h + grk(nhe3 + loff,i) he3h = he3h/(3.0d0 * gas * hsum) c..delta y, delta y / delta z and dydz deltay = grk(nhe4+loff,i)/gas - primhe4 delydelz = deltay/zsumf dydz = -100.0d0 c..log10 n(li7); abundance by number; h=12 if (nli7 .ne. 0) li72h = li72h + grk(nli7 + loff,i) li72h = li72h/(7.0d0 * gas) li72h = log10(max(tiny,li72h/hsum)) + 12.0d0 c..log10 n(be); abundance by number; h=12 if (nbe9 .ne. 0) be2h = be2h + grk(nbe9 + loff,i) be2h = be2h/(9.0d0 * gas) be2h = log10(max(tiny,be2h/hsum)) + 12.0d0 c..log10 n(b); abundance by number; h=12 if (nb10 .ne. 0) b2h = b2h + grk(nb10 + loff,i) if (nb11 .ne. 0) b2h = b2h + grk(nb11 + loff,i) b2h = b2h/(21.0d0 * gas) b2h = log10(max(tiny,b2h/hsum)) + 12.0d0 c..[fe/h] fesum = tiny if (nfe54 .ne. 0) fesum = fesum + grk(nfe54 + loff,i) if (nfe56 .ne. 0) fesum = fesum + grk(nfe56 + loff,i) if (nfe57 .ne. 0) fesum = fesum + grk(nfe57 + loff,i) if (nfe58 .ne. 0) fesum = fesum + grk(nfe58 + loff,i) fesum = fesum/gas fe2h = log10(max(tiny,fesum/fe_norm * h_norm/hsum)) dd = fe_norm/fesum c..[c/fe] if (nc12 .ne. 0) c2fe = c2fe + grk(nc12 + loff,i) if (nc13 .ne. 0) c2fe = c2fe + grk(nc13 + loff,i) c2fe = c2fe/gas c2fe = log10(max(tiny,c2fe/c_norm * dd)) c..[n/fe] if (nn14 .ne. 0) n2fe = n2fe + grk(nn14 + loff,i) if (nn15 .ne. 0) n2fe = n2fe + grk(nn15 + loff,i) n2fe = n2fe/gas n2fe = log10(max(tiny,n2fe/n_norm * dd)) c..[o/fe] if (no16 .ne. 0) o2fe = o2fe + grk(no16 + loff,i) if (no17 .ne. 0) o2fe = o2fe + grk(no17 + loff,i) if (no18 .ne. 0) o2fe = o2fe + grk(no18 + loff,i) o2fe = o2fe/gas o2fe = log10(max(tiny,o2fe/o_norm * dd)) c..[f/o] if (nf19 .ne. 0) f2o = f2o + grk(nf19 + loff,i) f2o = f2o/gas f2o = log10(max(tiny,f2o/f_norm * dd)) f2o = f2o - o2fe c..[ne/fe] if (nne20 .ne. 0) ne2fe = ne2fe + grk(nne20 + loff,i) if (nne21 .ne. 0) ne2fe = ne2fe + grk(nne21 + loff,i) if (nne22 .ne. 0) ne2fe = ne2fe + grk(nne22 + loff,i) ne2fe = ne2fe/gas ne2fe = log10(max(tiny,ne2fe/ne_norm * dd)) c..[na/fe] if (nna23 .ne. 0) na2fe = na2fe + grk(nna23 + loff,i) na2fe = na2fe/gas na2fe = log10(max(tiny,na2fe/na_norm * dd)) c..[al/fe] if (nal27 .ne. 0) al2fe = al2fe + grk(nal27 + loff,i) al2fe = al2fe/gas al2fe = log10(max(tiny,al2fe/al_norm * dd)) c..[mg/fe] if (nmg24 .ne. 0) mg2fe = mg2fe + grk(nmg24 + loff,i) if (nmg25 .ne. 0) mg2fe = mg2fe + grk(nmg25 + loff,i) if (nmg26 .ne. 0) mg2fe = mg2fe + grk(nmg26 + loff,i) mg2fe = mg2fe/gas mg2fe = log10(max(tiny,mg2fe/mg_norm * dd)) c..[na/mg],[al/mg] na2mg = na2fe - mg2fe al2mg = al2fe - mg2fe c..[si/fe] if (nsi28 .ne. 0) si2fe = si2fe + grk(nsi28 + loff,i) if (nsi29 .ne. 0) si2fe = si2fe + grk(nsi29 + loff,i) if (nsi30 .ne. 0) si2fe = si2fe + grk(nsi30 + loff,i) si2fe = si2fe/gas si2fe = log10(max(tiny,si2fe/si_norm * dd)) c..[p/fe] if (np31 .ne. 0) p2fe = p2fe + grk(np31 + loff,i) p2fe = p2fe/gas p2fe = log10(max(tiny,p2fe/p_norm * dd)) c..[s/fe] if (ns32 .ne. 0) s2fe = s2fe + grk(ns32 + loff,i) if (ns33 .ne. 0) s2fe = s2fe + grk(ns33 + loff,i) if (ns34 .ne. 0) s2fe = s2fe + grk(ns34 + loff,i) if (ns36 .ne. 0) s2fe = s2fe + grk(ns36 + loff,i) s2fe = s2fe/gas s2fe = log10(max(tiny,s2fe/su_norm * dd)) c..[cl/fe] if (ncl35 .ne. 0) cl2fe = cl2fe + grk(ncl35 + loff,i) if (ncl37 .ne. 0) cl2fe = cl2fe + grk(ncl37 + loff,i) cl2fe = cl2fe/gas cl2fe = log10(max(tiny,cl2fe/cl_norm * dd)) c..[ar/fe] if (nar36 .ne. 0) ar2fe = ar2fe + grk(nar36 + loff,i) if (nar38 .ne. 0) ar2fe = ar2fe + grk(nar38 + loff,i) if (nar40 .ne. 0) ar2fe = ar2fe + grk(nar40 + loff,i) ar2fe = ar2fe/gas ar2fe = log10(max(tiny,ar2fe/ar_norm * dd)) c..[k/fe] if (nk39 .ne. 0) k2fe = k2fe + grk(nk39 + loff,i) if (nk40 .ne. 0) k2fe = k2fe + grk(nk40 + loff,i) if (nk41 .ne. 0) k2fe = k2fe + grk(nk41 + loff,i) k2fe = k2fe/gas k2fe = log10(max(tiny,k2fe/k_norm * dd)) c..[ca/fe] if (nca40 .ne. 0) ca2fe = ca2fe + grk(nca40 + loff,i) if (nca42 .ne. 0) ca2fe = ca2fe + grk(nca42 + loff,i) if (nca43 .ne. 0) ca2fe = ca2fe + grk(nca43 + loff,i) if (nca44 .ne. 0) ca2fe = ca2fe + grk(nca44 + loff,i) if (nca46 .ne. 0) ca2fe = ca2fe + grk(nca46 + loff,i) if (nca48 .ne. 0) ca2fe = ca2fe + grk(nca48 + loff,i) ca2fe = ca2fe/gas ca2fe = log10(max(tiny,ca2fe/ca_norm * dd)) c..[sc/fe] if (nsc45 .ne. 0) sc2fe = sc2fe + grk(nsc45 + loff,i) sc2fe = sc2fe/gas sc2fe = log10(max(tiny,sc2fe/sc_norm * dd)) c..[ti/fe] if (nti46 .ne. 0) ti2fe = ti2fe + grk(nti46 + loff,i) if (nti47 .ne. 0) ti2fe = ti2fe + grk(nti47 + loff,i) if (nti48 .ne. 0) ti2fe = ti2fe + grk(nti48 + loff,i) if (nti49 .ne. 0) ti2fe = ti2fe + grk(nti49 + loff,i) if (nti50 .ne. 0) ti2fe = ti2fe + grk(nti50 + loff,i) ti2fe = ti2fe/gas ti2fe = log10(max(tiny,ti2fe/ti_norm * dd)) c..[v/fe] if (nv50 .ne. 0) v2fe = v2fe + grk(nv50 + loff,i) if (nv51 .ne. 0) v2fe = v2fe + grk(nv51 + loff,i) v2fe = v2fe/gas v2fe = log10(max(tiny,v2fe/v_norm * dd)) c..[cr/fe] if (ncr50 .ne. 0) cr2fe = cr2fe + grk(ncr50 + loff,i) if (ncr52 .ne. 0) cr2fe = cr2fe + grk(ncr52 + loff,i) if (ncr53 .ne. 0) cr2fe = cr2fe + grk(ncr53 + loff,i) if (ncr54 .ne. 0) cr2fe = cr2fe + grk(ncr54 + loff,i) cr2fe = cr2fe/gas cr2fe = log10(max(tiny,cr2fe/cr_norm * dd)) c..[mn/fe] if (nmn55 .ne. 0) mn2fe = mn2fe + grk(nmn55 + loff,i) mn2fe = mn2fe/gas mn2fe = log10(max(tiny,mn2fe/mn_norm * dd)) c..[co/fe] if (nco59 .ne. 0) co2fe = co2fe + grk(nco59 + loff,i) co2fe = co2fe/gas co2fe = log10(max(tiny,co2fe/co_norm * dd)) c..[ni/fe] if (nni58 .ne. 0) ni2fe = ni2fe + grk(nni58 + loff,i) if (nni60 .ne. 0) ni2fe = ni2fe + grk(nni60 + loff,i) if (nni61 .ne. 0) ni2fe = ni2fe + grk(nni61 + loff,i) if (nni62 .ne. 0) ni2fe = ni2fe + grk(nni62 + loff,i) if (nni64 .ne. 0) ni2fe = ni2fe + grk(nni64 + loff,i) ni2fe = ni2fe/gas ni2fe = log10(max(tiny,ni2fe/ni_norm * dd)) c..[cu/fe] if (ncu63 .ne. 0) cu2fe = cu2fe + grk(ncu63 + loff,i) if (ncu65 .ne. 0) cu2fe = cu2fe + grk(ncu65 + loff,i) cu2fe = cu2fe/gas cu2fe = log10(max(tiny,cu2fe/cu_norm * dd)) c..[zn/fe] if (nzn64 .ne. 0) zn2fe = zn2fe + grk(nzn64 + loff,i) if (nzn66 .ne. 0) zn2fe = zn2fe + grk(nzn66 + loff,i) if (nzn67 .ne. 0) zn2fe = zn2fe + grk(nzn67 + loff,i) if (nzn68 .ne. 0) zn2fe = zn2fe + grk(nzn68 + loff,i) if (nzn70 .ne. 0) zn2fe = zn2fe + grk(nzn70 + loff,i) zn2fe = zn2fe/gas zn2fe = log10(max(tiny,zn2fe/zn_norm * dd)) c..solution vectors if (zln .gt. -5.0) then aa = zsumf/zsol write(fil2,02) time(i),aa,fe2h,(eratio(k),k=1,7) write(fil3,02) time(i),aa,fe2h,(eratio(k),k=8,14) write(fil4,02) time(i),aa,fe2h,(eratio(k),k=15,21) write(fil8,02) time(i),aa,fe2h,(eratio(k),k=22,28) write(fil9,02) time(i),aa,fe2h,(eratio(k),k=29,neratio) end if c..back for another time point or close up the files enddo call bbb(clse,fil2,string,k) call bbb(clse,fil3,string,k) call bbb(clse,fil4,string,k) call bbb(clse,fil8,string,k) call bbb(clse,fil9,string,k) return end subroutine supernova_histories(nz) implicit none save include 'lunits.dek' include 'parze.dek' include 'chem3.dek' c..writes out the supernova histories for zone nz c..declare the pass integer nz c..local variables integer i,k double precision mlo,mtau,sn2r,sn1r,novr,gas,zsum,zsumf, 1 zln,aa,mtot,mdot,birthm,brate,x,tiny parameter (tiny = 1.0d-100) c..these functions contain the integrands external sn1rint,novrint,type2rn,imf double precision sn1rint,novrint,type2rn,imf c..quadrature storage double precision xxx(nqmax),www(nqmax),xa1(nqmax),wa1(nqmax), 1 xa2(nqmax),wa2(nqmax) c..format statements 01 format(' *',/, 1 ' ',t5,a,t17,a,t29,a,t41,a,t53,a,t65,a, 2 t77,a,t89,a,t101,a,t113,a,t125,a,t137,a) 02 format(1x,1p14e12.4) 08 format(a20,'_',a,'_z',i4,'.dat') 13 format(a,i2.2) c..get some of the quadrature weights call gauleg(mbmax,sn2mhi,xa1,wa1,nnhi) c..open the supernova rate files write(string,08) snfil,'type2',nz call sqeeze(string) call bbb(redwrt,fil3,string,i) call header(fil3) write(fil3,01) 'time','m(tau)','birthm','sn2rn' write(string,08) snfil,'type1',nz call sqeeze(string) call bbb(redwrt,fil8,string,i) call header(fil8) write(fil8,01) 'time','sn1rn',(t1a_comment(k),k=1,nt1a) write(string,08) snfil,'novae',nz call sqeeze(string) call bbb(redwrt,fil9,string,i) call header(fil9) write(fil9,01) 'time','novrn',(nova_comment(k),k=1,nnova) c..for each time point do i=1,kount timcom = time(i) mlo = mtau(time(i)) sn2r = 0.0d0 sn1r = 0.0d0 novr = 0.0d0 gas = gashist(i) zsum = zsumhist(i) zsumf = zsum/gas zln = log10(max(tiny,zsumf/zsol)) aa = max(mbmin,mlo) call gauleg(aa,mbmax,xa2,wa2,nnhi) c..get the birth rate so one can see the retardation call surfden(time(i),mtot,mdot) birthm = brate(mtot,gas) c..type 1a rate in this zone for type 1's if (mlo .lt. mbmax .and. sn1fac .ne. 0.0) then call qgaus(sn1rint,xa2,wa2,nnhi,sn1r) sn1r = sn1fac * sn1r end if c..nova rate for this zone if (mlo .lt. mbmax .and. novfac .ne. 0.0) then call qgaus(novrint,xa2,wa2,nnhi,novr) novr = trecur * novfac * novr end if c..type 2 rate for this zone if (mlo .lt. sn2mhi) then aa = max(sn2mlo,mlo) if (aa .lt. mbmax) then call gauleg(aa,mbmax,xxx,www,nnhi) call qgaus(type2rn,xxx,www,nnhi,sn2r) sn2r = (1.0d0 - sn1fac - novfac) * sn2r call qgaus(type2rn,xa1,wa1,nnhi,x) sn2r = sn2r + x else call gauleg(aa,sn2mhi,xxx,www,nnhi) call qgaus(type2rn,xxx,www,nnhi,sn2r) end if end if c..store for the radial integrals novtab(i,nz) = novr sn1tab(i,nz) = sn1r sn2tab(i,nz) = sn2r c..and write if (zln .gt. -5.0) then write(fil3,02) time(i),mlo,birthm,sn2r write(fil8,02) time(i),sn1r,(sn1r*fract1a(k),k=1,nt1a) write(fil9,02) time(i),novr,(sn1r*fracnova(k),k=1,nnova) end if c..back for another time point or close up enddo call bbb(clse,fil3,snfil,i) call bbb(clse,fil8,snfil,i) call bbb(clse,fil9,snfil,i) return end subroutine kinematics(nz) implicit none save include 'lunits.dek' include 'parze.dek' include 'chem3.dek' c..writes out the kinematic history for zone nz c..declare the pass integer nz c..local variables external imf integer i double precision mlo,mtau,gas,zsum,mtot,mdot,stars,gasf, 1 starf,zsumf,imf,imffac,xxx(nqmax),www(nqmax) c..format statements 01 format(' *',/, 1 ' ',t5,a,t17,a,t29,a,t41,a,t53,a,t65,a, 2 t77,a,t89,a,t101,a,t113,a,t125,a,t137,a) 02 format(1x,1p14e12.4) 07 format(a20,'_z',i4,'.dat') c..open the surface densities, accretion rates c..and global metallicity file write(string,07) densfil,nz call sqeeze(string) call bbb(redwrt,fil2,string,i) call header(fil2) write(fil2,01) 'time','m(tau)','gasmass','gas/tot', 1 'stars','star/tot','tot mass','mdot', 2 'z','z/tot','imf' c..for each time point do i = 1,kount mlo = mtau(time(i)) gas = gashist(i) zsum = zsumhist(i) c..surface density and accretion rate call surfden(time(i),mtot,mdot) c..total mass in stars and fractions of the total present mass stars = mtot - gas gasf = gas/mtot starf = stars/mtot zsumf = zsum/gas c..imf fraction imffac = 0.0d0 if (mlo .lt. imfmax) then call gauleg(mlo,imfmax,xxx,www,nnhi) call qgaus(imf,xxx,www,nnhi,imffac) imffac = imffac/ximf_norm end if c..global kinematics and supernova rates write(fil2,02) time(i),mlo,gas,gasf,stars,starf, 1 mtot,mdot,zsum,zsumf,imffac c..back for another time point or close up enddo call bbb(clse,fil2,densfil,i) return end subroutine gdwarf_distribution(nz) implicit none save include 'lunits.dek' include 'parze.dek' include 'chem3.dek' c..writes out the g-dwarf distribution c..declare the pass integer nz c..local variables external gdwfunc integer loff,i,j,igdwrf double precision gas,zsum,hsum,fesum,mtot,mdot,brate,z, 1 felo,fehi,femax,tlo,thi,gdwfunc,aa, 2 radeps,tiny parameter (tiny = 1.0d-30, radeps = 1.0e-6) c..for the gdwarf integrals integer ngdwmax parameter (ngdwmax = 30) double precision fe2htab(tphys),fecoef(tphys),stlo(ngdwmax), 1 sthi(ngdwmax),sfebin(ngdwmax),sgnum(ngdwmax) double precision bratab(tphys),timetab(tphys),bcoef(tphys) common /csexdw/ bratab,timetab,bcoef,igdwrf c..format statements 01 format(' *',/, 1 ' ',t5,a,t17,a,t29,a,t41,a,t53,a,t65,a, 2 t77,a,t89,a,t101,a,t113,a,t125,a,t137,a) 02 format(1x,1p14e12.4) 07 format(a20,'_z',i4,'.dat') c..set the zone offset loff = (nz - 1)*niso c..write out the gdwarf distribution; note that with the c..normalization to the total, the fraction of gdwarfs from c..the imf drops out. the expression being evaluated for the c..gdwarf distribution is from pardi & ferrini apj 421,491 1994 write(string,07) gdwfil,nz call sqeeze(string) call bbb(redwrt,fil2,string,i) call header(fil2) write(fil2,01) 'tlo','thi','[fe/h]','raw','normed' c..sweep over time igdwrf = 0 do i=1,kount gas = gashist(i) zsum = zsumhist(i) c..hydrogen gas fraction hsum = grk(1 + loff,i)/gas c..total iron fesum = tiny if (nfe54 .ne. 0) fesum = fesum + grk(nfe54 + loff,i) if (nfe56 .ne. 0) fesum = fesum + grk(nfe56 + loff,i) if (nfe57 .ne. 0) fesum = fesum + grk(nfe57 + loff,i) if (nfe58 .ne. 0) fesum = fesum + grk(nfe58 + loff,i) fesum = fesum/gas fe2h = log10(max(tiny,fesum/fesol * hsol/hsum)) if (fe2h .gt. -2.0) then igdwrf = igdwrf + 1 timetab(igdwrf) = time(i) fe2htab(igdwrf) = fe2h call surfden(time(i),mtot,mdot) bratab(igdwrf) = brate(mtot,gas) end if enddo c..spline fit the stored tables call spline(fe2htab,timetab,igdwrf,fecoef) call spline(timetab,bratab,igdwrf,bcoef) z = 0.0d0 felo = -1.20d0 femax = fe2htab(igdwrf) c..crank the integrals do i = 1,14 fehi = min(felo + 0.1d0,femax) j = i sfebin(i) = 0.5d0 * (felo + fehi) call splint(fe2htab,timetab,fecoef,igdwrf,felo,tlo) stlo(i) = tlo call splint(fe2htab,timetab,fecoef,igdwrf,fehi,thi) sthi(i) = thi call qromb(gdwfunc,tlo,thi,radeps,aa) sgnum(i) = aa z = z + aa felo = fehi if (fehi .eq. femax) go to 83 enddo 83 continue c..write it out z = 1.0d0/z do i=1,j write(fil2,02) stlo(i),sthi(i),sfebin(i),sgnum(i),sgnum(i)*z enddo call bbb(clse,fil2,gdwfil,i) return end subroutine solar_distribution(nz) implicit none save include 'lunits.dek' include 'parze.dek' include 'chem3.dek' c..writes out the solar composition comparison file c..declare the pass integer nz c..local variables integer i double precision gas,msol,mdot,birthm,sn2r,sn1r, 1 xsol,xnow,aa,zsumf,mtot,mdot,brate c..format statements 04 format(1x,a,t30,': ',1pe11.3,t45,a) 05 format(1x,2i4,' ',1p2e11.3,' ',a5,' ',1pe11.3) 06 format(1x,a,t30,': ',1pe11.3,t45,'at time = ',1pe11.3, 1 ' zsol =',1pe11.3) 07 format(a20,'_z',i4,'.dat') c..the mean ism at solar birth and anders & grevese values. write(string,07) agfil,nz call sqeeze(string) call bbb(redwrt,fil2,string,i) call header(fil2) call surfden(time(kount),mtot,mdot) birthm = brate(mtot,gashist(kount)) sn2r = sn2tab(kount,nz) sn1r = sn1tab(kount,nz) write(fil2,04) '* present gas density', 1 gashist(kount),'msol/pc**2' write(fil2,04) '* present accretion rate',mdot,'msol/pc**2/Gyr' write(fil2,04) '* present birthrate',birthm,'msol/pc**2/Gyr' write(fil2,04) '* present type II rate',sn2r,'num/pc**2/Gyr' write(fil2,04) '* present type Ia rate',sn1r,'num/pc**2/Gyr' write(fil2,04) '* present delta Y delta Z',delydelz,' ' c..the solar mass fractions at the desired time output point do i=1,niso if (i .ne. nal26 .and. i .ne. nfe60) then xsol = solnorm(i) nowiso = i call getgas(tagout,gas,zsumf,xnow) if (i .eq. 1) then write(fil2,06) '* gas metallicity',zsumf,tagout,zsol write(fil2,*) '* z a x x/xsol name xsol' end if aa = xnow/xsol write(fil2,05) int(zsolnorm(i)),int(asolnorm(i)), 1 xnow,aa,namnuc(i),xsol end if enddo call bbb(clse,fil2,agfil,i) return end subroutine radial_sums implicit none save include 'lunits.dek' include 'parze.dek' include 'chem3.dek' c..writes out the radial sums c..declare the pass integer nz c..local variables integer i,j,k,loff,ll double precision conv,pi,rlo,rhi,rtot,rgas,rt2rat,rt1rat,rnvrat, 1 ral26,ralt2,rfe60,rfet2,mtot,mdot,radialf, 2 gas parameter (pi = 3.1415926535897932384d0) double precision radeps,tiny,thalf1,tmean1,thalf2,tmean2 parameter (tiny = 1.0d-100, radeps = 1.0e-6, 1 thalf1 = 7.40d-4, 2 tmean1 = thalf1/0.693147181d0, 3 thalf2 = 1.50d-3, 4 tmean2 = thalf2/0.693147181d0) c..for the radial integrals double precision radum(nzmax),radcof(nzmax) common /chokey/ radum,radcof c..format statements 01 format(' *',/, 1 ' ',t5,a,t17,a,t29,a,t41,a,t53,a,t65,a, 2 t77,a,t89,a,t101,a,t113,a,t125,a,t137,a) 02 format(1x,1p14e12.4) 12 format(a20,a) 13 format(a,i2.2) c..open the files write(string,12) snfil,'_type2_rad.dat' call sqeeze(string) call bbb(redwrt,fil2,string,i) call header(fil2) write(fil2,01) 'time','mtot','mgas','snII', 1 'al26dot','al26','fe60dot','fe60' write(string,12) snfil,'_type1_rad.dat' call sqeeze(string) call bbb(redwrt,fil7,string,i) call header(fil7) write(fil7,01) 'time','snI',(t1a_comment(k),k=1,nt1a) write(string,12) snfil,'_novae_rad.dat' call sqeeze(string) call bbb(redwrt,fil8,string,i) call header(fil8) write(fil8,01) 'time','nova',(nova_comment(k),k=1,nnova) conv = 2.0d0 * pi rlo = rpoint(1) rhi = rpoint(nzone) c..sweep over time on the outer loop do i = 1,kount rtot = 0.0d0 rgas = 0.0d0 rt2rat = 0.0d0 rt1rat = 0.0d0 rnvrat = 0.0d0 ral26 = 0.0d0 ralt2 = 0.0d0 rfe60 = 0.0d0 rfet2 = 0.0d0 c..total mass do j=1,nzone nowzon = j call surfden(time(i),mtot,mdot) radum(j) = mtot * conv * rpoint(j) enddo call spline(rpoint,radum,nzone,radcof) call qromb(radialf,rlo,rhi,radeps,rtot) c..total gas do j=1,nzone loff = (j - 1)*niso gas = 0.0d0 do k=1,niso gas = gas + grk(k + loff,i) enddo radum(j) = gas * conv * rpoint(j) enddo call spline(rpoint,radum,nzone,radcof) call qromb(radialf,rlo,rhi,radeps,rgas) c..type II rate per century do j=1,nzone radum(j) = sn2tab(i,j) * conv * rpoint(j) * 1.0e-7 enddo call spline(rpoint,radum,nzone,radcof) call qromb(radialf,rlo,rhi,radeps,rt2rat) c..type I rate per century if (sn1fac .ne. 0.0) then do j=1,nzone radum(j) = sn1tab(i,j) * conv * rpoint(j) * 1.0e-7 enddo call spline(rpoint,radum,nzone,radcof) call qromb(radialf,rlo,rhi,radeps,rt1rat) end if c..nova rate per year if (novfac .ne. 0.0) then do j = 1,nzone radum(j) = novtab(i,j) * conv * rpoint(j) * 1.0e-9 enddo call spline(rpoint,radum,nzone,radcof) call qromb(radialf,rlo,rhi,radeps,rnvrat) end if c..aluminum 26 injection rate per million years and total if (nal26 .ne. 0) then do j = 1,nzone loff = (j - 1)*niso ll = nal26 + loff ral26 = abs(dgrk(ll,i)) if (modal26 .gt. 0) ral26 = ral26 + abs(grk(ll,i)/tmean1) radum(j) = ral26 * conv * rpoint(j) * 1.0e-3 enddo call spline(rpoint,radum,nzone,radcof) call qromb(radialf,rlo,rhi,radeps,ral26) do j = 1,nzone loff = (j - 1)*niso ll = nal26 + loff ralt2 = grk(ll,i) radum(j) = ralt2 * conv * rpoint(j) enddo call spline(rpoint,radum,nzone,radcof) call qromb(radialf,rlo,rhi,radeps,ralt2) end if c..iron 60 injection rate per million years and total if (nfe60 .ne. 0) then do j = 1,nzone loff = (j - 1)*niso ll = nfe60 + loff rfe60 = abs(dgrk(ll,i)) if (modal26 .gt. 0) rfe60 = rfe60 + abs(grk(ll,i)/tmean2) radum(j) = rfe60 * conv * rpoint(j) * 1.0e-3 enddo call spline(rpoint,radum,nzone,radcof) call qromb(radialf,rlo,rhi,radeps,rfe60) do j = 1,nzone loff = (j - 1)*niso ll = nfe60 + loff rfet2 = grk(ll,i) radum(j) = rfet2* conv * rpoint(j) enddo call spline(rpoint,radum,nzone,radcof) call qromb(radialf,rlo,rhi,radeps,rfet2) end if c..write it out write(fil2,02) time(i),rtot,rgas,rt2rat,ral26,ralt2,rfe60,rfet2 write(fil7,02) time(i),rt1rat,(rt1rat*fract1a(k),k=1,nt1a) write(fil8,02) time(i),rnvrat,(rnvrat*fracnova(k),k=1,nnova) c..back for another time point or close it up enddo call bbb(clse,fil2,string,i) call bbb(clse,fil7,string,i) call bbb(clse,fil8,string,i) return end double precision function radialf(x1) implicit none save include 'chem3.dek' c..for the radial integrals c.. c..declare double precision x1,y1 c..a small common for the radial integrals only double precision radum(nzmax),radcof(nzmax) common /chokey/ radum,radcof call splint(rpoint,radum,radcof,nzone,x1,y1) radialf = y1 return end double precision function gdwfunc(t1) implicit none save include 'chem3.dek' c..for the gdwarf integrals c..declare double precision t1,y1 c..small common for the gdwarf integrals only integer igdwrf double precision bratab(tphys),timetab(tphys),bcoef(tphys) common /csexdw/ bratab,timetab,bcoef,igdwrf call splint(timetab,bratab,bcoef,igdwrf,t1,y1) gdwfunc = y1 return end subroutine header(luz) implicit none save include 'chem3.dek' c..this routine writes out the date, time, and radius to logical unit luz c..declare the pass integer luz c..format statements 01 format(1x,'*',t4,a,t16,a,t27,a) 02 format(1x,'*',t4,a,t16,1pe11.3) write(luz,01) adat,atim if (luz .ne. 6) write(luz,02) 'radnow=',rpoint(nowzon) return end subroutine rnuc implicit none save include 'lunits.dek' include 'parze.dek' include 'chem3.dek' c..this routine reads in the stellar yields, c..sets up pointers, and c..checks the input. c..declare integer i,j,k,m,n,ii,kend,jj,mm,nn,kat,kfirst,jt1a,jnova double precision sum,mcheck,msave(masmax),zsave(metmax),zcheck,xm data kfirst/0/ c..formats 01 format(a5) 02 format(a80) 03 format(4(1x,a,' ',1pe10.2,' ')) 04 format(1x,a,1pe11.3,a,1pe11.3) 05 format(4x,a,a,1pe11.3,a,1pe11.3,a) 06 format(4x,a,1pe11.3,a,1pe11.3) c..open the nucleosynthesis file call bbb(redold,fil2,nukfil,i) if (i.eq. 1) stop 'could not open nukfil' c..start the read; get the name niso = 0 do i=1,gphys read(fil2,01,err=51,end=50) namnuc(i) c..make sure we dont count the remnants; they is special if (namnuc(i) .ne. 'rem ') niso = niso + 1 if (niso+1 .gt. isomax) then write(6,*) 'error on niso read',niso write(6,05) 'isotope name',namnuc(i) stop 'niso > isomax in routine rnuc' end if c..the bigbang mass fraction read(fil2,*,err=51) bigb(i) c..the number of type1a and their yields read(fil2,*,err=51) nt1a if (nt1a .gt. maxt1a) then write(6,*) 'number of type1a too big ',nt1a write(6,05) 'isotope name',namnuc(i) stop 'nt1a > maxt1a in routine rnuc' end if if (nt1a .ne. nfract1a) then write(6,*) 'nt1a .ne. nfract1a',nt1a,nfract1a write(6,05) 'isotope name',namnuc(i) stop 'nt1a .ne. nfract1a in routine rnuc' end if if (i .eq. 1) then jt1a = nt1a else if (nt1a .ne. jt1a) then write(6,*) 'unequal number of type1a ',jt1a,nt1a write(6,05) 'isotope name',namnuc(i) stop 'jt1a .ne. nt1a in routine rnuc' end if end if do j=1,nt1a read(fil2,*) mt1a(j,i) enddo c..the number of nova and their yields read(fil2,*,err=51) nnova if (nnova .gt. maxnova) then write(6,*) 'number of nova too big ',nnova write(6,05) 'isotope name',namnuc(i) stop 'nnova > maxnova in routine rnuc' end if if (nnova .ne. nfracnova) then write(6,*) 'nnova .ne. nfracnova',nt1a,nfracnova write(6,05) 'isotope name',namnuc(i) stop 'nnova .ne. nfracnova in routine rnuc' end if if (i .eq. 1) then jnova = nnova else if (nnova .ne. jnova) then write(6,*) 'unequal number of nova ',jnova,nnova write(6,05) 'isotope name',namnuc(i) stop 'jnova .ne. nnova in routine rnuc' end if end if do j=1,nnova read(fil2,*) mnova(j,i) enddo c..the number of metallicities read(fil2,*,err=51) metnuc(i) if (metnuc(i) .gt. metmax) then write(6,*) 'error around isotope number',niso write(6,05) 'isotope name',namnuc(i) stop ' > metmax points in rnuc' end if c..the metallicity read and loop do j=1,metnuc(i) read(fil2,*,err=51) zzznuc(j,i) c..number of mass points read(fil2,*,err=51) kend if (kend .gt. masmax) then write(6,*) 'error around isotope number',niso write(6,05) 'isotope name',namnuc(i) stop ' > masmax mass points in rnuc' end if numnuc(j,i) = kend c..stellar masses and yields do k=1,kend read(fil2,*,err=51) masnuc(j,k,i),isonuc(j,k,i) enddo c..end of metalicity loop enddo c..end of isotope loop enddo c..land here to close the file 50 call bbb(clse,fil2,nukfil,i) goto 52 c..a bad read somewhere in the above 51 write(6,*) 'error around isotope number',niso write(6,05) 'isotope name',namnuc(i) stop 'bad read in rnuc' c..checks on the input c..make sure the h1/remnant was first/last in the list 52 if (namnuc(1) .ne. 'h1 ') stop 'h1 not first in list' if (namnuc(niso+1) .ne. 'rem ') stop 'remnant is not last' c..set the metallicity start and end pointers nzend = niso nzbeg = 0 do i=1,niso if (.not. (namnuc(i) .eq. 'h1 ' .or. namnuc(i) .eq. 'h2 '.or. 1 namnuc(i) .eq. 'h3 ' .or. namnuc(i) .eq. 'he3 '.or. 2 namnuc(i) .eq. 'he4 ' )) then nzbeg = i goto 61 end if enddo if (nzbeg .eq. 0) stop 'could not set metal pointer nzbeg' 61 continue c..specific isotope pointers nh2 = 0 nhe3 = 0 nhe4 = 0 nli6 = 0 nli7 = 0 nbe9 = 0 nb10 = 0 nb11 = 0 nc12 = 0 nc13 = 0 nn14 = 0 nn15 = 0 no16 = 0 no17 = 0 no18 = 0 nf19 = 0 nne20 = 0 nne21 = 0 nne21 = 0 nna23 = 0 nmg24 = 0 nmg25 = 0 nmg26 = 0 nal26 = 0 nal27 = 0 nsi28 = 0 nsi29 = 0 nsi30 = 0 np31 = 0 ns32 = 0 ns33 = 0 ns34 = 0 ns36 = 0 ncl35 = 0 ncl37 = 0 nar36 = 0 nar38 = 0 nar40 = 0 nk39 = 0 nk40 = 0 nk41 = 0 nca40 = 0 nca42 = 0 nca43 = 0 nca44 = 0 nca46 = 0 nca48 = 0 nsc45 = 0 nti46 = 0 nti47 = 0 nti48 = 0 nti49 = 0 nti50 = 0 nv50 = 0 nv51 = 0 ncr50 = 0 ncr52 = 0 ncr53 = 0 ncr54 = 0 nmn55 = 0 nfe54 = 0 nfe56 = 0 nfe57 = 0 nfe58 = 0 nfe60 = 0 nco59 = 0 nni58 = 0 nni60 = 0 nni61 = 0 nni62 = 0 nni64 = 0 ncu63 = 0 ncu65 = 0 nzn64 = 0 nzn66 = 0 nzn67 = 0 nzn68 = 0 nzn70 = 0 do i=1,niso if (namnuc(i) .eq. 'h2 ') then nh2 = i else if (namnuc(i) .eq. 'he3 ') then nhe3 = i else if (namnuc(i) .eq. 'he4 ') then nhe4 = i else if (namnuc(i) .eq. 'li6 ') then nli6 = i else if (namnuc(i) .eq. 'li7 ') then nli7 = i else if (namnuc(i) .eq. 'be9 ') then nbe9 = i else if (namnuc(i) .eq. 'b10 ') then nb10 = i else if (namnuc(i) .eq. 'b11 ') then nb11 = i else if (namnuc(i) .eq. 'c12 ') then nc12 = i else if (namnuc(i) .eq. 'c13 ') then nc13 = i else if (namnuc(i) .eq. 'n14 ') then nn14 = i else if (namnuc(i) .eq. 'n15 ') then nn15 = i else if (namnuc(i) .eq. 'o16 ') then no16 = i else if (namnuc(i) .eq. 'o17 ') then no17 = i else if (namnuc(i) .eq. 'o18 ') then no18 = i else if (namnuc(i) .eq. 'f19 ') then nf19 = i else if (namnuc(i) .eq. 'ne20 ') then nne20 = i else if (namnuc(i) .eq. 'ne21 ') then nne21 = i else if (namnuc(i) .eq. 'ne22 ') then nne22 = i else if (namnuc(i) .eq. 'na23 ') then nna23 = i else if (namnuc(i) .eq. 'mg24 ') then nmg24 = i else if (namnuc(i) .eq. 'mg25 ') then nmg25 = i else if (namnuc(i) .eq. 'mg26 ') then nmg26 = i else if (namnuc(i) .eq. 'al26 ') then nal26 = i else if (namnuc(i) .eq. 'al27 ') then nal27 = i else if (namnuc(i) .eq. 'si28 ') then nsi28 = i else if (namnuc(i) .eq. 'si29 ') then nsi29 = i else if (namnuc(i) .eq. 'si30 ') then nsi30 = i else if (namnuc(i) .eq. 'p31 ') then np31 = i else if (namnuc(i) .eq. 's32 ') then ns32 = i else if (namnuc(i) .eq. 's33 ') then ns33 = i else if (namnuc(i) .eq. 's34 ') then ns34 = i else if (namnuc(i) .eq. 's36 ') then ns36 = i else if (namnuc(i) .eq. 'cl35 ') then ncl35 = i else if (namnuc(i) .eq. 'cl37 ') then ncl37 = i else if (namnuc(i) .eq. 'ar36 ') then nar36 = i else if (namnuc(i) .eq. 'ar38 ') then nar38 = i else if (namnuc(i) .eq. 'ar40 ') then nar40 = i else if (namnuc(i) .eq. 'k39 ') then nk39 = i else if (namnuc(i) .eq. 'k40 ') then nk40 = i else if (namnuc(i) .eq. 'ca40 ') then nca40 = i else if (namnuc(i) .eq. 'ca42 ') then nca42 = i else if (namnuc(i) .eq. 'ca43 ') then nca43 = i else if (namnuc(i) .eq. 'ca44 ') then nca44 = i else if (namnuc(i) .eq. 'ca46 ') then nca46 = i else if (namnuc(i) .eq. 'ca48 ') then nca48 = i else if (namnuc(i) .eq. 'sc45 ') then nsc45 = i else if (namnuc(i) .eq. 'ti46 ') then nti46 = i else if (namnuc(i) .eq. 'ti47 ') then nti47 = i else if (namnuc(i) .eq. 'ti48 ') then nti48 = i else if (namnuc(i) .eq. 'ti49 ') then nti49 = i else if (namnuc(i) .eq. 'ti50 ') then nti50 = i else if (namnuc(i) .eq. 'v50 ') then nv50 = i else if (namnuc(i) .eq. 'v51 ') then nv51 = i else if (namnuc(i) .eq. 'cr50 ') then ncr50 = i else if (namnuc(i) .eq. 'cr52 ') then ncr52 = i else if (namnuc(i) .eq. 'cr53 ') then ncr53 = i else if (namnuc(i) .eq. 'cr54 ') then ncr54 = i else if (namnuc(i) .eq. 'mn55 ') then nmn55 = i else if (namnuc(i) .eq. 'fe54 ') then nfe54 = i else if (namnuc(i) .eq. 'fe56 ') then nfe56 = i else if (namnuc(i) .eq. 'fe57 ') then nfe57 = i else if (namnuc(i) .eq. 'fe58 ') then nfe58 = i else if (namnuc(i) .eq. 'fe60 ') then nfe60 = i else if (namnuc(i) .eq. 'co59 ') then nco59 = i else if (namnuc(i) .eq. 'ni58 ') then nni58 = i else if (namnuc(i) .eq. 'ni60 ') then nni60 = i else if (namnuc(i) .eq. 'ni61 ') then nni61 = i else if (namnuc(i) .eq. 'ni62 ') then nni62 = i else if (namnuc(i) .eq. 'ni64 ') then nni64 = i else if (namnuc(i) .eq. 'cu63 ') then ncu63 = i else if (namnuc(i) .eq. 'cu65 ') then ncu65 = i else if (namnuc(i) .eq. 'zn64 ') then nzn64 = i else if (namnuc(i) .eq. 'zn66 ') then nzn66 = i else if (namnuc(i) .eq. 'zn67 ') then nzn67 = i else if (namnuc(i) .eq. 'zn68 ') then nzn68 = i else if (namnuc(i) .eq. 'zn70 ') then nzn70 = i end if enddo c..make sure we have big bang mass fraction conservation sum = 0.0d0 do i=1,niso sum = sum + bigb(i) enddo sum = 1.0d0 - sum if (sum .ne. 0.0d0) then write(6,04) 'bigbang 1- sum ',sum,' not zero; putting into h1' bigb(1) = bigb(1) + sum end if write(6,*) 'bigbang nucleosynthesis:' do i=1,niso if (bigb(i) .ne. 0.0) write(6,03) namnuc(i),bigb(i) enddo primhe4 = bigb(nhe4) c..sum the type1a masses do i=1,nt1a sum = 0.0d0 do j=1,niso sum = sum + mt1a(i,j) enddo write(6,*) ' ' write(6,04) 'type 1a mass sum ',sum write(6,*) 'type 1a nucleosynthesis:' write(6,03) (namnuc(j),mt1a(i,j), j=1,niso) enddo c..sum the nova masses do i=1,nnova sum = 0.0d0 do j=1,niso sum = sum + mnova(i,j) enddo write(6,*) ' ' write(6,04) 'nova mass sum ',sum write(6,*) 'nova nucleosynthesis:' write(6,03) (namnuc(j),mnova(i,j), j=1,niso) enddo c..check the nucleosynthesis sums. the data structure was chosen to make the c..computation have efficient do loops; this makes the error checking done c..here less than elegant or efficient write(6,*) ' ' write(6,*) 'checking nucleosynthesis sums:' nn = 0 do 170 jj=1,metmax mm = 1 95 continue zcheck = zzznuc(jj,mm) if (jj .gt. 1 .and. zcheck .eq. 0.0) goto 200 do ii=1,nn if (zsave(ii) .eq. zcheck) goto 170 enddo nn = nn + 1 zsave(nn) = zcheck write(6,*) ' ' write(6,04) 'metallicity ',zcheck n = 0 do 160 j=1,masmax m = 1 100 continue mcheck = masnuc(jj,j,m) if (mcheck .eq. 0.0) goto 165 c..check we have not done this mass before; if not increment do ii=1,n if (msave(ii) .eq. mcheck) goto 160 enddo n = n + 1 msave(n) = mcheck c..sum sum = 0.0d0 do i=1,niso+1 do k=1,numnuc(jj,i) if (masnuc(jj,k,i) .eq. mcheck) then if (namnuc(i) .eq. 'he4 ') kat=k c..make sure one ones this amount of output if (iprint .gt. 1) then if (mcheck .ne. isonuc(jj,k,niso+1)) then write(6,05) namnuc(i),' gives ',isonuc(jj,k,i), 1 ' solor masses; = ', 2 isonuc(jj,k,i)/(mcheck-isonuc(jj,k,niso+1)), 3 ' mass fraction' else write(6,05) namnuc(i),' gives ',isonuc(jj,k,i), 1 ' solor masses' end if end if c.. sum = sum + isonuc(jj,k,i) go to 130 end if enddo 130 continue enddo xm = mcheck - sum c..stuff nonconservation into oxygen (metals) usually for al26 only if (modal26 .gt. 0) then if (kfirst .eq. 0) then kfirst = 1 write(6,*) 'stuffing nonconservation into the oxygen' end if isonuc(jj,kat,no16) = isonuc(jj,kat,no16) + xm c..stuff nonconservation into hydrogen else if (kfirst .eq. 0) then kfirst = 1 write(6,*) 'stuffing nonconservation into hydrogen' end if isonuc(jj,kat,1) = isonuc(jj,kat,1) + xm end if c..stuff nonconservation into helium c.. else c.. if (kfirst .eq. 0) then c.. kfirst = 1 c.. write(6,*) 'stuffing nonconservation into the helium' c.. end if c.. isonuc(jj,kat,nhe4) = isonuc(jj,kat,nhe4) + xm c.. end if c..stuff nonconservation into the remnant c.. else c.. if (kfirst .eq. 0) then c.. kfirst = 1 c.. write(6,*) 'stuffing nonconservation into the remnant' c.. end if c.. isonuc(jj,kat,niso+1) = isonuc(jj,kat,niso+1) + xm c.. end if write(6,04) 'stellar mass',mcheck,' nonconservation ',xm c..update the mass pointer if (m+1 .le. numnuc(jj,j)) then m = m + 1 goto 100 end if 160 continue 165 continue c..update the metal pointer if (mm+1 .le. metnuc(jj)) then mm = mm + 1 goto 95 end if 170 continue 200 continue c..a few more blanks is pretty write(6,*) ' ' write(6,*) ' ' return end subroutine getyeld2(z,m,y) implicit none save include 'chem3.dek' c..this routine interpolates/extrapolates the stellar nucleosynthetic yields. c..input is the metallicity z, the mass m, and the isotope now being processed c..nowiso (which enters through common block). c..output is the stellar yield in solar masses. c..declare integer j,k,jat,kat,jend,kend double precision m,y,z,y1,y2,y3,y4,u,t,a,b data jat/1/, kat/1/ c..locate the lower bounding metallicity (assumes zero metals is at j=1) c..interpolate; extrapolate using the next to last metallicity as the lower c..bound; check if the last hit was good if (z .ge. zzznuc(jat,nowiso) .and. 1 z .lt. zzznuc(jat+1,nowiso)) goto 11 jend = metnuc(nowiso) do j=2,jend if (zzznuc(j,nowiso) .ge. z) then jat = j - 1 goto 11 end if enddo jat = jend - 1 11 continue c..locate the lower bounding mass at this lower bounding metallicity c..interpolate; extrapolate using the next to last mass point as the lower c..bound; check if the last hit was good if (m .ge. masnuc(jat,kat,nowiso) .and. 1 m .lt. masnuc(jat,kat+1,nowiso)) goto 21 kend = numnuc(jat,nowiso) do k=2,kend if (masnuc(jat,k,nowiso) .ge. m) then kat = k - 1 goto 21 end if enddo kat = kend - 1 21 continue c..bilinear interpolation; recipes page 96 y1 = isonuc(jat,kat,nowiso) y2 = isonuc(jat+1,kat,nowiso) y3 = isonuc(jat+1,kat+1,nowiso) y4 = isonuc(jat,kat+1,nowiso) a = zzznuc(jat,nowiso) t = (z - a) / (zzznuc(jat+1,nowiso) - a) a = masnuc(jat,kat,nowiso) u = (m - a) /(masnuc(jat,kat+1,nowiso) - a) a = 1.0d0 - t b = 1.0d0 - u y = a*b*y1 + t*b*y2 + t*u*y3 + a*u*y4 return end double precision function andgrev(nam,z,a,xelem) implicit none save c..anders and grevesse solar abundances. c.. c..input is the name the isotope, output is the mass fraction and grev, c..the proton number z, the real number of nucleons a, and the c..elemental mass fraction associated with this isotope. c.. c..warning message and x=1.0e100 if there is no a&g entry. c.. c..declare integer solsiz parameter (solsiz = 166) character*(*) nam character*5 namsol(solsiz) integer izsol(solsiz),iasol(solsiz),jcode(solsiz), 1 i,j,iz,ia,ifirst,jbeg,jend double precision agsol(solsiz),sum,z,a,x,xelem data ifirst/0/ data (namsol(j), j=1,120) / 1 'h1 ','h2 ','he3 ','he4 ','li6 ','li7 ','be9 ','b10 ', 2 'b11 ','c12 ','c13 ','n14 ','n15 ','o16 ','o17 ','o18 ', 3 'f19 ','ne20 ','ne21 ','ne22 ','na23 ','mg24 ','mg25 ','mg26 ', 4 'al27 ','si28 ','si29 ','si30 ','p31 ','s32 ','s33 ','s34 ', 5 's36 ','cl35 ','cl37 ','ar36 ','ar38 ','ar40 ','k39 ','k40 ', 6 'k41 ','ca40 ','ca42 ','ca43 ','ca44 ','ca46 ','ca48 ','sc45 ', 7 'ti46 ','ti47 ','ti48 ','ti49 ','ti50 ','v50 ','v51 ','cr50 ', 8 'cr52 ','cr53 ','cr54 ','mn55 ','fe54 ','fe56 ','fe57 ','fe58 ', 9 'co59 ','ni58 ','ni60 ','ni61 ','ni62 ','ni64 ','cu63 ','cu65 ', & 'zn64 ','zn66 ','zn67 ','zn68 ','zn70 ','ga69 ','ga71 ','ge70 ', 1 'ge72 ','ge73 ','ge74 ','ge76 ','as75 ','se74 ','se76 ','se77 ', 2 'se78 ','se80 ','se82 ','br79 ','br81 ','kr78 ','kr80 ','kr82 ', 3 'kr83 ','kr84 ','kr86 ','rb85 ','rb87 ','sr84 ','sr86 ','sr87 ', 4 'sr88 ','y89 ','zr90 ','zr91 ','zr92 ','zr94 ','zr96 ','nb93 ', 5 'mo92 ','mo94 ','mo95 ','mo96 ','mo97 ','mo98 ','mo100','ru96 '/ data (namsol(j), j=121,solsiz) / 1 'ru98 ','ru99 ','ru100','ru101','ru102','ru104','rh103','pd102', 2 'pd104','pd105','pd106','pd108','pd110','ag107','ag109','cd106', 3 'cd108','cd110','cd111','cd112','cd113','cd114','cd116','in113', 4 'in115','sn112','sn114','sn115','sn116','sn117','sn118','sn119', 5 'sn120','sn122','sn124','sb121','sb123','te120','te122','te123', 6 'te124','te125','te126','te128','te130','i127 '/ data (agsol(j), j=1,105) / 1 7.06e-01,4.80e-05,2.93e-05,2.75e-01,6.50e-10,9.35e-09,1.66e-10, 2 1.07e-09,4.73e-09,3.03e-03,3.65e-05,1.10e-03,4.36e-06,9.59e-03, 3 3.89e-06,2.17e-05,4.05e-07,1.62e-03,4.13e-06,1.30e-04,3.34e-05, 4 5.15e-04,6.77e-05,7.76e-05,5.80e-05,6.53e-04,3.43e-05,2.35e-05, 5 8.16e-06,3.96e-04,3.22e-06,1.87e-05,9.38e-08,2.53e-06,8.54e-07, 6 7.74e-05,1.54e-05,2.53e-08,3.47e-06,5.54e-09,2.63e-07,5.99e-05, 7 4.20e-07,8.97e-08,1.42e-06,2.79e-09,1.38e-07,3.89e-08,2.23e-07, 8 2.08e-07,2.15e-06,1.64e-07,1.64e-07,9.26e-10,3.77e-07,7.42e-07, 9 1.49e-05,1.72e-06,4.36e-07,1.33e-05,7.13e-05,1.17e-03,2.85e-05, & 3.70e-06,3.36e-06,4.94e-05,1.96e-05,8.59e-07,2.78e-06,7.27e-07, 1 5.75e-07,2.65e-07,9.92e-07,5.88e-07,8.76e-08,4.06e-07,1.38e-08, 2 3.96e-08,2.71e-08,4.32e-08,5.94e-08,1.71e-08,8.12e-08,1.78e-08, 3 1.24e-08,1.03e-09,1.08e-08,9.15e-09,2.90e-08,6.25e-08,1.18e-08, 4 1.20e-08,1.19e-08,3.02e-10,2.02e-09,1.07e-08,1.08e-08,5.46e-08, 5 1.71e-08,1.10e-08,4.64e-09,2.80e-10,5.05e-09,3.32e-09,4.32e-08/ data (agsol(j), j=106,solsiz) / 1 1.04e-08,1.34e-08,2.95e-09,4.56e-09,4.71e-09,7.77e-10,1.64e-09, 2 8.80e-10,5.61e-10,9.76e-10,1.03e-09,5.99e-10,1.52e-09,6.22e-10, 3 2.50e-10,8.68e-11,5.91e-10,5.92e-10,8.07e-10,1.52e-09,9.15e-10, 4 8.96e-10,3.66e-11,4.08e-10,8.23e-10,1.02e-09,1.01e-09,4.54e-10, 5 6.82e-10,6.45e-10,5.39e-11,3.91e-11,5.59e-10,5.78e-10,1.10e-09, 6 5.63e-10,1.34e-09,3.55e-10,2.26e-11,5.12e-10,1.05e-10,7.18e-11, 7 3.75e-11,1.63e-09,8.67e-10,2.76e-09,4.87e-10,3.79e-09,5.46e-10, 8 6.93e-10,5.42e-10,4.11e-10,1.31e-11,3.83e-10,1.33e-10,7.18e-10, 9 1.08e-09,2.90e-09,4.95e-09,5.36e-09,2.89e-09/ data izsol/1, 1, 2, 2, 3, 3, 4, 5, 5, 6, 6, 7, 7, 8, 1 8, 8, 9, 10, 10, 10, 11, 12, 12, 12, 13, 14, 14, 14, 2 15, 16, 16, 16, 16, 17, 17, 18, 18, 18, 19, 19, 19, 20, 3 20, 20, 20, 20, 20, 21, 22, 22, 22, 22, 22, 23, 23, 24, 4 24, 24, 24, 25, 26, 26, 26, 26, 27, 28, 28, 28, 28, 28, 5 29, 29, 30, 30, 30, 30, 30, 31, 31, 32, 32, 32, 32, 32, 6 33, 34, 34, 34, 34, 34, 34, 35, 35, 36, 36, 36, 36, 36, 7 36, 37, 37, 38, 38, 38, 38, 39, 40, 40, 40, 40, 40, 41, 8 42, 42, 42, 42, 42, 42, 42, 44, 44, 44, 44, 44, 44, 44, 9 45, 46, 46, 46, 46, 46, 46, 47, 47, 48, 48, 48, 48, 48, 1 48, 48, 48, 49, 49, 50, 50, 50, 50, 50, 50, 50, 50, 50, 2 50, 51, 51, 52, 52, 52, 52, 52, 52, 52, 52, 53/ data iasol / 1, 2, 3, 4, 6, 7, 9, 10, 11, 12, 13, 14, 15, 16, 1 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 2 31, 32, 33, 34, 36, 35, 37, 36, 38, 40, 39, 40, 41, 40, 3 42, 43, 44, 46, 48, 45, 46, 47, 48, 49, 50, 50, 51, 50, 4 52, 53, 54, 55, 54, 56, 57, 58, 59, 58, 60, 61, 62, 64, 5 63, 65, 64, 66, 67, 68, 70, 69, 71, 70, 72, 73, 74, 76, 6 75, 74, 76, 77, 78, 80, 82, 79, 81, 78, 80, 82, 83, 84, 7 86, 85, 87, 84, 86, 87, 88, 89, 90, 91, 92, 94, 96, 93, 8 92, 94, 95, 96, 97, 98,100, 96, 98, 99,100,101,102,104, 9103,102,104,105,106,108,110,107,109,106,108,110,111,112, 1113,114,116,113,115,112,114,115,116,117,118,119,120,122, 2124,121,123,120,122,123,124,125,126,128,130,127/ c..jcode tells the type progenitors each stable species can have. c..jcode = 0 if the species is the only stable one of that a c.. = 1 if the species can have proton-rich progenitors c.. = 2 if the species can have neutron-rich progenitors c.. = 3 if the species can only be made as itself (eg k40) data jcode / 0,0,2,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1 0,0,0,0,0,0,2,0,0,1,0,2,0,3,0,1,0,0,0,2,2,0,1,0,1,0,2,3,0, 2 1,0,0,2,0,1,0,0,2,0,1,0,0,0,2,0,0,1,0,0,0,2,0,0,1,0,0,2,2, 3 0,1,1,0,2,2,2,0,0,1,1,1,0,2,2,0,2,1,1,1,0,0,0,0,2,2,2,0,1, 4 1,0,3,0,2,2,1,1,0,1,0,2,2,0,1,1,0,2,2,2,0,0,1,1,1,0,2,2,2, 5 2,1,2,1,1,1,1,0,0,0,2,2,2,0,2,1,1,1,3,0,2,2,2,0/ c..sum; stuff residual into hydrogen if (ifirst .eq. 0) then ifirst = 1 sum = 0.0d0 do j=1,solsiz sum = sum + agsol(j) enddo sum = 1.0d0 - sum agsol(1) = agsol(1) + sum end if c..initialize; strait sweep andgrev = 1.0d30 z = 1.0d30 a = 1.0d30 if (len(nam) .lt. 5) stop 'nam < 5 characters in routine andgrev' do i=1,solsiz if ( namsol(i)(1:5) .eq. nam(1:5) ) then andgrev = agsol(i) z = float(izsol(i)) a = float(iasol(i)) c..now get the elmental mass fraction associated with this isotope xelem = 0.0d0 jbeg = max(1,i-12) jend = min(i+12,solsiz) do j=jbeg,jend if (izsol(j) .eq. z) xelem = xelem + agsol(j) enddo c..bail return end if enddo write(6,*) 'warning: no such entry ',nam(1:5) return end subroutine rchem implicit none save include 'lunits.dek' include 'parze.dek' include 'chem3.dek' c..reads the chem3 input deck c..declare integer i,j,ipos double precision taum c..communicate with the salpeter imf double precision salpmu common /imfcof/ salpmu c..formats 01 format(a) 02 format(1x,a,t26,'= ',a) 03 format(1x,a,t26,'=',1pe10.2) 04 format(1x,a,t26,'=',i3) 05 format(1x,a,' ',a,t26,'=',1pe10.2) 06 format(1x,a,' ',i2.2,t26,'=',1pe10.2) c..open and echo 10 write(6,*) 'give input deck =>' read(5,01) string call bbb(redold,fil2,string,i) if (i .eq. 1) goto 10 write(6,*) ' ' write(6,02) 'input file',string(1:index(string,' ')) c..read the input nucleosynthesis file read(fil2,01) string nukfil = string(1:index(string,' ')) c..read the output kinematics root file name read(fil2,01) string densfil = string(1:index(string,' ')) c..read the output supernova rate root file name read(fil2,01) string snfil = string(1:index(string,' ')) c..read the output mass fraction root file name read(fil2,01) string xmasfil = string(1:index(string,' ')) c..read the output mass ratio root file name read(fil2,01) string decfil = string(1:index(string,' ')) c..read the output solor comparison root file name read(fil2,01) string agfil = string(1:index(string,' ')) c..read the output g-dwarf root file name read(fil2,01) string gdwfil = string(1:index(string,' ')) c..read model parameters read(fil2,*) nu read(fil2,*) kk read(fil2,*) wimf read(fil2,*) salpmu read(fil2,*) radsol read(fil2,*) massol read(fil2,*) mcenter read(fil2,*) rbreak read(fil2,*) imfmin read(fil2,*) imfmax read(fil2,*) mmin read(fil2,*) mmax read(fil2,*) tend read(fil2,*) tgal read(fil2,*) atdisk read(fil2,*) btdisk read(fil2,*) tagout read(fil2,*) sigzer read(fil2,*) sn2mlo read(fil2,*) sn2mhi read(fil2,*) mbmin read(fil2,*) mbmax c..read the type1a supernova fractions read(fil2,*) sn1fac read(fil2,*) nfract1a if (nfract1a .gt. maxt1a) stop 'nfract1a > maxt1a in rchem' do i=1,nfract1a read(fil2,01) string ipos = 1 j = getnam(string,word,ipos) fract1a(i) = value(word) j = getnam(string,word,ipos) j = getnam(string,word,ipos) t1a_comment(i) = word enddo c..read the nova fractions read(fil2,*) novfac read(fil2,*) trecur read(fil2,*) nfracnova if (nfracnova .gt. maxnova) stop 'nfracnova > maxnova in rchem' do i=1,nfracnova read(fil2,01) string ipos = 1 j = getnam(string,word,ipos) fracnova(i) = value(word) j = getnam(string,word,ipos) j = getnam(string,word,ipos) nova_comment(i) = word enddo read(fil2,*) toler read(fil2,*) odescal read(fil2,*) stptry read(fil2,*) iprint read(fil2,*) nnhi read(fil2,*) nnmid read(fil2,*) nnlo if (nnhi .gt. nqmax) stop 'nnhi > nqmax in rchem' if (nnmid .gt. nqmax) stop 'nnmid > nqmax in rchem' if (nnlo .gt. nqmax) stop 'nnlo > nqmax in rchem' read(fil2,*) thalo read(fil2,*) nrcomp read(fil2,*) modal26 read(fil2,*) iznorm read(fil2,*) nzone if (nzone .gt. nzmax) stop 'nzone > nzmax in routine rchem' if (iznorm .gt. nzone) stop 'iznorm > nzone in routine rchem' nrcomp = min(nrcomp,nzone) do i=1,nzone read(fil2,*) rpoint(i) enddo c..close call bbb(clse,fil2,string,i) c..set some parameters that are based on the users input sn2mlo = max(sn2mlo,mmin) sn2mhi = min(sn2mhi,mmax) c..echo back the input for posterity write(6,*) ' ' write(6,02) 'nucleosynthesis file',nukfil(1:index(nukfil,' ')) write(6,02) 'kinematics file' ,densfil(1:index(densfil,' ')) write(6,02) 'supernova rate file' ,snfil(1:index(snfil,' ')) write(6,02) 'mass fraction file' ,xmasfil(1:index(xmasfil,' ')) write(6,02) 'element ratio file' ,decfil(1:index(decfil,' ')) write(6,02) 'a&g comparison file' ,agfil(1:index(agfil,' ')) write(6,02) 'g-dwarf file' ,gdwfil(1:index(gdwfil,' ')) write(6,*) ' ' write(6,03) 'birthrate constant nu',nu write(6,03) 'birthrate exponent kk',kk write(6,04) 'which imf to use',wimf write(6,03) 'salpeter exponent',salpmu write(6,03) 'distance to sun in pc',radsol write(6,03) 'mass/pc**2 at sun',massol write(6,03) 'mass/pc**2 at center',mcenter write(6,03) 'bulge-disk distance',rbreak write(6,03) 'imf min mass',imfmin write(6,03) 'imf max mass',imfmax write(6,03) 'min mass to evolve',mmin write(6,03) 'max mass to evolve',mmax write(6,03) 'time integration end',tend write(6,03) 'age of galaxy',tgal write(6,03) 'slope of timescale',atdisk write(6,03) 'intercept timescale',btdisk write(6,03) 'compare to a&g at',tagout write(6,03) 'initial mass/pc**2',sigzer write(6,03) 'sn2 rate min mass ',sn2mlo write(6,03) 'sn2 rate max mass',sn2mhi write(6,03) 'binary lo mass limit',mbmin write(6,03) 'binary hi mass limit',mbmax write(6,03) 'type1a factor',sn1fac write(6,04) 'different type1a types',nfract1a do i=1,nfract1a write(6,05) t1a_comment(i),'fraction',fract1a(i) enddo write(6,03) 'nova factor',novfac write(6,03) 'nova recur time',trecur write(6,04) 'different nova types',nfracnova do i=1,nfracnova write(6,05) nova_comment(i),'fraction',fracnova(i) enddo write(6,03) 'ode tolerance',toler write(6,03) 'ode scale',odescal write(6,03) 'first timestep to try',stptry write(6,04) 'iprint',iprint write(6,04) 'high mass quad points',nnhi write(6,04) 'mid mass quad points',nnmid write(6,04) 'lo mass quad points',nnlo write(6,03) 'switch infall comp',thalo write(6,04) 'zone to use new infall',nrcomp write(6,04) 'should we decay al26',modal26 write(6,04) 'zone for normalizations',iznorm write(6,04) 'number of radial zones',nzone do i=1,nzone write(6,06) 'radial zone',i,rpoint(i) enddo write(6,*) ' ' write(6,*) ' ' return end c.. c..some system and glue utility routines c.. c..routine getsym pulls a symbol from the operating system c..routine bbb opens and closed files c..routine today gets the date and clock time c..routine zsecond get the cpu time from the operating system c..routine lenstr finds the non-blank length of a string c..routine sqeeze compresses a string c..routine timlap converts total seconds into hours, minutes, seconds subroutine getsym(vblnam,vblval,lenval) implicit none save c..this routine gets a symbol from the operating system. c..input is the name of the symbol vblnam, and output is the value of the c..symbol vblval, of length lenval. c..declare character*(*) vblnam,vblval integer namlen,lenval,lent,lenstr c..initialize lenval = 0 lent = len(vblnam) namlen = lenstr(vblnam,lent) call getenv(vblnam(1:namlen),vblval) lent = len(vblval) lenval = lenstr(vblval,lent) return end subroutine bbb(id,lunit,luname,ierr) implicit none save c..this routine opens and closes files in various modes. c.. id function c.. 3 close file c.. 7 close with delete c.. 9 open old file c.. 10 open read/write unformatted new file c.. 11 open old unformatted file c.. 12 append old file c.. 13 open read/write new file c.. c..declare logical opened character*(*) luname integer id,lunit,ierr,i c..initialize ierr = 0 c..close a file if (id.eq.3) then if (lunit.ne.0) then inquire (lunit, opened=opened) if (opened) close (lunit) end if return c..close and delete the file else if (id.eq.7) then inquire (lunit, opened=opened) if (opened) close (lunit,status='delete') return c..open an old named file for reading, 1056 wide else if (id.eq.9) then i=index(luname,' ') -1 open(unit=lunit,file=luname(1:i),err=100,status='old') rewind lunit return c..open a binary file for reading and writing, 1056 wide else if (id.eq.10) then i=index(luname,' ') -1 open(unit=lunit,file=luname(1:i),form='unformatted', 1 err=100,status='unknown') rewind (lunit) return c..open an old binary file for reading, 1056 wide else if (id.eq.11) then i=index(luname,' ') -1 open(unit=lunit,file=luname(1:i),form='unformatted', 1 err=100,status='old') rewind lunit return c..open old files for writing (append), 1056 wide c..append not supported in fortran 90 c.. else if (id.eq.12) then c.. i=index(luname,' ') -1 c..cxcl vax c.. open(unit=lunit,file=luname(1:i), c.. 1 err=100, status='old', access='append') c..cbeg vax c.. open(unit=lunit,file=luname(1:i),recl=1056, c.. 1 err=100,status='old', access='append') c..cend c.. rewind (lunit) c.. do 1200 i=1,30000 c.. read(lunit,1199,end=1201) aline c..1199 format(a80) c..1200 continue c..1201 continue c.. return c..open a new file for reading and writing, 1056 wide else if (id.eq.13) then i=index(luname,' ') - 1 open(unit=lunit,file=luname(1:i),err=100,status='unknown') rewind (lunit) return end if c..error with the file 100 write(6,101) luname(1:20) 101 format(1x,'* error with file >',a,'<') ierr = 1 return end subroutine today(adat,atim) implicit none save c..this routine gets the date and time out of a machine. c..the output format is adat ddmmmyyyy c.. atim hh:mm:ss c..declare character*3 amon(12) character*8 atim character*9 adat integer idat(3),itim(3),idum data amon/ 'jan' , 'feb' , 'mar' , 'apr' , 'may' , 'jun' , 1 'jul' , 'aug' , 'sep' , 'oct' , 'nov' , 'dec' / c..popular format statements 113 format(i2.2,':',i2.2,':',i2.2) 114 format(i2.2,a3,i4.4) c..initialize adat=' ' atim=' ' call itime(itim) write(atim,113) itim call idate(idat) write(adat,114) idat(1),amon(idat(2)),idat(3) c..some older sgi boxes work off this scheme c.. call itime(itim) c.. write(atim,113) itim c.. call idate(idat(1),idat(2),idat(3)) c.. write(adat,114) idat(2),amon(idat(1)),idat(3) return end subroutine zsecond(time) implicit none save c.. c..this routine gets the elapsed time of a job from the machine c.. c..declare double precision time real tarray(2) call etime(tarray) time = tarray(1) + tarray(2) return end integer function lenstr(string,istrln) implicit none save c..lenstr returns the non blank length length of the string. c..declare integer istrln,i character*(*) string lenstr=0 do i=istrln,1,-1 if (string(i:i).ne. ' ') then if (ichar(string(i:i)).ne. 0 )then lenstr=i goto 20 end if end if enddo 20 return end subroutine sqeeze(line) implicit none save c..this routine takes line and removes all blanks, such as c..those from writing to string with fortran format statements c..declare character*(*) line character*1 achar integer l,n,k,lend,lsiz,lenstr c..find the end of the line lsiz = len(line) lend = lenstr(line,lsiz) n = 0 l = 0 c..do the compression in place 10 continue l = l + 1 achar = line(l:l) if (achar .eq. ' ') goto 10 n = n + 1 line(n:n) = achar if (l .lt. lend) goto 10 c..blank the rest of the line do k=n+1,lsiz line(k:k) = ' ' enddo return end subroutine timlap(tlap,hours,minut,sec,msec) implicit none save c..this routines converts seconds to hours, minutes, seconds and microseconds c.. c..declare integer hours,minut,sec,msec double precision tlap,x msec = 0 sec = 0 minut = 0 hours = 0 sec = int(tlap) msec = 1.0d6 * (tlap-sec) if (sec .ge. 60) then x = dble(sec)/60.0d0 minut = int(x) end if sec = sec - minut*60 if (minut .ge. 60) then x = dble(minut)/60.0d0 hours = int(x) end if minut = minut - hours*60 return end integer function getnam(string,word,ipos) implicit none save c..this routine is a basic string parser, only spaces and equal signs are c..used as delimiters. c.. c..input is the string to chop and where to start the parse from ipos. c.. c..output is the parsed word and the updated position on the string ipos. c..the length of the word is returned as getnam. if no word is found, c..meaning that there are no more words in the string, then getnam is c..returned as zero and the word is filled with spaces. c..declare character*(*) string,word integer kpos,ipos,k,lend,nbegin c..get the length of the input line, blank word and save ipos lend = len(string) word = ' ' kpos = ipos c..find begining of the word; if nothing found return getnam as zero do k = kpos,lend nbegin = k if (string(k:k).ne.' ' .and. string(k:k).ne.'=') goto 25 enddo getnam = 0 return c..find end of the word 25 continue do k = nbegin,lend ipos = k if( string(k:k).eq.' ' .or. string(k:k).eq.'=') goto 35 enddo c..build the word, set getnam and return 35 continue word = string(nbegin:ipos-1) getnam = ipos - nbegin return end double precision function value(string) implicit none save c..this routine takes a character string and converts it to a real number. c..on error during the conversion, a fortran stop is issued c..declare logical pflag character*(*) string character*1 plus,minus,decmal,blank,se,sd,se1,sd1 integer noblnk,long,ipoint,power,psign,iten,j,z,i double precision x,sign,factor,rten,temp parameter (plus = '+' , minus = '-' , decmal = '.' , 1 blank = ' ' , se = 'e' , sd = 'd' , 2 se1 = 'E' , sd1 = 'D' , rten = 10.0, 3 iten = 10 ) c..initialize x = 0.0d0 sign = 1.0d0 factor = rten pflag = .false. noblnk = 0 power = 0 psign = 1 long = len(string) c..remove any leading blanks and get the sign of the number do z = 1,7 noblnk = noblnk + 1 if ( string(noblnk:noblnk) .eq. blank) then if (noblnk .gt. 6 ) goto 30 else if (string(noblnk:noblnk) .eq. plus) then noblnk = noblnk + 1 else if (string(noblnk:noblnk) .eq. minus) then noblnk = noblnk + 1 sign = -1.0d0 end if goto 10 end if enddo c..main number conversion loop 10 continue do i = noblnk,long ipoint = i + 1 c..if a blank character then we are done if ( string(i:i) .eq. blank ) then x = x * sign value = x return c..if an exponent character, process the whole exponent, and return else if (string(i:i).eq.se .or. string(i:i).eq.sd .or. 1 string(i:i).eq.se1 .or. string(i:i).eq.sd1 ) then if (x .eq. 0.0 .and. ipoint.eq.2) x = 1.0d0 if (sign .eq. -1.0 .and. ipoint.eq.3) x = 1.0d0 if (string(ipoint:ipoint) .eq. plus) ipoint = ipoint + 1 if (string(ipoint:ipoint) .eq. minus) then ipoint = ipoint + 1 psign = -1 end if do z = ipoint,long if (string(z:z) .eq. blank) then x = sign * x * rten**(power*psign) value = x return else j = ichar(string(z:z)) - 48 if ( (j.lt.0) .or. (j.gt.9) ) goto 30 power= (power * iten) + j end if enddo c..if an ascii number character, process ie else if (string(i:i) .ne. decmal) then j = ichar(string(i:i)) - 48 if ( (j.lt.0) .or. (j.gt.9) ) goto 30 if (.not.(pflag) ) then x = (x*rten) + j else temp = j x = x + (temp/factor) factor = factor * rten goto 20 end if c..must be a decimal point if none of the above c..check that there are not two decimal points! else if (pflag) goto 30 pflag = .true. end if 20 continue end do c..if we got through the do loop ok, then we must be done x = x * sign value = x return c..error processing the number 30 write(6,40) long,string(1:long) 40 format(' error converting the ',i4,' characters ',/, 1 ' >',a,'< ',/, 2 ' into a real number in function value') stop ' error in routine value' end c..quadrature: c..routine qgaus does the gaussian quadrature summation c..routine gauleg finds the abscissas & weights for an n point gauss-legendre c..routine qgaus2 same as above, used to avoid recurrsive routines c..routine gauleg2 same as above, used to avoid recurrsive routines subroutine qgaus(func,x,w,n,ss) implicit none save c..returns as ss the quadrature summation of the function func as determined c..by the abcissas x and weights w. c.. c..declare external func integer i,n double precision func,ss,x(n),w(n) ss = 0.0d0 do i=1,n ss = ss + w(i)*func(x(i)) enddo return end subroutine gauleg(x1,x2,x,w,n) implicit none save c..given the lower and upper limits of integration x1 and x2, and given n, c..this routine returns arrays x and w of length n, containing the c..abscissas and weights of the gauss-legendre n-point quadrature formula c..declare integer i,m,j,n double precision x1,x2,x(n),w(n),eps,xm,xl,p1,p2,p3,pp,z,z1 parameter (eps=1.0e-14) c..roots are symmetric in the interval so we only have to find half of them m = (n+1)/2 xm = 0.5d0 * (x2 + x1) xl = 0.5d0 * (x2 - x1) c..loop over the desired roots and make a slick guess at each one do i=1,m z = cos(3.141592653589d0 * (i-0.25d0)/(n + 0.5d0)) c..newton loop 1 continue p1 = 1.0d0 p2 = 0.0d0 c..loop the recurrence relation to get the legendre polynomial at z do j=1,n p3 = p2 p2 = p1 p1 = ((2.0d0 * j - 1.0d0) * z * p2 - (j - 1.0d0)*p3)/j enddo c..p1 is now the desired legendre polynomial. pp is the derivative. pp = n * (z*p1 - p2)/(z*z - 1.0d0) z1 = z z = z1 - p1/pp if (abs(z-z1) .gt. eps) goto 1 c..scale to the users interval x(i) = xm - xl*z x(n+1-i) = xm + xl * z w(i) = 2.0d0 * xl/((1.0d0 - z*z)*pp*pp) w(n+1-i) = w(i) enddo return end subroutine qgaus2(func,x,w,n,ss) implicit none save c..returns as ss the quadrature summation of the function func as determined c..by the abcissas x and weights w. c.. c..declare external func integer i,n double precision func,ss,x(n),w(n) ss = 0.0d0 do i=1,n ss = ss + w(i)*func(x(i)) enddo return end subroutine gauleg2(x1,x2,x,w,n) implicit none save c..given the lower and upper limits of integration x1 and x2, and given n, c..this routine returns arrays x and w of length n, containing the c..abscissas and weights of the gauss-legendre n-point quadrature formula c..declare integer i,m,j,n double precision x1,x2,x(n),w(n),eps,xm,xl,p1,p2,p3,pp,z,z1 parameter (eps=1.0e-14) c..roots are symmetric in the interval so we only have to find half of them m = (n+1)/2 xm = 0.5d0 * (x2 + x1) xl = 0.5d0 * (x2 - x1) c..loop over the desired roots and make a slick guess at each one do i=1,m z = cos(3.141592653589d0 * (i-0.25d0)/(n + 0.5d0)) c..newton loop 1 continue p1 = 1.0d0 p2 = 0.0d0 c..loop the recurrence relation to get the legendre polynomial at z do j=1,n p3 = p2 p2 = p1 p1 = ((2.0d0 * j - 1.0d0) * z * p2 - (j - 1.0d0)*p3)/j enddo c..p1 is now the desired legendre polynomial. pp is the derivative. pp = n * (z*p1 - p2)/(z*z - 1.0d0) z1 = z z = z1 - p1/pp if (abs(z-z1) .gt. eps) goto 1 c..scale to the users interval x(i) = xm - xl*z x(n+1-i) = xm + xl * z w(i) = 2.0d0 * xl/((1.0d0 - z*z)*pp*pp) w(n+1-i) = w(i) enddo return end subroutine locate(xx,n,x,j) implicit none save c..given an array xx of length n, and a value of x, this routine returns c..a value j such that x is between xx(j) and xx(j+1). the array xx must be c..monotonic. j=0 or j=n indicates that x is out of range. bisection is used c..to find the entry c..declare integer n,j,jl,ju,jm double precision xx(n),x c..initialize jl = 0 ju = n+1 c..compute a midpoint, and replace either the upper or lower limit 10 if (ju-jl .gt. 1) then jm = (ju+jl)/2 if ( (xx(n) .ge. xx(1)) .eqv. (x .ge. xx(jm)) ) then jl = jm else ju = jm end if goto 10 end if if (x .eq. xx(1))then j = 1 else if(x .eq. xx(n))then j = n - 1 else j = jl end if return end subroutine hunt(xx,n,x,jlo) implicit none save c..given an array xx of length n and a value x, this routine returns a value c..jlo such that x is between xx(jlo) and xx(jlo+1). the array xx must be c..monotonic. j=0 or j=n indicates that x is out of range. on input jlo is a c..guess for the table entry, and should be used as input for the next hunt. c..one assumes the next table entry is relatively close to the old one. c..declare logical ascnd integer n,jlo,jhi,inc,jm double precision xx(n),x c..logical function; true if ascending table, false if descending ascnd = xx(n) .ge. xx(1) c..horrid initial guess; go to bisection immediatly if (jlo.le.0 .or. jlo.gt.n) then jlo = 0 jhi = n+1 goto 3 end if c..set the initial hunting increment inc = 1 c..this is the hunt up section if (x.ge.xx(jlo) .eqv. ascnd) then 1 jhi = jlo + inc c..if we are the end of the table then we are done hunting if (jhi .gt. n) then jhi = n + 1 c..otherwise replace the lower limit, double the increment and try again else if (x.ge.xx(jhi) .eqv. ascnd) then jlo = jhi inc = inc + inc goto 1 end if c..this is the hunt down section else jhi = jlo 2 jlo = jhi - inc c..if we are the end of the table then we are done hunting if (jlo.lt.1) then jlo = 0 c..otherwise replace the upper limit, double the increment and try again else if (x.lt.xx(jlo) .eqv. ascnd) then jhi = jlo inc = inc + inc goto 2 end if end if c..this is the final bisection phase 3 if (jhi-jlo .eq. 1) then if (x .eq. xx(n)) jlo = n - 1 if (x .eq. xx(1)) jlo = 1 return end if jm = (jhi+jlo)/2 if (x .ge. xx(jm) .eqv. ascnd) then jlo = jm else jhi = jm end if goto 3 end c..routine spline computes the natural cubic spline coefficients c..routine splint uses the spline coefficients and bisection to interpolate c..routine splinth uses the spline coefficients and hunting to interpolate c..routine splinth1 returns the spline value and first derivative c subroutine spline(x,y,n,yp1,ypn,y2) subroutine spline(x,y,n,y2) implicit none save c..given arrays x(n) and y(n) containing a tabulated function y=f(x), c..this routine returns the array y2(n) of spline coefficients. c..x must be a monotonically increasing c.. c..if yp1 and/or ypn are set larger than 1e30, the routine sets the boundary c..condtion for a natural spline (zero second derivative) at that boundary, c..otherwise yp1 and ypn are the first derivatives to use at the endpoints c.. c..this routine is only called once for any given x and y. c.. c..declare integer n,i,k,nmax parameter (nmax=1000) double precision x(n),y(n),y2(n),yp1,ypn,u(nmax),sig,p,qn,un parameter (yp1 = 2.0d30, ypn=2.0d30) c..the lower boundary condition is set to be either natural if (yp1 .gt. 1.0d30) then y2(1) = 0.0d0 u(1) = 0.0d0 c..or it has a specified first derivative else y2(1) = -0.5d0 u(1) = (3.0d0/(x(2)-x(1))) * ((y(2)-y(1))/(x(2)-x(1)) - yp1) end if c..this is the decomposition loop of the tridiagonal algorithm. y2 and u c..are used for temporary storage of the decomposition factors. do i=2,n-1 sig = (x(i) - x(i-1))/(x(i+1) - x(i-1)) p = sig*y2(i-1) + 2.0d0 y2(i) = (sig - 1.0d0)/p u(i) = (6.0d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) / 1 (x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p enddo c..the upper boundary condition is set to be either natural if (ypn .gt. 1.0d30) then qn = 0.0d0 un = 0.0d0 c..or it has a specified first derivative else qn = 0.5d0 un = (3.0d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) end if c..this is the backsubstitution loop of the tridiagonal algorithm y2(n) = (un-qn*u(n-1)) / (qn*y2(n-1) + 1.0d0) do k=n-1,1,-1 y2(k) = y2(k)*y2(k+1) + u(k) enddo return end subroutine splint(xa,ya,y2a,n,x,y) implicit none save c..given arrays xa and ya of length n, which tablulate a monotonic function, c..and given the array y2a, which is the output of spline (above), and given c..a value of x this routine returns a cubic spline interpolated value of y. c..declare integer n,klo,khi,k double precision xa(n),ya(n),y2a(n),x,y,h,a,b c..find the right place in the table by bisection. this is