instar.f

Go to the documentation of this file.

      SUBROUTINE instar
*
*
*       Initialization of stars.
*       ------------------------
*
      include 'common6.h'
      REAL*8 tscls(20),lums(10),gb(10),tm,tn
      REAL*8 m0,m1,rm,lum,age,mc,rcc
      REAL*8 menv,renv,k2,k3
      parameter(k3=0.21d0)
      REAL*8 jspin,ospin
      REAL*8 vrotf
      EXTERNAL vrotf
*
*
*       Initialize mass loss variables & counters.
      tphys = 0.d0
      zmrg = 0.d0
      zmhe = 0.d0
      zmrs = 0.d0
      zmwd = 0.d0
      zmsn = 0.d0
      zmdot = 0.d0
      nmdot = 0
      nrg = 0
      nhe = 0
      nrs = 0
      nwd = 0
      nsn = 0
      nbh = 0
      ntz = 0
      nbs = 0
      nkick = 0
      nbkick = 0
*
*NB4 NEW
      zmnh = 0.d0
      zmbh = 0.d0
      zmsy = 0.d0
      nroche = 0
      nchaos = 0
      iqcoll = 0
      iblue = 0
      ncha = 0
      nsp = 0
      nhg = 0
      nnh = 0
      nbr = 0
      nas = 0
      nro = 0
      ndd = 0
      nspir = 0
      ncirc = 0
      nslp = 0
      ncont = 0
      ncoal = 0
      nce = 0
      instab = 0
      neint = 0
      nemod = 0
      nhyp = 0
      nhypc = 0.0
      ngb = 0
      nms = n
*
      tmdot = 1.0d+10
*NB6  IF (KZ(27).EQ.0) THEN
*         TSTAR = TSCALE
*     END IF
*
      WRITE (6,9)  zpars(11), zpars(12), zmet
 9    FORMAT(//,12x,'ABUNDANCES: X =',f6.3,' Y =',f6.3,' Z =',f6.3)
*NB6  zpars(11)=xhyd
*     zpars(12)=yhel
*
*       Calculate scale factor for spin angular momentum.
      spnfac = zmbar*su**2/(1.0d+06*tstar)
*
      epoch1 = epoch0
      DO 10 i = 1,n
*
*       Obtain stellar parameters at current epoch.
          IF(kz(12).EQ.2)THEN
             READ(12,*)m1,kw,m0,epoch1,ospin
          ELSE
             m1 = body(i)*zmbar
             m0 = m1
             IF(kstar(i).GT.1)THEN
                kw = kstar(i)
             ELSE
*       Include optional pre-mainsequence evolution.
                IF (m1.LT.8.0.AND.kz(41).GT.0) THEN
                    kw = -1
                ELSE
                    kw = 1
                END IF
*               IF(M0.LE.0.01D0) KW = 10
*               IF(M0.GE.100.D0) KW = 14
                IF (kz(45).GT.0.AND.m0.GT.10.0.AND.i.LE.kz(24)) kw = 14
             ENDIF
          ENDIF
          mc = 0.d0
          age = time*tstar - epoch1
          CALL star(kw,m0,m1,tm,tn,tscls,lums,gb,zpars)
*       Define preMS EPOCH & AGE if relevant.
          IF (kw.EQ.-1) THEN
              epoch(i) = tscls(15)
              age = time*tstar - epoch(i)
          END IF
          CALL hrdiag(m0,age,m1,tm,tn,tscls,lums,gb,zpars,
     &                rm,lum,kw,mc,rcc,menv,renv,k2)
*
*       Assign initial spin angular momentum. 
          IF(kz(22).EQ.3)THEN
             jspin = (k2*(m1-mc)*rm**2 + k3*mc*rcc**2)*ospin
          ELSE
             ospin = 45.35d0*vrotf(m1)/rm
             jspin = k2*m1*rm**2*ospin
             IF(kw.GT.1) jspin = 0.d0
          ENDIF
*
*       Convert from solar radii to scaled units (assume Sun = 1/215 AU).
          IF (kz(27).EQ.-2) THEN
              radius(i) = 0.d0
              zlmsty(i) = 0.d0
              spin(i) = 0.d0
          ELSE
              radius(i) = rm/su
              zlmsty(i) = lum
              spin(i) = jspin/spnfac
          END IF
*
*       Check neutron star option for GR capture (avoid large masses).
          IF (kz(27).EQ.3.AND.kz(28).EQ.4) THEN
              radius(i) = 1.0d-12/(3.0*rbar)
              kw = 13
              tev(i) = 1000.0
              m0 = 25.0
              epoch(i) = -15.0/tstar
              IF (i.EQ.1) WRITE (6,5)  m1, kw, radius(i)*su
    5         FORMAT (/,12x,'FIRST STAR:    M K* R/SU ',f7.2,i4,1p,e9.1)
          END IF
*
*       Initialize the stellar classification type (KW = 0 - 15).
          kstar(i) = kw
*
*       Save the initial mass of all the stars in sequential order.
          body(i) = m1/zmbar
          body0(i) = m0/zmbar
*
*       Set initial look-up time.
          epoch(i) = time*tstar - age
          tev0(i) = time
          CALL trdot(i,dtm,m1)
          tev(i) = dtm
          IF (kw.GE.13) tev(i) = 0.0
*
*       Determine the time for next stellar evolution check.
          IF (tev(i).LT.tmdot) THEN
              tmdot = tev(i)
          END IF
          IF (i.EQ.1.AND.kz(41).GT.0) THEN
              WRITE (6,8)  m1, tev(i), radius(i)*su
    8         FORMAT (/,12x,'BEGIN STAR #1    M TEV r*  ',
     &                       f7.2,1p,2e10.2)
          END IF
*
   10 CONTINUE
*
*       Define first quantized step < 1000 yrs (minimum interval for MDOT).
*       (changed to 100 yr to be consistent with trdot: 2/8/02)
      dt = 1.0d-04/tscale
      CALL stepk(dt,dtn)
      IF(dtn*tscale.LT.100.0) dtn = 2.d0*dtn
      stepx = dtn
*
*        Ensure binary components will be updated at the same time.
      DO 15 i = 1,nbin0
         i1 = 2*i - 1
         i2 = i1 + 1
         tev(i1) = min(tev(i1),tev(i2))
         tev(i1) = min(tev(i1),10.d0*stepx)
         tmdot = min(tmdot,tev(i1))
         tev(i2) = tev(i1)
   15 CONTINUE
*
*       Initialize stellar collision matrix.
*
      ktype(0,0) = 1
      do 20 j = 1,6
         ktype(0,j) = j
         ktype(1,j) = j
 20   continue
      ktype(0,7) = 4
      ktype(1,7) = 4
      do 25 j = 8,12
         if(j.ne.10)then
            ktype(0,j) = 6
         else
            ktype(0,j) = 3
         endif
         ktype(1,j) = ktype(0,j)
 25   continue
      ktype(2,2) = 3
      do 30 i = 3,14
         ktype(i,i) = i
 30   continue
      ktype(5,5) = 4
      ktype(7,7) = 1
      ktype(10,10) = 15
*       Change for CO+CO owing to Tout theory (21/11/08).
      ktype(11,11) = 12
      ktype(13,13) = 14
      do 35 i = 2,5
         do 40 j = i+1,12
            ktype(i,j) = 4
 40      continue
 35   continue
      ktype(2,3) = 3
      ktype(2,6) = 5
      ktype(2,10) = 3
      ktype(2,11) = 5
      ktype(2,12) = 5
      ktype(3,6) = 5
      ktype(3,10) = 3
      ktype(3,11) = 5
      ktype(3,12) = 5
      ktype(6,7) = 4
      ktype(6,8) = 6
      ktype(6,9) = 6
      ktype(6,10) = 5
      ktype(6,11) = 6
      ktype(6,12) = 6
      ktype(7,8) = 8
      ktype(7,9) = 9
      ktype(7,10) = 7
      ktype(7,11) = 9
      ktype(7,12) = 9
      ktype(8,9) = 9
      ktype(8,10) = 7
      ktype(8,11) = 9
      ktype(8,12) = 9
      ktype(9,10) = 7
      ktype(9,11) = 9
      ktype(9,12) = 9
      ktype(10,11) = 9
      ktype(10,12) = 9
      ktype(11,12) = 12
      do 45 i = 0,12
         ktype(i,13) = 13
         ktype(i,14) = 14
 45   continue
      ktype(13,14) = 14
*
*       Increase common-envelope cases by 100.
      do 50 i = 0,9
         do 55 j = i,14
            if(i.le.1.or.i.eq.7)then
               if(j.ge.2.and.j.le.9.and.j.ne.7)then
                  ktype(i,j) = ktype(i,j) + 100
               endif
            else
               ktype(i,j) = ktype(i,j) + 100
            endif
 55      continue
 50   continue
*
*       Assign the remaining values by symmetry.
      do 60 i = 1,14
         do 65 j = 0,i-1
            ktype(i,j) = ktype(j,i)
 65      continue
 60   continue
*
      WRITE (6,75)  ktype
   75 FORMAT (/,11x,' KTYPE: ',15i4,14(/,19x,15i4))
*
      RETURN
      END