diff --git a/src/biosphere_biomee.mod.f90 b/src/biosphere_biomee.mod.f90
index 7edfbd43..9a2733d4 100644
--- a/src/biosphere_biomee.mod.f90
+++ b/src/biosphere_biomee.mod.f90
@@ -181,6 +181,8 @@ subroutine biosphere_annual(out_biosphere)
     ! Re-organize cohorts
     !---------------------------------------------
     call kill_lowdensity_cohorts( vegn )
+
+    call kill_old_grass( vegn ) 
     
     call relayer_cohorts( vegn )
     
diff --git a/src/datatypes.mod.f90 b/src/datatypes.mod.f90
index 36867655..0d1230d0 100644
--- a/src/datatypes.mod.f90
+++ b/src/datatypes.mod.f90
@@ -158,6 +158,7 @@ module datatypes
     integer :: layer      = 1.0          ! the layer of this cohort (numbered from top, top layer=1)
     integer :: firstlayer = 0.0          ! 0 = never been in the first layer; 1 = at least one year in first layer
     real    :: layerfrac  = 0.0          ! fraction of layer area occupied by this cohort
+    real    :: leaf_age   = 0.0          ! leaf age (years)
 
     !===== Population structure
     real :: nindivs   = 1.0          ! density of vegetation, individuals/m2
diff --git a/src/vegetation_biomee.mod.f90 b/src/vegetation_biomee.mod.f90
index 226d2973..b7204722 100755
--- a/src/vegetation_biomee.mod.f90
+++ b/src/vegetation_biomee.mod.f90
@@ -16,6 +16,7 @@ module md_vegetation_biomee
   public :: vegn_reproduction, vegn_annualLAImax_update !, annual_calls
   public :: vegn_nat_mortality, vegn_species_switch !, vegn_starvation
   public :: relayer_cohorts, vegn_mergecohorts, kill_lowdensity_cohorts
+  public :: kill_old_grass
   public :: vegn_annual_starvation, Zero_diagnostics, reset_vegn_initial
 
 contains
@@ -233,6 +234,9 @@ subroutine vegn_growth_EW( vegn )
 
       if (cc%status == LEAF_ON) then
         
+        !update leaf age 
+        cc%leaf_age = cc%leaf_age + 1.0/365.0
+
         ! Get carbon from NSC pool. This sets cc%C_growth
         call fetch_CN_for_growth( cc )
 
@@ -340,6 +344,7 @@ subroutine vegn_growth_EW( vegn )
         cc%psapw%c%c12    = cc%psapw%c%c12 + dBSW
         cc%pseed%c%c12    = cc%pseed%c%c12 + dSeed
         cc%plabl%c%c12    = cc%plabl%c%c12 - dBR - dBL - dSeed - dBSW
+        cc%leaf_age = (1.0 - dBL/cc%pleaf%c%c12) * cc%leaf_age !NEW
         cc%resg = 0.5 * (dBR + dBL + dSeed + dBSW) !  daily
 
         ! update nitrogen pools, Nitrogen allocation
@@ -490,7 +495,7 @@ subroutine vegn_phenology( vegn )
     integer :: i
     ! real    :: grassdensity   ! for grasses only
     ! real    :: BL_u,BL_c
-     integer :: GrassMaxL = 3 
+    integer :: GrassMaxL = 3 
     real    :: ccNSC, ccNSN
     logical :: cc_firstday != .false.
     logical :: TURN_ON_life, TURN_OFF_life
@@ -526,7 +531,7 @@ subroutine vegn_phenology( vegn )
         .and.(vegn%thetaS>sp%betaON)                                 &  ! Water
         .and.(.NOT.(sp%lifeform==0 .and. cc%layer > GrassMaxL))      &  ! If grasses, layer< 3
          )
-      
+
       cc_firstday = .false.
       if (TURN_ON_life) then
         cc%status = LEAF_ON ! Turn on a growing season
@@ -534,7 +539,8 @@ subroutine vegn_phenology( vegn )
       endif
 
       ! Reset grass density at the first day of a growing season
-      if (cc_firstday .and. sp%lifeform == 0 .and. cc%age > 2.0) then
+      !if (cc_firstday .and. sp%lifeform == 0 .and. cc%age > 2.0) then
+      if  (sp%lifeform ==0 .and. (cc_firstday .and. cc%age>0.5)) then
         
         ! reset grass density and size for perenials
         ccNSC   = (cc%plabl%c%c12 + cc%pleaf%c%c12 + cc%psapw%c%c12 + &
@@ -573,7 +579,8 @@ subroutine vegn_phenology( vegn )
 
     enddo cohortloop2
 
-    if (TURN_ON_life) call relayer_cohorts( vegn )
+    !if (TURN_ON_life) call relayer_cohorts( vegn )
+    if (cc_firstday) call relayer_cohorts(vegn)
 
     ! OFF of a growing season
     cohortloop3: do i = 1,vegn%n_cohorts
@@ -1873,6 +1880,54 @@ subroutine kill_lowdensity_cohorts( vegn )
     endif
   end subroutine kill_lowdensity_cohorts
 
+subroutine kill_old_grass(vegn)
+  ! kill old grass cohorts
+  ! Weng, 01/22/2023
+  type(vegn_tile_type), intent(inout) :: vegn
+
+  ! ---- local vars
+  type(cohort_type), pointer :: cx, cc(:) ! array to hold new cohorts
+  logical :: merged(vegn%n_cohorts)        ! mask to skip cohorts that were already merged
+  real, parameter :: mindensity = 0.25E-4
+  logical :: OldGrass
+  integer :: i,j,k
+
+ ! calculate the number of cohorts that are not old grass
+  k = 0
+  do i = 1, vegn%n_cohorts
+    cx =>vegn%cohorts(i)
+    associate(sp=>spdata(cx%species))
+      OldGrass = (sp%lifeform ==0 .and. cx%age > 3.0)
+      if (.not. OldGrass) k=k+1
+    end associate
+  enddo
+  if (k==0)then
+     write(*,*)'in kill_old_grass: All cohorts are old grass, No action!'
+     !stop
+  endif
+
+  ! exclude cohorts that are old grass
+  if (k>0 .and. k<vegn%n_cohorts)then
+     allocate(cc(k))
+     j=0
+     do i = 1,vegn%n_cohorts
+        cx =>vegn%cohorts(i)
+        associate(sp=>spdata(cx%species))
+        OldGrass = (sp%lifeform ==0 .and. cx%age > 3.0)
+        if (.not. OldGrass) then
+           j=j+1
+           cc(j) = cx
+        else
+           ! Carbon and Nitrogen from plants to soil pools
+           call plant2soil(vegn,cx,cx%nindivs)
+        endif
+        end associate
+     enddo
+     vegn%n_cohorts = j
+     deallocate (vegn%cohorts)
+     vegn%cohorts=>cc
+  endif
+  end subroutine kill_old_grass
 
   subroutine merge_cohorts(c1, c2)
     !////////////////////////////////////////////////////////////////
@@ -2230,6 +2285,7 @@ subroutine initialize_vegn_tile( vegn, nCohorts )
         cx => vegn%cohorts(i)
         cx%status  = LEAF_OFF ! ON=1, OFF=0 ! ON
         cx%layer   = 1
+        cx%age = 0
         cx%species = INT(myinterface%init_cohort(i)%init_cohort_species)
         cx%ccID    =  i
         cx%plabl%c%c12     = myinterface%init_cohort(i)%init_cohort_nsc