c  z<->time conversion for flat LCDM cosmologies.
      program LCDM
      implicit NONE
      real*8 Omega_m, H0, z, t, H, t0
      character Conv*1
      include 'const.f'


c      write (*,*) 'Omega_m:'
      read (*,*) Omega_m
c      write (*,*)  Omega_m
c      write (*,*) 'Hubble constant H0 (km/s/Mpc):'
      read (*,*) H0
c      write (*,*) H0
c      write (*,*) 'Conversion: t or z'
      read (*,'(A)') Conv
c      write (*,'(A)')  Conv

c Hublle units of Gyr^{-1}:
      H = H0 * 1d8 * yr / pc

      write (*,*)
      call LCDM_z_time (.true.,0.d0,t0, Omega_m,H)
      write (*,'(A,1PE13.5,A)') 'Age of the Universe: t0=', t0,' Gyr'
      write (*,*)

      if (Conv .eq. 't')  then         
c         write (*,*) 'z='
         read (*,*) z
         call LCDM_z_time (.true.,z,t, Omega_m,H)
         write (*,'(A,F9.5,A,2(1PE13.5,A))') 'z=',z,', t=',t, 
     #                                       ' Gyr, t_lb=',t0-t,' Gyr'
        elseif (Conv .eq. 'z')  then
c         write (*,*) 't(Gyr)='
         read (*,*) t
         call LCDM_z_time (.false.,z,t, Omega_m,H)         
         write (*,'(A,E13.5,A,F9.5)') 't=',t,' Gyr, z=',z
        else
         write (*,*) 'Wrong conversion!'
         stop
        end if

      write (*,*)

      stop
      end
c=====================================================================================

      subroutine LCDMztime (direct,z,t, OmegaM,H0)
c z->time (when direct=true) or time->z (direct=false) coversion for flat LCDM or dS cosmology.
c (Time is in H0^{-1} units). (From Longair 1998; Coles & Lucchin 1995).
      implicit NONE
      real*8 t,z,theta,C2,x,OmegaM, H0
      logical direct
      
      if (OmegaM .eq. 1.d0)  then
c  dS case:
         if (direct)  then
            t = 2.d0/3.d0/H0 / (1.d0+z)**1.5d0
            else
            z = (2.d0/3.d0/H0/t)**(2.d0/3.d0) - 1.d0
            end if
         else
c  Flat LCDM:
         if (direct)  then
            theta = datan(dsqrt(OmegaM/(1.d0-OmegaM)) * 
     #              (1.d0+z)**1.5d0)
            t = 2.d0/3.d0/H0/dsqrt(1.d0-OmegaM) * 
     #           dlog((1.d0+dcos(theta))/dsin(theta))
            else
            C2 = dexp(3d0*t*H0*dsqrt(1.d0-OmegaM))
            x = (C2-1.d0)/(C2+1.d0)
            z = ((1.d0-x**2)*(1.d0-OmegaM)/OmegaM)**(1.d0/3.d0) * 
     #           x**(-2.d0/3.d0) - 1.d0
            end if
         end if   

      return
      end
c=====================================================================================