Merge pull request #84 from noaa-oar-arl/develop

Merging Develop into Stable

README.md

@@ -19,12 +19,23 @@ Canopy-App requires NetCDF-Fortran Libraries (i.e., `-lnetcdf -lnetcdff`) when u
 See [the included Makefile](./src/Makefile), which detects NetCDF using `nf-config`, for an example (on GMU Hopper, you can use the `netcdf-c/4.7.4-vh` and `netcdf-fortran/4.5.3-ff` modules).
 
 Compilation options can be controlled with environment variables:
- - `DEBUG=0` (off; default) or `DEBUG=1` (on)
- - `NETCDF=0` (off) or `NETCDF=1` (on; default)
+
+- `FC=gfortran` (default) or compiler name/path (e.g. `FC=ifort`, `FC=gfortran-11`, `FC=/usr/bin/gfortran-11`)
+- `DEBUG=0` (off; default) or `DEBUG=1` (on)
+- `NC=0` (off) or `NC=1` (on; default)
 
 Example:
+a) with compiler set by `FC` environment variable (falling back to `gfortran` if unset), Debug flags ON and with NetCDF:
+```
+DEBUG=1 NC=1 make -C src
+```
+Note: Not supplying `FC` doesn't necessarily give `gfortran`, since `FC` might already be set in the environment (for example, `module load` situations may do this). In such case do:
+```
+DEBUG=1 NC=1 FC=gfortran make -C src
 ```
-DEBUG=1 NETCDF=1 make -C src
+b) with Intel Fortran (`ifort`), Debug flags ON and with NetCDF:
+```
+DEBUG=1 NC=1 FC=ifort make -C src
 ```
 
 ### Modify settings
@@ -38,7 +49,7 @@ which is read at runtime.
 ./canopy
 ```
 
-You can also [run with Python](./python/README.md).
+You can also [generate global inputs and run with Python](./python/README.md).
 
 ## Components
 
@@ -146,7 +157,7 @@ The Canopy-App input data in [Table 2](#table-2-canopy-app-required-input-variab
 | `frp`                            | Total Fire Radiative Power (MW/grid cell area) | [NOAA/NESDIS GBBEPx](https://www.ospo.noaa.gov/Products/land/gbbepx/) |
 | `csz`                            | Cosine of the solar zenith angle (dimensionless) | [Based on Python Pysolar](https://pysolar.readthedocs.io/en/latest/) |
 | `mol`                            | Monin-Obukhov Length (m)                    | Externally calculated using GFS `tmp2m`, `fricv`, and `shtfl`.  ([Essa, 1999](https://inis.iaea.org/collection/NCLCollectionStore/_Public/37/118/37118528.pdf)) |
-| `href`                           | Reference height above canopy (m) - 10 m    | Added as a constant reference height above surface (i.e., 10 m).  Can be taken from NL. |
+| `href`                           | Reference height above canopy (m) - 10 m    | Assumed constant (i.e., 10 m).  Can be taken from NL. |
 
 **More Information on Data Sources from [Table 2](#table-2-canopy-app-required-input-variables):**
 
@@ -168,6 +179,9 @@ Hourly gridded GFSv16 data is available on AWS from March 23, 2021 - Current Day
 https://nacc-in-the-cloud.s3.amazonaws.com/inputs/geo-files/gfs.canopy.t12z.2022MM01.sfcf000.nc
 ```
 
+You can also [generate global inputs using Python (see python/global_data_process.py)](./python/README.md).
+
+
 ### Table 3. Current User Namelist Options
 
 | Namelist Option | Namelist Description and Units                                                     |
@@ -182,7 +196,7 @@ https://nacc-in-the-cloud.s3.amazonaws.com/inputs/geo-files/gfs.canopy.t12z.2022
 | `ifcaneddy`     | logical canopy eddy Kz option (default: `.FALSE.`)                                 |
 | `ifcanphot`     | logical canopy photolysis option (default: `.FALSE.`)                              |
 | `ifcanbio`      | logical canopy biogenic emissions option (default: `.FALSE.`)                      |
-| `href_opt`      | integer for using href_set in namelist (= `0`, default) or array from file (= `1`) |
+| `href_opt`      | integer for using `href_set` in namelist (= `0`, default) or array from file (= `1`) |
 | `href_set`      | user-set real value of reference height above canopy associated with input wind speed (m) (only used if `href_opt=0`) **\*\*\*** |
 | `z0ghc`         | ratio of ground roughness length to canopy top height (Massman et al., 2017)       |
 | `rsl_opt`       | user-set option for either MOST or unified Roughness SubLayer (RSL) effects above and at canopy top (Uc).(= `0`, default: uses MOST and a constant lambdars factor only), (= `1`, under development: will use a more unified RSL approach from Bonan et al. (2018) and Abolafia-Rosenzweig et al., 2021)   |
@@ -197,6 +211,13 @@ https://nacc-in-the-cloud.s3.amazonaws.com/inputs/geo-files/gfs.canopy.t12z.2022
 | `lu_opt`        | integer for input model land use type (`0`: VIIRS 17 Cat (default) or `1`: MODIS-IGBP 20 Cat (valid LU types 1-10 and 12); input mapped to Massman et al.) |
 | `z0_opt`        | integer (`0`: use model input or `1`: vegtype dependent z0 for first estimate)     |
 | `bio_cce`       | user-set real value of MEGAN biogenic emissions "canopy environment coefficient" used to tune emissions to model inputs/calculations (default: `0.21`, based on Silva et al. 2020) |
+| `biovert_opt`   | user set biogenic vertical summing option (`0`: no sum, full leaf-level biogenic emissions, units=kg/m3/s; `1`: MEGANv3-like summing of LAD weighted activity coefficients using the canopy-app plant distributions, caution-- units=kg m-2 s-1 and puts the total emissions in the topmost canopy-app model layer only; `2`: Same as in option 1, but instead uses Gaussian/normally weighted activity coefficients acoss all sub-canopy layers -- also units of kg m-2 s-1 in topmost model layer; `3`: Same as in option 1, but instead uses evenly weighted activity coefficients acoss all sub-canopy layers -- also units of kg m-2 s-1 in topmost model layer          |
+| `ssg_opt`       | integer for using either input data  (= `0`, default) or user set shrub/savanna/grass (SSG) vegetation type heights from namelist (= `1`).  Currently, GEDI FCH input data only provides canopy heights for forests and not SSG.  Warning: use of ssg_opt=1 will overide typically higher resolution input data (e.g., GEDI) forest canopy heights where the lower resolution vegtype data indicates SSG  |
+| `ssg_set`       | user-set real value of constant SSG vegetation type heights (m) (only used if `ssg_opt=1`) |
+| `crop_opt`      | integer for using either input data  (= `0`, default) or user set crop vegetation type heights from namelist (= `1`).  Currently, GEDI FCH input data only provides canopy heights for forests and not crops.  Warning: use of crop_opt=1 will overide typically higher resolution input data (e.g., GEDI) forest canopy heights where the lower resolution vegtype data indicates crops  |
+| `crop_set`      | user-set real value of constant crop vegetation type heights (m) (only used if `crop_opt=1`) |
+| `co2_opt`       | user-set options for applying a CO2 inhibition factor for biogenic isoprene-only emissions using either the [Possell & Hewitt (2011)](https://doi.org/10.1111/j.1365-2486.2010.02306.x) (= `0`, default) or [Wilkinson et al. (2009)](https://doi.org/10.1111/j.1365-2486.2008.01803.x) method (= `1`). Use of option = `1` (Possell & Hewitt 2011) is especially recommended for sub-ambient CO2 concentrations.  To turn off co2 inhibition factor set `co2_opt=2`  |
+| `co2_set`       | user-set real value of atmospheric co2 concentration (ppmv) (only used if `co2_opt=0` or `co2_opt=1`) |
 | `lai_thresh`    | user-set real value of LAI threshold for contiguous canopy (m2/m2)                 |
 | `frt_thresh`    | user-set real value of forest fraction threshold for contiguous canopy             |
 | `fch_thresh`    | user-set real value of canopy height  threshold for contiguous canopy (m)          |

input/namelist.canopy

@@ -5,32 +5,39 @@
 /
 
 &USERDEFS
-  infmt_opt  =  0
-  nlat       =  43
-  nlon       =  86
-  modlays    =  100
-  modres     =  0.5
-  ifcanwind  = .TRUE.
-  ifcanwaf   = .TRUE.
-  ifcaneddy  = .TRUE.
-  ifcanphot  = .TRUE.
-  ifcanbio   = .TRUE.
-  href_opt   =  0
-  href_set   =  10.0
-  z0ghc      =  0.0025
-  rsl_opt    =  0
-  lambdars   =  1.25
-  dx_opt     =  0
-  dx_set     =  12000.0
-  flameh_opt =  2
-  flameh_set =  1.0
-  frp_fac    =  1.0
-  pai_opt    =  0
-  pai_set    =  4.0
-  lu_opt     =  1
-  z0_opt     =  0
-  bio_cce    =  0.21
-  lai_thresh =  0.1
-  frt_thresh =  0.1
-  fch_thresh =  0.5
+  infmt_opt   =  0
+  nlat        =  43
+  nlon        =  86
+  modlays     =  100
+  modres      =  0.5
+  ifcanwind   = .TRUE.
+  ifcanwaf    = .TRUE.
+  ifcaneddy   = .TRUE.
+  ifcanphot   = .TRUE.
+  ifcanbio    = .TRUE.
+  href_opt    =  0
+  href_set    =  10.0
+  z0ghc       =  0.0025
+  rsl_opt     =  0
+  lambdars    =  1.25
+  dx_opt      =  0
+  dx_set      =  12000.0
+  flameh_opt  =  2
+  flameh_set  =  1.0
+  frp_fac     =  1.0
+  pai_opt     =  0
+  pai_set     =  4.0
+  lu_opt      =  1
+  z0_opt      =  0
+  bio_cce     =  0.21
+  biovert_opt =  0
+  ssg_opt     =  0
+  ssg_set     =  1.0
+  crop_opt    =  0
+  crop_set    =  3.0
+  co2_opt     =  0
+  co2_set     =  400.0
+  lai_thresh  =  0.1
+  frt_thresh  =  0.1
+  fch_thresh  =  0.5
 /

python/README.md

@@ -51,3 +51,35 @@ cases = config_cases(
 ds = run_config_sens(cases)
 ```
 :point_up: We still get a single output dataset, but it has a `case` dimension.
+
+### Generating global data
+
+You can also download and generate global gridded canopy-app inputs using Python.
+
+1. Edit python script (`global_data_process.py`)
+
+2. Change user settings
+
+   ```python
+   '''User Options'''
+   path     = '/scratch/pcampbe8/canopy-app/input'  # work directory
+   year     = 2022   # year
+   month    = 7      # month
+   day      = 1      # day
+   houri    = 12     # gfs initialization hour in UTC (caution currently GFS files are initialized at 12 UTC only -- do not change)
+   hour     = 00     # gfs forecast hour in UTC
+   ref_lev  = 10     # reference height above the canopy (m)
+   frp_src  = 1      # frp data source (0: local source; 1: check local source first, switch to climatological file if no available data; 2: 12 month climatology; 3: all ones when ifcanwaf=.FALSE.)
+   ```
+
+2. Activate canopy-app Conda environment:
+
+   ```
+   conda activate canopy-app
+   ```
+
+3. Run Python script
+
+   ```
+   python global_data_process.py
+   ```

python/environment.yml

@@ -21,3 +21,5 @@ dependencies:
   #
   # Extra
   - ipython
+  - pysolar
+  - scipy

python/examples.ipynb

@@ -760,6 +760,64 @@
     "    ax.plot(clu, ds_c[vn])\n",
     "    ax.set(ylabel=vn, xlabel=\"CLU\")"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "13cc57aa-666c-4f27-ac01-400838ba94a0",
+   "metadata": {},
+   "source": [
+    "### Emissions calculation options"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "578f31ab-832e-4303-88c4-9cb444fe99f9",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "cases = config_cases(\n",
+    "    file_vars=\"../input/input_variables_point.txt\",\n",
+    "    infmt_opt=1,\n",
+    "    nlat=1,\n",
+    "    nlon=1,\n",
+    "    biovert_opt=[0, 1, 2, 3],\n",
+    ")\n",
+    "ds = run_config_sens(cases)\n",
+    "ds"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "affd849c-3464-4f27-abc3-1390a4395e31",
+   "metadata": {},
+   "source": [
+    "Currently, method 1 should give the same result as canopy-integrating the laywerwise method 0 results afterwards.\n",
+    "\n",
+    "The method 2 result is smaller since it places more weight on the lower--middle regions of the canopy, where light levels are lower. Method 3 is somewhere in between these (consider the shape of the LAD profile vs a Gaussian in height)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "86b197ac-9ecb-40b2-a7e5-1354fcc0a2d3",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "dz = ds.z.diff(\"z\")\n",
+    "\n",
+    "e0 = float((ds.emi_isop.isel(case=0) * dz).sum(\"z\"))\n",
+    "e1 = float(ds.emi_isop.isel(case=1, z=-1))\n",
+    "e2 = float(ds.emi_isop.isel(case=2, z=-1))\n",
+    "e3 = float(ds.emi_isop.isel(case=3, z=-1))\n",
+    "print(e0, e1, e2, e3, sep=\"\\n\")\n",
+    "\n",
+    "assert np.isclose(e0, e1)"
+   ]
   }
  ],
  "metadata": {