c Returns age of universe [Gyr] at redshift z for a universe with 
c H0 = h0 [km/s/Mpc], Omega_matter = om, Omega_lambda = ov
c
      FUNCTION tz(h0,om,ov,zz)
      IMPLICIT NONE
      INTEGER i,nstep,nstep1
      PARAMETER (nstep=100,nstep1=1000)
      REAL tz,h0,om,ov,z,dz,h_gyr,acosh,ee,x,zmid,dzi,
     .     s,s1,s2,t1,t2,a,b,zz
      PARAMETER (zmid=50.)
      acosh(x)=LOG(x+SQRT(x**2-1.))
      ee(x)=SQRT(om*(1.+x)**3+(1.-om-ov)*(1.+x)**2+ov)
c
      z=zz
c Force z=0 for negative redshifts to compute age.
      IF (zz.le.0.0) z=0.0 

      h_gyr=h0*1.022012E-3
      IF (om.eq.0..and.ov.eq.0.) THEN
         tz=1./(1.+z)
      ELSE IF (om.eq.1..and.ov.eq.0.) THEN
         tz=(2./3.)/(1.+z)**1.5
      ELSE IF (om.lt.1..and.ov.eq.0.) THEN
         tz=0.5*(om/(1.-om)**1.5)*
     .        ((2./om)*SQRT((1.-om)*(1.+om*z))/(1.+z)
     .        -acosh((2.-om*(1.-z))/(om*(1.+z))))
      ELSE IF (om.gt.1..and.ov.eq.0.) THEN
         tz=0.5*(om/(om-1.)**1.5)*
     .        (acos((2.-om*(1.-z))/(om*(1.+z)))
     .        -(2./om)*SQRT((om-1.)*(1.+om*z))/(1.+z))
      ELSE 
         IF (z.eq.0.) THEN
c first step: integral in [0,zmid]
            dz=zmid/FLOAT(nstep1)
            s=0.
            tz=0.
            t1=0.
            DO i=1,nstep1
               a=(i-1)*dz
               b=i*dz
               s1=1./((1.+a)*ee(a))
               s2=1./((1.+b)*ee(b))
               s=dz*(s1+s2)/2.
               t1=t1+s
            END DO
c second step: integral in [zmid,infty] with change of variable z -> 1/zi
            dzi=(1./zmid)/FLOAT(nstep)
            s=0.
            t2=0.
            DO i=2,nstep
               a=(i-1)*dzi
               b=i*dzi
               s1=1./((1.+1./a)*ee(1./a)*a*a)
               s2=1./((1.+1./b)*ee(1./b)*b*b)
               s=dzi*(s1+s2)/2.
               t2=t2+s
            END DO
            tz=t1+t2
         ELSE
c integral in [z,infty] with change of variable z -> 1/zi
            dzi=(1./z)/FLOAT(nstep1)
            s=0.
            tz=0.
            DO i=2,nstep1
               a=(i-1)*dzi
               b=i*dzi
               s1=1./((1.+1./a)*ee(1./a)*a*a)
               s2=1./((1.+1./b)*ee(1./b)*b*b)
               s=dzi*(s1+s2)/2.
               tz=tz+s
            END DO
         END IF
      END IF
      tz=tz/h_gyr
      RETURN
      END