c    Example code to call make_tree() and so set up a 
c    a series of merger trees.
c
c 
      program  example_tree
**************************************variables********************************
      implicit none
      integer ILEVMAX,NLOGM
      parameter(ILEVMAX=1000,NLOGM=7)
      integer itree,ntrees,iter,iseed,ierr,iseed_prev,iseed0,nlev,
     & nhalo,nhalomax,jlev,ilogm,ihalo
      real mphalo,ahalo,mres,alev(ILEVMAX),wlev(ILEVMAX),astart,dlga
      integer jphalo(ILEVMAX),nhalolev(ILEVMAX),ifraglev(ILEVMAX)
      integer jlevel(*),jpar(*),nchild(*),jchild(*)
      real     mhalo(*),wk1(*),wk2(*),wk3(*),wk4(*),wk5(*)
      pointer (pjlevel,jlevel),(pmhalo,mhalo),(pjpar,jpar),(pnchild,nchild),
     &        (pjchild,jchild),(pwk1,wk1),(pwk2,wk2),(pwk3,wk3),(pwk4,wk4),
     &        (pwk5,wk5)
      real     deltcrit,sigmacdm,split
      EXTERNAL deltcrit,sigmacdm,split
      real omega0,lambda0,h0,omegab,nspec,gamma,sigma8,eps1,eps2
      real logmin,dlogm
      integer igwave,itrans
      real mhist(NLOGM,ILEVMAX),logm(NLOGM)
      common /splitpars/ eps1,eps2
      common /cosmology/ omega0,lambda0,h0,omegab
      common /power_spectrum/ nspec,gamma,sigma8,igwave,itrans
      data nhalomax/0/
      data ierr/1/
*******************************************************************************
c     Cosmology and power spectrum parameters
      omega0 =1.0
      lambda0=0.0
      h0=0.5
      omegab=0.1
      nspec=1.0
      gamma=0.25
      sigma8=0.5
      igwave=0
      itrans=1
c     Numerical Parameter used in split()
      eps1=0.1
      eps2=0.1
c
      mphalo = 1.0e+12 ! Msol/h  (final or parent halo mass)
      ahalo= 1.0       ! a=1/(1+z) at redshift at which the parent halo 
c                        is specified
      astart= 1.0/(1.0+10.0) ! a=1/(1+z) at highest tabulated redshift 
      mres=1.0e+07     ! Msol/h  mass resolution
      ntrees=100         ! number of merger trees to generate. NB best to produce
c                        and process one a time rather than store many.
      nlev=50         ! number of time levels on which to store the tree
      iseed0 = -436    ! negative random number seed
c
c  Set up a grid of time steps space however you like
      dlga = log(ahalo/astart)/(nlev-1)
      do jlev = 1,nlev
         alev(jlev) = ahalo*exp(-real(jlev-1)*dlga)
      end do
c  Set up and initialize the mass histogram bins
      logmin=7.0
      dlogm= 6.0/real(NLOGM-1)
      do ilogm=1,NLOGM
         logm(ilogm)=logmin +dlogm*real(ilogm-1)
         do jlev = 1,nlev
            mhist(ilogm,jlev)=0.0
         end do
      end do
c
c
c Loop over ntrees
      do itree = 1,ntrees
         write(0,*) 'itree=',itree
         iter=1
         do while (ierr.ne.0 .or. iter.eq.1)
            if (iter.gt.1) then
               write(0,*) 'repeating 1 of ',ntrees,' halos of mass ',mphalo
               write(0,*) 'due to need to increase array sizes.'
               iseed=iseed_prev !reset random number seed
            else 
               iseed=-abs(iseed_prev)-abs(iseed0) !advance to next seed
               iseed_prev=iseed
            end if
            iter=iter+1
c     Allocate memory used for halo tree 
            write(0,*) 'alocating memory'
            call memory(nhalo,nhalomax,ierr,nlev,mphalo,mres,
     &           pjlevel,pmhalo,pjpar,pnchild,pjchild,pwk1,pwk2,pwk3,pwk4,pwk5)
c     make halo tree
            write(0,*) 'calling make_tree with iseed=',iseed
c     NB the last five arrays passed to make_tree are used
c     as temporary work space
            call make_tree(mphalo,ahalo,mres,alev,NLEV,iseed,
     &           split,sigmacdm,deltcrit,
     &           nhalomax,ierr,nhalo,nhalolev,jphalo,
     &           jlevel,mhalo,jpar,nchild,jchild,wlev,ifraglev,
     &           wk1 ,wk2  ,wk3  ,wk4 ,wk5   )
         end do
         
c         Accumulate the mass histogram
         do jlev=1,nlev
            do ihalo = jphalo(jlev),jphalo(jlev)+nhalolev(jlev)-1
               ilogm=1+nint((log10(mhalo(ihalo))-logmin)/dlogm)
               if (ilogm.ge.1 .and. ilogm.le.NLOGM) then
                  mhist(ilogm,jlev)=mhist(ilogm,jlev)+mhalo(ihalo)
               else
                  write(0,*) 'mass of fragment outside histogram range'
               end if   
            end do    
         end do   

      end do

c   Normalize and Output the histograms
      do jlev = 1,nlev
         do ilogm=1,NLOGM

            mhist(ilogm,jlev)=mhist(ilogm,jlev)/(real(ntrees)*mphalo)
         end do
         write(0,'(100e10.3)') alev(jlev),(mhist(ilogm,jlev),ilogm=1,NLOGM)
      end do      

      stop
      end