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

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,13 @@ 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).
+</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