!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! ! Copyright D.J. Jeffery, 2006jan01. ! ! This is a set of subroutines for reading in supernova models. ! They come in all different formats. ! ! References: ! ! \bibitem[Metcalf et al.(2004)]{metcalf2004} Metcalf, M., Reid, J., ! \&~Cohen, M. 2004, ! Fortran 95/2003 Explained ! (Oxford: Oxford University Press) (MRC) ! % Pretty good reference book. ! ! \bibitem[Nomoto et al.(1984)]{nomoto1984} ! Nomoto, K., Thielemann, F.-K., \&~Yokoi, \apj, 286, 644 ! % The original W7 paper. ! % Here the iron peak yield is .86 solar masses ! ! \bibitem[Press et al.(1992)]{press1992} Press, W. H., Teukolsky, S. A., ! Vetterling, W. T., \& Flannery, B. P. 1992, ! Numerical Recipes in Fortran ! (Cambridge: Cambridge University Press) ! % Bill should have given me a complimentary copy when ! % I co-taught with him in 1993! ! ! \bibitem[Thielemann et al.(1986)]{thielemann1986} Thielemann, F.-K., ! Nomoto, K., ! \&~Yokoi, A\&A, 158, 17 ! % The final W7 composition paper. ! % I sum the iron-peak to be .79 solar masses. ! % This is less than the old W7 which had .86 solar masses ! % of iron-peak ! ! \bibitem[Yoon \&~Langer(2005)]{yoon2005} ! Yoon, S.-C., \&~Langer, N. 2005, A\&A, 435, 985 ! % On the evolution of rapidly rotating massive white dwarfs toward ! % supernovae or collapses. ! % On p. 977 they give a figure of binding energy for ! % massive rotating WDs. ! % They also give fitting formula for those binding energies. ! !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! module model_mod use numerical use physical use astronomical use atomic_mod ! The module is in /aalib/constant.f implicit none ! integer, parameter :: ncell=400 real (kind=nprecision) :: abund(ncell,natom)=0.d0 integer :: iatom ! Model upper limit on atomic number. integer :: icell ! Number of cells in model. integer :: iradio integer :: ivel_center=0 ! 0 for none (default); 1 for simple averaging; ! ! 2 for tricky, but it's about the same as 1. integer :: ivel_kms=0 ! 0 for no conversion; 1 for conversion to km/s. real (kind=nprecision) :: radio(ncell,nradio_isobar,nradio)=0.d0 real (kind=nprecision) :: rho(ncell) real (kind=nprecision) :: rhoc real (kind=nprecision) :: time ! Time of model in seconds. real (kind=nprecision) :: vefold real (kind=nprecision) :: vel(ncell) real (kind=nprecision) :: xke_erg real (kind=nprecision) :: xmass_cell(ncell) real (kind=nprecision) :: xmass_cell_del(ncell) real (kind=nprecision) :: xmass_gram ! end module model_mod ! !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! ! ! subroutine model_diagnostic use numerical use physical use astronomical use atomic_mod ! The module is in /aalib/constant.f use model_mod implicit none ! integer :: i,j,k,l,m,n real (kind=nprecision) :: vel_ipe real (kind=nprecision) :: xmass_ipe real (kind=nprecision) :: xmass_nico ! !----------------------------------------------------------------------- write(*,*) write(*,*) 'Getting the exponential fitted vefold and rhoc.' !----------------------------------------------------------------------- ! vefold=sqrt(sixth*xke_erg/xmass_gram) ! Exponential model fit ! ! for vefold. rhoc=xmass_gram/(8.d0*pi*(vefold*daysec)**3) ! Exponential model ! ! fit for rhoc. write(*,*) 'Fitted vefold and rhoc for 1 day.' write(*,*) 'vefold,rhoc,xmass_gram,xke_erg' write(*,'(1p,4e25.15)') vefold,rhoc,xmass_gram,xke_erg ! !----------------------------------------------------------------------- write(*,*) write(*,*) 'Checking the Ni/Co-56, IPE mass, and v_ipe.' !----------------------------------------------------------------------- ! ! Sum up the Ni/Co-56 mass. xmass_nico=sum( ( radio(:icell,1,1)+radio(:icell,2,1) ) & & *xmass_cell_del(:icell) ) ! ! Sum up the IPE mass. xmass_ipe=0.d0 do410 : do i=1,icell xmass_ipe=xmass_ipe+sum(abund(i,21:32))*xmass_cell_del(i) end do do410 ! ! Determine the fiducial IPE mass core. ! Usually the velocities and masses are both at the cell walls I think. ! call interpolation(icell,xmass_cell(:icell), & & vel(:icell),1,1,xmass_ipe,k,vel_ipe) ! write(*,*) 'Ni/Co56 mass,IPE mass,g-factor,v_ipe,vel(k)' write(*,*) xmass_nico,xmass_ipe,xmass_nico/xmass_ipe, & & vel_ipe,vel(k) write(*,*) 'Si,S,IPE,Ni-56 at vel(k)' write(*,*) abund(k,14),abund(k,16),sum(abund(k,21:28)), & & radio(k,1,1),vel(k) write(*,*) 'Si,S,IPE at vel(k+1)' write(*,*) abund(k+1,14),abund(k+1,16),sum(abund(k+1,21:28)), & radio(k+1,1,1),vel(k+1) ! end subroutine model_diagnostic ! ! !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! ! Subroutine bravo is for reading the model files of Bravo & Garcia-Senz. ! subroutine bravo(afile) use numerical use physical use astronomical use atomic_mod use model_mod implicit none ! integer, parameter :: nbravo=15 ! character, intent(in) :: afile*(*) ! integer :: i,j,k,l,m,n integer :: ibravo(nbravo)=(/ 2, 6, 8,10, 12,14,16,18, & & 20,24,25,26, 27,28,30/) real (kind=nprecision) :: tmp real (kind=nprecision) :: tmp1,tmp2,tmp3 real (kind=nprecision) :: vel2(ncell) ! namelist/param/iatom,time ! open(unit=1,file='/home/jeffery/public_html/aalib/atomic.dat', & & status='old',action='read') read(1,atomic_data) close(unit=1) ! !----------------------------------------------------------------------- write(*,*) write(*,*) 'bravo: Before zeroing and reading in.' !----------------------------------------------------------------------- ! abund(:,:)=0.d0 radio(:,:,:)=0.d0 xke_erg=0.d0 xmass_cell_del(:)=0.d0 ! open(unit=1,file=afile,status='old',action='read') ! Open bravo data file. read(1,param) call readcopy(1,0,0,' ',icell) write(*,*) write(*,"('icell',3x,'vel (cm/s)',4x,'rho (cgs)', & & 2x,'m_i (m_sun)',5x,'ratio', & & 7x,'sum',3x,'X Ni-56')") do410 : do i=1,icell read(1,*) j,tmp,xmass_cell(i),vel(i), & & (abund(i,ibravo(j)),j=1,nbravo), & & (radio(i,1,j),radio(i,2,j),radio(i,3,j),j=1,2), & & radio(i,1,4),radio(i,2,4),radio(i,3,4), & & radio(i,2,3),radio(1,3,3) ! The 1st nuclide Co-55 is ! ! not in Bravo. ! ! -------------------------------------------------------------- ! Note in radio(i,j,k) that the last index is the decay chain and ! the 2nd to last is the radionuclide in the decay chain. ! Bravo doesn't have Co-55 abundances in his table. ! -------------------------------------------------------------- ! abund(i,28)=abund(i,28)+radio(i,1,1)+radio(i,1,2) abund(i,27)=abund(i,27)+radio(i,2,1)+radio(i,2,2) abund(i,26)=abund(i,26)+radio(i,3,1)+radio(i,3,2) ! abund(i,26)=abund(i,26)+radio(i,2,3) abund(i,25)=abund(i,25)+radio(i,3,3) ! abund(i,22)=abund(i,22)+radio(i,1,4) abund(i,21)=abund(i,21)+radio(i,2,4) abund(i,20)=abund(i,20)+radio(i,3,4) end do do410 ! ! The do420 loop is needed until I get the zone-edge velocities ! from Eduardo Bravo. ! ! vel2(:)=vel(:) ! do420 : do i=1,icell-1 ! vel(i)=.5d0*(vel(i)+vel(i+1)) ! end do do420 ! call interpolation(icell-1,xmass_cell,vel,1,1, & ! & xmass_cell(icell),i,vel(icell)) ! print*,'vel(icell)' ! print*,vel(icell) ! Gives 2.4318e+9 for PRD5.5. ! vel(icell)=tmp/time ! Gives a way too low velocity for PRD5.5: ! print*,vel(icell) ! i.e., 1.7696e+9. ! do430 : do i=1,icell ! if(i .gt. 1) then xmass_cell_del(i)=xmass_cell(i)-xmass_cell(i-1) tmp1=vel(i-1) tmp2=vel(i-1)**3 else xmass_cell_del(i)=xmass_cell(i) tmp1=0.d0 tmp2=0.d0 end if rho(i)=(xmass_cell_del(i)*xmass_sun) & & /(pi4L3*(vel(i)**3-tmp2)*daysec**3) xke_erg=xke_erg+.5d0*( .5d0*(vel(i)+tmp1) )**2 & & *xmass_cell_del(i)*xmass_sun tmp3=sum(abund(i,:iatom)) write(*,910) i,vel(i),rho(i), & & xmass_cell(i),tmp/(vel(i)*time),tmp3, & & abund(i,28) 910 format(i5,1p,3e13.4,0p,3f10.4) end do do430 close(unit=1) ! xmass_gram=xmass_cell(icell)*xmass_sun ! call model_diagnostic if(ivel_center .eq. 1) then call velocity_centering ! vel(:)=vel2(:) end if call velocity_cms_to_kms ! write(*,*) i=1 write(*,*) 'i,iatom,abund(i,8),awt(8)' write(*,*) i,iatom,abund(i,8),awt(8) ! end subroutine bravo ! !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! ! Subroutine hoeflich is for reading the old model files of Peter Hoeflich. ! subroutine hoeflich(afile) use numerical use physical use astronomical use atomic_mod use model_mod implicit none ! integer, parameter :: nnatom=30 ! character, intent(in) :: afile*(*) ! integer :: i,j,k,l,m,n integer :: iiatom(nnatom) character :: record*180 real (kind=nprecision) :: tmp real (kind=nprecision) :: tmp1,tmp2,tmp3 real (kind=nprecision) :: xmass_mult ! namelist/param/iatom,iiatom,icell,rhoc,time,xmass_mult ! open(unit=1,file='/home/jeffery/public_html/aalib/atomic.dat', & & status='old',action='read') read(1,atomic_data) close(unit=1) ! !----------------------------------------------------------------------- write(*,*) write(*,*) 'hoeflich: Before zeroing and reading in.' !----------------------------------------------------------------------- ! abund(:,:)=0.d0 radio(:,:,:)=0.d0 xke_erg=0.d0 xmass_cell_del(:)=0.d0 ! open(unit=1,file=afile,status='old',action='read') ! Open bravo data file. read(1,param) call readcopy(1,0,0,' ',k) read(1,*) read(1,*) record write(*,*) trim(record) write(*,"('icell',3x,'vel (cm/s)',4x,'rho (cgs)', & & 2x,'m_i (m_sun)',5x,'ratio', & & 7x,'sum',3x,'X Ni-56')") do410 : do i=1,icell read(1,*) j,xmass_cell(i),rho(i),vel(i),tmp2,tmp3,radio(i,1,1) end do do410 do420 : do i=1,11 read(1,*) end do do420 do430 : do i=1,icell read(1,*) j write(*,*) i,j read(1,*) (abund(i,iiatom(j)),j=1,7) read(1,*) (abund(i,iiatom(j)),j=8,14) read(1,*) (abund(i,iiatom(j)),j=15,21) abund(i,28)=abund(i,28)+radio(i,1,1) abund(i,26)=abund(i,26)-radio(i,1,1) write(*,*) i,abund(i,26),abund(i,28),radio(i,1,1) end do do430 ! ! Note in radio(i,j,k) that the last index is the decay chain, the 2nd to ! last is the radionuclide in the decay chain. ! rhoc=rhoc*(time/daysec)**3 rho(:)=rho(:)*rhoc xmass_cell(:)=xmass_cell(:)*xmass_mult ! do440 : do i=1,icell if(i .gt. 1) then xmass_cell_del(i)=xmass_cell(i)-xmass_cell(i-1) tmp1=vel(i-1) tmp2=vel(i-1)**3 else xmass_cell_del(i)=xmass_cell(i) tmp1=0.d0 tmp2=0.d0 end if tmp=(xmass_cell_del(i)*xmass_sun) & & /(pi4L3*(vel(i)**3-tmp2)*daysec**3) xke_erg=xke_erg+.5d0*( .5d0*(vel(i)+tmp1) )**2 & & *xmass_cell_del(i)*xmass_sun tmp3=sum(abund(i,:iatom)) write(*,910) i,vel(i),rho(i), & & xmass_cell(i),tmp/rho(i),tmp3, & & abund(i,28) 910 format(i5,1p,3e13.4,0p,3f10.4) end do do440 close(unit=1) ! xmass_gram=xmass_cell(icell)*xmass_sun ! call model_diagnostic call velocity_centering call velocity_cms_to_kms ! i=1 write(*,*) 'i,iatom,abund(i,8),awt(8)' write(*,*) i,iatom,abund(i,8),awt(8) ! end subroutine hoeflich ! !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! ! Subroutine w7 is for reading the w7.dat file. Not finished. ! subroutine w7(afile) use numerical use physical use astronomical use atomic_mod use model_mod implicit none ! character, intent(in) :: afile*(*) real (kind=nprecision) :: amass(nnuclide) real (kind=nprecision) :: frac(nnuclide) integer :: i,j,k,l,m,n integer :: inuclide(nnuclide) character :: record*80 real (kind=nprecision) :: rho_cell real (kind=nprecision) :: sum1,sum2,sum3 real (kind=nprecision) :: tmp real (kind=nprecision) :: tmp1,tmp2,tmp3 ! real (kind=nprecision) :: acarbon=.4875 real (kind=nprecision) :: aoxygen=.4875 real (kind=nprecision) :: aneon=.025 integer :: icell_comp integer :: iinuclide integer :: itmp real (kind=nprecision) :: xmass namelist/param/iatom,icell,icell_comp,iinuclide, & & iradio,time,xmass ! character :: anuclide(nnuclide)*4 real (kind=nprecision) :: radius(ncell) ! namelist/runs/anuclide,radius,rho,vel,xmass_cell ! open(unit=1,file='/home/jeffery/public_html/aalib/atomic.dat', & & status='old',action='read') read(1,atomic_data) close(unit=1) ! open(unit=1,file=afile,status='old',action='read') ! Open w7 data file. read(1,param) read(1,runs) ! ! Dealing with the first four nuclides. inuclide(1:4)=1 ! All hydrogen or they will be. amass(1:2)=1.d0 amass(3)=2.d0 amass(4)=3.d0 ! !----------------------------------------------------------------------- write(*,*) write(*,*) 'w7: Before identifiying the nuclides with '// & & 'their atoms.' !----------------------------------------------------------------------- ! write(*,"('nuclide',' Atomic No.',' Atomic mass')") do402 : do i=5,nnuclide ! The first four a N,P,D,T: dealt with above. do404 : do j=1,natom if(anuclide(i)(1:2) .eq. aname2(j)) then inuclide(i)=j read(anuclide(i)(3:4),'(f3.0)') amass(i) end if end do do404 write(*,'(i7,i11,f11.0)') i,inuclide(i),amass(i) end do do402 ! !----------------------------------------------------------------------- write(*,*) write(*,*) 'Before converting density and masses.' !----------------------------------------------------------------------- ! ! Convert to density at day 1. rho(:)=10.d0**rho(:)*(time/daysec)**3 ! Conversion to solar masses. Masses must be in log of fractional mass. xmass_cell(:icell)=( 10.d0**xmass_cell(:icell) )*xmass xmass_gram=0.d0 ! Mass from the mass in the cells. xke_erg=0.d0 ! Kinetic energy. abund(:,:)=0.d0 frac(:)=0.d0 radio(:,:,:)=0.d0 ! !----------------------------------------------------------------------- write(*,*) write(*,*) 'Before calculating mass and kinetic energy and' write(*,*) 'doing some consistency checks which passed so-so.' write(*,*) write(*,"(' zone',16x,'rho given',11x,'rho calculated', & & 1x,'rho ratio',1x,'cell mass')") ! !----------------------------------------------------------------------- ! do410 : do i=1,icell if(i .eq. 1) then tmp=0.d0 tmp1=0.d0 tmp2=0.d0 xmass_cell_del(i)=xmass_cell(i) else tmp=vel(i-1)**3 tmp1=xmass_cell(i-1) tmp2=vel(i-1) xmass_cell_del(i)=xmass_cell(i)-xmass_cell(i-1) end if ! rho_cell=xmass_cell_del(i)*xmass_sun & & /(pi4L3*(vel(i)**3-tmp)*daysec**3) write(*,'(i5,1p,2e25.15,0p,2f10.5)') i,rho(i),rho_cell, & & (rho(i)/rho_cell), & & xmass_cell(i) xmass_gram=xmass_gram+rho(i)*(pi4L3*(vel(i)**3-tmp)*daysec**3) ! xke_erg=xke_erg+.5d0*( .5d0*(vel(i)+tmp2) )**2 & & *rho(i)*(pi4L3*(vel(i)**3-tmp)*daysec**3) ! if(i .le. icell_comp) then ! w7 abundances only given up to icell_comp 110 continue read(1,'(a80)') record ! write(*,*) record if(record(2:4) .ne. 'ZO=') go to 110 do412 : do j=5,iinuclide+1,5 ! Reading the nuclide fractions. read(1,'(2x,5(4x,e10.3))') & & (frac(k),k=j-4,min(j,iinuclide)) ! reading the nuclide fractions end do do412 sum3=sum(frac(:)*amass(:)) frac(:)=frac(:)*amass(:)/sum3 do414 : do j=1,nnuclide abund(i,inuclide(j))=abund(i,inuclide(j))+frac(j) end do do414 radio(i,1,1)=frac(204) ! Ni-56 parent. radio(i,2,1)=frac(195) ! Co-56 daughter. else ! Set outer layer abundances. abund(i,6)=acarbon abund(i,8)=aoxygen abund(i,10)=aneon radio(i,1,1)=0.d0 radio(i,2,1)=0.d0 end if ! end do do410 ! close(unit=1) ! Close w7 data file. ! !----------------------------------------------------------------------- write(*,*) write(*,*) 'Checking the mass consistency and time' !----------------------------------------------------------------------- ! sum1=xmass_gram/xmass_sun tmp=(xmass/sum1)**third*time write(*,*) 'xmass,xmass cal.,time,time cal.,kinetic energy' write(*,'(2f25.15,1p,3e25.15)') xmass,sum1,time,tmp,xke_erg ! xmass_gram=xmass*xmass_sun ! This is just for w7 accuracy. ! call model_diagnostic call velocity_centering call velocity_cms_to_kms ! end subroutine w7 ! !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! ! Conversion to velocity centered in mass cell. ! ! The simple average agrees with the more elaborate result to ! better than 1 % almost always for any reasonably well resolved model ! based on one example test. But it seems that this should be generally ! true. One can always test if one is unsure. ! subroutine velocity_centering use numerical use physical use astronomical use atomic_mod use model_mod implicit none ! integer :: i,j,k,l,m,n real (kind=nprecision) :: tmp,tmp1 ! if(ivel_center .eq. 0) go to 200 ! tmp=0.d0 do410 : do i=1,icell tmp1=vel(i) if(ivel_center .eq. 1) then vel(i)=.5d0*(tmp+tmp1) ! print*,'vel(i)' ! print*,vel(i) else vel(i)=sqrt(xmass_cell_del(i)*xmass_sun & & /(pi4*rho(i)*daysec**3*(tmp1-tmp))) ! print*,vel(i) end if tmp=tmp1 end do do410 ! 200 continue ! end subroutine velocity_centering !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! ! Conversion of velocity from cm/s to km/s. ! subroutine velocity_cms_to_kms use numerical use physical use astronomical use atomic_mod use model_mod implicit none ! if(ivel_kms .eq. 0) return ! vel(:)=vel(:)*1.d-5 vefold=vefold*1.d-5 ! end subroutine velocity_cms_to_kms