program mkcloud
      implicit none

      double precision, parameter :: &
        pc     = 3.0856776e+18, & ! cm
        G_cgs  = 6.6743e-08,    & ! cm^3 / s^2 / g
        Msun   = 1.9884e+33,    & ! g
        l_unit = 1e3*pc,        & ! kpc
        m_unit = 1e10*Msun,     & ! 10^10 solar masses
        v_unit = 1e5,           & ! km/s
        t_unit = l_unit/v_unit, & ! resulting time unit
        G      = G_cgs * m_unit/l_unit/l_unit/l_unit*t_unit*t_unit, &
        pi     = 4.*atan(1.),   &
        H0     = 70./1e3,       & ! already in internal units
        rhoc   = 3.*H0**2/(8.*pi*G), &
        rhom   = 0.3*rhoc * 27. ! at z=2, the universe had 27 times the current density

      integer, parameter :: ngrid  = 20
      real,    parameter :: ngrid2 = ngrid/2.
      real, allocatable, dimension(:,:) :: pos,vel
      real, allocatable, dimension(:) :: mass
      integer, allocatable, dimension(:) :: id
      double precision :: redshift,time,masstab(6)

      real :: rcloud,Mtot,mpart,tff
      real :: x,y,z,r
      integer :: i,j,k,l,m,n

      redshift = 0.
      time     = 0.
      masstab  = 0.

      Mtot   = 1e5  ! = 1e15 Msun
!     rcloud = 10e3 ! = 10 Mpc
      rcloud = (3.*Mtot/(4.*pi*rhom))**(1./3.)
      tff    = sqrt(3.*pi/(32.*G*rhom))

! Two passes: first only counts, second also assigns.
      do l=1,2
        n = 0
        do i=1,ngrid
          x = (i-0.5-ngrid2)*rcloud/ngrid2
          do j=1,ngrid
            y = (j-0.5-ngrid2)*rcloud/ngrid2
            do k=1,ngrid
              z = (k-0.5-ngrid2)*rcloud/ngrid2
              r = sqrt(x**2+y**2+z**2)
              if(r.le.rcloud) then
                n = n+1
                if(l.eq.2) then
                  id(n)    = n
                  mass(n)  = mpart
                  pos(:,n) = [x,y,z]
                  vel(:,n) = [0.0,0.0,0.0]
                endif
              endif
            enddo
          enddo
        enddo
        if(l.eq.1) then
          m = n
          allocate(id(m),pos(3,m),vel(3,m),mass(m))
          mpart = Mtot/n
        endif
      enddo

      write(*,*) 'rcloud = ',rcloud
      write(*,*) 'mpart  = ',mpart
      write(*,*) 'tff    = ',tff

!     stop

! Write stuff to file.
      open(1,file='cloud.ic',form='unformatted')
      write(1) 'HEAD',8+256
      write(1) [0,m,0,0,0,0],masstab,redshift,time, &
               (0,i=1,(256-6*4-6*8-2*8)/4) ! pad to 256 bytes
      write(1) 'POS ',8+3*4*m
      write(1) pos
      write(1) 'VEL ',8+3*4*m
      write(1) vel
      write(1) 'ID  ',8+4*m
      write(1) id
      write(1) 'MASS',8+4*m
      write(1) mass
      close(1)

      end