Merge pull request #6 from rgknox/fire-threshold

merge up to tag 1.31.1_8.1.0 and add file version of cg_strikes

biogeochem/EDCohortDynamicsMod.F90

@@ -82,6 +82,7 @@ module EDCohortDynamicsMod
   use PRTAllometricCarbonMod, only : ac_bc_in_id_pft
   use PRTAllometricCarbonMod, only : ac_bc_in_id_ctrim
   use PRTAllometricCarbonMod, only : ac_bc_inout_id_dbh
+  use PRTAllometricCarbonMod, only : ac_bc_in_id_lstat
 
   !  use PRTAllometricCNPMod,    only : cnp_allom_prt_vartypes
 
@@ -354,7 +355,8 @@ subroutine InitPRTBoundaryConditions(new_cohort)
        call new_cohort%prt%RegisterBCInOut(ac_bc_inout_id_netdc,bc_rval = new_cohort%npp_acc)
        call new_cohort%prt%RegisterBCIn(ac_bc_in_id_pft,bc_ival = new_cohort%pft)
        call new_cohort%prt%RegisterBCIn(ac_bc_in_id_ctrim,bc_rval = new_cohort%canopy_trim)
-
+       call new_cohort%prt%RegisterBCIn(ac_bc_in_id_lstat,bc_ival = new_cohort%status_coh)
+
     case (prt_cnp_flex_allom_hyp)
 
        write(fates_log(),*) 'You have not specified the boundary conditions for the'

biogeochem/EDPhysiologyMod.F90

@@ -661,10 +661,31 @@ subroutine phenology( currentSite, bc_in )
     endif
     !
     ! accumulate the GDD using daily mean temperatures
-    if (bc_in%t_veg24_si .gt. tfrz) then
+    ! Don't accumulate GDD during the growing season (that wouldn't make sense)
+    if (bc_in%t_veg24_si .gt. tfrz.and. currentSite%cstatus == phen_cstat_iscold) then
        currentSite%grow_deg_days = currentSite%grow_deg_days + bc_in%t_veg24_si - tfrz
     endif
 
+    !this logic is to prevent GDD accumulating after the leaves have fallen and before the 
+    ! beginnning of the accumulation period, to prevend erroneous autumn leaf flushing. 
+    if(model_day_int>365)then !only do this after the first year to prevent odd behaviour
+
+    if(currentSite%lat .gt. 0.0_r8)then !Northern Hemisphere                                           
+      ! In the north, don't accumulate when we are past the leaf fall date.
+      ! Accumulation starts on day 1 of year in NH.  
+      ! The 180 is to prevent going into an 'always off' state after initialization
+      if( model_day_int .gt. currentSite%cleafoffdate.and.hlm_day_of_year.gt.180)then !
+       currentSite%grow_deg_days = 0._r8
+      endif
+     else !Southern Hemisphere 
+      ! In the South, don't accumulate after the leaf off date, and before the start of
+      ! the accumulation phase (day 181).  
+      if(model_day_int .gt. currentSite%cleafoffdate.and.hlm_day_of_year.lt.gddstart) then! 
+        currentSite%grow_deg_days = 0._r8
+      endif
+    endif
+    endif !year1 
+
     ! Calculate the number of days since the leaves last came on 
     ! and off. If this is the beginning of the simulation, that day might
     ! not had occured yet, so set it to last year to get things rolling
@@ -694,9 +715,12 @@ subroutine phenology( currentSite, bc_in )
     if ( (currentSite%cstatus == phen_cstat_iscold .or. &
           currentSite%cstatus == phen_cstat_nevercold) .and. &
          (currentSite%grow_deg_days > gdd_threshold) .and. &
+         (dayssincecleafoff > ED_val_phen_mindayson) .and. &
          (currentSite%nchilldays >= 1)) then
        currentSite%cstatus = phen_cstat_notcold  ! Set to not-cold status (leaves can come on)
        currentSite%cleafondate = model_day_int  
+       dayssincecleafon = 0 
+       currentSite%grow_deg_days = 0._r8 ! zero GDD for the rest of the year until counting season begins. 
        if ( debug ) write(fates_log(),*) 'leaves on'
     endif !GDD
 
@@ -907,7 +931,7 @@ subroutine phenology_leafonoff(currentSite)
     real(r8) :: store_c_transfer_frac  ! Fraction of storage carbon used to flush leaves
     integer  :: ipft
     real(r8), parameter :: leaf_drop_fraction = 1.0_r8
-
+    real(r8), parameter :: carbon_store_buffer = 0.10_r8
     !------------------------------------------------------------------------
 
     currentPatch => CurrentSite%oldest_patch   
@@ -936,8 +960,11 @@ subroutine phenology_leafonoff(currentSite)
                                                                 ! stop flow of carbon out of bstore. 
 
                    if(store_c>nearzero) then
-                      store_c_transfer_frac = &
-                            min(EDPftvarcon_inst%phenflush_fraction(ipft)*currentCohort%laimemory, store_c)/store_c
+                   ! flush either the amount required from the laimemory, or -most- of the storage pool
+                   ! RF: added a criterion to stop the entire store pool emptying and triggering termination mortality
+                   ! n.b. this might not be necessary if we adopted a more gradual approach to leaf flushing... 
+                       store_c_transfer_frac =  min((EDPftvarcon_inst%phenflush_fraction(ipft)* &
+                            currentCohort%laimemory)/store_c,(1.0_r8-carbon_store_buffer))
                    else
                       store_c_transfer_frac = 0.0_r8
                    end if
@@ -1204,16 +1231,16 @@ subroutine SeedDecay( litt )
 
     ! default value from Liscke and Loffler 2006 ; making this a PFT-specific parameter
     ! decays the seed pool according to exponential model
-    ! seed_decay_turnover is in yr-1
+    ! seed_decay_rate is in yr-1
     ! seed_decay is kg/day
     ! Assume that decay rates are same for all chemical species
 
     do pft = 1,numpft 
        litt%seed_decay(pft) = litt%seed(pft) * &
-             EDPftvarcon_inst%seed_decay_turnover(pft)*years_per_day
+             EDPftvarcon_inst%seed_decay_rate(pft)*years_per_day
 
        litt%seed_germ_decay(pft) = litt%seed_germ(pft) * &
-             EDPftvarcon_inst%seed_decay_turnover(pft)*years_per_day
+             EDPftvarcon_inst%seed_decay_rate(pft)*years_per_day
 
     enddo
 
@@ -1247,14 +1274,14 @@ subroutine SeedGermination( litt, cold_stat, drought_stat )
 
     !----------------------------------------------------------------------
 
-    ! germination_timescale is being pulled to PFT parameter; units are 1/yr
+    ! germination_rate is being pulled to PFT parameter; units are 1/yr
     ! thus the mortality rate of seed -> recruit (in units of carbon) 
-    ! is seed_decay_turnover(p)/germination_timescale(p)
+    ! is seed_decay_rate(p)/germination_rate(p)
     ! and thus the mortality rate (in units of individuals) is the product of 
     ! that times the ratio of (hypothetical) seed mass to recruit biomass
 
     do pft = 1,numpft
-       litt%seed_germ_in(pft) =  min(litt%seed(pft) * EDPftvarcon_inst%germination_timescale(pft), &
+       litt%seed_germ_in(pft) =  min(litt%seed(pft) * EDPftvarcon_inst%germination_rate(pft), &
                                      max_germination)*years_per_day
 
        !set the germination only under the growing season...c.xu

fire/SFMainMod.F90

@@ -664,6 +664,7 @@ subroutine  area_burnt_intensity ( currentSite )
 
     use FatesInterfaceMod, only : hlm_use_spitfire
     use EDParamsMod,       only : ED_val_nignitions
+    use EDParamsMod,       only : cg_strikes    ! fraction of cloud-to-ground ligtning strikes
     use FatesConstantsMod, only : years_per_day
     use SFParamsMod,       only : SF_val_fdi_alpha,SF_val_fuel_energy, &
          SF_val_max_durat, SF_val_durat_slope
@@ -680,8 +681,6 @@ subroutine  area_burnt_intensity ( currentSite )
 
     real(r8) size_of_fire !in m2
     real(r8),parameter :: km2_to_m2 = 1000000.0_r8 !area conversion for square km to square m
-    real(r8),parameter :: CG_strikes = 0.20_r8     !cloud to ground lightning strikes
-                                                   !Latham and Williams (2001)
 
     !  ---initialize site parameters to zero--- 
     currentSite%frac_burnt = 0.0_r8  
@@ -692,7 +691,7 @@ subroutine  area_burnt_intensity ( currentSite )
     currentSite%FDI  = 1.0_r8 - exp(-SF_val_fdi_alpha*currentSite%acc_NI)
 
     !NF = number of lighting strikes per day per km2 scaled by cloud to ground strikes
-    currentSite%NF = ED_val_nignitions * years_per_day * CG_strikes
+    currentSite%NF = ED_val_nignitions * years_per_day * cg_strikes
 
     ! If there are 15  lightning strikes per year, per km2. (approx from NASA product for S.A.) 
     ! then there are 15 * 1/365 strikes/km2 each day 

main/EDParamsMod.F90

@@ -1,9 +1,11 @@
 module EDParamsMod
+
    !
    ! module that deals with reading the ED parameter file
    !
 
    use FatesConstantsMod, only : r8 => fates_r8
+   use FatesConstantsMod, only : nearzero
    use FatesParametersInterface, only : param_string_length
    use FatesGlobals        , only : fates_log
    use FatesGlobals        , only : endrun => fates_endrun
@@ -42,6 +44,13 @@ module EDParamsMod
    real(r8),protected, public :: ED_val_patch_fusion_tol
    real(r8),protected, public :: ED_val_canopy_closure_thresh ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry
 
+
+   logical,protected, public :: active_crown_fire        ! flag, 1=active crown fire 0=no active crown fire
+   character(len=param_string_length),parameter :: fates_name_active_crown_fire = "fates_fire_active_crown_fire"
+
+   real(r8), protected, public :: cg_strikes             ! fraction of cloud to ground lightning strikes (0-1)
+   character(len=param_string_length),parameter :: fates_name_cg_strikes="fates_fire_cg_strikes"
+
    real(r8),protected,public  :: q10_mr     ! Q10 for respiration rate (for soil fragmenation and plant respiration)    (unitless)
    real(r8),protected,public  :: q10_froz   ! Q10 for frozen-soil respiration rates (for soil fragmentation)            (unitless)
 
@@ -50,6 +59,7 @@ module EDParamsMod
    real(r8),protected,allocatable,public :: ED_val_history_ageclass_bin_edges(:)
    real(r8),protected,allocatable,public :: ED_val_history_height_bin_edges(:)
 
+
    character(len=param_string_length),parameter,public :: ED_name_mort_disturb_frac = "fates_mort_disturb_frac"
    character(len=param_string_length),parameter,public :: ED_name_comp_excln = "fates_comp_excln"
    character(len=param_string_length),parameter,public :: ED_name_init_litter = "fates_init_litter"
@@ -69,10 +79,13 @@ module EDParamsMod
    character(len=param_string_length),parameter,public :: ED_name_phen_mindayson= "fates_phen_mindayson"   
    character(len=param_string_length),parameter,public :: ED_name_phen_ncolddayslim= "fates_phen_ncolddayslim"   
    character(len=param_string_length),parameter,public :: ED_name_phen_coldtemp= "fates_phen_coldtemp"   
-   character(len=param_string_length),parameter,public :: ED_name_cohort_fusion_tol= "fates_cohort_fusion_tol"   
+   character(len=param_string_length),parameter,public :: ED_name_cohort_fusion_tol= "fates_cohort_size_fusion_tol"   
    character(len=param_string_length),parameter,public :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol"
    character(len=param_string_length),parameter,public :: ED_name_canopy_closure_thresh= "fates_canopy_closure_thresh"      
 
+   ! Resistance to active crown fire
+
+
    character(len=param_string_length),parameter :: fates_name_q10_mr="fates_q10_mr"
    character(len=param_string_length),parameter :: fates_name_q10_froz="fates_q10_froz"
 
@@ -170,7 +183,6 @@ subroutine FatesParamsInit()
     ED_val_cohort_fusion_tol              = nan
     ED_val_patch_fusion_tol               = nan
     ED_val_canopy_closure_thresh          = nan    
-
     hydr_kmax_rsurf1                      = nan
     hydr_kmax_rsurf2                      = nan
 
@@ -333,6 +345,12 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=ED_name_history_height_bin_edges, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names_height)
 
+    call fates_params%RegisterParameter(name=fates_name_active_crown_fire, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
+    call fates_params%RegisterParameter(name=fates_name_cg_strikes, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
   end subroutine FatesRegisterParams
 
 
@@ -345,6 +363,9 @@ subroutine FatesReceiveParams(fates_params)
 
     class(fates_parameters_type), intent(inout) :: fates_params
 
+    real(r8) :: active_crown_fire_real !Local temp to transfer real data in file
+
+
     call fates_params%RetreiveParameter(name=ED_name_mort_disturb_frac, &
           data=fates_mortality_disturbance_fraction)
 
@@ -456,6 +477,13 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetreiveParameter(name=fates_name_q10_froz, &
           data=q10_froz)
 
+    call fates_params%RetreiveParameter(name=fates_name_active_crown_fire, &
+          data=active_crown_fire_real)
+    active_crown_fire = (abs(active_crown_fire_real-1.0_r8)<nearzero)
+
+    call fates_params%RetreiveParameter(name=fates_name_cg_strikes, &
+          data=cg_strikes)
+
     ! parameters that are arrays of size defined within the params file and thus need allocating as well
     call fates_params%RetreiveParameterAllocate(name=ED_name_history_sizeclass_bin_edges, &
           data=ED_val_history_sizeclass_bin_edges)
@@ -465,7 +493,7 @@ subroutine FatesReceiveParams(fates_params)
 
     call fates_params%RetreiveParameterAllocate(name=ED_name_history_height_bin_edges, &
           data=ED_val_history_height_bin_edges)
-
+    
 
   end subroutine FatesReceiveParams
 
@@ -517,6 +545,8 @@ subroutine FatesReportParams(is_master)
         write(fates_log(),fmt0) 'logging_dbhmax_infra = ',logging_dbhmax_infra
         write(fates_log(),fmt0) 'q10_mr = ',q10_mr
         write(fates_log(),fmt0) 'q10_froz = ',q10_froz
+        write(fates_log(),fmt0) 'cg_strikes = ',cg_strikes
+        write(fates_log(),'(a,L)') 'active_crown_fire = ',active_crown_fire
         write(fates_log(),*) '------------------------------------------------------'
 
      end if