This commit introduces multiple changes in the phenology code aiming at making the code a

bit simpler. This is not intended to be bit-for-bit, because some of the changes impact
the logic. However, the results shouldn't be dramatically different.

1. Moved the updates for days since last flushing and abscission events to the cumulative
   time subroutine.
2. Changed the initialisation of the dates of the last flushing/abscission deciduous
   events, so they can start at negative values. This change should eliminate the need
   for a few "if" statements in the main phenology subroutine.
3. Added a fourth case for cold deciduous phenology, which forces leaf flushing if plants
   remain dormant due to cold for a very long time. This should cause these plants to be
   vulnerable to frost in climates that are too cold for them to survive.
4. Saved some phenology transition tests to logical variables to make the code a bit more
   readable.
5. Restructured the cold-deciduous phenology code so the different checks are inside an
   if/elseif/else block (so only one change can happen at any time).
6. Replaced some additional checks for growing degree days that should be taken care by
   change 2 plus the definition of a growing season

Author: mpaiao

biogeochem/EDPhysiologyMod.F90

@@ -11,6 +11,7 @@ module EDInitMod
   use FatesConstantsMod         , only : primaryland
   use FatesConstantsMod         , only : nearzero
   use FatesConstantsMod         , only : rsnbl_math_prec
+  use FatesConstantsMod         , only : ndays_per_year
   use EDTypesMod                , only : min_patch_area_forced
   use FatesConstantsMod         , only : n_landuse_cats
   use FatesConstantsMod         , only : is_crop
@@ -48,11 +49,11 @@ module EDInitMod
   use FatesConstantsMod         , only : isemi_stress_decid
   use PRTGenericMod             , only : num_elements
   use PRTGenericMod             , only : element_list
-  use EDTypesMod                , only : phen_cstat_nevercold
-  use EDTypesMod                , only : phen_cstat_iscold
+  use EDTypesMod                , only : phen_cstat_timeoff
+  use EDTypesMod                , only : phen_cstat_tempoff
   use EDTypesMod                , only : phen_dstat_timeoff
   use EDTypesMod                , only : phen_dstat_moistoff
-  use EDTypesMod                , only : phen_cstat_notcold
+  use EDTypesMod                , only : phen_cstat_tempon
   use EDTypesMod                , only : phen_dstat_moiston
   use EDTypesMod                , only : phen_estat_evergreen
   use FatesInterfaceTypesMod         , only : bc_in_type,bc_out_type
@@ -393,11 +394,15 @@ subroutine set_site_properties( nsites, sites,bc_in )
     !
     ! !LOCAL VARIABLES:
     integer  :: s
-    real(r8) :: GDD               ! First guess for growing degree days
-    real(r8) :: liqvolmem         ! First guess for soil water content memory
-    real(r8) :: smpmem            ! First guess for soil matric potential memory
-    real(r8) :: tempmem           ! First guess for temperature memory
-    real(r8) :: elong_factor      ! Elongation factor (0 - fully off; 1 - fully on)
+    real(r8) :: gdd_1st           ! First guess for growing degree days
+    real(r8) :: liqvolmem_1st     ! First guess for soil water content memory
+    real(r8) :: smpmem_1st        ! First guess for soil matric potential memory
+    real(r8) :: tempmem_1st       ! First guess for temperature memory
+    real(r8) :: elong_factor_1st  ! First guess for elongation factor (fraction of leaves
+                                  !    flushed relative to the maximum given plant size)
+                                  !    0   - completely leafless
+                                  !    0.5 - 50% of leaves flushed
+                                  !    1.0 - 100% of leaves flushed
     integer  :: leafon            ! DOY for last leaf flushing event, initial guess
     integer  :: leafoff           ! DOY for last leaf abscission event, initial guess
     integer  :: ndleafon          ! Elapsed days since last flushing, initial guess
@@ -422,26 +427,26 @@ subroutine set_site_properties( nsites, sites,bc_in )
 
        !    For now, we impose that deciduous plants start the simulation fully flushed,
        ! but with the possibility of switching status if needed.
-       GDD           = 30.0_r8 ! Assume sufficiently warm conditions
-       liqvolmem     = 0.5_r8  ! Assume well watered conditions
-       smpmem        = 0._r8   ! Assume well watered conditions
-       tmpmem        = 15.0_r8 ! Assume sufficiently warm conditions
-       elong_factor  = 1._r8   ! Assume leaves fully flushed.
+       gdd_1st           = 30.0_r8 ! Assume sufficiently warm conditions
+       liqvolmem_1st     = 0.5_r8  ! Assume well watered conditions
+       smpmem_1st        = 0._r8   ! Assume well watered conditions
+       tempmem_1st       = 15.0_r8 ! Assume sufficiently warm conditions
+       elong_factor_1st  = 1._r8   ! Assume leaves fully flushed.
 
        !   Set the model integration date since the beginning of time.
        model_day_int = floor(hlm_model_day)
 
        do s = 1,nsites
-
           !    Define the last flushing and abscission dates based on the hemisphere. This
-          ! is to avoid the dates to be too far from optimal in the Southern Hemisphere, 
-          ! which could make deciduous plants spend multiple years until synchronising with
-          ! the local climate. For now, both cold- and drought-deciduous default to leaf
-          ! flushing during spring, and to leaf abscission during autumn. This assumes winter
-          ! dry season, which is common in many tropical regions. This is still not ideal
-          ! for some tropical and Mediterranean climates, which peak greenness during autumn
-          ! or spring. In the future, we could consider implementing some map that allows
-          ! assigning the expected dates for each location.
+          ! is to avoid first guess dates that are too far from optimal in the Southern 
+          ! Hemisphere, which could make deciduous plants spend multiple years until 
+          ! synchronising with the local climate. For now, both cold- and drought-deciduous
+          ! dates default to leaf flushing during spring, and to leaf abscission during 
+          ! autumn. This assumes winter dry season, which is common in many tropical 
+          ! regions. This is still not ideal for some tropical areas and for Mediterranean
+          ! climates, in which peak greenness occurs during autumn or spring, respectively.
+          ! In the future, we could consider implementing some map that allows assigning
+          ! the expected dates for each location.
           if (sites(s)%lat > 0._r8) then
              leafon  = 100 - hlm_day_of_year
              leafoff = 300 - hlm_day_of_year
@@ -450,10 +455,11 @@ subroutine set_site_properties( nsites, sites,bc_in )
              leafoff = 120 - hlm_day_of_year
           end if
 
-          !    Make sure the initial dates for last leaf abscission and leaf flushing events
-          ! precede the the beginning of the simulation. This check used to occur inside the
-          ! integrator, but it is better to ensure this during initialisation, so we reduce
-          ! the number of logical tests during time steps.
+
+          !    Make sure the initial dates for last leaf abscission and leaf flushing
+          ! events precede the the beginning of the simulation. This check used to occur
+          ! inside the integrator, but it is better to ensure this during initialisation,
+          ! so we reduce the number of logical tests during time steps.
           if ( leafon  > model_day_int ) then
              leafon  = leafon  - ndays_per_year
           end if
@@ -469,60 +475,57 @@ subroutine set_site_properties( nsites, sites,bc_in )
           ndleafoff = model_day_int - leafoff
 
 
-          estat         = phen_estat_evergreen
-          cstat         = phen_cstat_notcold     ! Leaves are on
-          dstat         = phen_dstat_moiston     ! Leaves are on
-
-
-          ! We set the number of chilling days and cold days to zero. This will be 
-          sites(s)%nchilldays(1:numpft)    = 0
-          sites(s)%ncolddays(1:numpft)     = 0 ! recalculated in phenology
-                                               ! immediately, so yes this
-                                               ! is memory-less, but needed
-                                               ! for first value in history file
-          sites(s)%phen_model_date         = model_day_int
-          sites(s)%grow_deg_days(1:numpft) = GDD
-
-          sites(s)%liqvol_memory(1:numWaterMem,1:numpft) = liqvolmem
-          sites(s)%smp_memory(1:numWaterMem,1:numpft) = smpmem
-          sites(s)%vegtemp_memory(1:num_vegtemp_mem) = tempmem
-
-          do ft = 1, numpft
-             ! These must be updated once the new parameter is integrated to FATES.
-             select case (prt_params%evergreen(ft))
-             case (itrue)
+          !---~---
+          !   Initialise site-level variables that are not PFT specific or that can be
+          ! set for all PFTs without considering the leaf phenological habit. Some of 
+          ! the memory and cumulative variables are currently used only for leaf
+          ! phenology, but they may be useful for other modules, so they are set for
+          ! all PFTs.
+          !    Some variables such as ncolddays are recalculated every day (memoryless)
+          ! but they need an initial value so they are defined for the first history 
+          ! file.
+          !---~---
+          sites(s)%phen_model_date                        = model_day_int
+          sites(s)%vegtemp_memory(1:num_vegtemp_mem)      = tempmem_1st
+          sites(s)%liqvol_memory (1:numWaterMem,1:numpft) = liqvolmem_1st
+          sites(s)%smp_memory    (1:numWaterMem,1:numpft) = smpmem_1st
+          sites(s)%grow_deg_days (1:numpft)               = gdd_1st
+          sites(s)%nchilldays    (1:numpft)               = 0
+          sites(s)%ncolddays     (1:numpft)               = 0 
+          sites(s)%elong_factor  (1:numpft)               = elong_factor_1st
+
+
+
+          !---~---
+          !   Fill in variables that do require checks for leaf phenology habit.
+          !---~---
+          do_phen_pft_init: do ft = 1, numpft
+             case_leaf_habit: select case (prt_params%phen_leaf_habit(ft))
+             case (ievergreen)
                 ! Evergreens. Phenological status is not used, set it accordingly.
-                !    Assign dummy values to other values.
-                sites(s)%phen_status (ft) = estat
-                sites(s)%leafondate  (ft) = cleafon    - hlm_day_of_year
-                sites(s)%leafoffdate (ft) = cleafoff   - hlm_day_of_year
-                sites(s)%ndaysleafon (ft) = 0
-                sites(s)%ndaysleafoff(ft) = 0
-             end select
-
-             select case (prt_params%season_decid(ft))
-             case (itrue)
-                ! Cold deciduous
-                sites(s)%phen_status (ft) = cstat
-                sites(s)%leafondate  (ft) = cleafon    - hlm_day_of_year
-                sites(s)%leafoffdate (ft) = cleafoff   - hlm_day_of_year
-                sites(s)%ndaysleafon (ft) = cndleafon
-                sites(s)%ndaysleafoff(ft) = cndleafoff
-             end select
-
-
-             select case (prt_params%stress_decid(ft))
+                sites(s)%phen_status (ft) = phen_estat_evergreen
+             case (ihard_season_decid)
+                !    Cold deciduous. Initialise both status and time since last flushing
+                ! and abscission days. We assume the initial status to be with leaves
+                ! fully flushed.
+                sites(s)%phen_status (ft) = phen_cstat_tempon
+                sites(s)%leafondate  (ft) = leafon
+                sites(s)%leafoffdate (ft) = leafoff
+                sites(s)%ndaysleafon (ft) = ndleafon
+                sites(s)%ndaysleafoff(ft) = ndleafoff
              case (ihard_stress_decid,isemi_stress_decid)
-                ! Drought deciduous
-                sites(s)%phen_status (ft) = dstat
-                sites(s)%leafondate  (ft) = dleafon    - hlm_day_of_year
-                sites(s)%leafoffdate (ft) = dleafoff   - hlm_day_of_year
-                sites(s)%ndaysleafon (ft) = dndleafon
-                sites(s)%ndaysleafoff(ft) = dndleafoff
-             end select
-          end do
+                ! Drought deciduous. Initialise both status and time since last flushing
+                ! and abscission days. We assume the initial status to be with leaves
+                ! fully flushed.
+                sites(s)%phen_status (ft) = phen_dstat_moiston
+                sites(s)%leafondate  (ft) = leafon
+                sites(s)%leafoffdate (ft) = leafoff
+                sites(s)%ndaysleafon (ft) = ndleafon
+                sites(s)%ndaysleafoff(ft) = ndleafoff
+             end select case_leaf_habit
+          end do do_phen_pft_init
+          !---~---
 
-          sites(s)%elong_factor(1:numpft) = elong_factor
 
           sites(s)%NF         = 0.0_r8
           sites(s)%NF_successful  = 0.0_r8
@@ -1206,7 +1209,7 @@ subroutine init_cohorts(site_in, patch_in, bc_in)
                phen_select: select case (prt_params%phen_leaf_habit(pft))
                case (ihard_season_decid)
                   select case (site_in%phen_status(pft))
-                  case (phen_cstat_nevercold, phen_cstat_iscold)
+                  case (phen_cstat_tempoff, phen_cstat_timeoff)
                      ! Cold deciduous, off season, assume complete abscission
                      efleaf_coh = 0.0_r8
                      effnrt_coh = 1.0_r8 - fnrt_drop_fraction 

main/EDInitMod.F90

@@ -98,16 +98,19 @@ module EDTypesMod
   !       The default phen_estat_evergreen is a dummy value reserved for evergreen PFTs.
 
   integer, parameter, public :: phen_estat_evergreen = 0  ! This PFT is evergreen.
-  integer, parameter, public :: phen_cstat_nevercold = 1  ! This (location/plant) has not experienced a cold period over a large number
-                                                          !    of days, leaves are dropped and flagged as non-cold region
-  integer, parameter, public :: phen_cstat_iscold    = 2  ! This (location/plant) is in a cold-state where leaves should have fallen
-  integer, parameter, public :: phen_cstat_notcold   = 3  ! This site is in a warm-state where leaves are allowed to flush
-
-  integer, parameter, public :: phen_dstat_timeoff   = 4  ! Leaves off due to time exceedance (drought phenology)
-  integer, parameter, public :: phen_dstat_moistoff  = 5  ! Leaves off due to moisture avail  (drought phenology)
-  integer, parameter, public :: phen_dstat_moiston   = 6  ! Leaves on due to moisture avail   (drought phenology)
-  integer, parameter, public :: phen_dstat_timeon    = 7  ! Leaves on due to time exceedance  (drought phenology)
-  integer, parameter, public :: phen_dstat_pshed     = 8  ! Leaves partially abscissing       (drought phenology)
+
+                                                          ! Season (cold) deciduous flags
+  integer, parameter, public :: phen_cstat_timeoff   = 1  ! - Leaves off due to time, leaves exceeded life span without cold days
+  integer, parameter, public :: phen_cstat_tempoff   = 2  ! - Leaves off after reaching multiple cold days
+  integer, parameter, public :: phen_cstat_tempon    = 3  ! - Leaves on after exceeding the growing degree day threshold
+  integer, parameter, public :: phen_cstat_timeon    = 4  ! - Leaves on after extended cold period, forced flushing
+
+                                                          ! Stress (drought) deciduous flags
+  integer, parameter, public :: phen_dstat_timeoff   = 5  ! Leaves off due to time, leaves exceeded life span without desiccation
+  integer, parameter, public :: phen_dstat_moistoff  = 6  ! Leaves off after moisture availability fell below the threshold
+  integer, parameter, public :: phen_dstat_moiston   = 7  ! Leaves on after moisture availability increased above the threshold
+  integer, parameter, public :: phen_dstat_timeon    = 8  ! Leaves on after extended dry period, forced flushing
+  integer, parameter, public :: phen_dstat_pshed     = 9  ! Leaves partially abscissing due to drying conditions
 
   ! PATCH FUSION 
   real(r8), parameter, public :: force_patchfuse_min_biomass = 0.005_r8   ! min biomass (kg / m2 patch area) below which to force-fuse patches

main/EDTypesMod.F90

@@ -14,16 +14,20 @@ module FatesCumulativeMemoryMod
    use EDTypesMod            , only : ed_site_type
    use EDTypesMod            , only : num_vegtemp_mem
    use EDTypesMod            , only : numWaterMem
-   use EDTypesMod            , only : phen_cstat_iscold
    use FatesAllometryMod     , only : set_root_fraction
+   use FatesConstantsMod     , only : ievergreen
+   use FatesConstantsMod     , only : ihard_season_decid
+   use FatesConstantsMod     , only : ihard_stress_decid
+   use FatesConstantsMod     , only : isemi_stress_decid
    use FatesConstantsMod     , only : ndays_per_year
    use FatesConstantsMod     , only : nearzero
    use FatesConstantsMod     , only : r8   => fates_r8
    use FatesConstantsMod     , only : tfrz => t_water_freeze_k_1atm
    use FatesInterfaceTypesMod, only : bc_in_type
-   use FatesInterfaceTypesMod, only : hlm_day_of_year
+   use FatesInterfaceTypesMod, only : hlm_model_day
    use FatesInterfaceTypesMod, only : numpft
    use FatesPatchMod         , only : fates_patch_type
+   use PRTParametersMod      , only : prt_params
 
    implicit none
    private
@@ -66,7 +70,6 @@ subroutine UpdateCumulativeMemoryVars(currentSite,bc_in)
       ! Update moisture-related memory variables.
       call UpdateMemoryMoisture(currentSite,bc_in)
 
-      return
    end subroutine UpdateCumulativeMemoryVars
 
 
@@ -84,53 +87,37 @@ subroutine UpdatePhenologyDate(currentSite)
       ! Arguments
       type(ed_site_type), intent(inout), target :: currentSite
       ! Local variables
-      integer :: ipft
-      integer :: date_offset_off
-      integer :: date_offset_on
+      integer :: ipft          ! PFT index
 
       !---~---
-      !    Advance elapsed time. The only reason this is a site variable instead of a 
-      ! global variable is that we need to save this information to the restart file,
-      ! and we do not have global scalars in the restart file.
+      !    Advance elapsed time, by using the host land model variable. The only reason
+      ! this is a site variable instead of a global variable is that we need to save this
+      ! information to the restart file, and we do not have global scalars in the restart
+      ! file. This value starts at zero and increases indefinitely.
       !---~---
-      currentSite%phen_model_date = currentSite%phen_model_date + 1
+      currentSite%phen_model_date = floor(hlm_model_day)
 
 
-      !---~---
-      !   Update the number of days since last flushing and abscission events.
-      !---~---
+      !   Update the number of days since last flushing and abscission events. This is
+      ! done for deciduous PFTs only. 
       do ipft = 1, numpft
-
-         !  If this is the beginning of the simulation, the last leaf abscission and/or
-         ! the last flushing events may have occurred before 
-         ! might not had occured yet, so set it to last year to get things rolling.
-         if ( model_day_int < currentSite%leafoffdate(ipft) ) then
-            date_offset = ndays_per_year
-         else 
-            date_offset = 0
-         end if
-         
-
-         !    Calculate the number of days since the last leaf flushing and leaf
-         ! abscission events.
-
+         select case (prt_params%phen_leaf_habit(ipft))
+         case (ihard_season_decid,ihard_stress_decid,isemi_stress_decid)
+            !---~---
+            !    Update the number of days since last flushing and abscission events, by
+            ! using the current date and the dates of the last events.  Note that we no
+            ! longer need to check whether this is the beginning of the simulation, 
+            ! because the cold start initialisation already fixes the dates of the last
+            ! event to be prior to the beginning of the simulation.
+            !---~---
             currentSite%ndaysleafoff(ipft) = &
-               currentSite%phen_model_date - (currentSite%leafoffdate(ipft) - ndays_per_year)
-         else
-            currentSite%ndaysleafoff(ipft) = 
-            model_day_int - currentSite%leafoffdate(ipft)
-         end if
-
-         if (model_day_int < currentSite%leafondate(ipft)) then
-            currentSite%ndaysleafon(ipft) = model_day_int - (currentSite%leafondate(ipft) - ndays_per_year)
-         else
-            currentSite%ndaysleafon(ipft) = model_day_int - currentSite%leafondate(ipft)
-         end if
+               currentSite%phen_model_date - currentSite%leafoffdate(ipft)
+            currentSite%ndaysleafon (ipft) = &
+               currentSite%phen_model_date - currentSite%leafondate (ipft)
+         end select
       end do
       !---~---
 
-
-      return
    end subroutine UpdatePhenologyDate
 
 
@@ -181,7 +168,6 @@ subroutine UpdateCumulativeThermal(currentSite,bc_in)
       currentSite%vegtemp_memory(1) = temp_in_C
 
 
-      return
    end subroutine UpdateCumulativeThermal
 
 
@@ -262,6 +248,5 @@ subroutine UpdateMemoryMoisture(currentSite,bc_in)
       end do pft_memory_loop
 
 
-      return
    end subroutine UpdateMemoryMoisture
 end module FatesCumulativeMemoryMod