Merge pull request #1196 from GEOS-ESM/feature/borescan_peat_and_hwsd

BCS processing to support HWSDv2 soil preprocessing and new Peat GPA22

Author: sdrabenh

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/Shared/catch_constants.f90

@@ -141,10 +141,11 @@ module catch_constants
   ! Use of peat-specific hydrology (PEATCLSM) is triggered by a porosity threshold.
   ! Porosity of peat tiles depends on bcs version.
   !
-  !     bcs version           | source of peat info       | porosity
-  !     -----------------------------------------------------------------
-  !     NLv3, NLv4            | HWSD                      | poros=0.80
-  !     NLv5                  | PEATMAP                   | poros=0.93
+  !     bcs version           | source of peat info       | porosity   | PEATCLSM 
+  !     --------------------------------------------------------------------------
+  !     NLv3, NLv4            | HWSDv1.21                 | poros=0.80 |   no
+  !     NLv5, v11, v12, v13   | HWSDv1.21 + PEATMAP       | poros=0.93 |   yes
+  !     v14                   | Global Peatland Map 2.0   | poros=0.93 |   yes
   !
   ! - reichle, 26 Jan 2022
 

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/create_README.csh

@@ -550,14 +550,27 @@ APPENDIX I - mkCatchParam input options and log ................................
 	 across a given catchment land element to determine that catchment’s effective soil 
 	 depth.
 
+         PEAT_INFO=0:
+         Peat tiles are determined based solely on the organic carbon content provided 
+         by the HWSDv1.21 texture data. 
+
+         PEAT_INFO=1:
 	 For PEATMAP, sixteen ARC/GIS shapefiles were obtained from Xu et al. (2017). Canada and 
 	 Hokaido-Mongolia-NorthKorea provided peat fractions inside designated polygons. All 
-	 other regional shape files give exact perimeter of of the peatland. A global 30-arcsec
+	 other regional shape files give exact perimeter of the peatland. A global 30-arcsec
 	 raster array (43200x21600) of peatland fraction was constructed using those 16 shapefiles.
-	 Soil types derived on tiles are now further updated using the peatmap data. If the computed
+	 Soil types derived on tiles are now further updated using the PEATMAP data. If the computed
 	 fractional coverage of peatland based on PEATMAP data at a given catchment-tile excceeds 0.5(?), 
 	 we assume the dominant soil type is as peatland. Soil hydraulic parameters for all
-	 peatland tiles are now obtained from Bechtold et al. (2019).
+	 peatland tiles are now obtained from Bechtold et al. (2019).  Bugs in the processing logic
+         result in inconsistent parameter values for some tiles; the number of tiles with
+         inconsistent parameters depends on resolution. 
+
+         PEAT_INFO=2:
+         The Global Peatland Map 2.0 (GPM 2.0) from the Greifwald Mire Centre is used to 
+         determine peat/mineral on the 30-arcsec raster grid.  Tiles are considered peat
+         if the majority of the contributing raster grid cells is peat based on GPM 2.0, 
+         regardless of the HWSDv1.21 organic carbon content values. 
          
    2.2 Data files and images
 

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py

@@ -184,26 +184,33 @@ def ask_questions(default_grid="Cubed-Sphere"):
     fixes that were implemented since the archived BCs in the above-mentioned \n \
     directories were originally created.  The impact of these differences on \n \
     science is insignificant, and the parameter files produced by current \n \
-    code are scientifically equivalent to the corresponding archived BCs. \n",
+    code are scientifically equivalent to corresponding archived BCs. \n \
+    Current code always applies \"mean land elevation fix\" of v13.  For v12 and earlier, land elevation produced \n \
+    by current code thus differs from that in bcs_shared project disk. \n",
             "choices": [ \
    "NL3 : Icarus-NLv3   (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Icarus-NLv3/)", \
    "NL4 : NLv4 [SMAPL4] (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Icarus-NLv4/) \n\
           = NL3 + JPL veg height", \
    "NL5 : NLv5 [SMAPL4] (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Icarus-NLv5/)\n \
          = NL3 + JPL veg height + PEATMAP", \
-   "v06 : NL3 + JPL veg height + PEATMAP + MODIS snow alb", \
-   "v07 : NL3 + PEATMAP", \
-   "v08 : NL3 + MODIS snow alb", \
-   "v09 : NL3 + PEATMAP + MODIS snow alb", \
-   "v10 : NL3 + PEATMAP + MODIS snow alb v2", \
    "v11 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2", \
    "v12 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2 + Argentina peatland fix", \
    "v13 : As in v12 + mean land elevation fix + MOM6 v2 (OM4) ocean-seaice bathymetry", \
+   "v14 : As in v13 but with GPM 2.0 peat", \
    "ICA : Icarus        (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Icarus/)", \
    "GM4 : Ganymed-4_0   (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Ganymed-4_0/)", \
    "F25 : Fortuna-2_5   (archived*: n/a)"], 
         },
 
+# The following bcs versions were never used for major model versions or products and thus are not in
+# the bcs_shared project directory on Discover:
+#
+#  "v06 : NL3 + JPL veg height + PEATMAP + MODIS snow alb", \
+#  "v07 : NL3 + PEATMAP", \
+#  "v08 : NL3 + MODIS snow alb", \
+#  "v09 : NL3 + PEATMAP + MODIS snow alb", \
+#  "v10 : NL3 + PEATMAP + MODIS snow alb v2", \
+      
        {
             "type": "select",
             "name": "grid_type",

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_shared.py

@@ -26,13 +26,36 @@
     "v11": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
     "v12": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"},
     "v13": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v2"},
+    "v14": {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v2"},
 }
 
-_DEFAULTS = {"TOPO_VERSION": "v1", "MOM6_BATHY_VERSION": "v1"}
+_REQUIRED_VERSION_KEYS = ("TOPO_VERSION", "MOM6_BATHY_VERSION")
 
 def resolve_bcs_matrix(bcs_version: str):
     key = (bcs_version or "").strip()
-    return {**_DEFAULTS, **_VERSION_MATRIX.get(key, {})}
+
+    if not key:
+        valid = ", ".join(sorted(_VERSION_MATRIX))
+        raise ValueError(
+            f"BCS version is not defined. Must be one of: {valid}"
+        )
+
+    if key not in _VERSION_MATRIX:
+        valid = ", ".join(sorted(_VERSION_MATRIX))
+        raise ValueError(
+            f"Unknown BCS version '{key}'. Must be one of: {valid}"
+        )
+
+    entry = _VERSION_MATRIX[key]
+
+    missing = [name for name in _REQUIRED_VERSION_KEYS if name not in entry]
+    if missing:
+        raise ValueError(
+            f"BCS version '{key}' is missing required version setting(s): "
+            f"{', '.join(missing)}"
+        )
+
+    return entry
 
 def topo_version_for_bcs(bcs_version: str) -> str:
     return resolve_bcs_matrix(bcs_version)["TOPO_VERSION"]
@@ -198,7 +221,7 @@ def get_script_mv(grid_type):
 
 /bin/rm -r {TMP_DIR}
 
-# if necessary, copy resolution-independent CO2 file from MAKE_BCS_INPUT_DIR to bcs dir
+# if necessary, copy resolution-independent CO2 and route files from MAKE_BCS_INPUT_DIR to bcs dir
 
 if(-f land/shared/CO2_MonthlyMean_DiurnalCycle.nc4) then
     echo "CO2_MonthlyMean_DiurnalCycle.nc4 already present in bcs dir."
@@ -207,10 +230,10 @@ def get_script_mv(grid_type):
     echo "Successfully copied CO2_MonthlyMean_DiurnalCycle.nc4 to bcs dir."
 endif
 
-if(-f land/shared/route_parameters.nc ) then
+if(-f route/route_parameters.nc ) then
     echo "route_parameters.nc already present in bcs dir."
 else
-    /bin/cp -p {MAKE_BCS_INPUT_DIR}/route/routing_model/v1/route_parameters.nc land/shared/route_parameters.nc
+    /bin/cp -p {MAKE_BCS_INPUT_DIR}/route/routing_model/v1/route_parameters.nc route/route_parameters.nc
     echo "Successfully copied route_parameters.nc to bcs dir."
 endif
 

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90

@@ -44,7 +44,7 @@ PROGRAM mkCatchParam
   character*128        :: ARG, MaskFile
   character*256        :: CMD
   character*1          :: opt
-  character*7          :: PEATSOURCE   = 'GDLHWSD'
+  character*128        :: PEATSOURCE   = ''
   character*3          :: VEGZSOURCE   = 'D&S'
   character*2          :: DL ='DC'    
   integer              :: II, JJ, Type
@@ -207,8 +207,19 @@ PROGRAM mkCatchParam
       fnameRst='rst/'//trim(Gridname)  
       fnameTil='til/'//trim(Gridname)  
     endif 
-
-    if(use_PEATMAP)  PEATSOURCE   = 'PEATMAP'
+    
+    select case (PEAT_INFO)
+    case (0)
+       PEATSOURCE = 'HWSDv1.21 (formerly "GDLHWSD")'
+    case (1)
+       PEATSOURCE = 'HWSDv1.21 + PEATMAP (Xu et al 2017, doi:10.5518/252)'
+    case (2)
+       PEATSOURCE = 'Global Peatland Map 2.0 (Greifswald)'
+    case default
+       write (log_file,'(a)') 'Unknown PEAT_INFO.  Stopping!'
+       stop
+    end select
+    
     if(jpl_height)   VEGZSOURCE   = 'JPL'
 
     if (trim(SNOWALB)=='MODC061' .or. trim(SNOWALB) =='MODC061v2')  process_snow_albedo=.true.