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