Merge pull request #2434 from samsrabin/hillslopes-off-fsurdat

Move hillslope variables off surface dataset
ESCOMP · Oct 8, 2024 · f95d192 · f95d192
2 parents f80f4c4 + f2c5d37
commit f95d192
Show file tree

Hide file tree

Showing 14 changed files with 175 additions and 121 deletions.
diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm
@@ -2759,8 +2759,7 @@ SIMYR:    foreach my $sim_yr ( @sim_years ) {
 
   # this check has to be here and not earlier since use_init_interp is set here and hillslope is already set above in setup_logic_hillslope
   if ( &value_is_true($nl->get_value($useinitvar)) && value_is_true($nl->get_value("use_hillslope")) ) {
-     $log->warning("WARNING: You have set use_hillslope while $useinitvar is TRUE.\n This means all hillslope columns in a gridcell will read identical values" .
-                   " from initial conditions. If you are sure you want this behaviour:")
+     $log->warning("WARNING: You have set use_hillslope while $useinitvar is TRUE.\n This means all hillslope columns in a gridcell will read identical values from initial conditions, even if the initial conditions (finidat) file has hillslope information. If you are sure you want this behaviour, add -ignore_warnings to CLM_BLDNML_OPTS.")
   }
 
 } # end initial conditions
@@ -3637,6 +3636,10 @@ sub setup_logic_hillslope {
   if ( (! &value_is_true($use_hillslope)) && &value_is_true($use_hillslope_routing) ) {
       $log->fatal_error("Cannot turn on use_hillslope_routing when use_hillslope is off\n" );
   }
+  my $hillslope_file = $nl->get_value('hillslope_file');
+  if ( &value_is_true($use_hillslope) && ( ! defined($hillslope_file) ) ) {
+    $log->fatal_error("You must provide hillslope_file if use_hillslope is .true.\n" );
+  }
 }
 
 #-------------------------------------------------------------------------------

diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml
@@ -962,6 +962,11 @@ Full pathname datafile with fates parameters
 Full pathname of surface data file.
 </entry>
 
+<entry id="hillslope_file" type="char*256"     category="datasets"
+       input_pathname="abs" group="clm_inparm" valid_values="" >
+Full pathname of hillslope data file.
+</entry>
+
 <entry id="fsnowoptics" type="char*256" category="datasets" 
        input_pathname="abs" group="clm_inparm" valid_values="" >
 SNICAR (SNow, ICe, and Aerosol Radiative model) optical data file name

diff --git a/cime_config/testdefs/ExpectedTestFails.xml b/cime_config/testdefs/ExpectedTestFails.xml
@@ -245,13 +245,6 @@
     </phase>
   </test>
 
-  <test name="SMS_D_Ld5.5x5_amazon.I1850Clm51Bgc.derecho_gnu.clm-HillslopeC">
-    <phase name="RUN">
-      <status>FAIL</status>
-      <issue>#2423</issue>
-    </phase>
-  </test>
-
   <test name="ERS_D_Ld15.f45_f45_mg37.I2000Clm50FatesRs.derecho_intel.clm-FatesColdTwoStream">
     <phase name="COMPARE_base_rest">
       <status>FAIL</status>

diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml
@@ -1773,7 +1773,7 @@
       <machine name="derecho" compiler="intel" category="aux_clm"/>
     </machines>
     <options>
-      <option name="wallclock">01:40:00</option>
+      <option name="wallclock">02:30:00</option>
       <option name="comment"  >Repeat ERS_Ly5_P128x1...cropMonthOutput test with matrixcnOn</option>
     </options>
   </test>
@@ -3511,23 +3511,24 @@
     </options>
   </test>
 
-  <test name="SMS_D_Mmpi-serial_Ld5" grid="5x5_amazon" compset="I1850Clm60Bgc" testmods="clm/HillslopeC">
+  <test name="SMS_D_Mmpi-serial_Ld5" grid="5x5_amazon" compset="I2000Clm60Bgc" testmods="clm/HillslopeC">
     <machines>
       <machine name="derecho" compiler="gnu" category="aux_clm"/>
       <machine name="derecho" compiler="gnu" category="hillslope"/>
     </machines>
     <options>
       <option name="wallclock">00:20:00</option>
+      <option name="comment">Serial test of 5x5_amazon hillslopes</option>
     </options>
   </test>
 
-  <test name="SMS_D_Ld5" grid="5x5_amazon" compset="I1850Clm60Bgc" testmods="clm/HillslopeC">
+  <test name="SMS_D_Ld5" grid="5x5_amazon" compset="I2000Clm60Bgc" testmods="clm/HillslopeC">
     <machines>
-      <machine name="derecho" compiler="gnu" category="aux_clm"/>
+      <machine name="derecho" compiler="gnu" category="hillslope"/>
     </machines>
     <options>
       <option name="wallclock">00:20:00</option>
-      <option name="comment">Add a non-FATES test to aux_clm for 5x5_amazon that fails because of the ccs_config issue, can remove once this passes</option>
+      <option name="comment">Parallel test of 5x5_amazon hillslopes</option>
     </options>
   </test>
 

diff --git a/cime_config/testdefs/testmods_dirs/clm/Hillslope/shell_commands b/cime_config/testdefs/testmods_dirs/clm/Hillslope/shell_commands
@@ -1,7 +1,19 @@
 ./xmlchange CLM_BLDNML_OPTS="-bgc sp"
 DIN_LOC_ROOT=$(./xmlquery --value DIN_LOC_ROOT)
-meshfile=$DIN_LOC_ROOT/lnd/clm2/testdata/ESMFmesh_10x15_synthetic_cosphill_1.0.nc
-./xmlchange ATM_DOMAIN_MESH=${meshfile},LND_DOMAIN_MESH=${meshfile}
+
+# Set hillslope_file. Needed for any grids without default hillslope_file already set by CTSM.
+lnd_grid=$(./xmlquery --value LND_GRID)
+if [[ ${lnd_grid} == "10x15" ]]; then
+    # Synthetic data
+    hillslope_file='$DIN_LOC_ROOT/lnd/clm2/testdata/surfdata_10x15_hist_1850_78pfts_c240216.synth_hillslopes_241001.nc'
+elif [[ ${lnd_grid} == "5x5_amazon" ]]; then
+    # Real data
+    hillslope_file='/glade/derecho/scratch/samrabin/hillslopes_5x5_amazon/hand_analysis_global/combined/hilldata_5x5_amazon_hist_2000_78pfts_c240216.nc'
+else
+    echo "ERROR: Hillslope file not found for LND_GRID=${lnd_grid}" >&2
+    exit 1
+fi
+echo -e "hillslope_file = '${hillslope_file}'\n" >> user_nl_clm
 
 # -ignore_warnings is needed as long as we don't allow use_hillslope and use_init_interp together
 ./xmlchange --append CLM_BLDNML_OPTS=-ignore_warnings
diff --git a/cime_config/testdefs/testmods_dirs/clm/Hillslope/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/Hillslope/user_nl_clm
@@ -6,6 +6,4 @@ hillslope_transmissivity_method = 'LayerSum'
 hillslope_pft_distribution_method = 'PftLowlandUpland'
 hillslope_soil_profile_method = 'Uniform'
 
-fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_esmf/ctsm5.3.0/synthetic/surfdata_10x15_hist_2000_78pfts_c240905.synthetic_hillslopes3.nc'
-
 use_ssre = .false.
diff --git a/src/biogeophys/HillslopeHydrologyMod.F90 b/src/biogeophys/HillslopeHydrologyMod.F90
@@ -168,7 +168,7 @@ end subroutine check_aquifer_layer
 
   !-----------------------------------------------------------------------
 
-  subroutine InitHillslope(bounds,fsurdat)
+  subroutine InitHillslope(bounds, hillslope_file)
     !
     ! !DESCRIPTION:
     ! Initialize hillslope geomorphology from input dataset
@@ -187,7 +187,7 @@ subroutine InitHillslope(bounds,fsurdat)
     !
     ! !ARGUMENTS:
     type(bounds_type), intent(in) :: bounds
-    character(len=*) , intent(in) :: fsurdat    ! surface data file name
+    character(len=*) , intent(in) :: hillslope_file    ! hillslope data file name
     integer,  pointer     :: ihillslope_in(:,:) ! read in - integer
     integer,  pointer     :: ncolumns_hillslope_in(:) ! read in number of columns
     integer,  allocatable :: ncolumns_hillslope(:)    ! number of hillslope columns
@@ -224,9 +224,9 @@ subroutine InitHillslope(bounds,fsurdat)
     ! consistency check
     call check_aquifer_layer()
 
-    ! Open surface dataset to read in data below
+    ! Open hillslope dataset to read in data below
 
-    call getfil (fsurdat, locfn, 0)
+    call getfil (hillslope_file, locfn, 0)
     call ncd_pio_openfile (ncid, locfn, 0)
 
     allocate( &
@@ -248,7 +248,7 @@ subroutine InitHillslope(bounds,fsurdat)
 
     call ncd_io(ncid=ncid, varname='nhillcolumns', flag='read', data=ncolumns_hillslope_in, dim1name=grlnd, readvar=readvar)
     if (masterproc .and. .not. readvar) then
-       call endrun( 'ERROR:: nhillcolumns not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+       call endrun( 'ERROR:: nhillcolumns not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
     end if
     do l = bounds%begl,bounds%endl
        g = lun%gridcell(l)
@@ -266,7 +266,7 @@ subroutine InitHillslope(bounds,fsurdat)
 
     call ncd_io(ncid=ncid, varname='pct_hillslope', flag='read', data=fhillslope_in, dim1name=grlnd, readvar=readvar)
     if (masterproc .and. .not. readvar) then
-       call endrun( 'ERROR:: pct_hillslope not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+       call endrun( 'ERROR:: pct_hillslope not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
     end if
     do l = bounds%begl,bounds%endl
        g = lun%gridcell(l)
@@ -278,7 +278,7 @@ subroutine InitHillslope(bounds,fsurdat)
 
     call ncd_io(ncid=ncid, varname='hillslope_index', flag='read', data=ihillslope_in, dim1name=grlnd, readvar=readvar)
     if (masterproc .and. .not. readvar) then
-       call endrun( 'ERROR:: hillslope_index not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+       call endrun( 'ERROR:: hillslope_index not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
     end if
     do l = bounds%begl,bounds%endl
        g = lun%gridcell(l)
@@ -287,7 +287,7 @@ subroutine InitHillslope(bounds,fsurdat)
 
     call ncd_io(ncid=ncid, varname='column_index', flag='read', data=ihillslope_in, dim1name=grlnd, readvar=readvar)
     if (masterproc .and. .not. readvar) then
-       call endrun( 'ERROR:: column_index not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+       call endrun( 'ERROR:: column_index not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
     end if
     do l = bounds%begl,bounds%endl
        g = lun%gridcell(l)
@@ -296,7 +296,7 @@ subroutine InitHillslope(bounds,fsurdat)
 
     call ncd_io(ncid=ncid, varname='downhill_column_index', flag='read', data=ihillslope_in, dim1name=grlnd, readvar=readvar)
     if (masterproc .and. .not. readvar) then
-       call endrun( 'ERROR:: downhill_column_index not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+       call endrun( 'ERROR:: downhill_column_index not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
     end if
     do l = bounds%begl,bounds%endl
        g = lun%gridcell(l)
@@ -307,7 +307,7 @@ subroutine InitHillslope(bounds,fsurdat)
     allocate(fhillslope_in(bounds%begg:bounds%endg,max_columns_hillslope))
     call ncd_io(ncid=ncid, varname='hillslope_slope', flag='read', data=fhillslope_in, dim1name=grlnd, readvar=readvar)
     if (masterproc .and. .not. readvar) then
-       call endrun( 'ERROR:: hillslope_slope not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+       call endrun( 'ERROR:: hillslope_slope not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
     end if
 
     do l = bounds%begl,bounds%endl
@@ -317,7 +317,7 @@ subroutine InitHillslope(bounds,fsurdat)
 
     call ncd_io(ncid=ncid, varname='hillslope_aspect', flag='read', data=fhillslope_in, dim1name=grlnd, readvar=readvar)
     if (masterproc .and. .not. readvar) then
-       call endrun( 'ERROR:: hillslope_aspect not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+       call endrun( 'ERROR:: hillslope_aspect not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
     end if
 
     do l = bounds%begl,bounds%endl
@@ -327,15 +327,15 @@ subroutine InitHillslope(bounds,fsurdat)
 
     call ncd_io(ncid=ncid, varname='hillslope_area', flag='read', data=fhillslope_in, dim1name=grlnd, readvar=readvar)
     if (masterproc .and. .not. readvar) then
-       call endrun( 'ERROR:: hillslope_area not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+       call endrun( 'ERROR:: hillslope_area not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
     end if
     do l = bounds%begl,bounds%endl
        g = lun%gridcell(l)
        hill_area(l,:) = fhillslope_in(g,:)
     enddo
     call ncd_io(ncid=ncid, varname='hillslope_distance', flag='read', data=fhillslope_in, dim1name=grlnd, readvar=readvar)
     if (masterproc .and. .not. readvar) then
-       call endrun( 'ERROR:: hillslope_length not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+       call endrun( 'ERROR:: hillslope_distance not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
     end if
 
     do l = bounds%begl,bounds%endl
@@ -345,7 +345,7 @@ subroutine InitHillslope(bounds,fsurdat)
 
     call ncd_io(ncid=ncid, varname='hillslope_width', flag='read', data=fhillslope_in, dim1name=grlnd, readvar=readvar)
     if (masterproc .and. .not. readvar) then
-       call endrun( 'ERROR:: hillslope_width not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+       call endrun( 'ERROR:: hillslope_width not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
     end if
     do l = bounds%begl,bounds%endl
        g = lun%gridcell(l)
@@ -354,7 +354,7 @@ subroutine InitHillslope(bounds,fsurdat)
 
     call ncd_io(ncid=ncid, varname='hillslope_elevation', flag='read', data=fhillslope_in, dim1name=grlnd, readvar=readvar)
     if (masterproc .and. .not. readvar) then
-       call endrun( 'ERROR:: hillslope_height not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+       call endrun( 'ERROR:: hillslope_elevation not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
     end if
     do l = bounds%begl,bounds%endl
        g = lun%gridcell(l)
@@ -380,7 +380,7 @@ subroutine InitHillslope(bounds,fsurdat)
 
        call ncd_io(ncid=ncid, varname='hillslope_stream_depth', flag='read', data=fstream_in, dim1name=grlnd, readvar=readvar)
        if (masterproc .and. .not. readvar) then
-          call endrun( 'ERROR:: hillslope_stream_depth not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+          call endrun( 'ERROR:: hillslope_stream_depth not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
        end if
        do l = bounds%begl,bounds%endl
           g = lun%gridcell(l)
@@ -389,7 +389,7 @@ subroutine InitHillslope(bounds,fsurdat)
 
        call ncd_io(ncid=ncid, varname='hillslope_stream_width', flag='read', data=fstream_in, dim1name=grlnd, readvar=readvar)
        if (masterproc .and. .not. readvar) then
-          call endrun( 'ERROR:: hillslope_stream_width not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+          call endrun( 'ERROR:: hillslope_stream_width not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
        end if
        do l = bounds%begl,bounds%endl
           g = lun%gridcell(l)
@@ -398,7 +398,7 @@ subroutine InitHillslope(bounds,fsurdat)
 
        call ncd_io(ncid=ncid, varname='hillslope_stream_slope', flag='read', data=fstream_in, dim1name=grlnd, readvar=readvar)
        if (masterproc .and. .not. readvar) then
-          call endrun( 'ERROR:: hillslope_stream_slope not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+          call endrun( 'ERROR:: hillslope_stream_slope not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
        end if
        do l = bounds%begl,bounds%endl
           g = lun%gridcell(l)
@@ -504,11 +504,11 @@ subroutine InitHillslope(bounds,fsurdat)
              enddo
           end if
 
-          ! if missing hillslope information on surface dataset,
+          ! if missing hillslope information on dataset,
           ! call endrun
           if (ncolumns_hillslope(l) > 0 .and. sum(hillslope_area) == 0._r8 .and. masterproc) then
              write(iulog,*) 'Problem with input data: nhillcolumns is non-zero, but hillslope area is zero'
-             write(iulog,*) 'Check surface data for gridcell at (lon/lat): ', grc%londeg(g),grc%latdeg(g)
+             write(iulog,*) 'Check hillslope data for gridcell at (lon/lat): ', grc%londeg(g),grc%latdeg(g)
              call endrun( 'ERROR:: sum of hillslope areas is zero.'//errmsg(sourcefile, __LINE__) )
           end if
 
@@ -559,7 +559,7 @@ end subroutine InitHillslope
 
   !-----------------------------------------------------------------------
 
-  subroutine SetHillslopeSoilThickness(bounds,fsurdat,soil_depth_lowland_in,soil_depth_upland_in)
+  subroutine SetHillslopeSoilThickness(bounds, hillslope_file, soil_depth_lowland_in, soil_depth_upland_in)
     !
     ! !DESCRIPTION:
     ! Set hillslope column nbedrock values
@@ -578,7 +578,7 @@ subroutine SetHillslopeSoilThickness(bounds,fsurdat,soil_depth_lowland_in,soil_d
     !
     ! !ARGUMENTS:
     type(bounds_type), intent(in) :: bounds
-    character(len=*) , intent(in) :: fsurdat    ! surface data file name
+    character(len=*) , intent(in) :: hillslope_file    ! hillslope data file name
     real(r8), intent(in), optional  :: soil_depth_lowland_in
     real(r8), intent(in), optional  :: soil_depth_upland_in
     real(r8), pointer     :: fhillslope_in(:,:) ! read in - float
@@ -599,14 +599,14 @@ subroutine SetHillslopeSoilThickness(bounds,fsurdat,soil_depth_lowland_in,soil_d
 
     if (soil_profile_method==soil_profile_from_file) then
 
-       ! Open surface dataset to read in data below
-       call getfil (fsurdat, locfn, 0)
+       ! Open hillslope dataset to read in data below
+       call getfil (hillslope_file, locfn, 0)
        call ncd_pio_openfile (ncid, locfn, 0)
 
        allocate(fhillslope_in(bounds%begg:bounds%endg,max_columns_hillslope))
        call ncd_io(ncid=ncid, varname='hillslope_bedrock_depth', flag='read', data=fhillslope_in, dim1name=grlnd, readvar=readvar)
        if (masterproc .and. .not. readvar) then
-          call endrun( 'ERROR:: soil_profile_method = "FromFile", but hillslope_bedrock not found on surface data set.'//errmsg(sourcefile, __LINE__) )
+          call endrun( 'ERROR:: soil_profile_method = "FromFile", but hillslope_bedrock not found on hillslope data set.'//errmsg(sourcefile, __LINE__) )
        end if
        do l = bounds%begl,bounds%endl
           g = lun%gridcell(l)
@@ -896,7 +896,7 @@ end subroutine HillslopeDominantLowlandPft
   subroutine HillslopePftFromFile(bounds,col_pftndx)
     !
     ! !DESCRIPTION:
-    ! Reassign patch type using indices from surface data file
+    ! Reassign patch type using indices from data file
     ! Assumes one patch per hillslope column
     ! In preparation for this reassignment of patch type, only the
     ! first patch was given a non-zero weight in surfrd_hillslope.

diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90
@@ -136,8 +136,8 @@ subroutine initialize2(ni,nj)
     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 : fsurdat, hillslope_file
+    use clm_varctl                    , only : finidat, finidat_interp_source, finidat_interp_dest
     use clm_varctl                    , only : use_cn, use_fates, use_fates_luh
     use clm_varctl                    , only : use_crop, ndep_from_cpl, fates_spitfire_mode
     use clm_varctl                    , only : use_hillslope
@@ -254,7 +254,7 @@ subroutine initialize2(ni,nj)
     call pftcon%Init()
 
     ! Read surface dataset and set up subgrid weight arrays
-    call surfrd_get_data(begg, endg, ldomain, fsurdat, actual_numcft)
+    call surfrd_get_data(begg, endg, ldomain, fsurdat, hillslope_file, actual_numcft)
 
     if(use_fates) then
 
@@ -305,7 +305,7 @@ subroutine initialize2(ni,nj)
 
     if (use_hillslope) then
        ! Initialize hillslope properties
-       call InitHillslope(bounds_proc, fsurdat)
+       call InitHillslope(bounds_proc, hillslope_file)
     endif
 
     ! Set filters

diff --git a/src/main/clm_instMod.F90 b/src/main/clm_instMod.F90
@@ -188,7 +188,7 @@ subroutine clm_instInit(bounds)
     !
     ! !USES: 
     use clm_varpar                         , only : nlevsno
-    use controlMod                         , only : nlfilename, fsurdat
+    use controlMod                         , only : nlfilename, fsurdat, hillslope_file
     use domainMod                          , only : ldomain
     use SoilBiogeochemDecompCascadeMIMICSMod, only : init_decompcascade_mimics
     use SoilBiogeochemDecompCascadeBGCMod  , only : init_decompcascade_bgc
@@ -280,7 +280,7 @@ subroutine clm_instInit(bounds)
 
     ! Set hillslope column bedrock values
     if (use_hillslope) then
-       call SetHillslopeSoilThickness(bounds,fsurdat, &
+       call SetHillslopeSoilThickness(bounds, hillslope_file, &
             soil_depth_lowland_in=8.5_r8,&
             soil_depth_upland_in =2.0_r8)
        call setSoilLayerClass(bounds)

diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90
@@ -109,6 +109,7 @@ module clm_varctl
 
   character(len=fname_len), public :: finidat    = ' '        ! initial conditions file name
   character(len=fname_len), public :: fsurdat    = ' '        ! surface data file name
+  character(len=fname_len), public :: hillslope_file = ' '    ! hillslope data file name
   character(len=fname_len), public :: paramfile  = ' '        ! ASCII data file with PFT physiological constants
   character(len=fname_len), public :: nrevsn     = ' '        ! restart data file name for branch run
   character(len=fname_len), public :: fsnowoptics  = ' '      ! snow optical properties file name