Merge pull request #7 from glemieux/fates-landuse-v2-systest

Fix `PVT` system test and avoid getting `fluh_timeseries` when `use_fates_potentialveg = .true.`
olyson · May 31, 2024 · eb9f1f3 · eb9f1f3
2 parents 071d076 + 709c3da
commit eb9f1f3
Show file tree

Hide file tree

Showing 5 changed files with 82 additions and 53 deletions.
diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml
@@ -1557,7 +1557,7 @@ lnd/clm2/surfdata_esmf/NEON/surfdata_1x1_NEON_TOOL_hist_78pfts_CMIP6_simyr2000_c
 
 <!-- This one is just for testing until Issue #2304 is resolved -->
 <fluh_timeseries hgrid="4x5" sim_year_range="constant" use_fates=".true."
- >lnd/clm2/surfdata_map/fates-sci.1.68.3_api.31.0.0_tools.1.0.1/LUH2_states_transitions_management.timeseries_4x5_hist_simyr1850-2015_c231101.nc</fluh_timeseries>
+ >lnd/clm2/surfdata_map/fates-sci.1.68.3_api.31.0.0_tools.1.0.1/LUH2_states_transitions_management.timeseries_4x5_hist_simyr1650-2015_c240216.nc</fluh_timeseries>
 
 <!-- Fixation and Uptake of Nitrogen Model (FUN2.0) -->
 <use_fun                                                            >.false.</use_fun>

diff --git a/cime_config/SystemTests/pvt.py b/cime_config/SystemTests/pvt.py
@@ -8,12 +8,10 @@
     - use CLM_ACCELERATED_SPINUP?
 2) run a transient landuse case with use_fates_lupft
     - start from the restart file generated in (1)
-
 """
 from CIME.XML.standard_module_setup import *
 from CIME.SystemTests.system_tests_common import SystemTestsCommon
 from CIME.SystemTests.test_utils.user_nl_utils import append_to_user_nl_files
-from systemtest_utils import find_user_nl_option
 import shutil, glob, os
 
 logger = logging.getLogger(__name__)
@@ -51,6 +49,11 @@ def run_phase(self):
         orig_casevar = self._case.get_value("CASE")
         caseroot = self._case.get_value("CASEROOT")
 
+        # Set the run start date based on the desired starting reference case year
+        refcase_year = 1700
+        stop_n_pveg = 5
+        startyear_pveg = refcase_year - stop_n_pveg
+
         # clone the main case to create spinup case
         logger.info("PVT log:  cloning setup")
         clone_path = "{}.potveg".format(caseroot)
@@ -65,10 +68,11 @@ def run_phase(self):
         self._set_active_case(clone)
 
         # set the clone case values
-        stop_n_pveg = 1
         with clone:
-            # clone.set_value("CLM_ACCELERATED_SPINUP", "on")
+            clone.set_value("CLM_ACCELERATED_SPINUP", "off")
             clone.set_value("STOP_N", stop_n_pveg)
+            clone.set_value("STOP_OPTION", "nyears")
+            clone.set_value("RUN_STARTDATE", "{}-01-01".format(startyear_pveg))
 
         # Modify the spin up case to use the potential vegetation mode.
         # Checks for incompatible cases and necessary mapping files are
@@ -93,11 +97,33 @@ def run_phase(self):
         os.chdir(caseroot)
         self._set_active_case(orig_case)
 
+        # Copy restart files from spin up to the transient case run directory
+        # obtain rpointer files and necessary restart files from short term archiving directory
+        rundir = self._case.get_value("RUNDIR")
+
+        refdate = str(refcase_year) + '-01-01-00000'
+        rest_path = os.path.join(dout_sr, "rest", "{}".format(refdate))
+
+        for item in glob.glob("{}/*{}*".format(rest_path, refdate)):
+            link_name = os.path.join(rundir, os.path.basename(item))
+            if os.path.islink(link_name) and os.readlink(link_name) == item:
+                # Link is already set up correctly: do nothing
+                # (os.symlink raises an exception if you try to replace an
+                # existing file)
+                pass
+            else:
+                os.symlink(item, link_name)
+
+        for item in glob.glob("{}/*rpointer*".format(rest_path)):
+            shutil.copy(item, rundir)
+
+        # Update run case settings
         self._case.set_value("CLM_ACCELERATED_SPINUP", "off")
         self._case.set_value("RUN_TYPE", "hybrid")
         self._case.set_value("GET_REFCASE", False)
         self._case.set_value("RUN_REFCASE", "{}.potveg".format(orig_casevar))
-        self._case.set_value("RUN_REFDATE", "1700-01-01")
+        self._case.set_value("RUN_REFDATE", "{}-01-01".format(refcase_year))
+        self._case.set_value("RUN_STARTDATE", "{}-01-01".format(refcase_year))
         self._case.set_value("DOUT_S", False)
         self._case.flush()
 

diff --git a/cime_config/SystemTests/systemtest_utils.py b/cime_config/SystemTests/systemtest_utils.py
@@ -85,9 +85,8 @@ def run_python_script(caseroot, this_conda_env, command_in, tool_path):
     except:
         print(f"ERROR trying to run {tool_name}.")
         raise
-"""
-Read a user_nl file and return the namelist option if found
-"""
+
+# Read a user_nl file and return the namelist option if found
 def find_user_nl_option(caseroot, component, namelist_option):
 
     # This is a copy of the CIME _get_list_of_user_nl_files

diff --git a/src/dyn_subgrid/dynFATESLandUseChangeMod.F90 b/src/dyn_subgrid/dynFATESLandUseChangeMod.F90
@@ -107,6 +107,7 @@ subroutine dynFatesLandUseInit(bounds, landuse_filename)
 
     ! !USES:
     use clm_varctl            , only : use_cn, use_fates_luh, fates_harvest_mode
+    use clm_varctl            , only : use_fates_potentialveg
     use dynVarTimeUninterpMod , only : dyn_var_time_uninterp_type
     use dynTimeInfoMod        , only : YEAR_POSITION_START_OF_TIMESTEP
     use dynTimeInfoMod        , only : YEAR_POSITION_END_OF_TIMESTEP
@@ -148,55 +149,58 @@ subroutine dynFatesLandUseInit(bounds, landuse_filename)
     landuse_transitions = 0._r8
     landuse_harvest = 0._r8
 
-    if (use_fates_luh) then
+    ! Avoid initializing the landuse timeseries file if in fates potential vegetation mode
+    if (.not. use_fates_potentialveg) then
+       if (use_fates_luh) then
 
-       ! Generate the dyn_file_type object.
-       ! Start calls get_prev_date, whereas end calls get_curr_date
-       dynFatesLandUse_file = dyn_file_type(landuse_filename, YEAR_POSITION_END_OF_TIMESTEP)
+          ! Generate the dyn_file_type object.
+          ! Start calls get_prev_date, whereas end calls get_curr_date
+          dynFatesLandUse_file = dyn_file_type(landuse_filename, YEAR_POSITION_END_OF_TIMESTEP)
 
-       ! Get initial land use data from the fates luh2 timeseries dataset
-       num_points = (bounds%endg - bounds%begg + 1)
-       landuse_shape(1) = num_points ! Does this need an explicit array shape to be passed to the constructor?
-       do varnum = 1, num_landuse_transition_vars
-          landuse_transition_vars(varnum) = dyn_var_time_uninterp_type( &
-               dyn_file=dynFatesLandUse_file, varname=landuse_transition_varnames(varnum), &
-               dim1name=grlnd, conversion_factor=1.0_r8, &
-               do_check_sums_equal_1=.false., data_shape=landuse_shape)
-       end do
-       do varnum = 1, num_landuse_state_vars
-          landuse_state_vars(varnum) = dyn_var_time_uninterp_type( &
-               dyn_file=dynFatesLandUse_file, varname=landuse_state_varnames(varnum), &
-               dim1name=grlnd, conversion_factor=1.0_r8, &
-               do_check_sums_equal_1=.false., data_shape=landuse_shape)
-       end do
-
-       ! Get the harvest rate data from the fates luh2 timeseries dataset if enabled
-       if (trim(fates_harvest_mode) .eq. fates_harvest_luh_area .or. &
-           trim(fates_harvest_mode) .eq. fates_harvest_luh_mass) then
-
-          ! change the harvest varnames being used depending on the mode selected
-          if (trim(fates_harvest_mode) .eq. fates_harvest_luh_area ) then
-             landuse_harvest_varnames => landuse_harvest_area_varnames
-             landuse_harvest_units = landuse_harvest_area_units
-          elseif (trim(fates_harvest_mode) .eq. fates_harvest_luh_mass ) then
-             landuse_harvest_varnames => landuse_harvest_mass_varnames
-             landuse_harvest_units = landuse_harvest_mass_units
-          else
-             call endrun(msg=' undefined fates harvest mode selected'//errMsg(__FILE__, __LINE__))
-          end if
-
-          do varnum = 1, num_landuse_harvest_vars
-             landuse_harvest_vars(varnum) = dyn_var_time_uninterp_type( &
-                  dyn_file=dynFatesLandUse_file, varname=landuse_harvest_varnames(varnum), &
+          ! Get initial land use data from the fates luh2 timeseries dataset
+          num_points = (bounds%endg - bounds%begg + 1)
+          landuse_shape(1) = num_points ! Does this need an explicit array shape to be passed to the constructor?
+          do varnum = 1, num_landuse_transition_vars
+             landuse_transition_vars(varnum) = dyn_var_time_uninterp_type( &
+                  dyn_file=dynFatesLandUse_file, varname=landuse_transition_varnames(varnum), &
+                  dim1name=grlnd, conversion_factor=1.0_r8, &
+                  do_check_sums_equal_1=.false., data_shape=landuse_shape)
+          end do
+          do varnum = 1, num_landuse_state_vars
+             landuse_state_vars(varnum) = dyn_var_time_uninterp_type( &
+                  dyn_file=dynFatesLandUse_file, varname=landuse_state_varnames(varnum), &
                   dim1name=grlnd, conversion_factor=1.0_r8, &
                   do_check_sums_equal_1=.false., data_shape=landuse_shape)
           end do
+
+          ! Get the harvest rate data from the fates luh2 timeseries dataset if enabled
+          if (trim(fates_harvest_mode) .eq. fates_harvest_luh_area .or. &
+              trim(fates_harvest_mode) .eq. fates_harvest_luh_mass) then
+
+             ! change the harvest varnames being used depending on the mode selected
+             if (trim(fates_harvest_mode) .eq. fates_harvest_luh_area ) then
+                landuse_harvest_varnames => landuse_harvest_area_varnames
+                landuse_harvest_units = landuse_harvest_area_units
+             elseif (trim(fates_harvest_mode) .eq. fates_harvest_luh_mass ) then
+                landuse_harvest_varnames => landuse_harvest_mass_varnames
+                landuse_harvest_units = landuse_harvest_mass_units
+             else
+                call endrun(msg=' undefined fates harvest mode selected'//errMsg(__FILE__, __LINE__))
+             end if
+
+             do varnum = 1, num_landuse_harvest_vars
+                landuse_harvest_vars(varnum) = dyn_var_time_uninterp_type( &
+                     dyn_file=dynFatesLandUse_file, varname=landuse_harvest_varnames(varnum), &
+                     dim1name=grlnd, conversion_factor=1.0_r8, &
+                     do_check_sums_equal_1=.false., data_shape=landuse_shape)
+             end do
+          end if
        end if
-    end if
 
-    ! Since fates needs state data during initialization, make sure to call
-    ! the interpolation routine at the start
-    call dynFatesLandUseInterp(bounds,init_state=.true.)
+       ! Since fates needs state data during initialization, make sure to call
+       ! the interpolation routine at the start
+       call dynFatesLandUseInterp(bounds,init_state=.true.)
+    end if
 
   end subroutine dynFatesLandUseInit
 

diff --git a/src/dyn_subgrid/dynSubgridDriverMod.F90 b/src/dyn_subgrid/dynSubgridDriverMod.F90
@@ -206,7 +206,7 @@ subroutine dynSubgrid_driver(bounds_proc,
     ! OUTSIDE any loops over clumps in the driver.
     !
     ! !USES:
-    use clm_varctl              , only : use_cn, use_fates, use_fates_luh
+    use clm_varctl              , only : use_cn, use_fates, use_fates_luh, use_fates_potentialveg
     use dynInitColumnsMod       , only : initialize_new_columns
     use dynConsBiogeophysMod    , only : dyn_hwcontent_init, dyn_hwcontent_final
     use dynEDMod                , only : dyn_ED
@@ -295,7 +295,7 @@ subroutine dynSubgrid_driver(bounds_proc,
        call dynurban_interp(bounds_proc)
     end if
 
-    if (use_fates_luh) then
+    if (use_fates_luh .and. .not. use_fates_potentialveg) then
        call dynFatesLandUseInterp(bounds_proc)
     end if