!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
