c----------------------------------------
       IMPLICIT REAL*8 (A-H,O-Z)
       Common /OmegRat/ OmLOm0,hsmall

       hsmall = 0.70
       Om     = 0.3
       write(*,'(a,f8.3)') ' Cosmological model: h=',
     &         hsmall,'     Omega_matter=',Om
       Do i=0,10
          z= i+1.d-10
          aexp =1./(1.+z)
          aa =Age(aexp)
          ab =Age2(aexp,Om,z)
         write(*,'(5(a,f8.4))') 'z=',z,' a=',aexp,
     &            ' Age(h=0.7)=',Age(aexp),
     &             ' ',Age2(aexp,Om,z)
       EndDo 
         Stop
         end
c     -----------------------------
       REAL*8 Function Age(aexp)
c     -----------------------------
c                 age for LCDM model h=0.70 
       IMPLICIT REAL*8 (A-H,O-Z)
       PARAMETER (Om0 =0.3)
       PARAMETER (tconst =9.30)
       
       x = ((1.-Om0)/Om0)**0.333333*aexp
       Age = tconst /sqrt(1.-Om0)*(log(x**1.5+sqrt(1.+x**3)))

       Return
       End
c     -----------------------------
      REAL*8  Function Age2(aexp,Om,z)
c     -----------------------------
c                                  age for LCDM model 
      IMPLICIT REAL*8 (A-H,O-Z)
      Common /OmegRat/ OmLOm0,hsmall
      Real*8 INTG
      external Hage

      aexp = 1./(1.+z)       
      zero = 1.d-12
      Oml    =Max(Abs(1.-Om),zero)
      OmLOm0 =Oml/Om
       Age2     =9.766/hsmall/sqrt(Om)*INTG(Hage,zero,aexp)
       Return
       End
C---------------------------------------
      REAL*8 FUNCTION Hnorm(x)
C---------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      Common /OmegRat/ OmLOm0,hsmall
         Hnorm =sqrt(1.d+0 
     &                  +OmLOm0*x**3 )
      Return
      End
C---------------------------------------
      REAL*8 FUNCTION Hage(x)
C---------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
         Hage =sqrt(x)/Hnorm(x) 
      Return
      End
C-------------------------------------------- Simpson integration
      REAL*8 FUNCTION INTG(FUNC,A,B)
C---------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (EPS=5.0d-4, JMAX=22)
      EXTERNAL FUNC
      OST=-1.D30
      OS= -1.D30
      ST =0.
      DO 11 J=1,JMAX
        CALL TRAPZD(FUNC,A,B,ST,J)
        INTG=(4.0d0*ST-OST)/3.0d0
        IF (ABS(INTG-OS).Le.EPS*ABS(OS)) RETURN
        OS=INTG
        OST=ST
11    CONTINUE
      WRITE (*,*)'Integration did not converge'
      RETURN
      END
C----------------------------------------------
      SUBROUTINE TRAPZD(FUNCC,A,B,S,N)
C---------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
        SAVE IT
        EXTERNAL FUNCC
      IF (N.EQ.1) THEN
        S=0.5d0*(B-A)*(FUNCC(A)+FUNCC(B))
        IT=1
      ELSE
        TNM=IT
        DEL=(B-A)/TNM
        X=A+0.5D0*DEL
        SUM=0.0D0
        DO 11 J=1,IT
          SUM=SUM+FUNCC(X)
          X=X+DEL
11      CONTINUE
        S=0.5D0*(S+(B-A)*SUM/TNM)
        IT=2*IT
      ENDIF
      RETURN
      END

