Merge pull request #287 from rgknox/rgknox-slascaler-fix

removed slascaler and various fates parameter changes

Author: rgknox

biogeochem/EDGrowthFunctionsMod.F90

@@ -83,10 +83,10 @@ real(r8) function Hite( cohort_in )
     ! if the hite is larger than the maximum allowable height (set by dbhmax) then 
     ! set the height to the maximum value. 
     ! this could do with at least re-factoring and probably re-thinking. RF
-    if(cohort_in%dbh <= EDPftvarcon_inst%max_dbh(cohort_in%pft)) then
+    if(cohort_in%dbh <= EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft)) then
        h = (10.0_r8**(log10(cohort_in%dbh) * m + c))
     else 
-       h = (10.0_r8**(log10(EDPftvarcon_inst%max_dbh(cohort_in%pft))*m + c))
+       h = (10.0_r8**(log10(EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft))*m + c))
     endif
     Hite = h 
 
@@ -107,26 +107,22 @@ real(r8) function Bleaf( cohort_in )
     real(r8) :: dbh2bl_a 
     real(r8) :: dbh2bl_b
     real(r8) :: dbh2bl_c
-    real(r8) :: slascaler ! changes the target biomass according to the SLA
 
     dbh2bl_a  =  EDPftvarcon_inst%allom_d2bl1(cohort_in%pft)
     dbh2bl_b  =  EDPftvarcon_inst%allom_d2bl2(cohort_in%pft)
     dbh2bl_c  =  EDPftvarcon_inst%allom_d2bl3(cohort_in%pft)
-    slascaler = EDPftvarcon_inst%allom_d2bl_slascaler(cohort_in%pft)/EDPftvarcon_inst%slatop(cohort_in%pft)
 
     if(cohort_in%dbh < 0._r8.or.cohort_in%pft == 0.or.cohort_in%dbh > 1000.0_r8)then
        write(fates_log(),*) 'problems in bleaf',cohort_in%dbh,cohort_in%pft
     endif
 
-    if(cohort_in%dbh <= EDPftvarcon_inst%max_dbh(cohort_in%pft))then
+    if(cohort_in%dbh <= EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft))then
        bleaf = dbh2bl_a * (cohort_in%dbh**dbh2bl_b) * EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bl_c 
     else  
-       bleaf = dbh2bl_a * (EDPftvarcon_inst%max_dbh(cohort_in%pft)**dbh2bl_b) * &
+       bleaf = dbh2bl_a * (EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft)**dbh2bl_b) * &
             EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bl_c
     endif  
 
-    bleaf = bleaf * slascaler
-    
     !write(fates_log(),*) 'bleaf',bleaf, slascaler,cohort_in%pft
 
     !Adjust for canopies that have become so deep that their bottom layer is not producing any carbon... 
@@ -240,14 +236,14 @@ real(r8) function c_area( cohort_in )
 
     if (DEBUG_growth) then
        write(fates_log(),*) 'z_area 1',cohort_in%dbh,cohort_in%pft
-       write(fates_log(),*) 'z_area 2',EDPftvarcon_inst%max_dbh
+       write(fates_log(),*) 'z_area 2',EDPftvarcon_inst%allom_dbh_maxheight
        write(fates_log(),*) 'z_area 3',EDPftvarcon_inst%woody
        write(fates_log(),*) 'z_area 4',cohort_in%n
        write(fates_log(),*) 'z_area 5',cohort_in%siteptr%spread
        write(fates_log(),*) 'z_area 6',cohort_in%canopy_layer
     end if
 
-    dbh = min(cohort_in%dbh,EDPftvarcon_inst%max_dbh(cohort_in%pft))
+    dbh = min(cohort_in%dbh,EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft))
 
     ! ----------------------------------------------------------------------------------
     ! The function c_area is called during the process of canopy position demotion
@@ -360,7 +356,7 @@ real(r8) function dDbhdBd( cohort_in )
     dBD_dDBH =  dbh2bd_c*dbh2bd_a*(cohort_in%hite**dbh2bd_b)*(cohort_in%dbh**(dbh2bd_c-1.0_r8))* &
          (EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bd_d)  
 
-    if(cohort_in%dbh < EDPftvarcon_inst%max_dbh(cohort_in%pft))then
+    if(cohort_in%dbh < EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft))then
        dH_dDBH = (10.0_r8**c)*m*(cohort_in%dbh**(m-1.0_r8))          
 
        dBD_dDBH =  dBD_dDBH + dbh2bd_b*dbh2bd_a*(cohort_in%hite**(dbh2bd_b - 1.0_r8))* &
@@ -396,7 +392,7 @@ real(r8) function dDbhdBl( cohort_in )
     dblddbh = dbh2bl_b*dbh2bl_a*(cohort_in%dbh**dbh2bl_b)*(EDPftvarcon_inst%wood_density(cohort_in%pft)**dbh2bl_c)
     dblddbh = dblddbh*cohort_in%canopy_trim
 
-    if( cohort_in%dbh<EDPftvarcon_inst%max_dbh(cohort_in%pft) ) then
+    if( cohort_in%dbh<EDPftvarcon_inst%allom_dbh_maxheight(cohort_in%pft) ) then
         dDbhdBl = 1.0_r8/dblddbh
     else
         dDbhdBl = 1.0d15  ! At maximum size, the leaf biomass is saturated, dbl=0

biogeochem/EDPhysiologyMod.F90

@@ -970,7 +970,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in)
     !only if carbon balance is +ve
     if ((currentCohort%balive >= target_balive).AND.(currentCohort%carbon_balance >  0._r8))then 
        ! fraction of carbon going into active vs structural carbon        
-       if (currentCohort%dbh <= EDPftvarcon_inst%max_dbh(currentCohort%pft))then ! cap on leaf biomass
+       if (currentCohort%dbh <= EDPftvarcon_inst%allom_dbh_maxheight(currentCohort%pft))then ! cap on leaf biomass
           dbldbd = dDbhdBd(currentCohort)/dDbhdBl(currentCohort) 
           dbrdbd = EDPftvarcon_inst%allom_l2fr(currentCohort%pft) * dbldbd
           dhdbd_fn = dhdbd(currentCohort)

main/EDParamsMod.F90

@@ -88,6 +88,10 @@ module EDParamsMod
 
    real(r8),protected :: logging_collateral_frac     ! Ratio of collateral mortality to direct logging mortality
    character(len=param_string_length),parameter :: logging_name_collateral_frac = "fates_logging_collateral_frac"
+
+   real(r8),protected :: logging_coll_under_frac ! Fraction of understory plants that die when logging disturbance
+                                                 ! is generated
+   character(len=param_string_length),parameter :: logging_name_coll_under_frac = "fates_logging_coll_under_frac"
 
    real(r8),protected :: logging_direct_frac         ! Fraction of stems logged per event
    character(len=param_string_length),parameter :: logging_name_direct_frac = "fates_logging_direct_frac"
@@ -251,6 +255,9 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=logging_name_collateral_frac, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
+    call fates_params%RegisterParameter(name=logging_name_coll_under_frac, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names)
+
     call fates_params%RegisterParameter(name=logging_name_direct_frac, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
@@ -356,6 +363,9 @@ subroutine FatesReceiveParams(fates_params)
 
     call fates_params%RetreiveParameter(name=logging_name_collateral_frac, &
           data=logging_collateral_frac)
+    
+    call fates_params%RetreiveParameter(name=logging_name_coll_under_frac, &
+          data=logging_coll_under_frac)
 
     call fates_params%RetreiveParameter(name=logging_name_direct_frac, &
           data=logging_direct_frac)

main/EDPftvarcon.F90

@@ -26,10 +26,11 @@ module EDPftvarcon
   !ED specific variables. 
   type, public ::  EDPftvarcon_type
      real(r8), allocatable :: pft_used           (:) ! Switch to turn on and off PFTs
-     real(r8), allocatable :: max_dbh            (:) ! maximum dbh at which height growth ceases...
+    
      real(r8), allocatable :: freezetol          (:) ! minimum temperature tolerance (NOT CURRENTY USED)
      real(r8), allocatable :: wood_density       (:) ! wood density  g cm^-3  ...
      real(r8), allocatable :: hgt_min            (:) ! sapling height m
+     real(r8), allocatable :: dbh_repro_threshold(:) ! diameter at which mature plants shift allocation
      real(r8), allocatable :: dleaf              (:) ! leaf characteristic dimension length (m)
      real(r8), allocatable :: z0mr               (:) ! ratio of roughness length of vegetation to height (-) 
      real(r8), allocatable :: displar            (:) ! ratio of displacement height to canopy top height (-)
@@ -45,6 +46,7 @@ module EDPftvarcon
      real(r8), allocatable :: root_long          (:) ! root longevity (yrs)
      real(r8), allocatable :: clone_alloc        (:) ! fraction of carbon balance allocated to clonal reproduction.
      real(r8), allocatable :: seed_alloc         (:) ! fraction of carbon balance allocated to seeds.
+     real(r8), allocatable :: c2b                (:) ! Carbon to biomass multiplier [kg/kgC]
      real(r8), allocatable :: woody(:)
      real(r8), allocatable :: stress_decid(:)
      real(r8), allocatable :: season_decid(:)
@@ -98,7 +100,9 @@ module EDPftvarcon
                                                     ! Equation 16 Thonicke et al 2010
 
      ! Allometry Parameters
-     ! --------------------------------------------------------------------------------------------
+     ! -------------------------------------------------------------------------------------------- 
+     real(r8), allocatable :: allom_dbh_maxheight(:) ! dbh at which height growth ceases
+     
      real(r8), allocatable :: allom_hmode(:)        ! height allometry function type
      real(r8), allocatable :: allom_lmode(:)        ! maximum leaf allometry function type
      real(r8), allocatable :: allom_fmode(:)        ! maximum root allometry function type
@@ -117,7 +121,6 @@ module EDPftvarcon
      real(r8), allocatable :: allom_d2bl2(:)        ! Parameter 2 for d2bl allometry (slope)
      real(r8), allocatable :: allom_d2bl3(:)           ! Parameter 3 for d2bl allometry (optional)
      real(r8), allocatable :: allom_sai_scaler(:)      ! 
-     real(r8), allocatable :: allom_d2bl_slascaler(:)  ! 
      real(r8), allocatable :: allom_blca_expnt_diff(:) ! Any difference in the exponent between the leaf
                                                        ! biomass and crown area scaling
      real(r8), allocatable :: allom_d2ca_coefficient_max(:)  ! upper (savanna) value for crown area to dbh coefficient
@@ -252,7 +255,7 @@ subroutine Register_PFT(this, fates_params)
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
-    name = 'fates_max_dbh'
+    name = 'fates_dbh_repro_threshold'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
@@ -312,6 +315,10 @@ subroutine Register_PFT(this, fates_params)
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
+    name = 'fates_c2b'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
     name = 'fates_woody'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
@@ -428,6 +435,10 @@ subroutine Register_PFT(this, fates_params)
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
           dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
+    name = 'fates_allom_dbh_maxheight'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
     name = 'fates_allom_hmode'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
           dimension_names=dim_names, lower_bounds=dim_lower_bound)
@@ -500,10 +511,6 @@ subroutine Register_PFT(this, fates_params)
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
-    name = 'fates_allom_d2bl_slascaler'
-    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
-         dimension_names=dim_names, lower_bounds=dim_lower_bound)
-
     name = 'fates_allom_sai_scaler'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
@@ -648,9 +655,9 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=this%pft_used)
 
-    name = 'fates_max_dbh'
+    name = 'fates_dbh_repro_threshold'
     call fates_params%RetreiveParameterAllocate(name=name, &
-         data=this%max_dbh)
+         data=this%dbh_repro_threshold)
 
     name = 'fates_freezetol'
     call fates_params%RetreiveParameterAllocate(name=name, &
@@ -708,6 +715,10 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=this%seed_alloc)
 
+    name = 'fates_c2b'
+    call fates_params%RetreiveParameterAllocate(name=name, &
+         data=this%c2b)
+
     name = 'fates_woody'
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=this%woody)
@@ -820,6 +831,10 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=this%fire_alpha_SH)
 
+    name = 'fates_allom_dbh_maxheight'
+    call fates_params%RetreiveParameterAllocate(name=name, &
+          data=this%allom_dbh_maxheight)
+
     name = 'fates_allom_hmode'
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=this%allom_hmode)
@@ -896,10 +911,6 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=this%allom_d2ca_coefficient_min)
 
-    name = 'fates_allom_d2bl_slascaler'
-    call fates_params%RetreiveParameterAllocate(name=name, &
-         data=this%allom_d2bl_slascaler)
-
     name = 'fates_allom_sai_scaler'
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=this%allom_sai_scaler)
@@ -1378,7 +1389,8 @@ subroutine FatesReportPFTParams(is_master)
 
         write(fates_log(),*) '-----------  FATES PFT Parameters -----------------'
         write(fates_log(),fmt0) 'pft_used = ',EDPftvarcon_inst%pft_used
-        write(fates_log(),fmt0) 'max_dbh = ',EDPftvarcon_inst%max_dbh
+        write(fates_log(),fmt0) 'dbh max height = ',EDPftvarcon_inst%allom_dbh_maxheight
+        write(fates_log(),fmt0) 'dbh mature = ',EDPftvarcon_inst%dbh_repro_threshold
         write(fates_log(),fmt0) 'freezetol = ',EDPftvarcon_inst%freezetol
         write(fates_log(),fmt0) 'wood_density = ',EDPftvarcon_inst%wood_density
         write(fates_log(),fmt0) 'hgt_min = ',EDPftvarcon_inst%hgt_min
@@ -1396,6 +1408,7 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmt0) 'root_long = ',EDPftvarcon_inst%root_long
         write(fates_log(),fmt0) 'clone_alloc = ',EDPftvarcon_inst%clone_alloc
         write(fates_log(),fmt0) 'seed_alloc = ',EDPftvarcon_inst%seed_alloc
+        write(fates_log(),fmt0) 'C2B = ',EDPftvarcon_inst%c2b
         write(fates_log(),fmt0) 'woody = ',EDPftvarcon_inst%woody
         write(fates_log(),fmt0) 'stress_decid = ',EDPftvarcon_inst%stress_decid
         write(fates_log(),fmt0) 'season_decid = ',EDPftvarcon_inst%season_decid
@@ -1457,7 +1470,6 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmt0) 'allom_d2bl2 = ',EDPftvarcon_inst%allom_d2bl2
         write(fates_log(),fmt0) 'allom_d2bl3 = ',EDPftvarcon_inst%allom_d2bl3
         write(fates_log(),fmt0) 'allom_sai_scaler = ',EDPftvarcon_inst%allom_sai_scaler
-        write(fates_log(),fmt0) 'allom_d2bl_slascaler = ',EDPftvarcon_inst%allom_d2bl_slascaler
         write(fates_log(),fmt0) 'allom_blca_expnt_diff = ',EDPftvarcon_inst%allom_blca_expnt_diff
         write(fates_log(),fmt0) 'allom_d2ca_coefficient_max = ',EDPftvarcon_inst%allom_d2ca_coefficient_max
         write(fates_log(),fmt0) 'allom_d2ca_coefficient_min = ',EDPftvarcon_inst%allom_d2ca_coefficient_min        

main/FatesHistoryInterfaceMod.F90

@@ -3192,16 +3192,12 @@ subroutine define_history_vars(this, initialize_variables)
 
 
     !Logging
-    call this%set_history_var(vname='M7_SCPF', units = 'N/ha/event',          &
-          long='logging mortalities by pft/size',use_default='inactive', &
+    call this%set_history_var(vname='M7_SCPF', units = 'N/ha/event',               &
+          long='logging mortalities by pft/size',use_default='inactive',           &
           avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8,    &
           upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m7_si_scpf )
 
 
-
-
-
-
     call this%set_history_var(vname='MORTALITY_CANOPY_SCPF', units = 'N/ha/yr',          &
           long='total mortality of canopy plants by pft/size', use_default='inactive', &
           avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8,    &