--- mix.f.orig	2009-03-12 14:59:57.000000000 -0700
+++ mix.f	2009-03-12 15:22:44.000000000 -0700
@@ -18,6 +18,10 @@
       PARAMETER(MCH=1.44D0)
       REAL*8 NETA,BWIND,HEWIND,MXNS
       COMMON /VALUE1/ NETA,BWIND,HEWIND,MXNS
+      REAL*8 PHI, Q, QC1, QC2, QC3, ALPHA, Z
+      REAL*8, EXTERNAL :: GET_HYDROGEN_BURNING_FRACTION
+      REAL*8, EXTERNAL :: GET_MIXING_PARAMETER
+      LOGICAL, PARAMETER :: COMPLETE_MIXING = .FALSE.
 *
 *
 *       Define global indices with body #I1 being most evolved.
@@ -80,10 +84,30 @@
 *       Evaluate apparent age and other parameters.
 *
       IF(ICASE.EQ.1)THEN
+         IF (COMPLETE_MIXING) THEN
 *       Specify new age based on complete mixing.
          IF(K1.EQ.7) KW = 7
          CALL star(KW,M03,M3,TMS3,TN,TSCLS,LUMS,GB,ZPARS)
          AGE3 = 0.1d0*TMS3*(AGE1*M1/TMS1 + AGE2*M2/TMS2)/M3
+         ELSE
+*       Compute new age based on recipe of Glebbeek&Pols, 2008 A&A 488.
+*       Metallicity
+         Z = 1.0d0 - ZPARS(11) - ZPARS(12)
+*       Get mixing alpha
+         ALPHA = GET_MIXING_PARAMETER(Z)
+*       Estimate mass loss and correct mass, equation (3)
+         Q = MIN(M1, M2) / MAX(M1, M2)
+         PHI = 0.3D0 * Q / (1.0D0 + Q)**2
+         M3 = M3 * (1.0D0 - PHI)
+         M03 = M03 * (1.0D0 - PHI)
+*       Get effective core masses, equation (6)
+         QC1 = GET_HYDROGEN_BURNING_FRACTION(M1, Z)
+         QC2 = GET_HYDROGEN_BURNING_FRACTION(M2, Z)
+         QC3 = GET_HYDROGEN_BURNING_FRACTION(M3, Z)
+*       Get fake age, equation (12) 
+         CALL star(KW,M03,M3,TMS3,TN,TSCLS,LUMS,GB,ZPARS)
+         AGE3 = TMS3*(QC1*AGE1*M1/TMS1+QC2*AGE2*M2/TMS2)/(QC3*ALPHA*M3)
+         END IF
       ELSEIF(ICASE.EQ.3.OR.ICASE.EQ.6.OR.ICASE.EQ.9)THEN
          MC3 = M1
          CALL gntage(MC3,M3,KW,ZPARS,M03,AGE3)