Merge branch 'no_wetlands_ctsm52' into ctsm5.2.mksurfdata

Add option to convert wetlands to land; apply it by default for CLM51
physics

From #1890

This PR adds a new namelist option, "convert_ocean_to_land", which
converts all wetland area to natural veg (typically bare ground) at
runtime. This is set to true by default when using CLM51 physics.

In addition, this PR makes a minor change to the setting of soil
properties in mksurfdata: Previously, soil properties were set to some
default value outside of the pctlnd_pft-determined land mask. But I
can't see any reason why that needs to be done, and doing so could
remove some valid soil properties in the situation where there is a grid
cell that has valid soil properties even though the pctpft dataset
claims the area is ocean.

Note that this PR is into the ctsm5.2.mksurfdata branch. Most of the
changes could be applied directly to master, but there are some related
changes to mksurfdata that motivated me to make this PR into the
ctsm5.2.mksurfdata branch. However, standard f09 surface datasets don't
show any changes due to the code changes in mksurfdata here, and the
mksurfdata changes aren't critical anyway, so it would be reasonable to
rebase the non-mksurfdata changes onto master if we want them available
sooner.

Resolves #1878

Are answers expected to change (and if so in what way)? Yes: changes
answers for CLM51 cases in a few grid cells that used to be wetland but
now are classified as bare ground.

Any User Interface Changes (namelist or namelist defaults changes)? Adds
new namelist variable, "convert_ocean_to_land"

Testing performed, if any:
- Ran the following tests with comparison against the ctsm5.2mksurfdata branch:
  - `ERP_D_P36x2_Ld3.f10_f10_mg37.I1850Clm50BgcCrop.cheyenne_intel.clm-default`: passes and bit-for-bit
  - `ERP_Ld9.f45_g37.I2000Clm51Bgc.cheyenne_intel.clm-default`: passes but changes answers, as expected
  - `ERI_D_Ld9.f10_f10_mg37.I1850Clm51Bgc.cheyenne_gnu.clm-default`: passes (didn't do baseline comparisons, but expect baseline comparisons to fail)
- Also made new surface datasets with these changes in soil property mapping, 1850 f09:
  - Standard surface datasets show no differences
  - Also did a test where I temporarily set pctlnd_pft to 0 everywhere
  after the call to mkpft. I verified that this only leads to
  differences in various PCT fields: there are no longer any differences
  in soil properties arising from this change in pctlnd_pft. In
  contrast, setting pctlnd_pft to 0 before the changes in this PR leads
  to differences in many soil properties, because in this baseline
  version, soil properties were set to a constant value over grid cells
  outside the pctlnd_pft-based land mask.

bld/CLMBuildNamelist.pm

@@ -2034,6 +2034,7 @@ sub setup_logic_subgrid {
 
    my $var = 'run_zero_weight_urban';
    add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var);
+   add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'convert_ocean_to_land');
    add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'collapse_urban',
                'structure'=>$nl_flags->{'structure'});
    add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'n_dom_landunits',

bld/namelist_files/namelist_defaults_ctsm.xml

@@ -1754,6 +1754,9 @@ lnd/clm2/surfdata_map/release-clm5.0.30/surfdata_ne0np4.CONUS.ne30x8_hist_78pfts
 
 <run_zero_weight_urban>.false.</run_zero_weight_urban>
 
+<convert_ocean_to_land>.false.</convert_ocean_to_land>
+<convert_ocean_to_land phys="clm5_1">.true.</convert_ocean_to_land>
+
 <collapse_urban structure="standard" >.false.</collapse_urban>
 <collapse_urban structure="fast"     >.true.</collapse_urban>
 

bld/namelist_files/namelist_definition_ctsm.xml

@@ -1161,6 +1161,14 @@ Toggle for vancouver specific logic.
 Toggle for mexico city specific logic.
 </entry>
 
+<entry id="convert_ocean_to_land" type="logical" category="physics"
+       group="clm_inparm" value=".false.">
+If true, any ocean (i.e., wetland) points on the surface dataset are
+converted to bare ground (or whatever vegetation is given in that grid
+cell - but typically this will be bare ground due to lack of vegetation
+in grid cells with 100% ocean).
+</entry>
+
 <entry id="n_dom_pfts" type="integer" category="physics"
        group="clm_inparm"
        valid_values="0,1,2,3,4,5,6,7,8,9,10,11,12,13,14" value="0">

src/main/clm_varctl.F90

@@ -168,6 +168,11 @@ module clm_varctl
   ! true => make ALL patches, cols & landunits active (even if weight is 0)
   logical, public :: all_active = .false.          
 
+  ! true => any ocean (i.e., "wetland") points on the surface dataset are converted to
+  ! bare ground (or whatever vegetation is given in that grid cell... but typically this
+  ! will be bare ground)
+  logical, public :: convert_ocean_to_land = .false.
+
   logical, public :: collapse_urban = .false.  ! true => collapse urban landunits to the dominant urban landunit; default = .false. means "do nothing" i.e. keep all urban landunits as found in the input data
   integer, public :: n_dom_landunits = -1  ! # of dominant landunits; determines the number of active landunits; default = 0 (set in namelist_defaults_ctsm.xml) means "do nothing"
   integer, public :: n_dom_pfts = -1  ! # of dominant pfts; determines the number of active pfts; default = 0 (set in namelist_defaults_ctsm.xml) means "do nothing"

src/main/controlMod.F90

@@ -259,6 +259,8 @@ subroutine control_init(dtime)
 
     ! All old cpp-ifdefs are below and have been converted to namelist variables
 
+    namelist /clm_inparm/ convert_ocean_to_land
+
     ! Number of dominant pfts and landunits. Enhance ctsm performance by
     ! reducing the number of active pfts to n_dom_pfts and
     ! active landunits to n_dom_landunits.
@@ -668,6 +670,7 @@ subroutine control_spmd()
     ! Other subgrid logic
     call mpi_bcast(run_zero_weight_urban, 1, MPI_LOGICAL, 0, mpicom, ier)
     call mpi_bcast(all_active, 1, MPI_LOGICAL, 0, mpicom, ier)
+    call mpi_bcast(convert_ocean_to_land, 1, MPI_LOGICAL, 0, mpicom, ier)
 
     ! Number of dominant pfts and landunits. Enhance ctsm performance by
     ! reducing the number of active pfts to n_dom_pfts and
@@ -890,6 +893,7 @@ subroutine control_print ()
     else
        write(iulog,*) '   land frac data = ',trim(fatmlndfrc)
     end if
+    write(iulog,*) '   Convert ocean to land = ', convert_ocean_to_land
     write(iulog,*) '   Number of ACTIVE PFTS (0 means input pft data NOT collapsed to n_dom_pfts) =', n_dom_pfts
     write(iulog,*) '   Number of ACTIVE LANDUNITS (0 means input landunit data NOT collapsed to n_dom_landunits) =', n_dom_landunits
     write(iulog,*) '   Collapse urban landunits; done before collapsing all landunits to n_dom_landunits; .false. means do nothing i.e. keep all the urban landunits, though n_dom_landunits may still remove them =', collapse_urban

src/main/surfrdMod.F90

@@ -15,7 +15,7 @@ module surfrdMod
   use clm_varcon      , only : grlnd
   use clm_varctl      , only : iulog
   use clm_varctl      , only : use_cndv, use_crop, use_fates
-  use surfrdUtilsMod  , only : check_sums_equal_1, collapse_crop_types
+  use surfrdUtilsMod  , only : check_sums_equal_1, apply_convert_ocean_to_land, 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
   use ncdio_pio       , only : ncd_io, check_var, ncd_inqfdims, check_dim_size, ncd_inqdid, ncd_inqdlen
@@ -70,7 +70,7 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft)
     !    o real % abundance PFTs (as a percent of vegetated area)
     !
     ! !USES:
-    use clm_varctl  , only : create_crop_landunit, collapse_urban, &
+    use clm_varctl  , only : create_crop_landunit, convert_ocean_to_land, collapse_urban, &
                              toosmall_soil, toosmall_crop, toosmall_glacier, &
                              toosmall_lake, toosmall_wetland, toosmall_urban, &
                              n_dom_landunits
@@ -193,6 +193,10 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft)
 
     call check_sums_equal_1(wt_lunit, begg, 'wt_lunit', subname)
 
+    if (convert_ocean_to_land) then
+       call apply_convert_ocean_to_land(wt_lunit(begg:endg,:), begg, endg)
+    end if
+
     ! if collapse_urban = .true.
     ! collapse urban landunits to the dominant urban landunit
 

src/main/surfrdUtilsMod.F90

@@ -20,7 +20,8 @@ module surfrdUtilsMod
   ! !PUBLIC MEMBER FUNCTIONS:
   public :: check_sums_equal_1  ! Confirm that sum(arr(n,:)) == 1 for all n
   public :: renormalize         ! Renormalize an array
-  public :: convert_cft_to_pft  ! Conversion of crop CFT to natural veg PFT:w
+  public :: apply_convert_ocean_to_land ! Apply the conversion of ocean to land points
+  public :: convert_cft_to_pft  ! Conversion of crop CFT to natural veg PFT
   public :: collapse_crop_types ! Collapse unused crop types into types used in this run
   public :: collapse_individual_lunits  ! Collapse landunits by user-defined thresholds
   public :: collapse_to_dominant ! Collapse to dominant pfts or landunits
@@ -112,6 +113,45 @@ subroutine renormalize(arr, lb, normal)
 
   end subroutine renormalize
 
+  !-----------------------------------------------------------------------
+  subroutine apply_convert_ocean_to_land(wt_lunit, begg, endg)
+    !
+    ! !DESCRIPTION:
+    ! Apply the conversion of ocean points to land, by changing all "wetland" points to
+    ! natveg; typically this will result in these points becoming bare ground.
+    !
+    ! The motivation for doing this is to avoid the negative runoff that sometimes comes
+    ! from wetlands.
+    !
+    ! !USES:
+    use landunit_varcon, only : istsoil, istwet, max_lunit
+    !
+    ! !ARGUMENTS:
+    integer, intent(in) :: begg  ! Beginning grid cell index
+    integer, intent(in) :: endg  ! Ending grid cell index
+    ! This array is modified in-place:
+    real(r8), intent(inout) :: wt_lunit(begg:endg, max_lunit) ! Weights of landunits per grid cell
+    !
+    ! !LOCAL VARIABLES:
+    integer :: g
+
+    character(len=*), parameter :: subname = 'apply_convert_ocean_to_land'
+    !-----------------------------------------------------------------------
+
+    ! BUG(wjs, 2022-10-27, ESCOMP/CTSM#1886) Ideally we would distinguish between ocean
+    ! vs. true wetland points on the surface dataset; for now oceans are included in the
+    ! wetland area on the surface dataset, so we convert all wetlands to land. (Typically
+    ! there are no true/inland wetlands on the surface dataset, so this is currently okay,
+    ! but this would become a problem if we started having inland wetlands on the surface
+    ! dataset again.)
+    do g = begg, endg
+       wt_lunit(g,istsoil) = wt_lunit(g,istsoil) + wt_lunit(g,istwet)
+       wt_lunit(g,istwet) = 0._r8
+    end do
+
+  end subroutine apply_convert_ocean_to_land
+
+
 !-----------------------------------------------------------------------
   subroutine convert_cft_to_pft( begg, endg, cftsize, wt_cft )
     !

tools/mksurfdata_esmf/src/mksoilcolMod.F90

@@ -29,13 +29,12 @@ module mksoilcolMod
 contains
 !=================================================================================
 
-  subroutine mksoilcol(file_data_i, file_mesh_i, mesh_o, pctlnd_pft_o, pioid_o, rc)
+  subroutine mksoilcol(file_data_i, file_mesh_i, mesh_o, pioid_o, rc)
 
     ! input/output variables
     character(len=*)  , intent(in)    :: file_mesh_i     ! input mesh file name
     character(len=*)  , intent(in)    :: file_data_i     ! input data file name
     type(ESMF_Mesh)   , intent(in)    :: mesh_o          ! model mesho
-    real(r8)          , intent(in)    :: pctlnd_pft_o(:)
     type(file_desc_t) , intent(inout) :: pioid_o
     integer           , intent(out)   :: rc
 
@@ -185,13 +184,6 @@ subroutine mksoilcol(file_data_i, file_mesh_i, mesh_o, pctlnd_pft_o, pioid_o, rc
        end if
     end do
 
-    ! assume medium soil color (15) and loamy texture if pct_lndpft is small
-    do no = 1,ns_o
-       if (pctlnd_pft_o(no) < 1.e-6_r8) then
-          soil_color_o(no) = 15
-       end if
-    end do
-
     ! Write output data
     if (root_task)  write(ndiag, '(a)') trim(subname)//" writing out soil color"
     call mkfile_output(pioid_o,  mesh_o, 'SOIL_COLOR', soil_color_o,  rc=rc)
@@ -277,7 +269,7 @@ subroutine get_dominant_soilcol(dynamicMaskList, dynamicSrcMaskValue, dynamicDst
              soil_color_o = maxindex(1)
           end if
 
-          ! If land but no color, set color to 15 (in older dataset generic soil color 4)
+          ! If no color, set color to 15 (in older dataset generic soil color 4)
           if (num_soilcolors == 8) then
              if (soil_color_o == 0) then
                 soil_color_o = 4

tools/mksurfdata_esmf/src/mksoiltexMod.F90

@@ -33,7 +33,7 @@ module mksoiltexMod
 contains
 !=================================================================================
 
-  subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o, pctlnd_pft_o, rc)
+  subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o, rc)
     !
     ! make %sand, %clay, organic carbon content, coarse fragments, bulk density,
     ! and pH measured in H2O
@@ -44,7 +44,6 @@ subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o
     character(len=*)  , intent(in)    :: file_lookup_i   ! input data file name
     type(ESMF_Mesh)   , intent(in)    :: mesh_o          ! output mesh
     type(file_desc_t) , intent(inout) :: pioid_o
-    real(r8)          , intent(in)    :: pctlnd_pft_o(:) ! PFT data: % of gridcell for PFTs
     integer           , intent(out)   :: rc
 
     ! local variables
@@ -291,9 +290,9 @@ subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o
 
     do no = 1,ns_o
 
-       if (pctlnd_pft_o(no) < 1.e-6_r8 .or. mapunit_o(no) == 0) then
+       if (mapunit_o(no) == 0) then
 
-          ! Set sand and clay to loam if pctlnd_pft is < 1.e-6 or mapunit is 0
+          ! Set sand and clay to loam if mapunit is 0
           sand_o(no,:) = 43._r4
           clay_o(no,:) = 18._r4
           orgc_o(no,:) = 0._r4

tools/mksurfdata_esmf/src/mksurfdata.F90

@@ -458,7 +458,7 @@ program mksurfdata
   ! -----------------------------------
   if (fsurdat /= ' ') then
      call mksoiltex( mksrf_fsoitex_mesh, file_mapunit_i=mksrf_fsoitex, file_lookup_i=mksrf_fsoitex_lookup, &
-          mesh_o=mesh_model, pioid_o=pioid, pctlnd_pft_o=pctlnd_pft, rc=rc)
+          mesh_o=mesh_model, pioid_o=pioid, rc=rc)
      if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mksoiltex')
   end if
 
@@ -467,7 +467,7 @@ program mksurfdata
   ! -----------------------------------
   if (fsurdat /= ' ') then
      ! SOIL_COLOR and mxsoil_color is written out in the subroutine
-     call mksoilcol( mksrf_fsoicol, mksrf_fsoicol_mesh, mesh_model, pctlnd_pft, pioid, rc)
+     call mksoilcol( mksrf_fsoicol, mksrf_fsoicol_mesh, mesh_model, pioid, rc)
      if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mksoilcol')
   end if
 
@@ -1080,7 +1080,6 @@ subroutine normalize_and_check_landuse(ns_o)
          pctgla(n) = float(nint(pctgla(n)))
 
          ! Assume wetland, glacier and/or lake when dataset landmask implies ocean
-         ! (assume medium soil color (15) and loamy texture).
 
          if (pctlnd_pft(n) < 1.e-6_r8) then
             if (pctgla(n) < 1.e-6_r8) then