diff --git a/Externals_CLM.cfg b/Externals_CLM.cfg
index 01943849fa..631695f4bd 100644
--- a/Externals_CLM.cfg
+++ b/Externals_CLM.cfg
@@ -2,7 +2,7 @@
 local_path = src/fates
 protocol = git
 repo_url = https://github.com/NGEET/fates
-tag = sci.1.56.0_api.23.0.0
+tag = sci.1.57.4_api.24.0.0
 required = True
 
 [externals_description]
diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm
index 3676cc57aa..cbd5e3ce44 100755
--- a/bld/CLMBuildNamelist.pm
+++ b/bld/CLMBuildNamelist.pm
@@ -741,7 +741,7 @@ sub setup_cmdl_fates_mode {
 
     $var = "use_fates";
     if ( &value_is_true($nl_flags->{$var}) ) {
-      # This section is a place-holder to test for modules that are not allowed with ED
+      # This section is a place-holder to test for modules that are not allowed with FATES
       # the defaults which are set in the logic section of the namelist builder will
       # automatically set these correctly (well that is the assumption), but here we
       # want to set a catch to fail and warn users if they explicitly set incompatible user namelist
@@ -774,7 +774,7 @@ sub setup_cmdl_fates_mode {
        # dis-allow fates specific namelist items with non-fates runs
        my @list  = (  "fates_spitfire_mode", "use_fates_planthydro", "use_fates_ed_st3", "use_fates_ed_prescribed_phys",
                       "use_fates_cohort_age_tracking",
-                      "use_fates_inventory_init","use_fates_fixed_biogeog","use_fates_nocomp","use_fates_sp","fates_inventory_ctrl_filename","use_fates_logging","fates_parteh_mode" );
+                      "use_fates_inventory_init","use_fates_fixed_biogeog","use_fates_nocomp","use_fates_sp","fates_inventory_ctrl_filename","use_fates_logging","fates_parteh_mode","use_fates_tree_damage" );
        # dis-allow fates specific namelist items with non-fates runs
        foreach my $var ( @list ) {
           if ( defined($nl->get_value($var)) ) {
@@ -1167,7 +1167,7 @@ sub setup_cmdl_spinup {
      $log->fatal_error("$var has an invalid value ($val). Valid values are: @valid_values");
   }
   if ( $nl_flags->{'bgc_spinup'} eq "on" && (not &value_is_true( $nl_flags->{'use_cn'} ))  && (not &value_is_true($nl_flags->{'use_fates'})) ) {
-     $log->fatal_error("$var can not be '$nl_flags->{'bgc_spinup'}' if neither CN nor ED is turned on (use_cn=$nl_flags->{'use_cn'}, use_fates=$nl_flags->{'use_fates'}).");
+     $log->fatal_error("$var can not be '$nl_flags->{'bgc_spinup'}' if neither CN nor FATES is turned on (use_cn=$nl_flags->{'use_cn'}, use_fates=$nl_flags->{'use_fates'}).");
   }
   if ( $nl->get_value("spinup_state") eq 0 && $nl_flags->{'bgc_spinup'} eq "on" ) {
      $log->fatal_error("Namelist spinup_state contradicts the command line option bgc_spinup" );
@@ -2046,6 +2046,10 @@ sub setup_logic_subgrid {
    add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'toosmall_lake');
    add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'toosmall_wetland');
    add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'toosmall_urban');
+
+   if ( &value_is_true($nl_flags->{'use_fates'}) && $nl->get_value('n_dom_pfts') != 0 ) {
+      $log->fatal_error( "FATES and n_dom_pfts can NOT be set at the same time" );
+   }
 }
 
 #-------------------------------------------------------------------------------
@@ -2611,7 +2615,7 @@ sub setup_logic_do_transient_crops {
       # In principle, use_fates should be compatible with
       # do_transient_crops. However, this hasn't been tested, so to be safe,
       # we are not allowing this combination for now.
-      $cannot_be_true = "$var has not been tested with ED, so for now these two options cannot be combined";
+      $cannot_be_true = "$var has not been tested with FATES, so for now these two options cannot be combined";
    }
 
    if ($cannot_be_true) {
@@ -3660,6 +3664,10 @@ sub setup_logic_dry_deposition {
   my ($opts, $nl_flags, $definition, $defaults, $nl) = @_;
 
   if ($opts->{'drydep'} ) {
+    if ( &value_is_true( $nl_flags->{'use_fates'}) && not &value_is_true( $nl_flags->{'use_fates_sp'}) ) {
+       $log->warning("DryDeposition can NOT be on when FATES is also on unless FATES-SP mode is on.\n" .
+                     "   Use the '--no-drydep' option when '-bgc fates' is activated");
+    }
     add_default($opts,  $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'drydep_list');
     add_default($opts,  $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'drydep_method');
   } else {
@@ -3705,7 +3713,7 @@ sub setup_logic_megan {
 
   if ($nl_flags->{'megan'} ) {
     if ( &value_is_true( $nl_flags->{'use_fates'} ) ) {
-       $log->fatal_error("MEGAN can NOT be on when ED is also on.\n" .
+       $log->warning("MEGAN can NOT be on when FATES is also on.\n" .
                    "   Use the '-no-megan' option when '-bgc fates' is activated");
     }
     add_default($opts,  $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'megan_specifier');
@@ -4107,7 +4115,7 @@ sub setup_logic_fates {
         add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'fates_paramfile', 'phys'=>$nl_flags->{'phys'});
         my @list  = (  "fates_spitfire_mode", "use_fates_planthydro", "use_fates_ed_st3", "use_fates_ed_prescribed_phys",
                        "use_fates_inventory_init","use_fates_fixed_biogeog","use_fates_nocomp",
-                       "use_fates_logging","fates_parteh_mode", "use_fates_cohort_age_tracking" );
+                       "use_fates_logging","fates_parteh_mode", "use_fates_cohort_age_tracking","use_fates_tree_damage" );
         foreach my $var ( @list ) {
  	  add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var, 'use_fates'=>$nl_flags->{'use_fates'},
                       'use_fates_sp'=>$nl_flags->{'use_fates_sp'} );
diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml
index 05898ccef4..2686d62b9a 100644
--- a/bld/namelist_files/namelist_defaults_ctsm.xml
+++ b/bld/namelist_files/namelist_defaults_ctsm.xml
@@ -491,7 +491,7 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
 <!-- FATES default parameter file                                       -->
 <!-- ================================================================== -->
 
-<fates_paramfile>lnd/clm2/paramdata/fates_params_api.22.1.0_12pft_c220307.nc</fates_paramfile>
+<fates_paramfile>lnd/clm2/paramdata/fates_params_api.24.0.0_12pft_c220608.nc</fates_paramfile>
 
 <!-- ========================================================================================  -->
 <!-- clm 5.0 BGC nitrogen model                                                                -->
@@ -2541,6 +2541,7 @@ lnd/clm2/surfdata_map/release-clm5.0.30/surfdata_ne0np4.CONUS.ne30x8_hist_78pfts
 <!-- ===== FATES DEFAULTS =========== -->
 <fates_spitfire_mode           use_fates=".true.">0</fates_spitfire_mode>
 <use_fates_planthydro          use_fates=".true.">.false.</use_fates_planthydro>
+<use_fates_tree_damage         use_fates=".true.">.false.</use_fates_tree_damage>
 <use_fates_cohort_age_tracking use_fates=".true.">.false.</use_fates_cohort_age_tracking>
 <use_fates_ed_st3              use_fates=".true.">.false.</use_fates_ed_st3>
 <use_fates_ed_prescribed_phys  use_fates=".true.">.false.</use_fates_ed_prescribed_phys>
diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml
index 4ad4b191d1..4691706d78 100644
--- a/bld/namelist_files/namelist_definition_ctsm.xml
+++ b/bld/namelist_files/namelist_definition_ctsm.xml
@@ -648,6 +648,12 @@ Switch deciding which nutrient model to use in FATES.
 (Only relevant if FATES is on)
 </entry>
 
+<entry id="use_fates_tree_damage" type="logical" category="physics"
+        group="clm_inparm" valid_values="" value=".false."> 
+Toggle to turn on the tree damage module in FATES
+(Only relevant if FATES is on)
+</entry>
+
 <entry id="fates_spitfire_mode" type="integer" category="physics"
         group="clm_inparm" valid_values="0,1,2,3,4,5" value=".false.">
 Turn on spitfire module to simulate fire by setting fates_spitfire_mode > 0.
@@ -2124,7 +2130,7 @@ Command line arguement for biogeochemistry mode for CLM4.5
         Nitrification/De-nitrification
         Methane model
         Vertically resolved Carbon
-   fates = FATES/ED ecosystem demography model with below ground BGC:
+   fates = FATES Functionally Assembled Terrestrial Ecosystem Simulator (ecosystem demography model with below ground BGC):
 </entry>
 
 <!-- ========================================================================================  -->
diff --git a/bld/unit_testers/build-namelist_test.pl b/bld/unit_testers/build-namelist_test.pl
index 824826d578..f6c8d75937 100755
--- a/bld/unit_testers/build-namelist_test.pl
+++ b/bld/unit_testers/build-namelist_test.pl
@@ -163,7 +163,7 @@ sub cat_and_create_namelistinfile {
 #
 # Figure out number of tests that will run
 #
-my $ntests = 1844;
+my $ntests = 1846;
 if ( defined($opts{'compare'}) ) {
    $ntests += 1254;
 }
@@ -956,6 +956,11 @@ sub cat_and_create_namelistinfile {
                                      GLC_TWO_WAY_COUPLING=>"FALSE",
                                      phys=>"clm5_0",
                                    },
+     "useFATESWn_dom_pft"        =>{ options=>"-bgc fates -envxml_dir . -no-megan",
+                                     namelst=>"n_dom_pfts = 1",
+                                     GLC_TWO_WAY_COUPLING=>"FALSE",
+                                     phys=>"clm5_0",
+                                   },
      "useFATESWbMH"              =>{ options=>"-bgc fates -envxml_dir . -no-megan",
                                      namelst=>"use_biomass_heat_storage=.true.",
                                      GLC_TWO_WAY_COUPLING=>"FALSE",
@@ -1001,6 +1006,11 @@ sub cat_and_create_namelistinfile {
                                      GLC_TWO_WAY_COUPLING=>"FALSE",
                                      phys=>"clm4_5",
                                    },
+     "useDRYDEPwithFATES"        =>{ options=>"--bgc fates --envxml_dir . --no-megan --drydep",
+                                     namelst=>"",
+                                     GLC_TWO_WAY_COUPLING=>"FALSE",
+                                     phys=>"clm4_5",
+                                   },
      "useFATESSPWONOCOMP"        =>{ options=>"-bgc fates -envxml_dir . -no-megan",
                                      namelst=>"use_fates_sp=T,use_fates_nocomp=F",
                                      GLC_TWO_WAY_COUPLING=>"FALSE",
diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml
index 7e460f4c43..e1e5764170 100644
--- a/cime_config/testdefs/testlist_clm.xml
+++ b/cime_config/testdefs/testlist_clm.xml
@@ -1804,6 +1804,24 @@
       <option name="wallclock">00:20:00</option>
     </options>
   </test>
+  <test name="SMS_D" grid="f45_f45_mg37" compset="I2000Clm50FatesRs" testmods="clm/FatesColdDefDryDep">
+    <machines>
+      <machine name="cheyenne" compiler="gnu" category="aux_clm"/>
+    </machines>
+    <options>
+      <option name="wallclock">00:20:00</option>
+      <option name="comment"  >Run a test with dry deposition on and regular FATES</option>
+    </options>
+  </test>
+  <test name="ERS_D_Ld15" grid="f45_f45_mg37" compset="I2000Clm50FatesRs" testmods="clm/FatesColdDefTreeDamage">
+    <machines>
+      <machine name="cheyenne" compiler="intel" category="fates"/>
+    </machines>
+    <options>
+      <option name="wallclock">00:20:00</option>
+      <option name="comment"  >Ensure functionality of the tree damage option in FATES</option>
+    </options>
+  </test>
   <test name="ERS_D_Ld3_PS" grid="f09_g17" compset="I2000Clm50FatesRs" testmods="clm/FatesColdDef">
     <machines>
       <machine name="cheyenne" compiler="intel" category="aux_clm"/>
@@ -2060,7 +2078,7 @@
       <option name="comment">Short ERP debug FATES test for f19_g17 grid.</option>
     </options>
   </test>
-  <test name="ERP_D_P32x2_Ld3" grid="f19_g17" compset="I2000Clm50FatesCru" testmods="clm/FatesColdDef">
+  <test name="ERP_D_P36x2_Ld3" grid="f19_g17" compset="I2000Clm50FatesCru" testmods="clm/FatesColdDef">
     <machines>
       <machine name="cheyenne" compiler="intel" category="fates"/>
     </machines>
@@ -2251,6 +2269,24 @@
       <option name="comment">Test with DEBUG on for FatesSP mode, NOTE: FatesSp has the largest difference in CTSM code for any FATES mode</option>
     </options>
   </test>
+  <test name="SMS_D" grid="1x1_brazil" compset="I2000Clm51FatesSpCruRsGs" testmods="clm/FatesColdDefDryDepReducedComplexSatPhen">
+    <machines>
+      <machine name="cheyenne" compiler="gnu" category="aux_clm"/>
+    </machines>
+    <options>
+      <option name="wallclock">00:20:00</option>
+      <option name="comment"  >Run a test with dry deposition on and FATES-SP</option>
+    </options>
+  </test>
+  <test name="SMS_D" grid="1x1_brazil" compset="I2000Clm51FatesSpCruRsGs" testmods="clm/FatesColdDefMeganReducedComplexSatPhen">
+    <machines>
+      <machine name="cheyenne" compiler="gnu" category="aux_clm"/>
+    </machines>
+    <options>
+      <option name="wallclock">00:20:00</option>
+      <option name="comment"  >Run a test with MEGAN on and FATES-SP</option>
+    </options>
+  </test>
 <!--    Don't active this test until #1722 is fixed
   <test name="SMS_Ld3" grid="f09_f09_mg17" compset="I2000Clm51FatesCruRsGs" testmods="clm/FatesColdDefReducedComplexSatPhen_prescribed">
     <machines>
diff --git a/cime_config/testdefs/testmods_dirs/clm/FatesColdDefDryDep/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefDryDep/include_user_mods
new file mode 100644
index 0000000000..1f2c03aebe
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefDryDep/include_user_mods
@@ -0,0 +1 @@
+../FatesColdDef
diff --git a/cime_config/testdefs/testmods_dirs/clm/FatesColdDefDryDep/shell_commands b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefDryDep/shell_commands
new file mode 100644
index 0000000000..870cd99bf5
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefDryDep/shell_commands
@@ -0,0 +1,3 @@
+./xmlchange CLM_BLDNML_OPTS='--drydep' --append
+# The following would be required if you want to run with DryDep and FATES without Fates-SP mode
+#./xmlchange CLM_BLDNML_OPTS="--ignore_warnings" --append
diff --git a/cime_config/testdefs/testmods_dirs/clm/FatesColdDefDryDepReducedComplexSatPhen/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefDryDepReducedComplexSatPhen/include_user_mods
new file mode 100644
index 0000000000..f754dacdab
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefDryDepReducedComplexSatPhen/include_user_mods
@@ -0,0 +1,2 @@
+../FatesColdDefDryDep
+../FatesColdDefReducedComplexSatPhen
diff --git a/cime_config/testdefs/testmods_dirs/clm/FatesColdDefMegan/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefMegan/include_user_mods
new file mode 100644
index 0000000000..1f2c03aebe
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefMegan/include_user_mods
@@ -0,0 +1 @@
+../FatesColdDef
diff --git a/cime_config/testdefs/testmods_dirs/clm/FatesColdDefMegan/shell_commands b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefMegan/shell_commands
new file mode 100644
index 0000000000..749af3486f
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefMegan/shell_commands
@@ -0,0 +1,2 @@
+./xmlchange CLM_BLDNML_OPTS='--megan' --append
+./xmlchange CLM_BLDNML_OPTS="--ignore_warnings" --append
diff --git a/cime_config/testdefs/testmods_dirs/clm/FatesColdDefMeganReducedComplexSatPhen/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefMeganReducedComplexSatPhen/include_user_mods
new file mode 100644
index 0000000000..5dd78011ec
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefMeganReducedComplexSatPhen/include_user_mods
@@ -0,0 +1,2 @@
+../FatesColdDefMegan
+../FatesColdDefReducedComplexSatPhen
diff --git a/cime_config/testdefs/testmods_dirs/clm/FatesColdDefTreeDamage/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefTreeDamage/include_user_mods
new file mode 100644
index 0000000000..e269feea6f
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefTreeDamage/include_user_mods
@@ -0,0 +1 @@
+../FatesColdDef
\ No newline at end of file
diff --git a/cime_config/testdefs/testmods_dirs/clm/FatesColdDefTreeDamage/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefTreeDamage/user_nl_clm
new file mode 100644
index 0000000000..426df087ac
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/clm/FatesColdDefTreeDamage/user_nl_clm
@@ -0,0 +1 @@
+use_fates_tree_damage=.true.
\ No newline at end of file
diff --git a/doc/ChangeLog b/doc/ChangeLog
index 15c36ca5f5..cf9f18c0d4 100644
--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -1,4 +1,91 @@
 ===============================================================
+Tag name: ctsm5.1.dev099
+Originator(s): rgknox (Ryan Knox)
+Date: Tue Jun 21 03:20:42 PM EDT 2022
+One-line Summary: Enabling FATES control over the number of patches on the natural land unit
+
+Purpose and description of changes
+----------------------------------
+
+This set of changes 1) allows FATES to dictate the number of patches on the naturaly vegetated land unit and 2) is compatible with a version of FATES that updates the parameter file format (changes in parameter naming and parameter set).  Efforts were made to adhere to the system of how patches are looped and loop bounds are defined. There is interaction with the reading of the surface file, which previously dictated all patch counts based on the number of PFTs in the file.  FATES-SP and biogeography preserves the ability to use surface file weighting factors associated with PFTs and CFTs to drive area fractions in FATES.  This change is motivated by FATES-side updates to the parameter file, which includes new settings where the user can control the number of primary "fates_maxpatch_primary" and secondary "fates_maxpatch_secondary" patches. Users should also be aware of a new namelist setting "use_fates_tree_damage" which defaults to false (inoperative).
+
+
+Significant changes to scientifically-supported configurations
+--------------------------------------------------------------
+
+Does this tag change answers significantly for any of the following physics configurations?
+
+None
+
+
+[ ] clm5_1
+
+[ ] clm5_0
+
+[ ] ctsm5_0-nwp
+
+[ ] clm4_5
+
+
+Bugs fixed or introduced
+------------------------
+
+CTSM Issues Fixed:
+Fixes #1764 "fates test fails to run on dev095"
+Fixes #1762 "Change threaded fates test from P32x2 to P36x2"
+Fixes #1784 "Get dry-deposition working with FATES-SP mode"
+
+Known bugs introduced in this tag (include issue #): None
+
+Known bugs found since the previous tag (include issue #): None
+
+
+Notes of particular relevance for users
+---------------------------------------
+
+Relevant issue #1785 "Refactor surfrd for FATES"
+Relevant issue #1780 "Some collapse options should be disabled with FATES"
+
+In the process of this PR, we learned that FATES is fundamentally incompatible with the dry-deposition code. FATES-SP may be compatibile, but still needs investigation.
+
+Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): There is a new namelist setting for FATES simulations that turns on the crown damage module: "use_fates_tree_damage". Settings for this module can be found in the FATES parameter file.  This namelist variable is a binary true/false, and it defaults to false.
+
+Changes to the datasets (e.g., parameter, surface or initial files): Users should also be aware that FATES has a new parameter file format and API. The default file on the ucar dataserver has been updated.
+
+Substantial timing or memory changes: none that we're aware of, but with more refined control of patch counts and allocation should come improvements in memory usage.
+
+
+Notes of particular relevance for developers:
+---------------------------------------------
+
+Caveats for developers (e.g., code that is duplicated that requires double maintenance): None.  Some new code was added to surfrdMod.F90, but this is not redundant. Efforts may be made in the future to have more unification there when we have FATES compatable with crops on the crop land-unit.
+
+Changes to tests or testing:  Tests were added to cover the new crown damage module, and combinations of dry-deposition and megan with FATES-SP.  The build system is also aware that these features will not and should not be used with FATES in a non-SP framework.
+
+
+Testing summary:
+----------------
+
+  regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing):
+
+    cheyenne ---- OK
+    izumi ------- OK
+
+  fates tests :
+    cheyenne ---- OK
+    izumi ------- OK
+
+If the tag used for baseline comparisons was NOT the previous tag, note that here: dev098 was updated to point to the moost recent FATES tag: sci.1.57.4_api.23.0.0.
+
+
+Answer changes
+--------------
+
+Changes answers relative to baseline:  NLCOMP changes (differences in namelist fields) were the only differences and they were expected.  All other tests generated B4B results.
+
+
+===============================================================
+===============================================================
 Tag name: ctsm5.1.dev098
 Originator(s): swensosc (Sean Swenson)
 Date: Fri May 20 18:15:53 MDT 2022
diff --git a/doc/ChangeSum b/doc/ChangeSum
index cf4df1cb81..9d5386422b 100644
--- a/doc/ChangeSum
+++ b/doc/ChangeSum
@@ -1,5 +1,6 @@
 Tag                   Who      Date  Summary
 ============================================================================================================================
+    ctsm5.1.dev099   rgknox 06/21/2022 Enabling FATES control over the number of patches on the natural land unit
     ctsm5.1.dev098 swensosc 05/20/2022 Correct perched water table calculation
     ctsm5.1.dev097    sacks 05/20/2022 Update ccs_config to fix issue on izumi and maybe elsewhere
     ctsm5.1.dev096     erik 05/18/2022 Fix a few glitches from the last tag
diff --git a/src/biogeochem/DryDepVelocity.F90 b/src/biogeochem/DryDepVelocity.F90
index 7646699ffa..d8d39ea7db 100644
--- a/src/biogeochem/DryDepVelocity.F90
+++ b/src/biogeochem/DryDepVelocity.F90
@@ -103,9 +103,14 @@ Module DryDepVelocity
   !------------------------------------------------------------------------
   subroutine Init(this, bounds)
 
+    use clm_varctl     , only : use_fates, use_fates_sp
     class(drydepvel_type) :: this
     type(bounds_type), intent(in) :: bounds  
 
+    if ( (.not. use_fates_sp) .and. use_fates .and. (n_drydep > 0 .and. drydep_method == DD_XLND) )then
+       call endrun( msg='ERROR: Dry-deposition currently does NOT work with FATES outside of FATES-SP mode (see github issue #1044)'//&
+                    errMsg(sourcefile, __LINE__))
+    end if
     call this%InitAllocate(bounds)
     call this%InitHistory(bounds)
 
diff --git a/src/biogeochem/EDBGCDynMod.F90 b/src/biogeochem/EDBGCDynMod.F90
index eaf6a750cd..eb13932d13 100644
--- a/src/biogeochem/EDBGCDynMod.F90
+++ b/src/biogeochem/EDBGCDynMod.F90
@@ -260,6 +260,25 @@ subroutine EDBGCDyn(bounds,                                                    &
        call t_stopf('SoilBiogeochemLittVertTransp')
     end if
 
+    ! Wood product fluxes will eventually be added to FATES-CLM. However
+    ! it is likely this will be implemented during or after we break away from
+    ! using this module. This module and the current coupling stategy bypasses
+    ! a number of processes in CLM, which includes the wood product modules.
+    ! Therefore the following call is a placeholder so that the wood-product 
+    ! wrapper code can be copied from here and applied at the right place when the time comes.
+    ! RGK 06-2022
+    
+    !call FatesWrapWoodProducts(bounds, num_soilc, filter_soilc,c_products_inst)
+    !call t_startf('CNWoodProducts')
+    !call c_products_inst%UpdateProducts(bounds, &
+    !     num_soilp, filter_soilp, &
+    !     dwt_wood_product_gain_patch = cnveg_carbonflux_inst%dwt_wood_productc_gain_patch(begp:endp), &
+    !     wood_harvest_patch = cnveg_carbonflux_inst%wood_harvestc_patch(begp:endp), &
+    !     dwt_crop_product_gain_patch = cnveg_carbonflux_inst%dwt_crop_productc_gain_patch(begp:endp), &
+    !     crop_harvest_to_cropprod_patch = cnveg_carbonflux_inst%crop_harvestc_to_cropprodc_patch(begp:endp))
+    !call t_stopf('CNWoodProducts')
+
+    
     end associate
 
   end subroutine EDBGCDyn
diff --git a/src/biogeochem/VOCEmissionMod.F90 b/src/biogeochem/VOCEmissionMod.F90
index 88ae7c08cc..f1865af3b7 100644
--- a/src/biogeochem/VOCEmissionMod.F90
+++ b/src/biogeochem/VOCEmissionMod.F90
@@ -84,10 +84,15 @@ module VOCEmissionMod
   !------------------------------------------------------------------------
   subroutine Init(this, bounds)
 
+    use clm_varctl     , only : use_fates, use_fates_sp
     class(vocemis_type) :: this
     type(bounds_type), intent(in)    :: bounds  
 
     if ( shr_megan_mechcomps_n > 0) then
+       if ( use_fates .and. (.not. use_fates_sp) ) then
+           call endrun( msg='ERROR: MEGAN currently does NOT work with FATES outside of FATES-SP mode (see github issue #115)'//&
+                     errMsg(sourcefile, __LINE__))
+       end if
        call this%InitAllocate(bounds) 
        call this%InitHistory(bounds)
        call this%InitCold(bounds)
diff --git a/src/biogeophys/test/Photosynthesis_test/test_Photosynthesis.pf b/src/biogeophys/test/Photosynthesis_test/test_Photosynthesis.pf
index 78c82e1bf4..74ca9c9c4c 100644
--- a/src/biogeophys/test/Photosynthesis_test/test_Photosynthesis.pf
+++ b/src/biogeophys/test/Photosynthesis_test/test_Photosynthesis.pf
@@ -38,7 +38,7 @@ contains
     soil_layerstruct_predefined = '20SL_8.5m'
     
     call setup_ncells_single_veg_patch(ncells=1, pft_type=1)
-    call clm_varpar_init( actual_maxsoil_patches=17, actual_numcft=2 )
+    call clm_varpar_init( actual_maxsoil_patches=17, surf_numpft=15, surf_numcft=2 )
     call this%photo%Init( bounds )
     call this%photo%setParamsForTesting( )
 
diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90
index e02f4b853a..b9a7bc7f36 100644
--- a/src/main/clm_driver.F90
+++ b/src/main/clm_driver.F90
@@ -221,24 +221,41 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
 
     ! ============================================================================
     ! Specified phenology
+    ! Done in SP mode, FATES-SP mode and also when dry-deposition is active
     ! ============================================================================
 
     if (use_cn) then
        ! For dry-deposition need to call CLMSP so that mlaidiff is obtained
+       ! NOTE: This is also true of FATES below
        if ( n_drydep > 0 .and. drydep_method == DD_XLND ) then
           call t_startf('interpMonthlyVeg')
           call interpMonthlyVeg(bounds_proc, canopystate_inst)
           call t_stopf('interpMonthlyVeg')
        endif
 
+    elseif(use_fates) then
+
+       ! For FATES-Specified phenology mode interpolate the weights for
+       ! time-interpolation of monthly vegetation data (as in SP mode below)
+       ! Also for FATES with dry-deposition as above need to call CLMSP so that mlaidiff is obtained
+       !if ( use_fates_sp .or. (n_drydep > 0 .and. drydep_method == DD_XLND ) ) then    ! Replace with this when we have dry-deposition working
+       ! For now don't allow for dry-deposition because of issues in #1044 EBK Jun/17/2022
+       if ( use_fates_sp ) then
+          call t_startf('interpMonthlyVeg')
+          call interpMonthlyVeg(bounds_proc, canopystate_inst)
+          call t_stopf('interpMonthlyVeg')
+       end if
+
     else
+
        ! Determine weights for time interpolation of monthly vegetation data.
        ! This also determines whether it is time to read new monthly vegetation and
        ! obtain updated leaf area index [mlai1,mlai2], stem area index [msai1,msai2],
        ! vegetation top [mhvt1,mhvt2] and vegetation bottom [mhvb1,mhvb2]. The
        ! weights obtained here are used in subroutine SatellitePhenology to obtain time
        ! interpolated values.
-       if (doalb .or. ( n_drydep > 0 .and. drydep_method == DD_XLND ) .or. use_fates_sp) then
+       ! This is also done for FATES-SP mode above
+       if ( doalb .or. ( n_drydep > 0 .and. drydep_method == DD_XLND ) )then
           call t_startf('interpMonthlyVeg')
           call interpMonthlyVeg(bounds_proc, canopystate_inst)
           call t_stopf('interpMonthlyVeg')
diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90
index e5f6673b7c..44ac50b27c 100644
--- a/src/main/clm_initializeMod.F90
+++ b/src/main/clm_initializeMod.F90
@@ -26,7 +26,7 @@ module clm_initializeMod
   use PatchType             , only : patch         ! instance
   use reweightMod           , only : reweight_wrapup
   use filterMod             , only : allocFilters, filter, filter_inactive_and_active
-  use CLMFatesInterfaceMod  , only : CLMFatesGlobals
+  use CLMFatesInterfaceMod  , only : CLMFatesGlobals1,CLMFatesGlobals2
   use CLMFatesInterfaceMod  , only : CLMFatesTimesteps
   use dynSubgridControlMod  , only : dynSubgridControl_init, get_reset_dynbal_baselines
   use SelfTestDriver        , only : self_test_driver
@@ -40,6 +40,7 @@ module clm_initializeMod
   public :: initialize2  ! Phase two initialization
 
   integer :: actual_numcft  ! numcft from sfc dataset
+  integer :: actual_numpft  ! numpft from sfc dataset
 
 !-----------------------------------------------------------------------
 contains
@@ -96,8 +97,16 @@ subroutine initialize1(dtime)
 
     call control_init(dtime)
     call ncd_pio_init()
-    call surfrd_get_num_patches(fsurdat, actual_maxsoil_patches, actual_numcft)
-    call clm_varpar_init(actual_maxsoil_patches, actual_numcft)
+    call surfrd_get_num_patches(fsurdat, actual_maxsoil_patches, actual_numpft, actual_numcft)
+
+    ! If fates is on, we override actual_maxsoil_patches. FATES dictates the
+    ! number of patches per column.  We still use numcft from the surface
+    ! file though...
+    if(use_fates) then
+       call CLMFatesGlobals1(actual_numpft, actual_numcft, actual_maxsoil_patches)
+    end if
+
+    call clm_varpar_init(actual_maxsoil_patches, actual_numpft, actual_numcft)
     call decomp_cascade_par_init( NLFilename )
     call clm_varcon_init( IsSimpleBuildTemp() )
     call landunit_varcon_init()
@@ -118,7 +127,9 @@ subroutine initialize2(ni,nj)
     ! !USES:
     use clm_varcon                    , only : spval
     use clm_varpar                    , only : natpft_lb, natpft_ub, cft_lb, cft_ub, maxpatch_glc
+    use clm_varpar                    , only : surfpft_lb, surfpft_ub
     use clm_varpar                    , only : nlevsno
+    use clm_varpar                    , only : natpft_size,cft_size
     use clm_varctl                    , only : fsurdat
     use clm_varctl                    , only : finidat, finidat_interp_source, finidat_interp_dest, fsurdat
     use clm_varctl                    , only : use_cn, use_fates
@@ -210,7 +221,6 @@ subroutine initialize2(ni,nj)
     ! Allocate surface grid dynamic memory (just gridcell bounds dependent)
     allocate (wt_lunit     (begg:endg, max_lunit           ))
     allocate (urban_valid  (begg:endg                      ))
-    allocate (wt_nat_patch (begg:endg, natpft_lb:natpft_ub ))
     allocate (wt_cft       (begg:endg, cft_lb:cft_ub       ))
     allocate (fert_cft     (begg:endg, cft_lb:cft_ub       ))
     allocate (irrig_method (begg:endg, cft_lb:cft_ub       ))
@@ -218,7 +228,8 @@ subroutine initialize2(ni,nj)
     allocate (topo_glc_mec (begg:endg, maxpatch_glc     ))
     allocate (haslake      (begg:endg                      ))
     allocate (pct_urban_max(begg:endg, numurbl             ))
-
+    allocate (wt_nat_patch (begg:endg, surfpft_lb:surfpft_ub ))
+    
     ! Read list of Patches and their corresponding parameter values
     ! Independent of model resolution, Needs to stay before surfrd_get_data
     call pftcon%Init()
@@ -226,20 +237,24 @@ subroutine initialize2(ni,nj)
     ! Read surface dataset and set up subgrid weight arrays
     call surfrd_get_data(begg, endg, ldomain, fsurdat, actual_numcft)
 
-    ! Ask Fates to evaluate its own dimensioning needs.
-    ! This determines the total amount of space it requires in its largest
-    ! dimension.  We are currently calling that the "cohort" dimension, but
-    ! it is really a utility dimension that captures the models largest
-    ! size need.
-    ! Sets:
-    !   fates_maxElementsPerPatch
-    !   fates_maxElementsPerSite (where a site is roughly equivalent to a column)
-    ! (Note: fates_maxELementsPerSite is the critical variable used by CLM
-    ! to allocate space)
-    ! This also sets up various global constants in FATES
-    ! ------------------------------------------------------------------------
-    
-    call CLMFatesGlobals()
+    if(use_fates) then
+
+       ! Ask Fates to evaluate its own dimensioning needs.
+       ! This determines the total amount of space it requires in its largest
+       ! dimension.  We are currently calling that the "cohort" dimension, but
+       ! it is really a utility dimension that captures the models largest
+       ! size need.
+       ! Sets:
+       !   fates_maxElementsPerPatch
+       !   fates_maxElementsPerSite (where a site is roughly equivalent to a column)
+       ! (Note: fates_maxELementsPerSite is the critical variable used by CLM
+       ! to allocate space)
+       ! This also sets up various global constants in FATES
+       ! ------------------------------------------------------------------------
+
+       call CLMFatesGlobals2()
+
+    end if
 
     ! Determine decomposition of subgrid scale landunits, columns, patches
     call decompInit_clumps(ni, nj, glc_behavior)
@@ -417,6 +432,7 @@ subroutine initialize2(ni,nj)
        if (n_drydep > 0 .and. drydep_method == DD_XLND) then
           ! Must do this also when drydeposition is used so that estimates of monthly
           ! differences in LAI can be computed
+          ! Also do this for FATES see below
           call SatellitePhenologyInit(bounds_proc)
        end if
        if ( use_c14 .and. use_c14_bombspike ) then
@@ -425,12 +441,20 @@ subroutine initialize2(ni,nj)
        if ( use_c13 .and. use_c13_timeseries ) then
           call C13_init_TimeSeries()
        end if
-    else
-       call SatellitePhenologyInit(bounds_proc)
+
+    else ! FATES OR Satellite phenology
+
+       ! For SP FATES-SP Initialize SP
+       ! Also for FATES with Dry-Deposition on as well (see above)
+       !if(use_fates_sp .or. (.not.use_cn) .or. (n_drydep > 0 .and.  drydep_method == DD_XLND) )then  !  Replace with this when we have dry-deposition working
+       ! For now don't allow for dry-deposition because of issues in #1044 EBK Jun/17/2022
+       if( use_fates_sp .or. .not. use_fates )then
+          call SatellitePhenologyInit(bounds_proc)
+       end if
 
        ! fates_spitfire_mode is assigned an integer value in the namelist
        ! see bld/namelist_files/namelist_definition_clm4_5.xml for details
-       if (fates_spitfire_mode > scalar_lightning) then
+       if(use_fates .and. (fates_spitfire_mode > scalar_lightning)) then
           call clm_fates%Init2(bounds_proc, NLFilename)
        end if
     end if
@@ -589,19 +613,17 @@ subroutine initialize2(ni,nj)
     end if
 
     ! Read monthly vegetation
-    ! Even if CN is on, and dry-deposition is active, read CLMSP annual vegetation
+    ! Even if CN or FATES is on, and dry-deposition is active, read CLMSP annual vegetation
     ! to get estimates of monthly LAI
     if ( n_drydep > 0 .and. drydep_method == DD_XLND )then
        call readAnnualVegetation(bounds_proc, canopystate_inst)
-       if (nsrest == nsrStartup .and. finidat /= ' ') then
-          ! Call interpMonthlyVeg for dry-deposition so that mlaidiff will be calculated
-          ! This needs to be done even if CN or CNDV is on!
-          call interpMonthlyVeg(bounds_proc, canopystate_inst)
-       end if
+       ! Call interpMonthlyVeg for dry-deposition so that mlaidiff will be calculated
+       ! This needs to be done even if FATES, CN or CNDV is on!
+       call interpMonthlyVeg(bounds_proc, canopystate_inst)
     ! If fates has satellite phenology enabled, get the monthly veg values
     ! prior to the first call to SatellitePhenology()
     elseif ( use_fates_sp ) then
-          call interpMonthlyVeg(bounds_proc, canopystate_inst)
+       call interpMonthlyVeg(bounds_proc, canopystate_inst)
     end if
 
     ! Determine gridcell averaged properties to send to atm
diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90
index 40135e8575..c84c21ffe5 100644
--- a/src/main/clm_varctl.F90
+++ b/src/main/clm_varctl.F90
@@ -259,6 +259,7 @@ module clm_varctl
   integer, public            :: fates_spitfire_mode = 0                
   ! 0 for no fire; 1 for constant ignitions; > 1 for external data (lightning and/or anthropogenic ignitions)
   ! see bld/namelist_files/namelist_definition_clm4_5.xml for details
+  logical, public            :: use_fates_tree_damage = .false.        ! true => turn on tree damage module
   logical, public            :: use_fates_logging = .false.            ! true => turn on logging module
   logical, public            :: use_fates_planthydro = .false.         ! true => turn on fates hydro
   logical, public            :: use_fates_cohort_age_tracking = .false. ! true => turn on cohort age tracking
diff --git a/src/main/clm_varpar.F90 b/src/main/clm_varpar.F90
index 7f0d5faaec..23b237d6f3 100644
--- a/src/main/clm_varpar.F90
+++ b/src/main/clm_varpar.F90
@@ -60,7 +60,7 @@ module clm_varpar
   integer, public    :: nlayert               ! number of VIC soil layer + 3 lower thermal layers
   integer, public, parameter :: nvariants   =   2     ! number of variants of PFT constants
 
-  integer, public :: maxveg           ! # of pfts + cfts
+  integer, public :: maxveg                ! # of pfts + cfts
   integer, public :: maxpatch_urb= 5       ! max number of urban patches (columns) in urban landunit
 
   integer, public :: maxsoil_patches  ! # of pfts + cfts + bare ground; replaces maxpatch_pft, which is obsolete
@@ -90,6 +90,12 @@ module clm_varpar
   integer, public :: natpft_ub          ! In PATCH arrays, upper bound of Patches on the natural veg landunit
   integer, public :: natpft_size        ! Number of Patches on natural veg landunit (including bare ground)
 
+  integer, public :: surfpft_lb         ! Lower bound of PFTs in the surface file
+                                        ! synonymous with natpft_lb for non-fates and fates-sp
+  integer, public :: surfpft_ub         ! Upper bound of PFTs in the surface file
+                                        ! synonymous with natpft_ub for non-fates and fates-sp
+
+  
   ! The following variables pertain to arrays of all PFTs - e.g., those dimensioned (g,
   ! pft_index). These include unused CFTs that are merged into other CFTs. Thus, these
   ! variables do NOT give the actual number of CFTs on the crop landunit - that number
@@ -110,15 +116,19 @@ module clm_varpar
 contains
 
   !------------------------------------------------------------------------------
-  subroutine clm_varpar_init(actual_maxsoil_patches, actual_numcft)
+  subroutine clm_varpar_init(actual_maxsoil_patches, surf_numpft, surf_numcft)
     !
     ! !DESCRIPTION:
     ! Initialize module variables 
     !
     ! !ARGUMENTS:
     implicit none
-    integer, intent(in) :: actual_maxsoil_patches  ! value from surface dataset
-    integer, intent(in) :: actual_numcft  ! Actual number of crops
+    integer, intent(in) :: actual_maxsoil_patches  ! Number of soil patches to allocate
+                                                   ! This value comes either from the
+                                                   ! surface dataset (non-fates) or 
+                                                   ! from fates (via its parameter file)
+    integer, intent(in) :: surf_numpft             ! Number of PFTs in the surf dataset
+    integer, intent(in) :: surf_numcft             ! Number of CFTs in the surf dataset
     !
     ! !LOCAL VARIABLES:
     !
@@ -126,37 +136,64 @@ subroutine clm_varpar_init(actual_maxsoil_patches, actual_numcft)
     character(len=32) :: subname = 'clm_varpar_init'  ! subroutine name
     !------------------------------------------------------------------------------
 
-    ! actual_maxsoil_patches and actual_numcft were read directly from the
-    ! surface dataset
-    maxsoil_patches = actual_maxsoil_patches  ! # of patches with bare ground
+    ! actual_maxsoil_patches is either the total number of cfts+pfts in the surface
+    ! file (for non-fates), or the number of patches plus the bareground that
+    ! fates requests.  If this is a fates-sp run, this value will also be the number 
+    ! of cfts+pfts as in a nonfates run
+
+    maxsoil_patches = actual_maxsoil_patches
+    
     maxveg = maxsoil_patches - 1  ! # of patches without bare ground
 
     ! For arrays containing all Patches (natural veg & crop), determine lower and upper bounds
     ! for (1) Patches on the natural vegetation landunit (includes bare ground, and includes
     ! crops if create_crop_landunit=false), and (2) CFTs on the crop landunit (no elements
     ! if create_crop_landunit=false)
+    ! As for when we don't have a crop LU, which is currently when FATES is on...
+    ! These values are used to create the wt_nat_patch array that is used by fates_sp 
+    ! and fixed biogeog. Also, the pft and cft vectors are concatenated into
+    ! the natpft vector (wt_nat_patch), the wt_cft array is unused (size zero)
+    ! The following values should not be used for allocating patch structures
+    ! though.  That should be handled completely by maxoil_patches and maxveg
 
     if (create_crop_landunit) then
-       natpft_size = maxsoil_patches - actual_numcft  ! includes bare ground
-       cft_size    = actual_numcft
-    else
-       natpft_size = maxsoil_patches  ! includes bare ground
+       
+       natpft_size = surf_numpft    ! includes bare ground + natveg pfts
+       cft_size    = surf_numcft
+       natpft_lb   = 0
+       natpft_ub   = natpft_lb + natpft_size - 1
+       cft_lb      = natpft_ub + 1
+       cft_ub      = cft_lb + cft_size - 1
+       surfpft_lb  = natpft_lb
+       surfpft_ub  = natpft_ub
+       
+    else ! only true when FATES is active
+       
+       natpft_size = maxsoil_patches
        cft_size    = 0
+       natpft_lb   = 0
+       natpft_ub   = natpft_lb + natpft_size - 1
+       cft_lb      = 0
+       cft_ub      = 0
+       surfpft_lb  = 0
+       surfpft_ub  = surf_numpft+surf_numcft-1
+       
     end if
-
-    natpft_lb = 0
-    natpft_ub = natpft_lb + natpft_size - 1
-    cft_lb = natpft_ub + 1
-    cft_ub = cft_lb + cft_size - 1
-
+       
+    if(use_fates_sp .and. (natpft_ub .ne. maxveg) ) then
+       write(iulog,*) 'maxveg should match the upper bound for non-fates and fates-sp runs'
+       write(iulog,*) 'the surface dataset PFT+CFT indices (ie lsmft), yours: ',natpft_ub,maxveg
+       call shr_sys_abort(subname//' ERROR: conflict in maxveg and pft bounds')
+    end if
+    
     mxharvests = mxsowings + 1
 
-    ! TODO(wjs, 2015-10-04, bugz 2227) Using actual_numcft in this 'max' gives a significant
+    ! TODO(wjs, 2015-10-04, bugz 2227) Using surf_numcft in this 'max' gives a significant
     ! overestimate of max_patch_per_col when use_crop is true. This should be reworked -
     ! or, better, removed from the code entirely (because it is a maintenance problem, and
     ! I can't imagine that looping idioms that use it help performance that much, and
     ! likely they hurt performance.)
-    max_patch_per_col= max(maxsoil_patches, actual_numcft, maxpatch_urb)
+    max_patch_per_col= max(maxsoil_patches, surf_numcft, maxpatch_urb)
 
     nlevsoifl   =  10
     nlevurb     =  5
diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90
index b0f1dabaea..f5ed6ff986 100644
--- a/src/main/controlMod.F90
+++ b/src/main/controlMod.F90
@@ -229,7 +229,8 @@ subroutine control_init(dtime)
           use_fates_nocomp,                             &
           use_fates_sp,                                 &
           fates_inventory_ctrl_filename,                &
-          fates_parteh_mode
+          fates_parteh_mode,                            &
+          use_fates_tree_damage
     
    ! Ozone vegetation stress method 
    namelist / clm_inparam / o3_veg_stress_method
@@ -706,6 +707,7 @@ subroutine control_spmd()
     call mpi_bcast (fates_spitfire_mode, 1, MPI_INTEGER, 0, mpicom, ier)
     call mpi_bcast (use_fates_logging, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_fates_planthydro, 1, MPI_LOGICAL, 0, mpicom, ier)
+    call mpi_bcast (use_fates_tree_damage, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_fates_cohort_age_tracking, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_fates_ed_st3, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast (use_fates_ed_prescribed_phys,  1, MPI_LOGICAL, 0, mpicom, ier)
@@ -1048,6 +1050,7 @@ subroutine control_print ()
        write(iulog, *) '    fates_paramfile = ', fates_paramfile
        write(iulog, *) '    fates_parteh_mode = ', fates_parteh_mode
        write(iulog, *) '    use_fates_planthydro = ', use_fates_planthydro
+       write(iulog, *) '    use_fates_tree_damage = ', use_fates_tree_damage
        write(iulog, *) '    use_fates_cohort_age_tracking = ', use_fates_cohort_age_tracking
        write(iulog, *) '    use_fates_ed_st3 = ',use_fates_ed_st3
        write(iulog, *) '    use_fates_ed_prescribed_phys = ',use_fates_ed_prescribed_phys
diff --git a/src/main/histFileMod.F90 b/src/main/histFileMod.F90
index 5224f44ae2..a7753ac058 100644
--- a/src/main/histFileMod.F90
+++ b/src/main/histFileMod.F90
@@ -25,6 +25,7 @@ module histFileMod
   use EDTypesMod     , only : nlevleaf
   use FatesInterfaceTypesMod , only : nlevsclass, nlevage, nlevcoage
   use FatesInterfaceTypesMod , only : nlevheight
+  use FatesInterfaceTypesMod , only : nlevdamage
   use EDTypesMod        , only : nfsc
   use FatesLitterMod    , only : ncwd
   use PRTGenericMod     , only : num_elements_fates  => num_elements
@@ -2459,6 +2460,9 @@ subroutine htape_create (t, histrest)
        call ncd_defdim(lnfid, 'fates_levleaf', nlevleaf, dimid)
        call ncd_defdim(lnfid, 'fates_levcnlf', nlevleaf * nclmax, dimid)
        call ncd_defdim(lnfid, 'fates_levcnlfpf', nlevleaf * nclmax * numpft_fates, dimid)
+       call ncd_defdim(lnfid, 'fates_levcdsc', nlevdamage * nlevsclass, dimid)
+       call ncd_defdim(lnfid, 'fates_levcdpf', nlevdamage * nlevsclass * numpft_fates, dimid)
+       call ncd_defdim(lnfid, 'fates_levcdam', nlevdamage, dimid)
        call ncd_defdim(lnfid, 'fates_levelem', num_elements_fates, dimid)
        call ncd_defdim(lnfid, 'fates_levelpft', num_elements_fates * numpft_fates, dimid)
        call ncd_defdim(lnfid, 'fates_levelcwd', num_elements_fates * ncwd, dimid)
@@ -3010,6 +3014,11 @@ subroutine htape_timeconst(t, mode)
     use FatesInterfaceTypesMod, only : fates_hdim_agemap_levelage
     use FatesInterfaceTypesMod, only : fates_hdim_agmap_levagefuel
     use FatesInterfaceTypesMod, only : fates_hdim_fscmap_levagefuel
+    use FatesInterfaceTypesMod, only : fates_hdim_scmap_levcdsc
+    use FatesInterfaceTypesMod, only : fates_hdim_cdmap_levcdsc
+    use FatesInterfaceTypesMod, only : fates_hdim_scmap_levcdpf
+    use FatesInterfaceTypesMod, only : fates_hdim_cdmap_levcdpf
+    use FatesInterfaceTypesMod, only : fates_hdim_pftmap_levcdpf
 
 
     !
@@ -3126,8 +3135,18 @@ subroutine htape_timeconst(t, mode)
                    long_name='FATES age-class map into patch age x fuel size', units='-', ncid=nfid(t))
              call ncd_defvar(varname='fates_fscmap_levagefuel', xtype=ncd_int, dim1name='fates_levagefuel', &
                    long_name='FATES fuel size-class map into patch age x fuel size', units='-', ncid=nfid(t))
-
-
+             call ncd_defvar(varname='fates_cdmap_levcdsc',xtype=ncd_int, dim1name='fates_levcdsc', &
+                  long_name='FATES damage index of the combined damage-size dimension', ncid=nfid(t))
+             call ncd_defvar(varname='fates_scmap_levcdsc',xtype=ncd_int, dim1name='fates_levcdsc', &
+                  long_name='FATES size index of the combined damage-size dimension', ncid=nfid(t))
+             call ncd_defvar(varname='fates_cdmap_levcdpf',xtype=ncd_int, dim1name='fates_levcdpf', &
+                  long_name='FATES damage index of the combined damage-size-PFT dimension', ncid=nfid(t))
+             call ncd_defvar(varname='fates_scmap_levcdpf',xtype=ncd_int, dim1name='fates_levcdpf', &
+                  long_name='FATES size index of the combined damage-size-PFT dimension', ncid=nfid(t))
+             call ncd_defvar(varname='fates_pftmap_levcdpf',xtype=ncd_int, dim1name='fates_levcdpf', &
+                  long_name='FATES pft index of the combined damage-size-PFT dimension', ncid=nfid(t))
+             call ncd_defvar(varname='fates_levcdam', xtype=tape(t)%ncprec, dim1name='fates_levcdam', &
+                  long_name='FATES damage class lower bound', units='unitless', ncid=nfid(t))
 
           end if
 
@@ -3171,6 +3190,11 @@ subroutine htape_timeconst(t, mode)
              call ncd_io(varname='fates_agmap_levagepft',data=fates_hdim_agmap_levagepft, ncid=nfid(t), flag='write')
              call ncd_io(varname='fates_agmap_levagefuel',data=fates_hdim_agmap_levagefuel, ncid=nfid(t), flag='write')
              call ncd_io(varname='fates_fscmap_levagefuel',data=fates_hdim_fscmap_levagefuel, ncid=nfid(t), flag='write')
+             call ncd_io(varname='fates_scmap_levcdsc',data=fates_hdim_scmap_levcdsc, ncid=nfid(t), flag='write')
+             call ncd_io(varname='fates_cdmap_levcdsc',data=fates_hdim_cdmap_levcdsc, ncid=nfid(t), flag='write')
+             call ncd_io(varname='fates_scmap_levcdpf',data=fates_hdim_scmap_levcdpf, ncid=nfid(t), flag='write')
+             call ncd_io(varname='fates_cdmap_levcdpf',data=fates_hdim_cdmap_levcdpf, ncid=nfid(t), flag='write')
+             call ncd_io(varname='fates_pftmap_levcdpf',data=fates_hdim_pftmap_levcdpf, ncid=nfid(t), flag='write')
           end if
 
        endif
@@ -5466,6 +5490,12 @@ subroutine hist_addfld2d (fname, type2d, units, avgflag, long_name, type1d_out,
        num2d = nlevleaf * nclmax
     case ('fates_levcnlfpf')
        num2d = nlevleaf * nclmax * numpft_fates
+    case ('fates_levcdsc')
+       num2d = nlevdamage * nlevsclass
+    case ('fates_levcdpf')
+       num2d = nlevdamage * nlevsclass * numpft_fates
+    case ('fates_levcdam')
+       num2d = nlevdamage
     case ('ltype')
        num2d = max_lunit
     case ('natpft')
diff --git a/src/main/initGridCellsMod.F90 b/src/main/initGridCellsMod.F90
index 7a8a9404c0..a330658c1f 100644
--- a/src/main/initGridCellsMod.F90
+++ b/src/main/initGridCellsMod.F90
@@ -16,7 +16,7 @@ module initGridCellsMod
   use shr_log_mod    , only : errMsg => shr_log_errMsg
   use spmdMod        , only : masterproc,iam
   use abortutils     , only : endrun
-  use clm_varctl     , only : iulog
+  use clm_varctl     , only : iulog, use_fates, use_fates_sp
   use clm_varcon     , only : namep, namec, namel, nameg
   use decompMod      , only : bounds_type
   use GridcellType   , only : grc                
@@ -206,7 +206,7 @@ subroutine set_landunit_veg_compete (ltype, gi, li, ci, pi)
     ! !USES
     use clm_instur, only : wt_lunit, wt_nat_patch
     use subgridMod, only : subgrid_get_info_natveg, natveg_patch_exists
-    use clm_varpar, only : natpft_lb, natpft_ub
+    use clm_varpar, only : natpft_lb, natpft_ub, natpft_size
     !
     ! !ARGUMENTS:
     integer , intent(in)    :: ltype             ! landunit type
@@ -224,6 +224,7 @@ subroutine set_landunit_veg_compete (ltype, gi, li, ci, pi)
     integer  :: ncols_added                      ! number of columns actually added
     integer  :: nlunits_added                    ! number of landunits actually added
     real(r8) :: wtlunit2gcell                    ! landunit weight in gridcell
+    real(r8) :: p_wt                             ! patch weight (0-1)
     !------------------------------------------------------------------------
 
     ! Set decomposition properties
@@ -244,12 +245,26 @@ subroutine set_landunit_veg_compete (ltype, gi, li, ci, pi)
        call add_column(ci=ci, li=li, ctype=1, wtlunit=1.0_r8)
        ncols_added = ncols_added + 1
 
+       ! For FATES: the total number of patches may not match what is in the surface
+       ! file, and therefor the weighting can't be used. The weightings in
+       ! wt_nat_patch may be meaningful (like with fixed biogeography), but they
+       ! they need a mapping table to connect to the allocated patches (in fates)
+       ! so the wt_nat_patch array is not applicable to these area weights
+       ! A subsequent call, via the clmfates interface will update these weights
+       ! by using said mapping table
+       
        do m = natpft_lb,natpft_ub
           if (natveg_patch_exists(gi, m)) then
-             call add_patch(pi=pi, ci=ci, ptype=m, wtcol=wt_nat_patch(gi,m))
+             if(use_fates .and. .not.use_fates_sp)then
+                p_wt = 1.0_r8/real(natpft_size,r8)
+             else
+                p_wt = wt_nat_patch(gi,m)
+             end if
+             call add_patch(pi=pi, ci=ci, ptype=m, wtcol=p_wt)
              npatches_added = npatches_added + 1
           end if
        end do
+
     end if
 
     SHR_ASSERT_FL(nlunits_added == nlunits, sourcefile, __LINE__)
diff --git a/src/main/surfrdMod.F90 b/src/main/surfrdMod.F90
index 27da9e24df..78c9d8492e 100644
--- a/src/main/surfrdMod.F90
+++ b/src/main/surfrdMod.F90
@@ -14,7 +14,7 @@ module surfrdMod
   use landunit_varcon , only : numurbl
   use clm_varcon      , only : grlnd
   use clm_varctl      , only : iulog
-  use clm_varctl      , only : use_cndv, use_crop
+  use clm_varctl      , only : use_cndv, use_crop, use_fates
   use surfrdUtilsMod  , only : check_sums_equal_1, collapse_crop_types
   use surfrdUtilsMod  , only : collapse_to_dominant, collapse_crop_var, collapse_individual_lunits
   use ncdio_pio       , only : file_desc_t, var_desc_t, ncd_pio_openfile, ncd_pio_closefile
@@ -251,7 +251,7 @@ end subroutine surfrd_get_data
 
 !-----------------------------------------------------------------------
 
-  subroutine surfrd_get_num_patches (lfsurdat, actual_maxsoil_patches, actual_numcft)
+  subroutine surfrd_get_num_patches (lfsurdat, actual_maxsoil_patches, actual_numpft, actual_numcft)
     !
     ! !DESCRIPTION:
     ! Read maxsoil_patches and numcft from the surface dataset
@@ -263,17 +263,25 @@ subroutine surfrd_get_num_patches (lfsurdat, actual_maxsoil_patches, actual_numc
     character(len=*), intent(in) :: lfsurdat  ! surface dataset filename
     integer, intent(out) :: actual_maxsoil_patches  ! value from surface dataset
     integer, intent(out) :: actual_numcft           ! cft value from sfc dataset
+    integer, intent(out) :: actual_numpft           ! pft value from sfc dataset
+    
     !
     ! !LOCAL VARIABLES:
     character(len=256):: locfn                ! local file name
     type(file_desc_t) :: ncid                 ! netcdf file id
     integer :: dimid                          ! netCDF dimension id
     logical :: cft_dim_exists                 ! dimension exists on dataset
+    integer :: check_numpft                   ! Surface dataset count of numpft, should
+                                              ! match maxsoil_patches - actual_numcft
     character(len=32) :: subname = 'surfrd_get_num_patches'  ! subroutine name
     !-----------------------------------------------------------------------
 
     if (masterproc) then
-       write(iulog,*) 'Attempting to read maxsoil_patches and numcft from the surface data .....'
+       if(use_fates)then
+          write(iulog,*) 'Attempting to read numcft from the surface data .....'
+       else
+          write(iulog,*) 'Attempting to read maxsoil_patches and numcft from the surface data .....'
+       end if
        if (lfsurdat == ' ') then
           write(iulog,*)'lfsurdat must be specified'
           call endrun(msg=errMsg(sourcefile, __LINE__))
@@ -284,15 +292,28 @@ subroutine surfrd_get_num_patches (lfsurdat, actual_maxsoil_patches, actual_numc
     call getfil( lfsurdat, locfn, 0 )
     call ncd_pio_openfile (ncid, trim(locfn), 0)
 
-    ! Read maxsoil_patches and numcft
-    call ncd_inqdlen(ncid, dimid, actual_maxsoil_patches, 'lsmpft')
+    ! Read numcft
     call ncd_inqdid(ncid, 'cft', dimid, cft_dim_exists)
     if ( cft_dim_exists ) then
        call ncd_inqdlen(ncid, dimid, actual_numcft, 'cft')
     else
        actual_numcft = 0
     end if
+    
+    ! Read maxsoil_patches
+    call ncd_inqdlen(ncid, dimid, actual_maxsoil_patches, 'lsmpft')
+    actual_numpft = actual_maxsoil_patches - actual_numcft
+
+    call ncd_inqdlen(ncid, dimid, check_numpft, 'natpft')
 
+    if(check_numpft.ne.actual_numpft)then
+       write(iulog,*)'the sum of the cftdim and the natpft dim should match the lsmpft dim in the surface file'
+       write(iulog,*)'natpft: ',check_numpft
+       write(iulog,*)'lsmpft: ',actual_maxsoil_patches
+       write(iulog,*)'cft: ',actual_numcft
+       call endrun(msg=errMsg(sourcefile, __LINE__))
+    end if
+    
     if ( masterproc )then
        write(iulog,*) 'Successfully read maxsoil_patches and numcft from the surface data'
        write(iulog,*)
@@ -455,7 +476,7 @@ subroutine surfrd_special(begg, endg, ncid, ns)
 
   end subroutine surfrd_special
 
-!-----------------------------------------------------------------------
+  !-----------------------------------------------------------------------
   subroutine surfrd_cftformat( ncid, begg, endg, wt_cft, fert_cft, cftsize, natpft_size )
     !
     ! !DESCRIPTION:
@@ -463,7 +484,7 @@ subroutine surfrd_cftformat( ncid, begg, endg, wt_cft, fert_cft, cftsize, natpft
     !     crop landunit and read in as Crop Function Types.
     ! !USES:
     use clm_instur      , only : wt_nat_patch, irrig_method
-    use clm_varpar      , only : cft_size, cft_lb, natpft_lb
+    use clm_varpar      , only : cft_size, cft_lb, natpft_lb, cft_ub, natpft_ub
     use IrrigationMod   , only : irrig_method_unset
     ! !ARGUMENTS:
     implicit none
@@ -489,7 +510,7 @@ subroutine surfrd_cftformat( ncid, begg, endg, wt_cft, fert_cft, cftsize, natpft
 
     call check_dim_size(ncid, 'cft',    cftsize)
     call check_dim_size(ncid, 'natpft', natpft_size)
-
+    
     call ncd_io(ncid=ncid, varname='PCT_CFT', flag='read', data=wt_cft, &
             dim1name=grlnd, readvar=readvar)
     if (.not. readvar) call endrun( msg=' ERROR: PCT_CFT NOT on surfdata file'//errMsg(sourcefile, __LINE__))
@@ -524,10 +545,113 @@ subroutine surfrd_cftformat( ncid, begg, endg, wt_cft, fert_cft, cftsize, natpft
     if (.not. readvar) call endrun( msg=' ERROR: PCT_NAT_PFT NOT on surfdata file'//errMsg(sourcefile, __LINE__))
     wt_nat_patch(begg:,natpft_lb:natpft_size-1+natpft_lb) = array2D(begg:,:)
     deallocate( array2D )
+ 
 
   end subroutine surfrd_cftformat
 
-!-----------------------------------------------------------------------
+  subroutine surfrd_wtfates( ncid, begg, endg )
+
+    !--------------------------------------------------------------------------
+    !     This routine evaluates the natural and crop functional
+    !     type fractions in the surface file and returns them to
+    !     a single, concatenated vector.  These weights
+    !     are only used for a satellite phenology run.
+    !     Note that FATES will actually allocate a different number of patches
+    !     and will use a mapping table to connect its own pft and cft
+    !     definitions to those it finds in the surface file.
+    !--------------------------------------------------------------------------
+    
+    ! !USES:
+    use clm_instur      , only : wt_nat_patch, wt_lunit
+    use clm_varpar      , only : cft_size, surfpft_lb, surfpft_ub
+    use landunit_varcon , only : istsoil, istcrop
+    
+    ! !ARGUMENTS:
+    implicit none
+    type(file_desc_t), intent(inout) :: ncid         ! netcdf id
+    integer          , intent(in)    :: begg, endg
+
+    !
+    ! !LOCAL VARIABLES:
+    logical  :: readvar                              ! is variable on dataset
+    real(r8),pointer :: array2d_pft(:,:)                 ! local array
+    real(r8),pointer :: array2d_cft(:,:)                 ! local array
+    integer :: g,p
+    integer :: cft_dimlen,natpft_dimlen,dimid
+    
+    character(len=32) :: subname = 'surfrd_fates'! subroutine name
+    
+    call ncd_inqdlen(ncid, dimid, cft_dimlen, 'cft')
+    call ncd_inqdlen(ncid, dimid, natpft_dimlen, 'natpft')
+
+    ! double check that cft_dimlen+natpft_dimlen = natpft_size
+    if((cft_dimlen+natpft_dimlen).ne.(surfpft_ub-surfpft_lb+1))then
+       call endrun( msg=' ERROR: PCT+CFT dimlen does not match array size for wt_nat_patch when fates is on'//errMsg(sourcefile, __LINE__))
+    end if
+    
+    allocate( array2d_cft(begg:endg,1:cft_dimlen) )
+    allocate( array2d_pft(begg:endg,1:natpft_dimlen) )
+    
+    call ncd_io(ncid=ncid, varname='PCT_CFT', flag='read', data=array2d_cft, &
+         dim1name=grlnd, readvar=readvar)
+    if (.not. readvar) call endrun( msg=' ERROR: PCT_CFT NOT on surfdata file'//errMsg(sourcefile, __LINE__))
+    
+    call ncd_io(ncid=ncid, varname='PCT_NAT_PFT', flag='read', data=array2d_pft, &
+         dim1name=grlnd, readvar=readvar)
+    if (.not. readvar) call endrun( msg=' ERROR: PCT_NAT_PFT NOT on surfdata file'//errMsg(sourcefile, __LINE__))
+
+    ! In fates, all the weights in both the cft and pfts go into this array
+    ! It is only used by SP mode, and it can choose what PFTs to align with
+    
+    wt_nat_patch(begg:,0:natpft_dimlen-1) = array2d_pft(begg:,:)
+    wt_nat_patch(begg:,natpft_dimlen:natpft_dimlen+cft_dimlen-1) = array2d_cft(begg:,:)
+
+    if(.false.)then
+       ! Scale the weights by the lu weights from the dataset
+       do g = begg, endg
+          
+          if((wt_lunit(g,istsoil)+wt_lunit(g,istcrop))>0._r8)then
+             do p = 0,natpft_dimlen-1
+                wt_nat_patch(g,p) = wt_nat_patch(g,p) * wt_lunit(g,istsoil)/(wt_lunit(g,istsoil)+wt_lunit(g,istcrop))
+             end do
+             do p = natpft_dimlen,natpft_dimlen+cft_dimlen-1
+                wt_nat_patch(g,p) = wt_nat_patch(g,p) * wt_lunit(g,istcrop)/(wt_lunit(g,istsoil)+wt_lunit(g,istcrop))
+             end do
+          else
+             wt_nat_patch(g,:) = 0._r8
+             wt_nat_patch(g,0) = 100._r8
+          end if
+       end do
+       ! Add the crop weight to the natveg weight and zero the crop weight
+       wt_lunit(begg:,istsoil) = wt_lunit(begg:,istsoil) + wt_lunit(begg:,istcrop)
+       wt_lunit(begg:,istcrop) = 0._r8
+       
+    else
+       ! Legacy method
+       do g = begg, endg
+          if ( wt_lunit(g,istcrop) > 0.0_r8 )then
+             ! Move CFT over to PFT and do weighted average of the crop and soil parts
+             wt_nat_patch(g,0:natpft_dimlen-1) = wt_nat_patch(g,0:natpft_dimlen-1) * wt_lunit(g,istsoil)
+             wt_nat_patch(g,natpft_dimlen:natpft_dimlen+cft_dimlen-1)       = &
+                  wt_nat_patch(g,natpft_dimlen:natpft_dimlen+cft_dimlen-1) * wt_lunit(g,istcrop)
+             wt_lunit(g,istsoil) = (wt_lunit(g,istsoil) + wt_lunit(g,istcrop)) ! Add crop landunit to soil landunit
+             wt_nat_patch(g,:)   =  wt_nat_patch(g,:) / wt_lunit(g,istsoil)
+             wt_lunit(g,istcrop) = 0.0_r8                ! Zero out crop CFT's
+          else
+             wt_nat_patch(g,natpft_dimlen:natpft_dimlen+cft_dimlen-1) = 0.0_r8    ! Make sure generic crops are zeroed out
+          end if
+       end do
+       
+    end if
+    
+    deallocate(array2d_cft,array2d_pft)
+    
+
+  end subroutine surfrd_wtfates
+
+  
+  !-----------------------------------------------------------------------
+  
   subroutine surfrd_pftformat( begg, endg, ncid )
     !
     ! !DESCRIPTION:
@@ -551,6 +675,7 @@ subroutine surfrd_pftformat( begg, endg, ncid )
     SHR_ASSERT_ALL_FL((ubound(wt_nat_patch) == (/endg, natpft_size-1+natpft_lb/)), sourcefile, __LINE__)
 
     call check_dim_size(ncid, 'natpft', natpft_size)
+    
     ! If cft_size == 0, then we expect to be running with a surface dataset
     ! that does
     ! NOT have a PCT_CFT array (or CONST_FERTNITRO_CFT array), and thus does not have a 'cft' dimension.
@@ -586,7 +711,7 @@ subroutine surfrd_pftformat( begg, endg, ncid )
     call ncd_io(ncid=ncid, varname='PCT_NAT_PFT', flag='read', data=wt_nat_patch, &
          dim1name=grlnd, readvar=readvar)
     if (.not. readvar) call endrun( msg=' ERROR: PCT_NAT_PFT NOT on surfdata file'//errMsg(sourcefile, __LINE__))
-
+       
   end subroutine surfrd_pftformat
 
 !-----------------------------------------------------------------------
@@ -601,6 +726,7 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns, actual_numcft)
     use clm_instur      , only : wt_lunit, wt_nat_patch, wt_cft, fert_cft
     use landunit_varcon , only : istsoil, istcrop
     use surfrdUtilsMod  , only : convert_cft_to_pft
+
     !
     ! !ARGUMENTS:
     implicit none
@@ -617,7 +743,8 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns, actual_numcft)
     real(r8),pointer :: array2DCFT(:,:)              ! local 2D array for CFTs
     real(r8),pointer :: array2DFERT(:,:)             ! local 2D array for fertilizer
     character(len=32) :: subname = 'surfrd_veg_all'  ! subroutine name
-!-----------------------------------------------------------------------
+
+    !-----------------------------------------------------------------------
     !
     ! Read in variables that are handled the same for all formats
     !
@@ -638,6 +765,7 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns, actual_numcft)
 
     deallocate(arrayl)
 
+    
     ! Check the file format for CFT's and handle accordingly
     if ( actual_numcft > 0 ) then
        if ( create_crop_landunit )then
@@ -658,12 +786,12 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns, actual_numcft)
           if ( use_fates ) then
              if ( masterproc ) write(iulog,*) "WARNING: When fates is on we allow new CFT based surface datasets ", &
                                               "to be used with create_crop_land FALSE"
-             cftsize = 2
-             allocate(array2DCFT (begg:endg,cft_lb:cftsize-1+cft_lb))
-             allocate(array2DFERT(begg:endg,cft_lb:cftsize-1+cft_lb))
-             call surfrd_cftformat( ncid, begg, endg, array2DCFT, array2DFERT, cftsize, natpft_size-cftsize ) ! Read crops in as CFT's
-             call convert_cft_to_pft( begg, endg, cftsize, array2DCFT )                          ! Convert from CFT to natural veg. landunit
-             fert_cft(begg:,cft_lb:) = 0.0_r8
+
+             call surfrd_wtfates( ncid, begg, endg )
+             ! Set the weighting on the crop patches to zero
+             fert_cft(begg:,cft_lb:cft_ub) = 0.0_r8
+             wt_cft(begg:,cft_lb:cft_ub) = 0.0_r8
+
           else
              call endrun( msg=' ERROR: New format surface datasets require create_crop_landunit TRUE'//errMsg(sourcefile, __LINE__))
           end if
@@ -686,27 +814,30 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns, actual_numcft)
                ' must also have a separate crop landunit, and vice versa)'//&
                errMsg(sourcefile, __LINE__))
     end if
+    
     ! Convert from percent to fraction
     wt_lunit(begg:endg,istsoil) = wt_lunit(begg:endg,istsoil) / 100._r8
     wt_lunit(begg:endg,istcrop) = wt_lunit(begg:endg,istcrop) / 100._r8
     wt_nat_patch(begg:endg,:)   = wt_nat_patch(begg:endg,:) / 100._r8
-    wt_cft(begg:endg,:) = wt_cft(begg:endg,:) / 100._r8
+    wt_cft(begg:endg,:)         = wt_cft(begg:endg,:) / 100._r8
 
     ! Check sum of vegetation adds to 1
     call check_sums_equal_1(wt_nat_patch, begg, 'wt_nat_patch', subname)
+
     ! if ( use_fates ) wt_cft = 0 because called convert_cft_to_pft, else...
     if ( .not. use_fates ) then
        ! Check sum of vegetation adds to 1
        call check_sums_equal_1(wt_cft, begg, 'wt_cft', subname)
     end if
+
     ! Call collapse_crop_types: allows need to maintain only 78-pft input data
     ! For use_crop = .false. collapsing 78->16 pfts or 16->16 or some new
-    !    configuration
+       !    configuration
     ! For use_crop = .true. most likely collapsing 78 to the list of crops for
     !    which the CLM includes parameterizations
     ! The call collapse_crop_types also appears in subroutine dyncrop_interp
     call collapse_crop_types(wt_cft(begg:endg,:), fert_cft(begg:endg,:), cft_size, begg, endg, verbose=.true.)
-
+    
     ! Collapse crop variables as needed
     ! The call to collapse_crop_var also appears in subroutine dyncrop_interp
     ! - fert_cft TODO Is this call redundant because it simply sets the crop
@@ -719,8 +850,10 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns, actual_numcft)
     ! - Pfts could be up to 14 before collapsing if create_crop_landunit = .T.
     ! - Pfts could be up to 16 before collapsing if create_crop_landunit = .F.
     ! TODO Add the same call to subroutine dynpft_interp for transient runs
+    
     call collapse_to_dominant(wt_nat_patch(begg:endg,:), natpft_lb, natpft_ub, &
-                              begg, endg, n_dom_pfts)
+         begg, endg, n_dom_pfts)
+    
   end subroutine surfrd_veg_all
 
   !-----------------------------------------------------------------------
diff --git a/src/main/surfrdUtilsMod.F90 b/src/main/surfrdUtilsMod.F90
index b48d8cbf74..0763d43a16 100644
--- a/src/main/surfrdUtilsMod.F90
+++ b/src/main/surfrdUtilsMod.F90
@@ -7,7 +7,7 @@ module surfrdUtilsMod
   ! !USES:
 #include "shr_assert.h"
   use shr_kind_mod , only : r8 => shr_kind_r8
-  use clm_varctl   , only : iulog
+  use clm_varctl   , only : iulog,use_fates
   use abortutils   , only : endrun
   use shr_log_mod  , only : errMsg => shr_log_errMsg
   use spmdMod      , only : masterproc
@@ -131,10 +131,11 @@ subroutine convert_cft_to_pft( begg, endg, cftsize, wt_cft )
     !
     ! !LOCAL VARIABLES:
     integer :: g    ! index
-!-----------------------------------------------------------------------
+    !-----------------------------------------------------------------------
+
     SHR_ASSERT_ALL_FL((ubound(wt_cft) == (/endg, cftsize/)), sourcefile, __LINE__)
     SHR_ASSERT_ALL_FL((ubound(wt_nat_patch) == (/endg, nc3crop+cftsize-1/)), sourcefile, __LINE__)
-
+    
     do g = begg, endg
        if ( wt_lunit(g,istcrop) > 0.0_r8 )then
           ! Move CFT over to PFT and do weighted average of the crop and soil parts
diff --git a/src/utils/clmfates_interfaceMod.F90 b/src/utils/clmfates_interfaceMod.F90
index f4a6bd8c6c..9a04a9d66f 100644
--- a/src/utils/clmfates_interfaceMod.F90
+++ b/src/utils/clmfates_interfaceMod.F90
@@ -45,10 +45,12 @@ module CLMFatesInterfaceMod
    use TemperatureType   , only : temperature_type
    use EnergyFluxType    , only : energyflux_type
    use SoilStateType     , only : soilstate_type
+   use CNProductsMod     , only : cn_products_type
    use clm_varctl        , only : iulog
    use clm_varctl        , only : fates_parteh_mode
    use clm_varctl        , only : use_fates
    use clm_varctl        , only : fates_spitfire_mode
+   use clm_varctl        , only : use_fates_tree_damage
    use clm_varctl        , only : use_fates_planthydro
    use clm_varctl        , only : use_fates_cohort_age_tracking
    use clm_varctl        , only : use_fates_ed_st3
@@ -65,8 +67,7 @@ module CLMFatesInterfaceMod
    use clm_varcon        , only : spval
    use clm_varcon        , only : denice
    use clm_varcon        , only : ispval
-
-   use clm_varpar        , only : natpft_size, natpft_ub, natpft_lb
+   use clm_varpar        , only : surfpft_lb,surfpft_ub
    use clm_varpar        , only : numrad
    use clm_varpar        , only : ivis
    use clm_varpar        , only : inir
@@ -109,7 +110,8 @@ module CLMFatesInterfaceMod
    ! Used FATES Modules
    use FatesInterfaceMod , only : fates_interface_type
    use FatesInterfaceMod, only : FatesInterfaceInit, FatesReportParameters
-   use FatesInterfaceMod, only : SetFatesGlobalElements
+   use FatesInterfaceMod, only : SetFatesGlobalElements1
+   use FatesInterfaceMod, only : SetFatesGlobalElements2
    use FatesInterfaceMod     , only : allocate_bcin
    use FatesInterfaceMod     , only : allocate_bcout
    use FatesInterfaceMod     , only : allocate_bcpconst
@@ -119,12 +121,14 @@ module CLMFatesInterfaceMod
    use FatesInterfaceMod     , only : set_fates_ctrlparms
    use FatesInterfaceMod     , only : UpdateFatesRMeansTStep
    use FatesInterfaceMod     , only : InitTimeAveragingGlobals
+   
    use FatesHistoryInterfaceMod, only : fates_hist
    use FatesRestartInterfaceMod, only : fates_restart_interface_type
 
    use EDTypesMod            , only : ed_patch_type
    use PRTGenericMod         , only : num_elements
    use FatesInterfaceTypesMod, only : hlm_stepsize
+   use FatesInterfaceTypesMod, only : fates_maxPatchesPerSite
    use EDMainMod             , only : ed_ecosystem_dynamics
    use EDMainMod             , only : ed_update_site
    use EDInitMod             , only : zero_site
@@ -219,6 +223,7 @@ module CLMFatesInterfaceMod
       procedure, public  :: ComputeRootSoilFlux
       procedure, public  :: wrap_hydraulics_drive
       procedure, public  :: WrapUpdateFatesRmean
+      procedure, public  :: wrap_WoodProducts
       
    end type hlm_fates_interface_type
 
@@ -239,12 +244,91 @@ module CLMFatesInterfaceMod
    character(len=*), parameter, private :: sourcefile = &
         __FILE__
 
-   public  :: CLMFatesGlobals
+   public  :: CLMFatesGlobals1
+   public  :: CLMFatesGlobals2
 
  contains
 
+   subroutine CLMFatesGlobals1(surf_numpft,surf_numcft,maxsoil_patches)
+
+     ! --------------------------------------------------------------------------------
+     ! This is the first call to fates
+     ! We open the fates parameter file. And use that and some info on
+     ! namelist variables to determine how many patches need to be allocated
+     ! in CTSM
+     ! --------------------------------------------------------------------------------
+
+     integer,intent(in)                             :: surf_numpft
+     integer,intent(in)                             :: surf_numcft
+     integer,intent(out)                            :: maxsoil_patches
+     integer                                        :: pass_biogeog
+     integer                                        :: pass_nocomp
+     integer                                        :: pass_sp
+     integer                                        :: pass_masterproc
+     logical                                        :: verbose_output
+     
+     call t_startf('fates_globals1')
+
+     if (use_fates) then
+
+        verbose_output = .false.
+        call FatesInterfaceInit(iulog, verbose_output)
+
+        ! Force FATES parameters that are recieve type, to the unset value
+        call set_fates_ctrlparms('flush_to_unset')
+
+        ! Send parameters individually
+        
+        if(use_fates_fixed_biogeog)then
+           pass_biogeog = 1
+        else
+           pass_biogeog = 0
+        end if
+        call set_fates_ctrlparms('use_fixed_biogeog',ival=pass_biogeog)
+        
+        if(use_fates_nocomp)then
+           pass_nocomp = 1
+        else
+           pass_nocomp = 0
+        end if
+        call set_fates_ctrlparms('use_nocomp',ival=pass_nocomp)
+        
+        if(use_fates_sp)then
+           pass_sp = 1
+        else
+           pass_sp = 0
+        end if
+        call set_fates_ctrlparms('use_sp',ival=pass_sp)
+        
+        if(masterproc)then
+           pass_masterproc = 1
+        else
+           pass_masterproc = 0
+        end if
+        call set_fates_ctrlparms('masterproc',ival=pass_masterproc)
+
+     end if
+
+     ! The following call reads in the parameter file
+     ! and then uses that to determine the number of patches
+     ! FATES requires. We pass that to CLM here
+     ! so that it can perform some of its allocations.
+     ! During init 2, we will perform more size checks
+     ! and allocations on the FATES side, which require
+     ! some allocations from CLM (like soil layering)
+
+     call SetFatesGlobalElements1(use_fates,surf_numpft,surf_numcft)
+
+     maxsoil_patches = fates_maxPatchesPerSite
+     
+     call t_stopf('fates_globals1')
+
+     return
+   end subroutine CLMFatesGlobals1
+
+   ! ===================================================================================
 
-   subroutine CLMFatesGlobals()
+   subroutine CLMFatesGlobals2()
 
      ! --------------------------------------------------------------------------------
      ! This is one of the first calls to fates
@@ -255,9 +339,7 @@ subroutine CLMFatesGlobals()
      ! is used in the history file, we also transfer
      ! over the NL variables to FATES global settings.
      ! --------------------------------------------------------------------------------
-
-     logical                                        :: verbose_output
-     integer                                        :: pass_masterproc
+       
      integer                                        :: pass_vertsoilc
      integer                                        :: pass_ch4
      integer                                        :: pass_spitfire
@@ -270,19 +352,13 @@ subroutine CLMFatesGlobals()
      integer                                        :: pass_inventory_init
      integer                                        :: pass_is_restart
      integer                                        :: pass_cohort_age_tracking
-     integer                                        :: pass_biogeog
-     integer                                        :: pass_nocomp
-     integer                                        :: pass_sp
+     integer                                        :: pass_tree_damage
 
-     call t_startf('fates_globals')
+     call t_startf('fates_globals2')
 
      if (use_fates) then
 
-        verbose_output = .false.
-        call FatesInterfaceInit(iulog, verbose_output)
-
-        ! Force FATES parameters that are recieve type, to the unset value
-        call set_fates_ctrlparms('flush_to_unset')
+        
 
         ! Send parameters individually
         call set_fates_ctrlparms('num_sw_bbands',ival=numrad)
@@ -293,8 +369,7 @@ subroutine CLMFatesGlobals()
         call set_fates_ctrlparms('hlm_name',cval='CLM')
         call set_fates_ctrlparms('hio_ignore_val',rval=spval)
         call set_fates_ctrlparms('soilwater_ipedof',ival=get_ipedof(0))
-        call set_fates_ctrlparms('max_patch_per_site',ival=(natpft_size-1))
-
+        
         call set_fates_ctrlparms('parteh_mode',ival=fates_parteh_mode)
 
         ! CTSM-FATES is not fully coupled (yet)
@@ -309,6 +384,13 @@ subroutine CLMFatesGlobals()
         elseif(decomp_method == no_soil_decomp ) then
            call set_fates_ctrlparms('decomp_method',cval='NONE')
         end if
+
+        if(use_fates_tree_damage)then
+           pass_tree_damage = 1
+        else
+           pass_tree_damage = 0
+        end if
+        call set_fates_ctrlparms('use_tree_damage',ival=pass_tree_damage)
         
         ! These may be in a non-limiting status (ie when supplements)
         ! are added, but they are always allocated and cycled non-the less
@@ -332,7 +414,9 @@ subroutine CLMFatesGlobals()
         call set_fates_ctrlparms('sf_scalar_lightning_def',ival=scalar_lightning)
         call set_fates_ctrlparms('sf_successful_ignitions_def',ival=successful_ignitions)
         call set_fates_ctrlparms('sf_anthro_ignitions_def',ival=anthro_ignitions)
-        call set_fates_ctrlparms('sf_anthro_suppression_def',ival=anthro_suppression)
+
+        ! This has no variable on the FATES side yet (RGK)
+        !call set_fates_ctrlparms('sf_anthro_suppression_def',ival=anthro_suppression)
 
         if(is_restart()) then
            pass_is_restart = 1
@@ -352,28 +436,6 @@ subroutine CLMFatesGlobals()
         pass_vertsoilc = 1
         call set_fates_ctrlparms('use_vertsoilc',ival=pass_vertsoilc)
 
-        if(use_fates_fixed_biogeog)then
-           pass_biogeog = 1
-        else
-           pass_biogeog = 0
-        end if
-        call set_fates_ctrlparms('use_fixed_biogeog',ival=pass_biogeog)
-
-        if(use_fates_nocomp)then
-           pass_nocomp = 1
-	   else
-           pass_nocomp = 0
-	   end if
-        call set_fates_ctrlparms('use_nocomp',ival=pass_nocomp)
-
-        if(use_fates_sp)then
-           pass_sp = 1
-              else
-           pass_sp = 0
-              end if
-        call set_fates_ctrlparms('use_sp',ival=pass_sp)
-
-
         if(use_fates_ed_st3) then
            pass_ed_st3 = 1
         else
@@ -438,12 +500,6 @@ subroutine CLMFatesGlobals()
 
         call set_fates_ctrlparms('inventory_ctrl_file',cval=fates_inventory_ctrl_filename)
 
-        if(masterproc)then
-           pass_masterproc = 1
-        else
-           pass_masterproc = 0
-        end if
-        call set_fates_ctrlparms('masterproc',ival=pass_masterproc)
 
         ! Check through FATES parameters to see if all have been set
         call set_fates_ctrlparms('check_allset')
@@ -461,13 +517,12 @@ subroutine CLMFatesGlobals()
      ! (Note: this needs to be called when use_fates=.false. as well, becuase
      ! it will return some nominal dimension sizes of 1
 
-     call SetFatesGlobalElements(use_fates)
+     call SetFatesGlobalElements2(use_fates)
 
-     call t_stopf('fates_globals')
+     call t_stopf('fates_globals2')
 
      return
-   end subroutine CLMFatesGlobals
-
+   end subroutine CLMFatesGlobals2
 
    ! ===================================================================================
   
@@ -638,7 +693,7 @@ subroutine init(this, bounds_proc )
 
             ndecomp = col%nbedrock(c)
 
-            call allocate_bcin(this%fates(nc)%bc_in(s),col%nbedrock(c),ndecomp, num_harvest_inst)
+            call allocate_bcin(this%fates(nc)%bc_in(s),col%nbedrock(c),ndecomp, num_harvest_inst,surfpft_lb,surfpft_ub)
             call allocate_bcout(this%fates(nc)%bc_out(s),col%nbedrock(c),ndecomp)
             call zero_bcs(this%fates(nc),s)
 
@@ -652,12 +707,12 @@ subroutine init(this, bounds_proc )
             this%fates(nc)%bc_in(s)%pft_areafrac(:)=0._r8
             ! initialize static layers for reduced complexity FATES versions from HLM
             ! maybe make this into a subroutine of it's own later.
-            do m = natpft_lb,natpft_ub
-               ft = m-natpft_lb
+            do m = surfpft_lb,surfpft_ub
+               ft = m - surfpft_lb
                this%fates(nc)%bc_in(s)%pft_areafrac(ft)=wt_nat_patch(g,m)
             end do
 
-            if(abs(sum(this%fates(nc)%bc_in(s)%pft_areafrac(natpft_lb:natpft_ub))-1.0_r8).gt.1.0e-9)then
+            if(abs(sum(this%fates(nc)%bc_in(s)%pft_areafrac(surfpft_lb:surfpft_ub))-1.0_r8).gt.1.0e-9)then
                write(iulog,*) 'pft_area error in interfc ',s, sum(this%fates(nc)%bc_in(s)%pft_areafrac(:))-1.0_r8
                call endrun(msg=errMsg(sourcefile, __LINE__))
               endif
@@ -898,7 +953,7 @@ subroutine dynamics_driv(this, nc, bounds_clump,      &
          ! in FATES.
          ! N.B. Fow now these are fixed values pending HLM updates.
          if(use_fates_sp)then
-           do ft = natpft_lb,natpft_ub !set of pfts in HLM
+           do ft = surfpft_lb,surfpft_ub
                ! here we are mapping from P space in the HLM to FT space in the sp_input arrays.
                p = ft + col%patchi(c) ! for an FT of 1 we want to use
                this%fates(nc)%bc_in(s)%hlm_sp_tlai(ft) = canopystate_inst%tlai_patch(p)
@@ -1537,7 +1592,7 @@ subroutine restart( this, bounds_proc, ncid, flag, waterdiagnosticbulk_inst, &
                if(use_fates_sp)then
                   do s = 1,this%fates(nc)%nsites
                      c = this%f2hmap(nc)%fcolumn(s)
-                     do ft = natpft_lb,natpft_ub !set of pfts in HLM
+                     do ft = surfpft_lb,surfpft_ub  !set of pfts in HLM
                         ! here we are mapping from P space in the HLM to FT space in the sp_input arrays.
                         p = ft + col%patchi(c) ! for an FT of 1 we want to use
                         this%fates(nc)%bc_in(s)%hlm_sp_tlai(ft) = canopystate_inst%tlai_patch(p)
@@ -1684,7 +1739,7 @@ subroutine init_coldstart(this, waterstatebulk_inst, waterdiagnosticbulk_inst, &
             if(use_fates_sp)then
                do s = 1,this%fates(nc)%nsites
                   c = this%f2hmap(nc)%fcolumn(s)
-                  do ft = natpft_lb,natpft_ub !set of pfts in HLM
+                  do ft = surfpft_lb,surfpft_ub
                      ! here we are mapping from P space in the HLM to FT space in the sp_input arrays.
                      p = ft + col%patchi(c) ! for an FT of 1 we want to use
                      this%fates(nc)%bc_in(s)%hlm_sp_tlai(ft) = canopystate_inst%tlai_patch(p)
@@ -2268,6 +2323,54 @@ end subroutine wrap_accumulatefluxes
 
  ! ======================================================================================
 
+ subroutine wrap_WoodProducts(this, bounds_clump, fc, filterc, c_products_inst)
+
+   ! !ARGUMENTS:
+   class(hlm_fates_interface_type), intent(inout) :: this
+   type(bounds_type)              , intent(in)    :: bounds_clump
+   integer                        , intent(in)    :: fc                   ! size of column filter
+   integer                        , intent(in)    :: filterc(fc)          ! column filter
+   type(cn_products_type)         , intent(inout) :: c_products_inst
+   
+   ! Locals
+   integer                                        :: s,c,icc,g
+   integer                                        :: nc
+
+   ! This wrapper is not active.  This is just place-holder code until
+   ! harvest-product flux is fully implemented.  RGK-05-2022
+   
+   !associate( & 
+   !     prod10c => c_products_inst%hrv_deadstem_to_prod10_grc, & 
+   !     prod100c => c_products_inst%hrv_deadstem_to_prod100_grc)
+     
+   !  nc = bounds_clump%clump_index
+     ! Loop over columns
+     do icc = 1,fc
+        c = filterc(icc)
+        g = col%gridcell(c)
+        s = this%f2hmap(nc)%hsites(c)
+        
+        ! Shijie: Pass harvested wood products to ELM variable
+        !     prod10c(g) = prod10c(g) + &
+        !          this%fates(nc)%bc_out(s)%hrv_deadstemc_to_prod10c
+        !     prod100c(g) = prod100c(g) + &
+        !          this%fates(nc)%bc_out(s)%hrv_deadstemc_to_prod100c
+        
+        ! RGK: THere is also a patch level variable
+        !do ifp = 1,this%fates(nc)%sites(s)%youngest_patch%patchno
+        !   p = ifp+col%patchi(c)
+        !   hrv_deadstemc_to_prod10c(p)  = 
+        !   hrv_deadstemc_to_prod100c(p)
+        !end do
+          
+    end do
+    
+    !  end associate
+  return
+ end subroutine wrap_WoodProducts
+ 
+ ! ======================================================================================
+
  subroutine wrap_canopy_radiation(this, bounds_clump, nc, &
          num_vegsol, filter_vegsol, coszen, fcansno,  surfalb_inst)
 
@@ -2642,6 +2745,7 @@ subroutine init_history_io(this,bounds_proc)
    use FatesIOVariableKindMod, only : site_can_r8, site_cnlf_r8, site_cnlfpft_r8
    use FatesIOVariableKindMod, only : site_height_r8, site_elem_r8, site_elpft_r8
    use FatesIOVariableKindMod, only : site_elcwd_r8, site_elage_r8, site_agefuel_r8
+   use FatesIOVariableKindMod, only : site_cdpf_r8, site_cdsc_r8
    use FatesIODimensionsMod, only : fates_bounds_type
 
 
@@ -2743,7 +2847,8 @@ subroutine init_history_io(this,bounds_proc)
              site_fuel_r8, site_cwdsc_r8, &
              site_can_r8,site_cnlf_r8, site_cnlfpft_r8, site_scag_r8, &
              site_scagpft_r8, site_agepft_r8, site_elem_r8, site_elpft_r8, &
-             site_elcwd_r8, site_elage_r8, site_agefuel_r8)
+             site_elcwd_r8, site_elage_r8, site_agefuel_r8, &
+             site_cdsc_r8, site_cdpf_r8)
 
 
            d_index = fates_hist%dim_kinds(dk_index)%dim2_index
@@ -3031,6 +3136,7 @@ subroutine hlm_bounds_to_fates_bounds(hlm, fates)
    use FatesIODimensionsMod, only : fates_bounds_type
    use FatesInterfaceTypesMod, only : nlevsclass, nlevage, nlevcoage
    use FatesInterfaceTypesMod, only : nlevheight
+   use FatesInterfaceTypesMod, only : nlevdamage
    use EDtypesMod,        only : nfsc
    use FatesLitterMod,    only : ncwd
    use EDtypesMod,        only : nlevleaf, nclmax
@@ -3113,7 +3219,15 @@ subroutine hlm_bounds_to_fates_bounds(hlm, fates)
    fates%agefuel_begin = 1
    fates%agefuel_end   = nlevage * nfsc
 
+   fates%cdpf_begin = 1
+   fates%cdpf_end = nlevdamage * numpft_fates * nlevsclass
+
+   fates%cdsc_begin = 1
+   fates%cdsc_end = nlevdamage * nlevsclass
 
+   fates%cdam_begin = 1
+   fates%cdam_end = nlevdamage
+   
    call t_stopf('fates_hlm2fatesbnds')
 
  end subroutine hlm_bounds_to_fates_bounds