REAL FUNCTION deltcrit(a)
c     
*********************************************************************
c     subroutine to calc critical overdensity for collapse at time t,
c     for density field normalized at reference epoch a0=1.
c     NB delta is value extrapolated from collapse epoch to a0=1
c
c     eg For Omega=1 deltcrit = 1.686 * (1+z)
c     
c     Notation:
c     a = expansion factor rel to a0=1
c     t = time relative to t0=1
c     omega0 = omega at a0=1
c     delc = critical delta at time t
c     dldelcdlt = dln(delc)/dln(t)
*********************************************************************
c     

      IMPLICIT none
c
      REAL a,dldelcdlt
c
      REAL EPSOM,PI,delc0,eta0,sh0,ch0,tomega,d0,ch,sh,eta,t,acosh
      REAL omega0,lambda0,h0,omegab      
c
      INTEGER i,io,is,NTABLE,NSUM,NV
      REAL AMIN
      PARAMETER(AMIN=0.1,NTABLE=200,NSUM=2000,NV=1000)
      REAL omega0_save,lambda0_save
      REAL density,omflat(NV),delflat(NV),aflat(NTABLE),delta_flat(NTABLE),
     & sum,dlin,dlin0,x,x0,xp,dxp,h,aa,lambda,omega
      save omega0_save,lambda0_save,aflat,delta_flat
c
      PARAMETER(EPSOM=1e-5)
      PARAMETER(PI=3.141592654)
      external acosh
      common /cosmology/ omega0,lambda0,h0,omegab
      data omega0_save/0.0/
      data lambda0_save/0.0/
c     
      if(abs(1-omega0).le.EPSOM) then ! OMEGA=1
         delc0 = 3*(12*PI)**(2.0/3.0)/20
         deltcrit =  delc0/a 
         dldelcdlt = -2.0/3.0
c     
      else if((1-omega0).gt.EPSOM .and. lambda0.lt.EPSOM) then!OMEGA<1 LAMBDA=0
c     calc props at t=t0
         eta0 = acosh(2/omega0-1)
         sh0 = sinh(eta0)
         ch0 = cosh(eta0)
         tomega = 2*PI/(sh0-eta0)
         d0 = 3*sh0*(sh0-eta0)/(ch0-1)**2 -2 !linear growth factor
c     
c     calc props at expansion factor a
         ch = a*(ch0-1) + 1
         eta = acosh(ch)
         sh = sinh(eta)
         t = (sh-eta)/(sh0-eta0)
         deltcrit = 1.5*d0 *(1 + (tomega/t)**(2.0/3.0))
         dldelcdlt = -2.0/3.0 /(1 + (t/tomega)**(2.0/3.0))
c     
      else                      ! OMEGA>1
         if (omega0.ne.omega0_save .or. lambda0.ne.lambda0_save ) then
            write(*,*) 'Constructing look-up table for deltacrit(a)'
            write(*,*) 'for a flat Omega+Lambda=1 cosmology.'
            write(*,*) 'Note: derivative  dldelcdlt not implemented'
            omega0_save=omega0
            lambda0_save=lambda0
c          
c           Read Vince's file that tabulates deltcrit0 against omega0
            open(33,file='flat.data',status='old')
            read(33,*) ! skip header
            do i=1,NV
               read(33,*) omflat(i),density,delflat(i)
d              write(0,*) i,omflat(i),density,delflat(i)
            end do
            close(33)
c
c           Evaluate constant required to normalize the linear growth
c           factor
            x0=(2.0*(1.0/omega0-1.0))**0.333333
            sum=0.0
            dxp=x0/real(NSUM)
            do is=1,NSUM  
               xp=x0*(real(is)-0.5)/real(NSUM)
               sum=sum+((xp/(xp**3+2))**1.5 ) *dxp
            end do
            dlin0=sum*sqrt(x0**3+2.0)/sqrt(x0**3)
c           Tabulate deltcrit versus a for the specified values of
c           omega0 lambda0. Spacing in a is linear in order to enable
c           quick look up.
            do i=1,NTABLE
               aa=AMIN+(1.0-AMIN)*real(i-1)/real(NTABLE-1)
               aflat(i)=aa
               lambda=lambda0/(lambda0+(1.0-omega0-lambda0)/aa**2+omega0/aa**3)
               omega=omega0*lambda/(lambda0*aa**3)
               x=x0*aa
               sum=0.0
               dxp=x/real(NSUM)
               do is=1,NSUM  
                  xp=x*(real(is)-0.5)/real(NSUM)
                  sum=sum+((xp/(xp**3+2))**1.5 ) *dxp
               end do
               dlin=(sum*sqrt(x**3+2.0)/sqrt(x**3) )/dlin0
               call locate(omflat,NV,omega,io)
               if (io.lt.NV) then
                   h=(omflat(io+1)-omega)/(omflat(io+1)-omflat(io))
                   delta_flat(i)=(delflat(io)*h+delflat(io+1)*(1.0-h))/dlin
               else
                   delta_flat(i)=delflat(NV)/dlin
               end if
           end do
         end if
c
c        Evaluate deltcrit using look-up table
c     allow for a slightly greater than 1
         if (a.gt.AMIN .and. a.le.1.001) then
           i=1+int((a-AMIN)*real(NTABLE-1)/(1.0-AMIN))
           i=min(NTABLE-1,i)
           h=(aflat(i+1)-a)/(aflat(i+1)-aflat(i))
           deltcrit=delta_flat(i)*h+delta_flat(i+1)*(1.0-h)
           dldelcdlt=1.0  !not implemented
         else if (a.le.AMIN) then
c          extrapolate to higher redshift by approximating to
c          omega=1 at high redshift 
           deltcrit=delta_flat(1)*AMIN/a
           dldelcdlt=1.0  !not implemented           
         else
           write(*,*) 'a=',a
           stop 'deltcrit() look-up table only for a<1'
         end if

      end if
c     
      return
      END
c
***************************************************************************