data.f

Go to the documentation of this file.

      SUBROUTINE data
*
*
*       Initial conditions.
*       -------------------
*
      include 'common6.h'
      REAL*8  ran2
*
*
*       Initialize the portable random number generator (range: 0 to 1).
      kdum = -1
      rn1 = ran2(kdum)
*       Skip the first random numbers (IDUM1 specified at input).
      DO 1 k = 1,idum1
          rn1 = ran2(kdum)
    1 CONTINUE
*
*       Save random number sequence in COMMON for future use.
      idum1 = kdum
*
*       Read IMF parameters, # primordials, Z-abundance, epoch & HR interval.
      READ (5,*) alphas, body1, bodyn, nbin0, nhi0, zmet, epoch0, dtplot
      IF (n + 2*nbin0 + 2*nhi0.GE.nmax - 2) THEN
          WRITE (6,2)  n, nbin0, nhi0
    2     FORMAT (' FATAL ERROR!    BAD INPUT    N NBIN0 NHI0 ',i7,2i6)
          stop
      END IF
      IF (zmet.LE.0.0d0) zmet = 1.0d-04
      IF (zmet.GT.0.03) zmet = 0.03
      IF (kz(12).GT.0) dtplot = max(dtplot,deltat)
*
      IF (kz(19).GE.3) THEN
          CALL zcnsts(zmet,zpars)
          WRITE (6,4)  zpars(11), zpars(12), zmet
    4     FORMAT (//,12x,'ABUNDANCES:  X =',f6.3,'  Y =',f6.3,
     &                                           '  Z =',f7.4)
      END IF
*
*       Check options for reading initial conditions from input file.
      IF (kz(22).GE.2.OR.kz(22).EQ.-1) THEN
          zmass = 0.0
          DO 5 i = 1,n
              READ (10,*)  body(i), (x(k,i),k=1,3), (xdot(k,i),k=1,3)
              zmass = zmass + body(i)
    5     CONTINUE
*       Include possibility of a new IMF via option #20.
          IF (kz(22).GT.2.OR.kz(22).EQ.-1) go to 50
      END IF
*
*       Include optional initial conditions on #10 in astrophysical units.
*     IF (KZ(22).EQ.-1) THEN
*         CALL SETUP2
*         GO TO 30
*     END IF
*
*       Include the case of equal masses (ALPHAS = 1 or BODY1 = BODYN).
      IF (alphas.EQ.1.0.OR.body1.EQ.bodyn) THEN
          DO 10 i = 1,n
              body(i) = 1.0
   10     CONTINUE
*       Set provisional total mass (rescaled in routine SCALE).
          zmass = float(n)
          go to 40
      END IF
*
*       Choose between two realistic IMF's and standard Salpeter function.
      IF (kz(20).EQ.1) THEN
          CALL imf(body1,bodyn)
          go to 30
      ELSE IF (kz(20).GE.2) THEN
          CALL imf2(body1,bodyn)
          go to 30
      END IF
*       Assign #20 = 0 for no new IMF and no scaling with input on fort.10.
      IF (kz(22).EQ.2.AND.kz(20).EQ.0) go to 50
*
      WRITE (6,15)  alphas, body1, bodyn
   15 FORMAT (/,12x,'STANDARD IMF    ALPHAS =',f5.2,
     &              '  BODY1 =',f5.1,'  BODYN =',f5.2)
*
*       Generate a power-law mass function with exponent ALPHAS.
      alpha1 = alphas - 1.0
      fm1 = 1.0/body1**alpha1
      fmn = (fm1 - 1.0/bodyn**alpha1)/(float(n) - 1.0)
      zmass = 0.0d0
      const = 1.0/alpha1
*
*       Assign individual masses sequentially.
      DO 20 i = 1,n
          fmi = fm1 - float(i - 1)*fmn
          body(i) = 1.0/fmi**const
          zmass = zmass + body(i)
   20 CONTINUE
*
*       Scale the masses to <M> = 1 for now and set consistent total mass.
   30 zmbar1 = zmass/float(n)
      DO 35 i = 1,n
          body(i) = body(i)/zmbar1
   35 CONTINUE
      zmass = float(n)
*
*       Set up initial coordinates & velocities (uniform or Plummer model).
   40 IF (kz(22).EQ.0.OR.kz(22).EQ.1) THEN
          CALL setup
      END IF
*
   50 RETURN
*
      END