Merge pull request #462 from rgknox/cxu-rgknox-gpp-season-cycle

leaf age classifications

Author: rgknox

biogeochem/EDCanopyStructureMod.F90

@@ -1242,7 +1242,7 @@ subroutine canopy_summarization( nsites, sites, bc_in )
 
              currentCohort%treelai = tree_lai(leaf_c,             &
                   currentCohort%pft, currentCohort%c_area, currentCohort%n, &
-                  currentCohort%canopy_layer, currentPatch%canopy_layer_tlai )
+                  currentCohort%canopy_layer, currentPatch%canopy_layer_tlai,currentCohort%vcmax25top )
 
              canopy_leaf_area = canopy_leaf_area + currentCohort%treelai *currentCohort%c_area
 
@@ -1417,11 +1417,11 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
 
           currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, &
                                            currentCohort%n, currentCohort%canopy_layer,               &
-                                           currentPatch%canopy_layer_tlai )    
+                                           currentPatch%canopy_layer_tlai,currentCohort%vcmax25top )    
 
           currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, &
                                            currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, &
-                                           currentPatch%canopy_layer_tlai, currentCohort%treelai )  
+                                           currentPatch%canopy_layer_tlai, currentCohort%treelai ,currentCohort%vcmax25top)  
 
           currentCohort%lai =  currentCohort%treelai *currentCohort%c_area/currentPatch%total_canopy_area 
           currentCohort%sai =  currentCohort%treesai *currentCohort%c_area/currentPatch%total_canopy_area  

biogeochem/EDCohortDynamicsMod.F90

@@ -407,7 +407,7 @@ subroutine spawn_patches( currentSite, bc_in)
           currentCohort => currentPatch%shortest
           do while(associated(currentCohort))       
 
-             allocate(nc)             
+             allocate(nc)
              if(hlm_use_planthydro.eq.itrue) call InitHydrCohort(CurrentSite,nc)
              call InitPRTCohort(nc)
              call zero_cohort(nc)

biogeochem/EDPatchDynamicsMod.F90

@@ -550,7 +550,7 @@ end subroutine storage_fraction_of_target
 
   ! =====================================================================================
 
-  real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai)
+  real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top)
 
     ! -----------------------------------------------------------------------------------
     ! LAI of individual trees is a function of the total leaf area and the total 
@@ -565,6 +565,8 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai)
     integer, intent(in)  :: cl                        ! canopy layer index
     real(r8), intent(in) :: canopy_lai(nclmax)        ! total leaf area index of 
                                                       ! each canopy layer
+    real(r8), intent(in) :: vcmax25top                ! maximum carboxylation rate at canopy
+                                                      ! top, ref 25C
 
     ! !LOCAL VARIABLES:
     real(r8) :: leafc_per_unitarea ! KgC of leaf per m2 area of ground.
@@ -601,7 +603,7 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai)
        end if
 
        ! Coefficient for exponential decay of 1/sla with canopy depth:
-       kn = decay_coeff_kn(pft)
+       kn = decay_coeff_kn(pft,vcmax25top)
 
        ! take PFT-level maximum SLA value, even if under a thick canopy (which has units of m2/gC),
        ! and put into units of m2/kgC
@@ -677,7 +679,7 @@ end function tree_lai
 
   ! ============================================================================
 
-  real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, canopy_lai, treelai )
+  real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, canopy_lai, treelai, vcmax25top )
 
     ! ============================================================================
     !  SAI of individual trees is a function of the LAI of individual trees
@@ -692,13 +694,15 @@ real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, canopy_la
     real(r8), intent(in) :: canopy_lai(nclmax) ! total leaf area index of 
                                                ! each canopy layer
     real(r8), intent(in) :: treelai            ! tree LAI for checking purposes only
+    real(r8), intent(in) :: vcmax25top         ! maximum carboxylation rate at top of crown
+                                               
 
     real(r8)             :: target_bleaf
     real(r8)             :: target_lai
 
     call bleaf(dbh,pft,canopy_trim,target_bleaf)
 
-    target_lai = tree_lai( target_bleaf, pft, c_area, nplant, cl, canopy_lai) 
+    target_lai = tree_lai( target_bleaf, pft, c_area, nplant, cl, canopy_lai, vcmax25top) 
 
     tree_sai   =  EDPftvarcon_inst%allom_sai_scaler(pft) * target_lai
 
@@ -2109,7 +2113,7 @@ end subroutine jackson_beta_root_profile
   ! =====================================================================================
 
 
-  real(r8) function decay_coeff_kn(pft)
+  real(r8) function decay_coeff_kn(pft,vcmax25top)
 
     ! ---------------------------------------------------------------------------------
     ! This function estimates the decay coefficient used to estimate vertical
@@ -2123,6 +2127,8 @@ real(r8) function decay_coeff_kn(pft)
 
     !ARGUMENTS
     integer, intent(in) :: pft
+    real(r8),intent(in) :: vcmax25top
+
 
     !LOCAL VARIABLES
     ! -----------------------------------------------------------------------------------
@@ -2131,7 +2137,7 @@ real(r8) function decay_coeff_kn(pft)
     ! kn = 0.11. Here, we derive kn from vcmax25 as in Lloyd et al 
     ! (2010) Biogeosciences, 7, 1833-1859
 
-    decay_coeff_kn = exp(0.00963_r8 * EDPftvarcon_inst%vcmax25top(pft) - 2.43_r8)
+    decay_coeff_kn = exp(0.00963_r8 * vcmax25top - 2.43_r8)
 
     return
   end function decay_coeff_kn

biogeochem/EDPhysiologyMod.F90

@@ -24,12 +24,15 @@ module FATESPlantRespPhotosynthMod
    use FatesGlobals,      only : fates_log
    use FatesConstantsMod, only : r8 => fates_r8
    use FatesConstantsMod, only : itrue
+   use FatesConstantsMod, only : nearzero
    use FatesInterfaceMod, only : hlm_use_planthydro
    use FatesInterfaceMod, only : hlm_parteh_mode
    use FatesInterfaceMod, only : numpft
+   use FatesInterfaceMod, only : nleafage
    use EDTypesMod,        only : maxpft
    use EDTypesMod,        only : nlevleaf
    use EDTypesMod,        only : nclmax
+   use EDTypesMod,        only : max_nleafage
    use EDTypesMod,        only : do_fates_salinity 
    use PRTGenericMod,     only : prt_carbon_allom_hyp
    use PRTGenericMod,     only : prt_cnp_flex_allom_hyp 
@@ -156,7 +159,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
     real(r8) :: co2_cpoint         ! CO2 compensation point (Pa)
     real(r8) :: btran_eff          ! effective transpiration wetness factor (0 to 1) 
     real(r8) :: bbb                ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s)
-    real(r8) :: kn(maxpft)         ! leaf nitrogen decay coefficient
+    real(r8) :: kn                 ! leaf nitrogen decay coefficient
     real(r8) :: cf                 ! s m**2/umol -> s/m (ideal gas conversion) [umol/m3]
     real(r8) :: gb_mol             ! leaf boundary layer conductance (molar form: [umol /m**2/s])
     real(r8) :: ceair              ! vapor pressure of air, constrained (Pa)
@@ -202,8 +205,6 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
     real(r8) :: lai_current        ! the LAI in the current leaf layer
     real(r8) :: cumulative_lai     ! the cumulative LAI, top down, to the leaf layer of interest
 
-
-    
     ! -----------------------------------------------------------------------------------
     ! Keeping these two definitions in case they need to be added later
     !
@@ -214,6 +215,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
     integer  :: cl,s,iv,j,ps,ft,ifp ! indices
     integer  :: nv                  ! number of leaf layers
     integer  :: NCL_p               ! number of canopy layers in patch
+    integer  :: iage                ! loop counter for leaf age classes
 
     ! Parameters
     ! -----------------------------------------------------------------------
@@ -313,11 +315,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
 
                do ft = 1,numpft
 
-                  ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used
-                  ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al 
-                  ! (2010) Biogeosciences, 7, 1833-1859
 
-                  kn(ft) = decay_coeff_kn(ft)
 
 
                   ! This is probably unnecessary and already calculated
@@ -373,7 +371,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
 
                      ! are there any leaves of this pft in this layer?
                      if(currentPatch%canopy_mask(cl,ft) == 1)then 
-                        
+  
                         ! Loop over leaf-layers
                         do iv = 1,currentCohort%nv
 
@@ -385,14 +383,17 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                            ! cohort-layer combo of interest.  
                            ! but in the vanilla case, we only re-calculate if it has
                            ! not been done yet.
+                           ! Other cases where we need to solve for every cohort
+                           ! in every leaf layer:  nutrient dynamic mode, multiple leaf
+                           ! age classes
                            ! ------------------------------------------------------------
 
                            if ( .not.rate_mask_z(iv,ft,cl) .or. &
                                  (hlm_use_planthydro.eq.itrue) .or. &
+                                 (nleafage > 1) .or. &
                                  (hlm_parteh_mode .ne. prt_carbon_allom_hyp )   ) then
 
-                              if (hlm_use_planthydro.eq.itrue) then
-
+                              if (hlm_use_planthydro.eq.itrue ) then
 
                                  bbb       = max (bbbopt(nint(c3psn(ft)))*currentCohort%co_hydr%btran(1), 1._r8)
                                  btran_eff = currentCohort%co_hydr%btran(1) 
@@ -427,10 +428,16 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                               if(do_fates_salinity)then
                                 btran_eff = btran_eff*currentPatch%bstress_sal_ft(ft)
                               endif 
-
 
+                              
+                              ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used
+                              ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al 
+                              ! (2010) Biogeosciences, 7, 1833-1859
+                              
+                              kn = decay_coeff_kn(ft,currentCohort%vcmax25top)
+
                               ! Scale for leaf nitrogen profile
-                              nscaler = exp(-kn(ft) * cumulative_lai)
+                              nscaler = exp(-kn * cumulative_lai)
 
                               ! Leaf maintenance respiration to match the base rate used in CN
                               ! but with the new temperature functions for C3 and C4 plants.
@@ -473,13 +480,15 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                               ! These rates are the specific rates used in the actual photosynthesis
                               ! calculations that take localized environmental effects (temperature)
                               ! into consideration.
+                              
+                              
 
                               call LeafLayerBiophysicalRates(currentPatch%ed_parsun_z(cl,ft,iv), &  ! in
                                                              ft,                                 &  ! in
-                                                             EDPftvarcon_inst%vcmax25top(ft),    &  ! in
-                                                             param_derived%jmax25top(ft),        &  ! in
-                                                             param_derived%tpu25top(ft),         &  ! in
-                                                             param_derived%kp25top(ft),          &  ! in
+                                                             currentCohort%vcmax25top,           &  ! in
+                                                             currentCohort%jmax25top,            &  ! in
+                                                             currentCohort%tpu25top,             &  ! in
+                                                             currentCohort%kp25top,              &  ! in
                                                              nscaler,                            &  ! in
                                                              bc_in(s)%t_veg_pa(ifp),             &  ! in
                                                              btran_eff,                          &  ! in
@@ -1742,7 +1751,7 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, &
                                               ! for this pft (umol electrons/m**2/s)
       real(r8), intent(in) :: tpu25top_ft     ! canopy top triose phosphate utilization rate at 25C 
                                               ! for this pft (umol CO2/m**2/s)
-      real(r8), intent(in) :: co2_rcurve_islope25top_ft      ! initial slope of CO2 response curve
+      real(r8), intent(in) :: co2_rcurve_islope25top_ft ! initial slope of CO2 response curve
                                               ! (C4 plants) at 25C, canopy top, this pft
       real(r8), intent(in) :: veg_tempk           ! vegetation temperature
       real(r8), intent(in) :: btran           ! transpiration wetness factor (0 to 1) 
@@ -1803,6 +1812,8 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, &
          tpu               = 0._r8
          co2_rcurve_islope = 0._r8
       else                                     ! day time
+
+         ! Vcmax25top was already calculated to derive the nscaler function
          vcmax25 = vcmax25top_ft * nscaler
          jmax25  = jmax25top_ft * nscaler
          tpu25   = tpu25top_ft * nscaler

biogeochem/FatesAllometryMod.F90

@@ -21,10 +21,14 @@ module EDInitMod
   use EDTypesMod                , only : AREA
   use EDTypesMod                , only : init_spread_near_bare_ground
   use EDTypesMod                , only : init_spread_inventory
+  use EDTypesMod                , only : first_leaf_aclass
+  use EDTypesMod                , only : leaves_on
+  use EDTypesMod                , only : leaves_off
   use FatesInterfaceMod         , only : bc_in_type
   use FatesInterfaceMod         , only : hlm_use_planthydro
   use FatesInterfaceMod         , only : hlm_use_inventory_init
   use FatesInterfaceMod         , only : numpft
+  use FatesInterfaceMod         , only : nleafage
   use ChecksBalancesMod         , only : SiteCarbonStock
   use FatesInterfaceMod         , only : nlevsclass
   use FatesAllometryMod         , only : h2d_allom
@@ -105,8 +109,8 @@ subroutine zero_site( site_in )
     site_in%total_burn_flux_to_atm = 0._r8
 
     ! PHENOLOGY 
-    site_in%status           = 0    ! are leaves in this pixel on or off?
-    site_in%dstatus          = 0
+    site_in%is_cold          = .false.    ! Is cold deciduous leaf-off triggered?
+    site_in%is_drought       = .false.    ! Is drought deciduous leaf-off triggered?
     site_in%ED_GDD_site      = nan  ! growing degree days
     site_in%ncd              = nan  ! no chilling days
     site_in%last_n_days(:)   = 999  ! record of last 10 days temperature for senescence model.
@@ -182,10 +186,10 @@ subroutine set_site_properties( nsites, sites)
     integer  :: s
     real(r8) :: leafon
     real(r8) :: leafoff
-    real(r8) :: stat
+    logical  :: stat
     real(r8) :: NCD
     real(r8) :: GDD
-    real(r8) :: dstat
+    logical  :: dstat
     real(r8) :: acc_NI
     real(r8) :: watermem
     integer  :: dleafoff
@@ -198,9 +202,9 @@ subroutine set_site_properties( nsites, sites)
        GDD      = 30.0_r8
        leafon   = 100.0_r8
        leafoff  = 300.0_r8
-       stat     = 2
+       stat     = .false.
        acc_NI   = 0.0_r8
-       dstat    = 2
+       dstat    = .false.
        dleafoff = 300
        dleafon  = 100
        watermem = 0.5_r8
@@ -211,9 +215,9 @@ subroutine set_site_properties( nsites, sites)
        GDD      = 0.0_r8
        leafon   = 0.0_r8
        leafoff  = 0.0_r8
-       stat     = 1
+       stat     = .false.
        acc_NI   = 0.0_r8
-       dstat    = 2
+       dstat    = .false.
        dleafoff = 300
        dleafon  = 100
        watermem = 0.5_r8
@@ -232,9 +236,8 @@ subroutine set_site_properties( nsites, sites)
           sites(s)%water_memory(1:numWaterMem) = watermem
        end if
 
-       sites(s)%status = stat
-       !start off with leaves off to initialise
-       sites(s)%dstatus= dstat
+       sites(s)%is_cold    = stat
+       sites(s)%is_drought = dstat
 
        sites(s)%acc_NI     = acc_NI
        sites(s)%frac_burnt = 0.0_r8
@@ -432,37 +435,27 @@ subroutine init_cohorts( site_in, patch_in, bc_in)
 
        call bstore_allom(temp_cohort%dbh, pft, temp_cohort%canopy_trim, b_store)
 
-
-       if( EDPftvarcon_inst%evergreen(pft) == 1) then
-          temp_cohort%laimemory = 0._r8
-          cstatus = 2
-       endif
-
-       if( EDPftvarcon_inst%season_decid(pft) == 1 ) then !for dorment places
-          if(site_in%status == 2)then 
-             temp_cohort%laimemory = 0.0_r8
-          else
-             temp_cohort%laimemory = b_leaf
-          endif
-          ! reduce biomass according to size of store, this will be recovered when elaves com on.
-          cstatus = site_in%status
+       temp_cohort%laimemory = 0._r8
+       cstatus = leaves_on
+       
+       if( EDPftvarcon_inst%season_decid(pft) == itrue .and. site_in%is_cold ) then
+          temp_cohort%laimemory = b_leaf
+          b_leaf = 0._r8
+          cstatus = leaves_off
        endif
-
-       if ( EDPftvarcon_inst%stress_decid(pft) == 1 ) then
-          if(site_in%dstatus == 2)then 
-             temp_cohort%laimemory = 0.0_r8
-          else
-             temp_cohort%laimemory = b_leaf
-          endif
-          cstatus = site_in%dstatus
+       
+       if ( EDPftvarcon_inst%stress_decid(pft) == itrue .and. site_in%is_drought ) then
+          temp_cohort%laimemory = b_leaf
+          b_leaf = 0._r8
+          cstatus = leaves_off
        endif
 
        if ( debug ) write(fates_log(),*) 'EDInitMod.F90 call create_cohort '
 
        call create_cohort(site_in, patch_in, pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, &
             b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & 
-            temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, 1, site_in%spread, bc_in)
-
+            temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, 1, &
+            site_in%spread, first_leaf_aclass, bc_in)
 
        deallocate(temp_cohort) ! get rid of temporary cohort