documentation and style for replace_surfdata.py script

Author: jjokella

README.md

@@ -232,15 +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.
 
-Additional surface file processing is implemented in:
+
+## 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  # For BGC mode, 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.
+

mksurfdata/replace_surfdata.py

@@ -1,6 +1,6 @@
 import xarray as xr
 import numpy as np
-import sys
+# 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"
@@ -14,47 +14,54 @@
 # 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_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_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_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_patch = (
+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_patch:", np.nanmin(sum_patch))
-print("Max sum_patch:", np.nanmax(sum_patch))
+    + 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_patch > 0, 100 / sum_patch, 1.0)
+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_patch (min, max):", np.nanmin(sum_patch * scale), np.nanmax(sum_patch * 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)
+    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  
+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) 
+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))
@@ -67,14 +74,15 @@
 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)))
+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()
-