Landcover from GLC2000 and soils from SOILGRIDS

README.md

@@ -140,19 +140,30 @@ For the creation of the mapping files of CLM inputdata to our grid use `mkmapdat
 Adjust the Slurm directives to your compute time project and partition.
 Below the Slurm directives, modify `GRIDNAME` and `GRIDFILE` to your grid and previously created SCRIP file.
 The script can be used on JURECA and JUWELS but it is advisable to use the large memory partitions for larger domains.
+
 The files are at JSC available below `$CSMDATA/lnd/clm2/mappingdata/grids/`.
 If you don't have access to the CLM mappingdata you have to download it.
-Use:
 
+There are two possible ways to download the grids
+
+1. Use:
 ```
 wget --no-check-certificate -i clm_mappingfiles.txt
 ```
-Then start the creation of the weights with
+
+2. Run `download_grids.sh` after adapting inputs (in particular the path `myraw`).
+   You need Subversion (`svn`) for this.
+```
+./download_grids.sh
+```
+
+When all grids are downloaded, start the creation of the weights with
 ```
 sbatch runscript_mkmapdata.sh
 ```
 
-This will create regridding and netCDF mapping files in the current directory.
+This will create regridding and netCDF mapping files in the current
+directory.
 
 ### Icosahedral grid
 
@@ -221,9 +232,25 @@ In a later step we are anyway going to replace variables in this file that are s
 ./mksurfdata.pl -r usrspec -usr_gname $GRIDNAME -usr_gdate $CDATE -l $CSMDATA -allownofile -y 2005 -hirespft -usr_mapdir="../mkmapdata/" -no-crop -pft_idx 13 -pft_frc 100 -soil_cly 60 -soil_col 10 -soil_fmx 0.5 -soil_snd 40
 ```
 
+
 ## Modification of the surface and domain file
 
 The created surface and domain file have negative longitudes that CLM5 does not accept and inherently has no landmask. To modify the longitudes and to add a landmask use `mod_domain.sh` after inserting the paths to your files.
 
 At least for TSMP2, further modification of the surface file is needed and not yet included in this (tested) workflow.
 The necessary replacement routines can be found in the `dev_replace_tsmp2` branch in the [`mksurfdata/`](https://github.com/HPSCTerrSys/eCLM_static-file-generator/tree/dev_replace_tsmp2/mksurfdata) directory.
+
+
+## Surface File: Use landcover GLC2000 and soil texture from SOILGRIDS
+
+Only for BGC mode!
+
+The following script updates landcover using GLC2000 and soils using SOILGRIDS.
+
+```
+./replace_surfdata.py
+```
+
+Additionally, this script checks CLM gridcells to make sure the
+percentages per landunit in each gridcell sum up to one.
+

mkmapdata/download_grids.sh

@@ -0,0 +1,42 @@
+#!/usr/bin/env bash
+set -e
+
+### URL base for the grids
+SVNBASE="https://svn-ccsm-inputdata.cgd.ucar.edu/trunk/inputdata/lnd/clm2/mappingdata/grids"
+
+### Local directory
+myraw="$PWD/download_grids"
+mkdir -p "$myraw"
+
+### List of files to download
+files=(
+  SCRIPgrid_0.5x0.5_AVHRR_c110228.nc
+  SCRIPgrid_0.25x0.25_MODIS_c170321.nc
+  SCRIPgrid_3minx3min_LandScan2004_c120517.nc
+  SCRIPgrid_3minx3min_MODISv2_c190503.nc
+  SCRIPgrid_3minx3min_MODISwcspsea_c151020.nc
+  SCRIPgrid_3x3_USGS_c120912.nc
+  SCRIPgrid_5x5min_nomask_c110530.nc
+  SCRIPgrid_5x5min_IGBP-GSDP_c110228.nc
+  SCRIPgrid_5x5min_ISRIC-WISE_c111114.nc
+  SCRIPgrid_5x5min_ORNL-Soil_c170630.nc
+  SCRIPgrid_10x10min_nomask_c110228.nc
+  SCRIPgrid_10x10min_IGBPmergeICESatGIS_c110818.nc
+  SCRIPgrid_3minx3min_GLOBE-Gardner_c120922.nc
+  SCRIPgrid_3minx3min_GLOBE-Gardner-mergeGIS_c120922.nc
+  SCRIPgrid_0.9x1.25_GRDC_c130307.nc
+  SCRIPgrid_360x720_cruncep_c120830.nc
+  UGRID_1km-merge-10min_HYDRO1K-merge-nomask_c130402.nc
+)
+
+echo "Downloading grid files to $myraw ..."
+for file in "${files[@]}"; do
+  localf="${myraw}/${file}"
+  url="${SVNBASE}/${file}"
+  if [ ! -f "$localf" ]; then
+    echo "Downloading $file ..."
+    svn export "$url" "$localf"
+  else
+    echo "File $file already exists — skipping download."
+  fi
+done

mkmapdata/runscript_mkmapdata.sh

@@ -42,7 +42,20 @@ do
         ;;
       *)
         echo "Creating conservative remapping files for your grid ${GRIDNAME}..."
-        srun $ESMFBIN_PATH/ESMF_RegridWeightGen --ignore_unmapped -s ${rawpath}/${file} -d $(realpath $GRIDFILE) -m conserve -w ${OUTPUT}/${OUTFILE} --dst_regional --netcdf4
+
+	# Conservative remapping: Specify UGRID source type
+	#
+	# In general automatic detection should work:
+	# https://earthsystemmodeling.org/docs/release/latest/ESMF_refdoc/node3.html#SECTION03020000000000000000
+	#
+	# This overrides the automatic detection.
+	if [[ "$basename" == UGRID* ]]; then
+	  SRCTYPE="--src_type UGRID"
+	else
+	  SRCTYPE=""
+	fi
+
+        srun $ESMFBIN_PATH/ESMF_RegridWeightGen --ignore_unmapped -s ${rawpath}/${file} -d $(realpath $GRIDFILE) -m conserve $SRCTYPE -w ${OUTPUT}/${OUTFILE} --dst_regional --netcdf4
         ;;
     esac
 done

mksurfdata/replace_surfdata.py

@@ -0,0 +1,89 @@
+#!/usr/bin/env python3
+import xarray as xr
+import numpy as np
+# import sys
+
+file_clm5 = "/p/project1/cjibg36/jibg3674/eCLM_static-file-generator/surfdata_EUR-11_hist_78pfts_CMIP6_simyr2005_c251022.nc"
+file_glc2000 = "/p/project1/cjibg36/jibg3674/TSMP_EUR-11/static.resource/08_LandCover/GLC2000_PFT_urban.nc"
+file_clm = "/p/project/cjibg36/jibg3674/shared_DA/setup_eclm_cordex_444x432_v9/input_clm/surfdata_EUR-11_hist_16pfts_Irrig_CMIP6_simyr2000_c230808_GLC2000.nc"
+file_dest = "/p/project/cjibg36/jibg3674/shared_DA/setup_eclm_cordex_444x432_v9/input_clm/surfdata_EUR-11_hist_78pfts_CMIP6_simyr2005_c251022_GLC2000.nc"
+
+open_clm5 = xr.open_dataset(file_clm5)
+open_glc2000 = xr.open_dataset(file_glc2000)
+open_clm = xr.open_dataset(file_clm)
+
+# print(open_clm5["PCT_CLAY"].dims, open_clm["PCT_CLAY"].dims)
+# print(open_clm5["PCT_CLAY"].shape, open_clm["PCT_CLAY"].shape)
+
+open_clm5["PCT_CLAY"][:] = open_clm["PCT_CLAY"].values
+open_clm5["PCT_SAND"][:] = open_clm["PCT_SAND"].values
+open_clm5["ORGANIC"][:] = open_clm["ORGANIC"].values
+
+open_clm5["PCT_WETLAND"][:] = 0
+open_clm5["PCT_URBAN"][2, :, :] = open_glc2000["PCT_PFT"][19, :, :].values
+open_clm5["PCT_URBAN"][1, :, :] = 0
+open_clm5["PCT_GLACIER"][:, :] = open_glc2000["PCT_PFT"][18, :, :].values
+open_clm5["PCT_LAKE"][:, :] = open_glc2000["PCT_PFT"][17, :, :].values
+
+open_clm5["PCT_CROP"][:, :] = open_glc2000["PCT_PFT"][15, :, :].values
+open_clm5["PCT_NATVEG"][:, :] = np.sum(
+    open_glc2000["PCT_PFT"][0:15, :, :].values, axis=0
+)
+open_clm5["PCT_NAT_PFT"][0:15, :, :] = open_glc2000["PCT_PFT"][0:15, :, :].values
+
+sum_gridcell = (
+    open_clm5["PCT_NATVEG"].values
+    + open_clm5["PCT_CROP"].values
+    + np.sum(open_clm5["PCT_URBAN"].values, axis=0)
+    + open_clm5["PCT_LAKE"].values
+    + open_clm5["PCT_GLACIER"].values
+)
+print("Min sum_gridcell:", np.nanmin(sum_gridcell))
+print("Max sum_gridcell:", np.nanmax(sum_gridcell))
+
+# Compute scaling factor, with safety margin
+scale = np.where(sum_gridcell > 0, 100 / sum_gridcell, 1.0)
+for v in ["PCT_NATVEG", "PCT_CROP", "PCT_LAKE", "PCT_GLACIER", "PCT_URBAN"]:
+    open_clm5[v].values *= scale
+print(
+    "New sum_gridcell (min, max):",
+    np.nanmin(sum_gridcell * scale),
+    np.nanmax(sum_gridcell * scale),
+)
+
+
+pctspec = (
+    open_clm5["PCT_LAKE"].values
+    + open_clm5["PCT_WETLAND"].values
+    + open_clm5["PCT_GLACIER"].values
+    + np.sum(open_clm5["PCT_URBAN"].values, axis=0)
+)
+print("Max pctspec:", np.nanmax(pctspec))
+
+natpft = open_clm5["PCT_NAT_PFT"].values
+sum_natpft = np.sum(natpft, axis=0)
+sum_crop = np.sum(open_clm5["PCT_CFT"].values, axis=0)
+print("Min sum_natpft:", np.nanmin(sum_natpft))
+print("Max sum_natpft:", np.nanmax(sum_natpft))
+print("Min sum_crop:", np.nanmin(sum_crop))
+print("Max sum_crop:", np.nanmax(sum_crop))
+
+mask = sum_natpft > 0
+natpft[:, mask] = natpft[:, mask] * (100.0 / sum_natpft[mask])
+open_clm5["PCT_NAT_PFT"].values = natpft
+
+mask_zero = sum_natpft == 0
+open_clm5["PCT_NAT_PFT"].values[0, mask_zero] = 100.0
+
+print(
+    "New sum_natpft (min, max):",
+    np.nanmin(np.sum(open_clm5["PCT_NAT_PFT"].values, axis=0)),
+    np.nanmax(np.sum(open_clm5["PCT_NAT_PFT"].values, axis=0)),
+)
+
+open_clm5.to_netcdf(file_dest)
+print(f"Saved modified surfdata to:\n{file_dest}")
+
+open_clm5.close()
+open_glc2000.close()
+open_clm.close()