merged master and resolved conflict

README.md

@@ -28,3 +28,14 @@ https://github.com/E3SM-Project/E3SM
 https://github.com/ESCOMP/cesm
 https://github.com/ESCOMP/ctsm
 
+
+## Important Note About Host-Models and Compatible Branches:
+------------------------------------------------------------
+
+The FATES and E3SM teams maintain compatability of the NGEET/FATES master branch with the **E3SM master** branch. When changes to the FATES API force compatability updates with E3SM, there may be some modest lag time.
+
+The FATES team maintains compatability of the NGEET/FATES master branch with the **CTSM fates_next_api** branch.  Since the FATES team uses this branch for its internal testing, this compatability is tightly (immediately) maintained and these two should always be in sync.  However, CTSM master may become out of sync with FATES master for large periods (months) of time.
+
+
+
+

biogeochem/EDPatchDynamicsMod.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/FatesAllometryMod.F90

@@ -583,7 +583,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 
@@ -598,6 +598,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.
@@ -634,7 +636,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
@@ -710,7 +712,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
@@ -725,13 +727,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
 
@@ -2147,7 +2151,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
@@ -2161,6 +2165,8 @@ real(r8) function decay_coeff_kn(pft)
 
     !ARGUMENTS
     integer, intent(in) :: pft
+    real(r8),intent(in) :: vcmax25top
+
 
     !LOCAL VARIABLES
     ! -----------------------------------------------------------------------------------
@@ -2169,7 +2175,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

biogeophys/FatesPlantRespPhotosynthMod.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

main/EDInitMod.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