diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90
index 4e86d22ae7..20cd97b64d 100644
--- a/biogeophys/FatesPlantHydraulicsMod.F90
+++ b/biogeophys/FatesPlantHydraulicsMod.F90
@@ -47,6 +47,7 @@ module FatesPlantHydraulicsMod
   use EDParamsMod       , only : hydr_kmax_rsurf2
   use EDParamsMod       , only : hydr_psi0
   use EDParamsMod       , only : hydr_psicap
+  use EDParamsMod       , only : hydr_htftype_node
   
   use EDTypesMod        , only : ed_site_type
   use EDTypesMod        , only : ed_patch_type
@@ -190,12 +191,24 @@ module FatesPlantHydraulicsMod
        __FILE__
 
 
-  integer, public, parameter :: van_genuchten_type      = 1
-  integer, public, parameter :: campbell_type           = 2
-  integer, public, parameter :: tfs_type                = 3
+  ! These index flags specify which pressure-volumen and pressure
+  ! conductivity relationship are available.
+  ! For plants: Users can option between useing tfs and van_genuchten
+  ! by specifying their choice in the parameter file,
+  ! with the model parameter hydr_htftype_node,
+  ! the value should be 1 for TFS or 2 for VG (as shown below).
+  ! Campbell, could technically be used, but the parameters for
+  ! that hypothesis are not in the parameter file, so it not currently available.
+  ! For soil: The soil hypothesis should follow the hypothesis for water transfer
+  ! in the Host Land Model. At this time campbell is the default for both
+  ! ELM and ALM. However, if alternatives arise (like VG), we still need to write
+  ! interface routines to transfer over parameters.  Right now we just hard-code
+  ! the use of campbell_type for the soil (see a few lines below).
+  
+  integer, public, parameter :: van_genuchten_type      = 2
+  integer, public, parameter :: campbell_type           = 3
+  integer, public, parameter :: tfs_type                = 1
   
-  integer, parameter :: plant_wrf_type = tfs_type
-  integer, parameter :: plant_wkf_type = tfs_type
   integer, parameter :: soil_wrf_type  = campbell_type
   integer, parameter :: soil_wkf_type  = campbell_type
   
@@ -5312,81 +5325,67 @@ subroutine InitHydroGlobals()
     ! Initialize the Water Retention Functions
     ! -----------------------------------------------------------------------------------
 
-    select case(plant_wrf_type)
-    case(van_genuchten_type)
-       do ft = 1,numpft
-            do pm = 1, n_plant_media
-                allocate(wrf_vg)
-                wrf_plant(pm,ft)%p => wrf_vg
-                call wrf_vg%set_wrf_param([alpha_vg, psd_vg, th_sat_vg, th_res_vg])
-            end do
-        end do
-     case(campbell_type)
-        do ft = 1,numpft
-           do pm = 1,n_plant_media
-              allocate(wrf_cch)
-              wrf_plant(pm,ft)%p => wrf_cch
-              call wrf_cch%set_wrf_param([EDPftvarcon_inst%hydr_thetas_node(ft,pm), &
-                                          EDPftvarcon_inst%hydr_pinot_node(ft,pm), &
-                                          9._r8])
-           end do
-        end do
-     case(tfs_type)
-        do ft = 1,numpft
-           do pm = 1,n_plant_media
-              allocate(wrf_tfs)
-              wrf_plant(pm,ft)%p => wrf_tfs
-
-              if (pm.eq.leaf_p_media) then   ! Leaf tissue
-                 cap_slp    = 0.0_r8
-                 cap_int    = 0.0_r8
-                 cap_corr   = 1.0_r8
-              else               ! Non leaf tissues
-                 cap_slp    = (hydr_psi0 - hydr_psicap )/(1.0_r8 - rwccap(pm))  
-                 cap_int    = -cap_slp + hydr_psi0    
-                 cap_corr   = -cap_int/cap_slp
-              end if
-              
-              call wrf_tfs%set_wrf_param([EDPftvarcon_inst%hydr_thetas_node(ft,pm), &
-                                          EDPftvarcon_inst%hydr_resid_node(ft,pm), &
-                                          EDPftvarcon_inst%hydr_pinot_node(ft,pm), &
-                                          EDPftvarcon_inst%hydr_epsil_node(ft,pm), &
-                                          rwcft(pm), & 
-                                          cap_corr, &
-                                          cap_int, &
-                                          cap_slp,real(pm,r8)])
+    do pm = 1, n_plant_media
+       select case(hydr_htftype_node(pm))
+       case(van_genuchten_type)
+          do ft = 1,numpft
+             allocate(wrf_vg)
+             wrf_plant(pm,ft)%p => wrf_vg
+             call wrf_vg%set_wrf_param([EDPftvarcon_inst%hydr_vg_alpha_node(ft,pm), &
+                                        EDPftvarcon_inst%hydr_vg_m_node(ft,pm), &
+                                        EDPftvarcon_inst%hydr_thetas_node(ft,pm), &
+                                        EDPftvarcon_inst%hydr_resid_node(ft,pm)])
+          end do
+       case(tfs_type)
+          do ft = 1,numpft
+             allocate(wrf_tfs)
+             wrf_plant(pm,ft)%p => wrf_tfs
+             if (pm.eq.leaf_p_media) then   ! Leaf tissue
+                cap_slp    = 0.0_r8
+                cap_int    = 0.0_r8
+                cap_corr   = 1.0_r8
+             else               ! Non leaf tissues
+                cap_slp    = (hydr_psi0 - hydr_psicap )/(1.0_r8 - rwccap(pm))  
+                cap_int    = -cap_slp + hydr_psi0    
+                cap_corr   = -cap_int/cap_slp
+             end if
+             call wrf_tfs%set_wrf_param([EDPftvarcon_inst%hydr_thetas_node(ft,pm), &
+                                         EDPftvarcon_inst%hydr_resid_node(ft,pm), &
+                                         EDPftvarcon_inst%hydr_pinot_node(ft,pm), &
+                                         EDPftvarcon_inst%hydr_epsil_node(ft,pm), &
+                                         rwcft(pm), & 
+                                         cap_corr, &
+                                         cap_int, &
+                                         cap_slp,real(pm,r8)])
            end do
-        end do
-
-    end select
+        end select
+     end do
 
     ! -----------------------------------------------------------------------------------
     ! Initialize the Water Conductance (K) Functions
     ! -----------------------------------------------------------------------------------
-
-    select case(plant_wkf_type)
-    case(van_genuchten_type)
-        do ft = 1,numpft
-            do pm = 1, n_plant_media
-                allocate(wkf_vg)
-                wkf_plant(pm,ft)%p => wkf_vg
-                call wkf_vg%set_wkf_param([alpha_vg, psd_vg, th_sat_vg, th_res_vg, tort_vg])
-            end do
-           
-        end do
-    case(campbell_type)
-        write(fates_log(),*) 'campbell/clapp-hornberger conductance not used in plants'
-        call endrun(msg=errMsg(sourcefile, __LINE__))
-    case(tfs_type)
-        do ft = 1,numpft
-            do pm = 1, n_plant_media
-               allocate(wkf_tfs)
-               wkf_plant(pm,ft)%p => wkf_tfs
-               call wkf_tfs%set_wkf_param([EDPftvarcon_inst%hydr_p50_node(ft,pm), &
-                                       EDPftvarcon_inst%hydr_avuln_node(ft,pm)])
-            end do
-        end do
-    end select
+    do pm = 1, n_plant_media
+       select case(hydr_htftype_node(pm))
+    
+       case(van_genuchten_type)
+          do ft = 1,numpft
+             allocate(wkf_vg)
+             wkf_plant(pm,ft)%p => wkf_vg
+             call wkf_vg%set_wkf_param([EDPftvarcon_inst%hydr_vg_alpha_node(ft,pm), &
+                                        EDPftvarcon_inst%hydr_vg_m_node(ft,pm), &
+                                        EDPftvarcon_inst%hydr_thetas_node(ft,pm), &
+                                        EDPftvarcon_inst%hydr_resid_node(ft,pm), &
+                                        tort_vg])
+          end do
+       case(tfs_type)
+          do ft = 1,numpft
+             allocate(wkf_tfs)
+             wkf_plant(pm,ft)%p => wkf_tfs
+             call wkf_tfs%set_wkf_param([EDPftvarcon_inst%hydr_p50_node(ft,pm), &
+                  EDPftvarcon_inst%hydr_avuln_node(ft,pm)])
+          end do
+       end select
+    end do
 
     ! There is only 1 stomata conductance hypothesis which uses the p50 and 
     ! vulnerability parameters
diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90
index 8b919f52c1..349f6473bf 100644
--- a/biogeophys/FatesPlantRespPhotosynthMod.F90
+++ b/biogeophys/FatesPlantRespPhotosynthMod.F90
@@ -865,6 +865,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
     ! ------------------------------------------------------------------------------------
     
     use EDPftvarcon       , only : EDPftvarcon_inst
+    use EDParamsMod       , only : theta_cj_c3, theta_cj_c4
 
     
     ! Arguments
@@ -962,11 +963,6 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
    ! quantum efficiency, used only for C4 (mol CO2 / mol photons) (index 0)
    real(r8),parameter,dimension(0:1) :: quant_eff = [0.05_r8,0.0_r8]
 
-   ! empirical curvature parameter for ac, aj photosynthesis co-limitation. 
-   ! Changed theta_cj and theta_ip to 0.999 to effectively remove smoothing logic 
-   ! following Anthony Walker's findings from MAAT. 
-   real(r8),parameter,dimension(0:1) :: theta_cj  = [0.999_r8,0.999_r8]
-
    ! empirical curvature parameter for ap photosynthesis co-limitation
    real(r8),parameter :: theta_ip = 0.999_r8
 
@@ -1063,7 +1059,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
                           (4._r8*co2_inter_c+8._r8*co2_cpoint)
 
                     ! Gross photosynthesis smoothing calculations. Co-limit ac and aj.                                                                                             
-                    aquad = theta_cj(c3c4_path_index)
+                    aquad = theta_cj_c3
                     bquad = -(ac + aj)
                     cquad = ac * aj
                     call quadratic_f (aquad, bquad, cquad, r1, r2)
@@ -1094,7 +1090,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
 
                     ! Gross photosynthesis smoothing calculations. First co-limit ac and aj. Then co-limit ap                                                                     
 
-                    aquad = theta_cj(c3c4_path_index)
+                    aquad = theta_cj_c4
                     bquad = -(ac + aj)
                     cquad = ac * aj
                     call quadratic_f (aquad, bquad, cquad, r1, r2)
diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90
index 8162939bc3..1f10aa2c7f 100644
--- a/main/EDParamsMod.F90
+++ b/main/EDParamsMod.F90
@@ -20,6 +20,22 @@ module EDParamsMod
    !
    ! this is what the user can use for the actual values
    !
+
+   real(r8),protected, public :: vai_top_bin_width           ! width in VAI units of uppermost leaf+stem
+                                                             ! layer scattering element in each canopy layer [m2/m2]
+                                                             ! (NOT YET IMPLEMENTED)
+   real(r8),protected, public :: vai_width_increase_factor   ! factor by which each leaf+stem scattering element
+                                                             ! increases in VAI width (1 = uniform spacing)
+                                                             ! (NOT YET IMPLEMENTED)
+   real(r8),protected, public :: photo_temp_acclim_timescale ! Length of the window for the exponential moving average (ema)
+                                                             ! of vegetation temperature used in photosynthesis
+                                                             ! temperature acclimation (NOT YET IMPLEMENTED)
+
+   integer,protected, public :: maintresp_model       ! switch for choosing between leaf maintenance
+                                                      ! respiration model. 1=Ryan (1991) (NOT YET IMPLEMENTED)
+   integer,protected, public :: photo_tempsens_model  ! switch for choosing the model that defines the temperature
+                                                      ! sensitivity of photosynthetic parameters (vcmax, jmax).
+                                                      ! 1=non-acclimating (NOT YET IMPLEMENTED)
    
    real(r8),protected, public :: fates_mortality_disturbance_fraction ! the fraction of canopy mortality that results in disturbance
    real(r8),protected, public :: ED_val_comp_excln
@@ -50,16 +66,37 @@ module EDParamsMod
 
    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"
+
+   ! empirical curvature parameters for ac, aj photosynthesis co-limitation, c3 and c4 plants respectively
+   real(r8),protected,public  :: theta_cj_c3
+   real(r8),protected,public  :: theta_cj_c4
    
    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)
 
-   ! two special parameters whose size is defined in the parameter file
+   ! Unassociated pft dimensioned free parameter that developers can use for testing arbitrary new hypotheses
+   ! (THIS PARAMETER IS UNUSED, FEEL FREE TO USE IT FOR WHATEVER PURPOSE YOU LIKE. WE CAN
+   !  HELP MIGRATE YOUR USAGE OF THE PARMETER TO A PERMANENT HOME LATER)
+   real(r8),protected,public  :: dev_arbitrary 
+   character(len=param_string_length),parameter,public :: name_dev_arbitrary = "fates_dev_arbitrary"
+   
+   ! parameters whose size is defined in the parameter file
    real(r8),protected,allocatable,public :: ED_val_history_sizeclass_bin_edges(:)
    real(r8),protected,allocatable,public :: ED_val_history_ageclass_bin_edges(:)
    real(r8),protected,allocatable,public :: ED_val_history_height_bin_edges(:)
    real(r8),protected,allocatable,public :: ED_val_history_coageclass_bin_edges(:)
+
+   ! Switch that defines the current pressure-volume and pressure-conductivity model
+   ! to be used at each node (compartment/organ)
+   ! 1  = Christofferson et al. 2016 (TFS),   2 = Van Genuchten 1980
+   integer, protected,allocatable,public :: hydr_htftype_node(:) 
    
+   character(len=param_string_length),parameter,public :: ED_name_vai_top_bin_width = "fates_vai_top_bin_width"
+   character(len=param_string_length),parameter,public :: ED_name_vai_width_increase_factor = "fates_vai_width_increase_factor"
+   character(len=param_string_length),parameter,public :: ED_name_photo_temp_acclim_timescale = "fates_photo_temp_acclim_timescale"
+   character(len=param_string_length),parameter,public :: name_photo_tempsens_model = "fates_photo_tempsens_model"
+   character(len=param_string_length),parameter,public :: name_maintresp_model = "fates_maintresp_model"
+   character(len=param_string_length),parameter,public :: ED_name_hydr_htftype_node = "fates_hydr_htftype_node"
    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"
@@ -83,13 +120,12 @@ module EDParamsMod
    character(len=param_string_length),parameter,public :: ED_name_canopy_closure_thresh= "fates_canopy_closure_thresh"      
    character(len=param_string_length),parameter,public :: ED_name_stomatal_model= "fates_leaf_stomatal_model"
 
-   ! Resistance to active crown fire
-  
-
+   character(len=param_string_length),parameter,public :: name_theta_cj_c3 = "fates_theta_cj_c3"
+   character(len=param_string_length),parameter,public :: name_theta_cj_c4 = "fates_theta_cj_c4"
+   
    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"
 
-
    ! non-scalar parameter names
    character(len=param_string_length),parameter,public :: ED_name_history_sizeclass_bin_edges= "fates_history_sizeclass_bin_edges"      
    character(len=param_string_length),parameter,public :: ED_name_history_ageclass_bin_edges= "fates_history_ageclass_bin_edges"      
@@ -173,6 +209,11 @@ subroutine FatesParamsInit()
 
     implicit none
 
+    vai_top_bin_width                     = nan
+    vai_width_increase_factor             = nan
+    photo_temp_acclim_timescale           = nan
+    photo_tempsens_model                  = -9
+    maintresp_model                       = -9
     fates_mortality_disturbance_fraction  = nan
     ED_val_comp_excln                     = nan
     ED_val_init_litter                    = nan
@@ -211,7 +252,9 @@ subroutine FatesParamsInit()
     eca_plant_escalar                     = nan
     q10_mr                                = nan
     q10_froz                              = nan
-
+    theta_cj_c3                           = nan
+    theta_cj_c4                           = nan
+    dev_arbitrary                         = nan
   end subroutine FatesParamsInit
 
   !-----------------------------------------------------------------------
@@ -222,7 +265,7 @@ subroutine FatesRegisterParams(fates_params)
 
     use FatesParametersInterface, only : fates_parameters_type, dimension_name_scalar, dimension_shape_1d
     use FatesParametersInterface, only : dimension_name_history_size_bins, dimension_name_history_age_bins
-    use FatesParametersInterface, only : dimension_name_history_height_bins
+    use FatesParametersInterface, only : dimension_name_history_height_bins, dimension_name_hydr_organs
     use FatesParametersInterface, only : dimension_name_history_coage_bins
     use FatesParametersInterface, only : dimension_shape_scalar
 
@@ -236,10 +279,31 @@ subroutine FatesRegisterParams(fates_params)
     character(len=param_string_length), parameter :: dim_names_ageclass(1) = (/dimension_name_history_age_bins/)
     character(len=param_string_length), parameter :: dim_names_height(1) = (/dimension_name_history_height_bins/)
     character(len=param_string_length), parameter :: dim_names_coageclass(1) = (/dimension_name_history_coage_bins/)
-
+    character(len=param_string_length), parameter :: dim_names_hydro_organs(1) = (/dimension_name_hydr_organs/)
        
     call FatesParamsInit()
 
+    call fates_params%RegisterParameter(name=ED_name_vai_top_bin_width, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
+    call fates_params%RegisterParameter(name=ED_name_vai_width_increase_factor, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+    
+    call fates_params%RegisterParameter(name=ED_name_photo_temp_acclim_timescale, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
+    call fates_params%RegisterParameter(name=name_photo_tempsens_model,dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
+    call fates_params%RegisterParameter(name=name_maintresp_model,dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+    
+    call fates_params%RegisterParameter(name=name_theta_cj_c3, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+    
+    call fates_params%RegisterParameter(name=name_theta_cj_c4, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+    
     call fates_params%RegisterParameter(name=ED_name_mort_disturb_frac, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
@@ -357,7 +421,14 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=fates_name_q10_froz, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
+    call fates_params%RegisterParameter(name=name_dev_arbitrary, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+    
     ! non-scalar parameters
+
+    call fates_params%RegisterParameter(name=ED_name_hydr_htftype_node, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names_hydro_organs)
+    
     call fates_params%RegisterParameter(name=ED_name_history_sizeclass_bin_edges, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names_sizeclass)
 
@@ -390,6 +461,24 @@ subroutine FatesReceiveParams(fates_params)
     class(fates_parameters_type), intent(inout) :: fates_params
 
     real(r8) :: tmpreal ! local real variable for changing type on read
+    real(r8), allocatable :: hydr_htftype_real(:)
+    
+    call fates_params%RetreiveParameter(name=ED_name_vai_top_bin_width, &
+         data=vai_top_bin_width)
+
+    call fates_params%RetreiveParameter(name=ED_name_vai_width_increase_factor, &
+         data=vai_width_increase_factor)
+
+    call fates_params%RetreiveParameter(name=ED_name_photo_temp_acclim_timescale, &
+         data=photo_temp_acclim_timescale)
+
+    call fates_params%RetreiveParameter(name=name_photo_tempsens_model, &
+         data=tmpreal)
+    photo_tempsens_model = nint(tmpreal)
+
+    call fates_params%RetreiveParameter(name=name_maintresp_model, &
+         data=tmpreal)
+    maintresp_model = nint(tmpreal)
     
     call fates_params%RetreiveParameter(name=ED_name_mort_disturb_frac, &
           data=fates_mortality_disturbance_fraction)
@@ -503,11 +592,20 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetreiveParameter(name=eca_name_plant_escalar, &
           data=eca_plant_escalar)
 
+    call fates_params%RetreiveParameter(name=name_theta_cj_c3, &
+          data=theta_cj_c3)
+
+     call fates_params%RetreiveParameter(name=name_theta_cj_c4, &
+          data=theta_cj_c4)
+     
     call fates_params%RetreiveParameter(name=fates_name_q10_mr, &
           data=q10_mr)
     
     call fates_params%RetreiveParameter(name=fates_name_q10_froz, &
-          data=q10_froz)
+         data=q10_froz)
+    
+    call fates_params%RetreiveParameter(name=name_dev_arbitrary, &
+         data=dev_arbitrary)
 
     call fates_params%RetreiveParameter(name=fates_name_active_crown_fire, & 
           data=tmpreal)
@@ -529,6 +627,11 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetreiveParameterAllocate(name=ED_name_history_coageclass_bin_edges, &
          data=ED_val_history_coageclass_bin_edges)
 
+    call fates_params%RetreiveParameterAllocate(name=ED_name_hydr_htftype_node, &
+         data=hydr_htftype_real)
+    allocate(hydr_htftype_node(size(hydr_htftype_real)))
+    hydr_htftype_node(:) = nint(hydr_htftype_real(:))
+    deallocate(hydr_htftype_real)
 
   end subroutine FatesReceiveParams
   
@@ -539,11 +642,16 @@ subroutine FatesReportParams(is_master)
      logical,intent(in) :: is_master
 
      character(len=32),parameter :: fmt0 = '(a,(F12.4))'
+     character(len=32),parameter :: fmti = '(a,(I4))'
      logical, parameter :: debug_report = .false.
      
      if(debug_report .and. is_master) then
         
         write(fates_log(),*) '-----------  FATES Scalar Parameters -----------------'
+        write(fates_log(),fmt0) 'vai_top_bin_width = ',vai_top_bin_width
+        write(fates_log(),fmt0) 'vai_width_increase_factor = ',vai_width_increase_factor
+        write(fates_log(),fmt0) 'photo_temp_acclim_timescale = ',photo_temp_acclim_timescale
+        write(fates_log(),fmti) 'hydr_htftype_node = ',hydr_htftype_node
         write(fates_log(),fmt0) 'fates_mortality_disturbance_fraction = ',fates_mortality_disturbance_fraction
         write(fates_log(),fmt0) 'ED_val_comp_excln = ',ED_val_comp_excln
         write(fates_log(),fmt0) 'ED_val_init_litter = ',ED_val_init_litter
diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90
index 61f095a758..b6f17bdd70 100644
--- a/main/EDPftvarcon.F90
+++ b/main/EDPftvarcon.F90
@@ -15,7 +15,6 @@ module EDPftvarcon
   use FatesGlobals,   only : fates_log
   use FatesGlobals,   only : endrun => fates_endrun
   use FatesLitterMod, only : ilabile,icellulose,ilignin
-  use PRTGenericMod,  only : num_organ_types
   use PRTGenericMod,  only : leaf_organ, fnrt_organ, store_organ
   use PRTGenericMod,  only : sapw_organ, struct_organ, repro_organ
   use PRTGenericMod,  only : prt_cnp_flex_allom_hyp,prt_carbon_allom_hyp
@@ -42,7 +41,7 @@ module EDPftvarcon
 
   !ED specific variables. 
   type, public ::  EDPftvarcon_type
-
+      
      real(r8), allocatable :: freezetol(:)           ! minimum temperature tolerance
      real(r8), allocatable :: hgt_min(:)             ! sapling height m
      real(r8), allocatable :: dleaf(:)               ! leaf characteristic dimension length (m)
@@ -180,6 +179,11 @@ module EDPftvarcon
      real(r8), allocatable :: prescribed_puptake(:)   ! If there is no soil BGC model active,
                                                       ! prescribe an uptake rate for phosphorus
                                                       ! This is the fraction of plant demand
+
+
+     ! Unassociated pft dimensioned free parameter that
+     ! developers can use for testing arbitrary new hypothese
+     real(r8), allocatable :: dev_arbitrary_pft(:) 
      
      ! Parameters dimensioned by PFT and leaf age
      real(r8), allocatable :: vcmax25top(:,:)             ! maximum carboxylation rate of Rub. at 25C, 
@@ -189,23 +193,41 @@ module EDPftvarcon
      ! ---------------------------------------------------------------------------------------------
 
      ! PFT Dimension
-     real(r8), allocatable :: hydr_p_taper(:)       ! xylem taper exponent
-     real(r8), allocatable :: hydr_rs2(:)           ! absorbing root radius (m)
-     real(r8), allocatable :: hydr_srl(:)           ! specific root length (m g-1)
-     real(r8), allocatable :: hydr_rfrac_stem(:)    ! fraction of total tree resistance from troot to canopy
-     real(r8), allocatable :: hydr_avuln_gs(:)      ! shape parameter for stomatal control of water vapor exiting leaf 
-     real(r8), allocatable :: hydr_p50_gs(:)        ! water potential at 50% loss of stomatal conductance
-
+     real(r8), allocatable :: hydr_p_taper(:)    ! xylem taper exponent
+     real(r8), allocatable :: hydr_rs2(:)        ! absorbing root radius (m)
+     real(r8), allocatable :: hydr_srl(:)        ! specific root length (m g-1)
+     real(r8), allocatable :: hydr_rfrac_stem(:) ! fraction of total tree resistance from troot to canopy
+     real(r8), allocatable :: hydr_avuln_gs(:)   ! shape parameter for stomatal control of water vapor exiting leaf 
+     real(r8), allocatable :: hydr_p50_gs(:)     ! water potential at 50% loss of stomatal conductance
+     real(r8), allocatable :: hydr_k_lwp(:)      ! inner leaf humidity scaling coefficient
+     
      ! PFT x Organ Dimension  (organs are: 1=leaf, 2=stem, 3=transporting root, 4=absorbing root)
+     ! ----------------------------------------------------------------------------------
+
+     ! Van Genuchten PV PK curves  (NOT IMPLEMENTED)
+     real(r8), allocatable :: hydr_vg_alpha_node(:,:)   ! capilary length parameter in van Genuchten model
+     real(r8), allocatable :: hydr_vg_m_node(:,:)       ! pore size distribution, m in van Genuchten 1980 model, range (0,1)
+     real(r8), allocatable :: hydr_vg_n_node(:,:)       ! pore size distribution, n in van Genuchten 1980 model, range >2
+
+     ! TFS PV-PK curves
      real(r8), allocatable :: hydr_avuln_node(:,:)  ! xylem vulernability curve shape parameter 
      real(r8), allocatable :: hydr_p50_node(:,:)    ! xylem water potential at 50% conductivity loss (MPa)
-     real(r8), allocatable :: hydr_thetas_node(:,:) ! saturated water content (cm3/cm3)
      real(r8), allocatable :: hydr_epsil_node(:,:)  ! bulk elastic modulus (MPa)
      real(r8), allocatable :: hydr_pitlp_node(:,:)  ! turgor loss point (MPa)
-     real(r8), allocatable :: hydr_resid_node(:,:)  ! residual fraction (fraction)
      real(r8), allocatable :: hydr_fcap_node(:,:)   ! fraction of (1-resid_node) that is capillary in source
      real(r8), allocatable :: hydr_pinot_node(:,:)  ! osmotic potential at full turgor
      real(r8), allocatable :: hydr_kmax_node(:,:)   ! maximum xylem conductivity per unit conducting xylem area
+
+     ! Parameters for both VG and TFS PV-PK curves
+     real(r8), allocatable :: hydr_resid_node(:,:)  ! residual fraction (fraction)
+     real(r8), allocatable :: hydr_thetas_node(:,:) ! saturated water content (cm3/cm3)
+
+
+     ! Table that maps HLM pfts to FATES pfts for fixed biogeography mode
+     ! The values are area fractions (NOT IMPLEMENTED)
+     real(r8), allocatable :: hlm_pft_map(:,:)
+
+     
      
    contains
      procedure, public :: Init => EDpftconInit
@@ -284,16 +306,19 @@ subroutine Register_PFT(this, fates_params)
 
     use FatesParametersInterface, only : fates_parameters_type, param_string_length
     use FatesParametersInterface, only : dimension_name_pft, dimension_shape_1d
-
+    use FatesParametersInterface, only : dimension_name_hlm_pftno, dimension_shape_2d
+    
     implicit none
 
     class(EDPftvarcon_type), intent(inout) :: this
     class(fates_parameters_type), intent(inout) :: fates_params
 
     character(len=param_string_length), parameter :: dim_names(1) = (/dimension_name_pft/)
-
+    character(len=param_string_length) :: pftmap_dim_names(2)
+    
     integer, parameter :: dim_lower_bound(1) = (/ lower_bound_pft /)
-
+    
+    
     character(len=param_string_length) :: name
 
     !X!    name = ''
@@ -444,6 +469,10 @@ subroutine Register_PFT(this, fates_params)
     name = 'fates_hydr_p50_gs'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
           dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
+    name = 'fates_hydr_k_lwp'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
     
     name = 'fates_mort_bmort'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
@@ -603,7 +632,19 @@ subroutine Register_PFT(this, fates_params)
     name = 'fates_prescribed_puptake'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
           dimension_names=dim_names, lower_bounds=dim_lower_bound)    
-  
+
+    name = 'fates_dev_arbitrary_pft'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+          dimension_names=dim_names, lower_bounds=dim_lower_bound)
+    
+    ! adding the hlm_pft_map variable with two dimensions - FATES PFTno and HLM PFTno
+    pftmap_dim_names(1) = dimension_name_pft 
+    pftmap_dim_names(2) = dimension_name_hlm_pftno 
+
+    name = 'fates_hlm_pft_map'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, &
+         dimension_names=pftmap_dim_names, lower_bounds=dim_lower_bound)
+    
   end subroutine Register_PFT
 
   !-----------------------------------------------------------------------
@@ -766,6 +807,10 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetreiveParameterAllocate(name=name, &
           data=this%hydr_p50_gs)
 
+    name = 'fates_hydr_k_lwp'
+    call fates_params%RetreiveParameterAllocate(name=name, &
+         data=this%hydr_k_lwp)
+    
     name = 'fates_mort_bmort'
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=this%bmort)
@@ -887,6 +932,10 @@ subroutine Receive_PFT(this, fates_params)
     name = 'fates_prescribed_puptake'
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=this%prescribed_puptake)
+
+    name = 'fates_dev_arbitrary_pft'
+    call fates_params%RetreiveParameterAllocate(name=name, &
+         data=this%dev_arbitrary_pft)
     
     name = 'fates_eca_decompmicc'
     call fates_params%RetreiveParameterAllocate(name=name, &
@@ -932,6 +981,10 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=this%eca_lambda_ptase)
 
+    name = 'fates_hlm_pft_map' 
+    call fates_params%RetreiveParameterAllocate(name=name, &
+         data=this%hlm_pft_map)
+    
   end subroutine Receive_PFT
 
   !-----------------------------------------------------------------------
@@ -1172,6 +1225,18 @@ subroutine Register_PFT_hydr_organs(this, fates_params)
     dim_names(1) = dimension_name_pft
     dim_names(2) = dimension_name_hydr_organs
 
+    name = 'fates_hydr_vg_alpha_node'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+    
+    name = 'fates_hydr_vg_m_node'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+    
+    name = 'fates_hydr_vg_n_node'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, &
+          dimension_names=dim_names, lower_bounds=dim_lower_bound)
+    
     name = 'fates_hydr_avuln_node'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, &
           dimension_names=dim_names, lower_bounds=dim_lower_bound)
@@ -1207,8 +1272,19 @@ subroutine Register_PFT_hydr_organs(this, fates_params)
     name = 'fates_hydr_kmax_node'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, &
           dimension_names=dim_names, lower_bounds=dim_lower_bound)
-    
 
+    name = 'fates_hydr_vg_alpha_node'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, &
+          dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
+    name = 'fates_hydr_vg_m_node'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, &
+          dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
+    name = 'fates_hydr_vg_n_node'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+    
   end subroutine Register_PFT_hydr_organs
 
   !-----------------------------------------------------------------------
@@ -1224,6 +1300,19 @@ subroutine Receive_PFT_hydr_organs(this, fates_params)
      class(fates_parameters_type), intent(inout) :: fates_params
      
      character(len=param_string_length) :: name
+
+
+     name = 'fates_hydr_vg_alpha_node'
+     call fates_params%RetreiveParameterAllocate(name=name, &
+          data=this%hydr_vg_alpha_node)
+
+     name = 'fates_hydr_vg_m_node'
+     call fates_params%RetreiveParameterAllocate(name=name, &
+           data=this%hydr_vg_m_node)
+
+     name = 'fates_hydr_vg_n_node'
+     call fates_params%RetreiveParameterAllocate(name=name, &
+          data=this%hydr_vg_n_node)
      
      name = 'fates_hydr_avuln_node'
      call fates_params%RetreiveParameterAllocate(name=name, &
@@ -1261,6 +1350,18 @@ subroutine Receive_PFT_hydr_organs(this, fates_params)
      call fates_params%RetreiveParameterAllocate(name=name, &
            data=this%hydr_kmax_node)
 
+     name = 'fates_hydr_vg_alpha_node'
+     call fates_params%RetreiveParameterAllocate(name=name, &
+          data=this%hydr_vg_alpha_node)
+     
+     name = 'fates_hydr_vg_m_node'
+     call fates_params%RetreiveParameterAllocate(name=name, &
+          data=this%hydr_vg_m_node)
+     
+     name = 'fates_hydr_vg_n_node'
+     call fates_params%RetreiveParameterAllocate(name=name, &
+          data=this%hydr_vg_n_node)
+     
   end subroutine Receive_PFT_hydr_organs
 
   ! ===============================================================================================
@@ -1347,6 +1448,7 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmt0) 'hydr_rfrac_stem = ',EDPftvarcon_inst%hydr_rfrac_stem
         write(fates_log(),fmt0) 'hydr_avuln_gs = ',EDPftvarcon_inst%hydr_avuln_gs
         write(fates_log(),fmt0) 'hydr_p50_gs = ',EDPftvarcon_inst%hydr_p50_gs
+        write(fates_log(),fmt0) 'hydr_k_lwp = ',EDPftvarcon_inst%hydr_k_lwp
         write(fates_log(),fmt0) 'hydr_avuln_node = ',EDPftvarcon_inst%hydr_avuln_node
         write(fates_log(),fmt0) 'hydr_p50_node = ',EDPftvarcon_inst%hydr_p50_node
         write(fates_log(),fmt0) 'hydr_thetas_node = ',EDPftvarcon_inst%hydr_thetas_node 
@@ -1356,6 +1458,9 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmt0) 'hydr_fcap_node = ',EDPftvarcon_inst%hydr_fcap_node
         write(fates_log(),fmt0) 'hydr_pinot_node = ',EDPftvarcon_inst%hydr_pinot_node
         write(fates_log(),fmt0) 'hydr_kmax_node = ',EDPftvarcon_inst%hydr_kmax_node
+        write(fates_log(),fmt0) 'hydr_vg_alpha_node  = ',EDPftvarcon_inst%hydr_vg_alpha_node
+        write(fates_log(),fmt0) 'hydr_vg_m_node  = ',EDPftvarcon_inst%hydr_vg_m_node
+        write(fates_log(),fmt0) 'hydr_vg_n_node  = ',EDPftvarcon_inst%hydr_vg_n_node
         write(fates_log(),*) '-------------------------------------------------'
 
      end if
diff --git a/main/FatesParametersInterface.F90 b/main/FatesParametersInterface.F90
index ebaad3fa7c..f69d4ef5bf 100644
--- a/main/FatesParametersInterface.F90
+++ b/main/FatesParametersInterface.F90
@@ -35,6 +35,7 @@ module FatesParametersInterface
   character(len=*), parameter, public :: dimension_name_history_age_bins = 'fates_history_age_bins'
   character(len=*), parameter, public :: dimension_name_history_height_bins = 'fates_history_height_bins'
   character(len=*), parameter, public :: dimension_name_history_coage_bins = 'fates_history_coage_bins'
+  character(len=*), parameter, public :: dimension_name_hlm_pftno = 'fates_hlm_pftno'
   
   ! Dimensions in the host namespace:
   character(len=*), parameter, public :: dimension_name_host_allpfts = 'allpfts'
diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl
index 86d46710da..674c0a5ee7 100644
--- a/parameter_files/fates_params_default.cdl
+++ b/parameter_files/fates_params_default.cdl
@@ -1,4 +1,4 @@
-netcdf fates_params_default {
+netcdf fates_params_default.c210629_sorted {
 dimensions:
 	fates_NCWD = 4 ;
 	fates_history_age_bins = 7 ;
@@ -11,6 +11,7 @@ dimensions:
 	fates_pft = 12 ;
 	fates_prt_organs = 4 ;
 	fates_string_length = 60 ;
+	fates_hlm_pftno = 14 ;
 variables:
 	double fates_history_ageclass_bin_edges(fates_history_age_bins) ;
 		fates_history_ageclass_bin_edges:units = "yr" ;
@@ -24,15 +25,25 @@ variables:
 	double fates_history_sizeclass_bin_edges(fates_history_size_bins) ;
 		fates_history_sizeclass_bin_edges:units = "cm" ;
 		fates_history_sizeclass_bin_edges:long_name = "Lower edges for DBH size class bins used in size-resolved cohort history output" ;
+	double fates_hydr_htftype_node(fates_hydr_organs) ;
+		fates_hydr_htftype_node:units = "unitless" ;
+		fates_hydr_htftype_node:long_name = "Switch that defines the hydraulic transfer functions for each organ." ;
+		fates_hydr_htftype_node:possible_values = "1: Christofferson et al. 2016 (TFS); 2: Van Genuchten 1980" ;
+	double fates_prt_organ_id(fates_prt_organs) ;
+		fates_prt_organ_id:units = "index, unitless" ;
+		fates_prt_organ_id:long_name = "This is the global index the organ in this file is associated with in PRTGenericMod.F90" ;
 	char fates_pftname(fates_pft, fates_string_length) ;
 		fates_pftname:units = "unitless - string" ;
 		fates_pftname:long_name = "Description of plant type" ;
+	char fates_hydr_organname_node(fates_hydr_organs, fates_string_length) ;
+		fates_hydr_organname_node:units = "unitless - string" ;
+		fates_hydr_organname_node:long_name = "Name of plant hydraulics organs (DONT CHANGE, order matches media list in FatesHydraulicsMemMod.F90)" ;
+	char fates_litterclass_name(fates_litterclass, fates_string_length) ;
+		fates_litterclass_name:units = "unitless - string" ;
+		fates_litterclass_name:long_name = "Name of the litter classes, for variables associated with dimension fates_litterclass" ;
 	char fates_prt_organ_name(fates_prt_organs, fates_string_length) ;
 		fates_prt_organ_name:units = "unitless - string" ;
 		fates_prt_organ_name:long_name = "Name of plant organs (order must match PRTGenericMod.F90)" ;
-   double fates_prt_organ_id(fates_prt_organs) ;
-	   fates_prt_organ_id:units = "index, unitless" ;
-		fates_prt_organ_id:long_name = "This is the global index the organ in this file is associated with in PRTGenericMod.F90" ;
 	double fates_alloc_storage_cushion(fates_pft) ;
 		fates_alloc_storage_cushion:units = "fraction" ;
 		fates_alloc_storage_cushion:long_name = "maximum size of storage C pool, relative to maximum size of leaf C pool" ;
@@ -124,12 +135,30 @@ variables:
 		fates_allom_stmode:units = "index" ;
 		fates_allom_stmode:long_name = "storage allometry function index." ;
 		fates_allom_stmode:possible_values = "1: target storage proportional to trimmed maximum leaf biomass." ;
+	double fates_allom_zroot_k(fates_pft) ;
+		fates_allom_zroot_k:units = "unitless" ;
+		fates_allom_zroot_k:long_name = "scale coefficient of logistic rooting depth model (NOT USED)" ;
+	double fates_allom_zroot_max_dbh(fates_pft) ;
+		fates_allom_zroot_max_dbh:units = "cm" ;
+		fates_allom_zroot_max_dbh:long_name = "dbh at which a plant reaches the maximum value for its maximum rooting depth (NOT USED)" ;
+	double fates_allom_zroot_max_z(fates_pft) ;
+		fates_allom_zroot_max_z:units = "m" ;
+		fates_allom_zroot_max_z:long_name = "the maximum rooting depth defined at dbh = fates_allom_zroot_max_dbh (NOT USED). note: max_z=min_z=large, sets rooting depth to soil depth" ;
+	double fates_allom_zroot_min_dbh(fates_pft) ;
+		fates_allom_zroot_min_dbh:units = "cm" ;
+		fates_allom_zroot_min_dbh:long_name = "dbh at which the maximum rooting depth for a recruit is defined (NOT USED)" ;
+	double fates_allom_zroot_min_z(fates_pft) ;
+		fates_allom_zroot_min_z:units = "m" ;
+		fates_allom_zroot_min_z:long_name = "the maximum rooting depth defined at dbh = fates_allom_zroot_min_dbh (NOT USED) note: max_z=min_z=large, sets rooting depth to soil depth" ;
 	double fates_branch_turnover(fates_pft) ;
 		fates_branch_turnover:units = "yr" ;
 		fates_branch_turnover:long_name = "turnover time of branches" ;
 	double fates_c2b(fates_pft) ;
 		fates_c2b:units = "ratio" ;
 		fates_c2b:long_name = "Carbon to biomass multiplier of bulk structural tissues" ;
+	double fates_dev_arbitrary_pft(fates_pft) ;
+		fates_dev_arbitrary_pft:units = "unknown" ;
+		fates_dev_arbitrary_pft:long_name = "Unassociated pft dimensioned free parameter that developers can use for testing arbitrary new hypotheses" ;
 	double fates_displar(fates_pft) ;
 		fates_displar:units = "unitless" ;
 		fates_displar:long_name = "Ratio of displacement height to canopy top height" ;
@@ -211,6 +240,10 @@ variables:
 	double fates_hydr_fcap_node(fates_hydr_organs, fates_pft) ;
 		fates_hydr_fcap_node:units = "unitless" ;
 		fates_hydr_fcap_node:long_name = "fraction of non-residual water that is capillary in source" ;
+	double fates_hydr_k_lwp(fates_pft) ;
+		fates_hydr_k_lwp:units = "unitless" ;
+		fates_hydr_k_lwp:long_name = "inner leaf humidity scaling coefficient" ;
+		fates_hydr_k_lwp:possible_values = "0: turns off leaf humidity effects on conductance. 1-10 activates humidity effects" ;
 	double fates_hydr_kmax_node(fates_hydr_organs, fates_pft) ;
 		fates_hydr_kmax_node:units = "kg/MPa/m/s" ;
 		fates_hydr_kmax_node:long_name = "maximum xylem conductivity per unit conducting xylem area" ;
@@ -244,6 +277,15 @@ variables:
 	double fates_hydr_thetas_node(fates_hydr_organs, fates_pft) ;
 		fates_hydr_thetas_node:units = "cm3/cm3" ;
 		fates_hydr_thetas_node:long_name = "saturated water content" ;
+	double fates_hydr_vg_alpha_node(fates_hydr_organs, fates_pft) ;
+		fates_hydr_vg_alpha_node:units = "MPa-1" ;
+		fates_hydr_vg_alpha_node:long_name = "(used if hydr_htftype_node = 2), capillary length parameter in van Genuchten model" ;
+	double fates_hydr_vg_m_node(fates_hydr_organs, fates_pft) ;
+		fates_hydr_vg_m_node:units = "unitless" ;
+		fates_hydr_vg_m_node:long_name = "(used if hydr_htftype_node = 2),m in van Genuchten 1980 model, 2nd pore size distribution parameter" ;
+	double fates_hydr_vg_n_node(fates_hydr_organs, fates_pft) ;
+		fates_hydr_vg_n_node:units = "unitless" ;
+		fates_hydr_vg_n_node:long_name = "(used if hydr_htftype_node = 2),n in van Genuchten 1980 model, pore size distribution parameter" ;
 	double fates_leaf_c3psn(fates_pft) ;
 		fates_leaf_c3psn:units = "flag" ;
 		fates_leaf_c3psn:long_name = "Photosynthetic pathway (1=c3, 0=c4)" ;
@@ -352,6 +394,9 @@ variables:
 	double fates_nfix2(fates_pft) ;
 		fates_nfix2:units = "NA" ;
 		fates_nfix2:long_name = "place-holder for future n-fixation parameter (NOT IMPLEMENTED)" ;
+	double fates_nitr_store_ratio(fates_pft) ;
+		fates_nitr_store_ratio:units = "(gN/gN)" ;
+		fates_nitr_store_ratio:long_name = "ratio of storeable N, to functional N bound in cell structures of leaf,root,sap" ;
 	double fates_phen_cold_size_threshold(fates_pft) ;
 		fates_phen_cold_size_threshold:units = "cm" ;
 		fates_phen_cold_size_threshold:long_name = "the dbh size above which will lead to phenology-related stem and leaf drop" ;
@@ -370,6 +415,9 @@ variables:
 	double fates_phenflush_fraction(fates_pft) ;
 		fates_phenflush_fraction:units = "fraction" ;
 		fates_phenflush_fraction:long_name = "Upon bud-burst, the maximum fraction of storage carbon used for flushing leaves" ;
+	double fates_phos_store_ratio(fates_pft) ;
+		fates_phos_store_ratio:units = "(gP/gP)" ;
+		fates_phos_store_ratio:long_name = "ratio of storeable P, to functional P bound in cell structures of leaf,root,sap" ;
 	double fates_prescribed_mortality_canopy(fates_pft) ;
 		fates_prescribed_mortality_canopy:units = "1/yr" ;
 		fates_prescribed_mortality_canopy:long_name = "mortality rate of canopy trees for prescribed physiology mode" ;
@@ -383,10 +431,10 @@ variables:
 		fates_prescribed_npp_understory:units = "kgC / m^2 / yr" ;
 		fates_prescribed_npp_understory:long_name = "NPP per unit crown area of understory trees for prescribed physiology mode" ;
 	double fates_prescribed_nuptake(fates_pft) ;
-	   fates_prescribed_nuptake:units = "fraction" ;
+		fates_prescribed_nuptake:units = "fraction" ;
 		fates_prescribed_nuptake:long_name = "Prescribed N uptake flux. 0=fully coupled simulation >0=prescribed (experimental)" ;
 	double fates_prescribed_puptake(fates_pft) ;
-    	fates_prescribed_puptake:units = "fraction" ;
+		fates_prescribed_puptake:units = "fraction" ;
 		fates_prescribed_puptake:long_name = "Prescribed P uptake flux. 0=fully coupled simulation, >0=prescribed (experimental)" ;
 	double fates_prescribed_recruitment(fates_pft) ;
 		fates_prescribed_recruitment:units = "n/yr" ;
@@ -406,13 +454,6 @@ variables:
 	double fates_prt_phos_stoich_p2(fates_prt_organs, fates_pft) ;
 		fates_prt_phos_stoich_p2:units = "(gP/gC)" ;
 		fates_prt_phos_stoich_p2:long_name = "phosphorous stoichiometry, parameter 2" ;
-	double fates_nitr_store_ratio(fates_pft) ;	  
-		fates_nitr_store_ratio:units = "(gN/gN)" ;
-		fates_nitr_store_ratio:long_name = "ratio of storeable N, to functional N bound in cell structures of leaf,root,sap" ;
-	double fates_phos_store_ratio(fates_pft) ;	  
-		fates_phos_store_ratio:units = "(gP/gP)" ;
-		fates_phos_store_ratio:long_name = "ratio of storeable P, to functional P bound in cell structures of leaf,root,sap" ;
-		
 	double fates_recruit_hgt_min(fates_pft) ;
 		fates_recruit_hgt_min:units = "m" ;
 		fates_recruit_hgt_min:long_name = "the minimum height (ie starting height) of a newly recruited plant" ;
@@ -501,6 +542,9 @@ variables:
 	double fates_z0mr(fates_pft) ;
 		fates_z0mr:units = "unitless" ;
 		fates_z0mr:long_name = "Ratio of momentum roughness length to canopy top height" ;
+	double fates_hlm_pft_map(fates_hlm_pftno, fates_pft) ;
+		fates_hlm_pft_map:units = "area fraction" ;
+		fates_hlm_pft_map:long_name = "In fixed biogeog mode, fraction of HLM area associated with each FATES PFT" ;
 	double fates_fire_FBD(fates_litterclass) ;
 		fates_fire_FBD:units = "kg Biomass/m3" ;
 		fates_fire_FBD:long_name = "fuel bulk density" ;
@@ -552,6 +596,9 @@ variables:
 	double fates_cwd_flig ;
 		fates_cwd_flig:units = "unitless" ;
 		fates_cwd_flig:long_name = "Lignin fraction of coarse woody debris" ;
+	double fates_dev_arbitrary ;
+		fates_dev_arbitrary:units = "unknown" ;
+		fates_dev_arbitrary:long_name = "Unassociated free parameter that developers can use for testing arbitrary new hypotheses" ;
 	double fates_eca_plant_escalar ;
 		fates_eca_plant_escalar:units = "" ;
 		fates_eca_plant_escalar:long_name = "scaling factor for plant fine root biomass to calculate nutrient carrier enzyme abundance (ECA)" ;
@@ -642,6 +689,9 @@ variables:
 	double fates_logging_mechanical_frac ;
 		fates_logging_mechanical_frac:units = "fraction" ;
 		fates_logging_mechanical_frac:long_name = "Fraction of stems killed due infrastructure an other mechanical means" ;
+	double fates_maintresp_model ;
+		fates_maintresp_model:units = "unitless" ;
+		fates_maintresp_model:long_name = "switch for choosing between maintenance respiration models. 1=Ryan (1991) (NOT USED)" ;
 	double fates_mort_disturb_frac ;
 		fates_mort_disturb_frac:units = "fraction" ;
 		fates_mort_disturb_frac:long_name = "fraction of canopy mortality that results in disturbance (i.e. transfer of area from new to old patch)" ;
@@ -678,6 +728,12 @@ variables:
 	double fates_phen_ncolddayslim ;
 		fates_phen_ncolddayslim:units = "days" ;
 		fates_phen_ncolddayslim:long_name = "day threshold exceedance for temperature leaf-drop" ;
+	double fates_photo_temp_acclim_timescale ;
+		fates_photo_temp_acclim_timescale:units = "days" ;
+		fates_photo_temp_acclim_timescale:long_name = "Length of the window for the exponential moving average (ema) of vegetation temperature used in photosynthesis temperature acclimation (NOT USED)" ;
+	double fates_photo_tempsens_model ;
+		fates_photo_tempsens_model:units = "unitless" ;
+		fates_photo_tempsens_model:long_name = "switch for choosing the model that defines the temperature sensitivity of photosynthetic parameters (vcmax, jmax). 1=non-acclimating (NOT USED)" ;
 	double fates_q10_froz ;
 		fates_q10_froz:units = "unitless" ;
 		fates_q10_froz:long_name = "Q10 for frozen-soil respiration rates" ;
@@ -687,6 +743,18 @@ variables:
 	double fates_soil_salinity ;
 		fates_soil_salinity:units = "ppt" ;
 		fates_soil_salinity:long_name = "soil salinity used for model when not coupled to dynamic soil salinity" ;
+	double fates_theta_cj_c3 ;
+		fates_theta_cj_c3:units = "unitless" ;
+		fates_theta_cj_c3:long_name = "Empirical curvature parameter for ac, aj photosynthesis co-limitation in c3 plants" ;
+	double fates_theta_cj_c4 ;
+		fates_theta_cj_c4:units = "unitless" ;
+		fates_theta_cj_c4:long_name = "Empirical curvature parameter for ac, aj photosynthesis co-limitation in c4 plants" ;
+	double fates_vai_top_bin_width ;
+		fates_vai_top_bin_width:units = "m2/m2" ;
+		fates_vai_top_bin_width:long_name = "width in VAI units of uppermost leaf+stem layer scattering element in each canopy layer (NOT USED)" ;
+	double fates_vai_width_increase_factor ;
+		fates_vai_width_increase_factor:units = "unitless" ;
+		fates_vai_width_increase_factor:long_name = "factor by which each leaf+stem scattering element increases in VAI width (1 = uniform spacing) (NOT USED)" ;
 
 // global attributes:
 		:history = "This parameter file is maintained in version control\nSee https://github.com/NGEET/fates/blob/master/parameter_files/fates_params_default.cdl \nFor changes, use git blame \n" ;
@@ -701,6 +769,10 @@ data:
  fates_history_sizeclass_bin_edges = 0, 5, 10, 15, 20, 30, 40, 50, 60, 70, 
     80, 90, 100 ;
 
+ fates_hydr_htftype_node = 1, 1, 1, 1 ;
+
+ fates_prt_organ_id = 1, 2, 3, 6 ;
+
  fates_pftname =
   "broadleaf_evergreen_tropical_tree          ",
   "needleleaf_evergreen_extratrop_tree        ",
@@ -715,14 +787,26 @@ data:
   "cool_c3_grass                              ",
   "c4_grass                                   " ;
 
+ fates_hydr_organname_node =
+  "leaf                                             ",
+  "stem                                             ",
+  "transporting root                                ",
+  "absorbing root                                   " ;
+
+ fates_litterclass_name =
+  "twig                                             ",
+  "small branch                                     ",
+  "large branch                                     ",
+  "trunk                                            ",
+  "dead leaves                                      ",
+  "live grass                                       " ;
+
  fates_prt_organ_name =
   "leaf                                             ",
   "fine root                                        ",
   "sapwood                                          ",
   "structure                                        " ;
 
- fates_prt_organ_id = 1, 2, 3, 6 ;
-
  fates_alloc_storage_cushion = 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 
     1.2, 1.2, 1.2 ;
 
@@ -797,10 +881,25 @@ data:
 
  fates_allom_stmode = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;
 
+ fates_allom_zroot_k = 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10 ;
+
+ fates_allom_zroot_max_dbh = 100, 100, 100, 100, 100, 100, 2, 2, 2, 2, 2, 2 ;
+
+ fates_allom_zroot_max_z = 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
+    100, 100 ;
+
+ fates_allom_zroot_min_dbh = 1, 1, 1, 2.5, 2.5, 2.5, 0.1, 0.1, 0.1, 0.1, 0.1, 
+    0.1 ;
+
+ fates_allom_zroot_min_z = 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
+    100, 100 ;
+
  fates_branch_turnover = 150, 150, 150, 150, 150, 150, 150, 150, 150, 0, 0, 0 ;
 
  fates_c2b = 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 ;
 
+ fates_dev_arbitrary_pft = _, _, _, _, _, _, _, _, _, _, _, _ ;
+
  fates_displar = 0.67, 0.67, 0.67, 0.67, 0.67, 0.67, 0.67, 0.67, 0.67, 0.67, 
     0.67, 0.67 ;
 
@@ -884,6 +983,8 @@ data:
   0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08, 0.08,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
 
+ fates_hydr_k_lwp = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
+
  fates_hydr_kmax_node =
   -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999, -999,
   3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
@@ -943,6 +1044,28 @@ data:
   0.65, 0.65, 0.65, 0.65, 0.65, 0.65, 0.65, 0.65, 0.65, 0.65, 0.65, 0.65,
   0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75, 0.75 ;
 
+ fates_hydr_vg_alpha_node =
+  0.0005, 0.0005, 0.0005, 0.0005, 0.0005, 0.005, 0.005, 0.005, 0.005, 0.005, 
+    0.005, 0.005,
+  0.0005, 0.0005, 0.0005, 0.0005, 0.0005, 0.005, 0.005, 0.005, 0.005, 0.005, 
+    0.005, 0.005,
+  0.0005, 0.0005, 0.0005, 0.0005, 0.0005, 0.005, 0.005, 0.005, 0.005, 0.005, 
+    0.005, 0.005,
+  0.0005, 0.0005, 0.0005, 0.0005, 0.0005, 0.005, 0.005, 0.005, 0.005, 0.005, 
+    0.005, 0.005 ;
+
+ fates_hydr_vg_m_node =
+  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
+  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
+  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
+  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5 ;
+
+ fates_hydr_vg_n_node =
+  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
+  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
+  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
+  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 ;
+
  fates_leaf_c3psn = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0 ;
 
  fates_leaf_clumping_index = 0.85, 0.85, 0.8, 0.85, 0.85, 0.9, 0.85, 0.9, 
@@ -1006,7 +1129,7 @@ data:
  fates_maintresp_reduction_curvature = 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 
     0.01, 0.01, 0.01, 0.01, 0.01, 0.01 ;
 
- fates_maintresp_reduction_intercept = 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1 ;
+ fates_maintresp_reduction_intercept = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;
 
  fates_mort_bmort = 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 
     0.014, 0.014, 0.014, 0.014 ;
@@ -1040,6 +1163,9 @@ data:
 
  fates_nfix2 = _, _, _, _, _, _, _, _, _, _, _, _ ;
 
+ fates_nitr_store_ratio = 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 
+    1.5, 1.5 ;
+
  fates_phen_cold_size_threshold = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
 
  fates_phen_evergreen = 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0 ;
@@ -1052,6 +1178,9 @@ data:
 
  fates_phenflush_fraction = _, _, 0.5, _, 0.5, 0.5, _, 0.5, 0.5, 0.5, 0.5, 0.5 ;
 
+ fates_phos_store_ratio = 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 
+    1.5, 1.5 ;
+
  fates_prescribed_mortality_canopy = 0.0194, 0.0194, 0.0194, 0.0194, 0.0194, 
     0.0194, 0.0194, 0.0194, 0.0194, 0.0194, 0.0194, 0.0194 ;
 
@@ -1064,9 +1193,9 @@ data:
  fates_prescribed_npp_understory = 0.03125, 0.03125, 0.03125, 0.03125, 
     0.03125, 0.03125, 0.03125, 0.03125, 0.03125, 0.03125, 0.03125, 0.03125 ;
 
- fates_prescribed_nuptake = 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 ;
+ fates_prescribed_nuptake = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;
 
- fates_prescribed_puptake = 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 ;
+ fates_prescribed_puptake = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;
 
  fates_prescribed_recruitment = 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 
     0.02, 0.02, 0.02, 0.02, 0.02 ;
@@ -1115,10 +1244,6 @@ data:
   0.00047, 0.00047, 0.00047, 0.00047, 0.00047, 0.00047, 0.00047, 0.00047, 
     0.00047, 0.00047, 0.00047, 0.00047 ;
 
- fates_nitr_store_ratio = 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5;
-
- fates_phos_store_ratio = 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5;
-
  fates_recruit_hgt_min = 1.3, 1.3, 1.3, 1.3, 1.3, 1.3, 0.75, 0.75, 0.75, 
     0.125, 0.125, 0.125 ;
 
@@ -1206,6 +1331,22 @@ data:
  fates_z0mr = 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 
     0.055, 0.055, 0.055 ;
 
+ fates_hlm_pft_map =
+  0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
+  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 ;
+
  fates_fire_FBD = 15.4, 16.8, 19.6, 999, 4, 4 ;
 
  fates_fire_low_moisture_Coeff = 1.12, 1.09, 0.98, 0.8, 1.15, 1.15 ;
@@ -1240,6 +1381,8 @@ data:
 
  fates_cwd_flig = 0.24 ;
 
+ fates_dev_arbitrary = _ ;
+
  fates_eca_plant_escalar = 1.25e-05 ;
 
  fates_fire_active_crown_fire = 0 ;
@@ -1300,6 +1443,8 @@ data:
 
  fates_logging_mechanical_frac = 0.05 ;
 
+ fates_maintresp_model = 1 ;
+
  fates_mort_disturb_frac = 1 ;
 
  fates_mort_understorey_death = 0.55983 ;
@@ -1324,9 +1469,21 @@ data:
 
  fates_phen_ncolddayslim = 5 ;
 
+ fates_photo_temp_acclim_timescale = 30 ;
+
+ fates_photo_tempsens_model = 1 ;
+
  fates_q10_froz = 1.5 ;
 
  fates_q10_mr = 1.5 ;
 
  fates_soil_salinity = 0.4 ;
+
+ fates_theta_cj_c3 = 0.999 ;
+
+ fates_theta_cj_c4 = 0.999 ;
+
+ fates_vai_top_bin_width = 1 ;
+
+ fates_vai_width_increase_factor = 1 ;
 }
diff --git a/parteh/PRTParametersMod.F90 b/parteh/PRTParametersMod.F90
index 6e2c17ac66..dcf20dbd14 100644
--- a/parteh/PRTParametersMod.F90
+++ b/parteh/PRTParametersMod.F90
@@ -137,6 +137,13 @@ module PRTParametersMod
      real(r8), allocatable :: allom_agb3(:)                 ! Parameter 3 for agb allometry
      real(r8), allocatable :: allom_agb4(:)                 ! Parameter 3 for agb allometry
 
+     ! ------------------------ (NOT YET IMPLEMENTED) -------------------------
+     real(r8), allocatable :: allom_zroot_max_dbh(:)        ! dbh at which maximum rooting depth saturates (largest possible) [cm]
+     real(r8), allocatable :: allom_zroot_max_z(:)          ! the maximum rooting depth defined at dbh = fates_allom_zroot_max_dbh [m]
+     real(r8), allocatable :: allom_zroot_min_dbh(:)        ! dbh at which the maximum rooting depth for a recruit is defined [cm]
+     real(r8), allocatable :: allom_zroot_min_z(:)          ! the maximum rooting depth defined at dbh = fates_allom_zroot_min_dbh [m]
+     real(r8), allocatable :: allom_zroot_k(:)              ! scale coefficient of logistic rooting depth model
+     
      
   end type prt_param_type
 
diff --git a/parteh/PRTParamsFATESMod.F90 b/parteh/PRTParamsFATESMod.F90
index 4442c090e8..dce172d47d 100644
--- a/parteh/PRTParamsFATESMod.F90
+++ b/parteh/PRTParamsFATESMod.F90
@@ -331,6 +331,26 @@ subroutine PRTRegisterPFT(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_zroot_max_dbh'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+    
+    name = 'fates_allom_zroot_max_z'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+    
+    name = 'fates_allom_zroot_min_dbh'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+    
+    name = 'fates_allom_zroot_min_z'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+    
+    name = 'fates_allom_zroot_k'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+        
     name = 'fates_turnover_retrans_mode'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
           dimension_names=dim_names, lower_bounds=dim_lower_bound)
@@ -554,7 +574,27 @@ subroutine PRTReceivePFT(fates_params)
     name = 'fates_allom_agb4'
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=prt_params%allom_agb4)
-	 
+
+    name = 'fates_allom_zroot_max_dbh'
+    call fates_params%RetreiveParameterAllocate(name=name, &
+         data=prt_params%allom_zroot_max_dbh)
+    
+    name = 'fates_allom_zroot_max_z'
+    call fates_params%RetreiveParameterAllocate(name=name, &
+         data=prt_params%allom_zroot_max_z)
+    
+    name = 'fates_allom_zroot_min_dbh'
+    call fates_params%RetreiveParameterAllocate(name=name, &
+         data=prt_params%allom_zroot_min_dbh)
+    
+    name = 'fates_allom_zroot_min_z'
+    call fates_params%RetreiveParameterAllocate(name=name, &
+         data=prt_params%allom_zroot_min_z)
+    
+    name = 'fates_allom_zroot_k'
+    call fates_params%RetreiveParameterAllocate(name=name, &
+         data=prt_params%allom_zroot_k)
+    
     name = 'fates_branch_turnover'
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=prt_params%branch_long)
@@ -850,6 +890,13 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmt0) 'allom_agb2 = ',prt_params%allom_agb2
         write(fates_log(),fmt0) 'allom_agb3 = ',prt_params%allom_agb3
         write(fates_log(),fmt0) 'allom_agb4 = ',prt_params%allom_agb4
+
+        write(fates_log(),fmt0) 'allom_zroot_max_dbh = ',prt_params%allom_zroot_max_dbh
+        write(fates_log(),fmt0) 'allom_zroot_max_z = ',prt_params%allom_zroot_max_z
+        write(fates_log(),fmt0) 'allom_zroot_min_dbh = ',prt_params%allom_zroot_min_dbh
+        write(fates_log(),fmt0) 'allom_zroot_min_z = ',prt_params%allom_zroot_min_z
+        write(fates_log(),fmt0) 'allom_zroot_k = ',prt_params%allom_zroot_k
+        
         write(fates_log(),fmt0) 'prt_nitr_stoich_p1 = ',prt_params%nitr_stoich_p1
         write(fates_log(),fmt0) 'prt_nitr_stoich_p2 = ',prt_params%nitr_stoich_p2
         write(fates_log(),fmt0) 'prt_phos_stoich_p1 = ',prt_params%phos_stoich_p1
diff --git a/tools/ncvarsort.py b/tools/ncvarsort.py
index bd6587378d..75d80c3799 100755
--- a/tools/ncvarsort.py
+++ b/tools/ncvarsort.py
@@ -29,7 +29,7 @@ def main():
     # make empty lists to hold the variable names in. the first of these is a list of sub-lists,
     # one for each type of variable (based on dimensionality).
     # the second is the master list that will contain all variables.
-    varnames_list = [[],[],[],[],[],[],[],[],[],[]]
+    varnames_list = [[],[],[],[],[],[],[],[],[],[],[],[],[]]
     varnames_list_sorted = []
     #
     # sort the variables by dimensionality, but mix the PFT x other dimension in with the regular PFT-indexed variables
@@ -38,15 +38,20 @@ def main():
     (u'fates_history_coage_bins',):1,
     (u'fates_history_height_bins',):2,
     (u'fates_history_size_bins',):3,
-    (u'fates_pft', u'fates_string_length'):4,
-    (u'fates_prt_organs', u'fates_string_length'):5,
-    (u'fates_pft',):6,
-    (u'fates_hydr_organs', u'fates_pft'):6,
-    (u'fates_leafage_class', u'fates_pft'):6,
-    (u'fates_prt_organs', u'fates_pft'):6,
-    (u'fates_litterclass',):7,
-    (u'fates_NCWD',):8,
-    ():9}
+    (u'fates_hydr_organs',):4,
+    (u'fates_prt_organs',):4,
+    (u'fates_pft', u'fates_string_length'):5,
+    (u'fates_hydr_organs', u'fates_string_length'):6,
+    (u'fates_prt_organs', u'fates_string_length'):7,
+    (u'fates_litterclass', u'fates_string_length'):7,
+    (u'fates_pft',):8,
+    (u'fates_hydr_organs', u'fates_pft'):8,
+    (u'fates_leafage_class', u'fates_pft'):8,
+    (u'fates_prt_organs', u'fates_pft'):8,
+    (u'fates_hlm_pftno', u'fates_pft'):9,
+    (u'fates_litterclass',):10,
+    (u'fates_NCWD',):11,
+    ():12}
     #
     # go through each of the variables and assign it to one of the sub-lists based on its dimensionality
     for v_name, varin in dsin.variables.items():