diff --git a/README.md b/README.md index 409f0e3e..92388d32 100644 --- a/README.md +++ b/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) | diff --git a/input/namelist.canopy b/input/namelist.canopy index a07e3417..04781010 100755 --- a/input/namelist.canopy +++ b/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 / diff --git a/python/README.md b/python/README.md index c348d048..c7dfe85d 100644 --- a/python/README.md +++ b/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 + ``` diff --git a/python/environment.yml b/python/environment.yml index c8421258..0a80ee97 100644 --- a/python/environment.yml +++ b/python/environment.yml @@ -21,3 +21,5 @@ dependencies: # # Extra - ipython + - pysolar + - scipy diff --git a/python/examples.ipynb b/python/examples.ipynb index d9a75e91..62310959 100644 --- a/python/examples.ipynb +++ b/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": { diff --git a/python/global_data_process.py b/python/global_data_process.py new file mode 100755 index 00000000..9f64347b --- /dev/null +++ b/python/global_data_process.py @@ -0,0 +1,434 @@ +""" +Created on Sun Jun 4 15:45:29 2023 + +Author: Wei-Ting Hung +""" + +import os +import subprocess +from datetime import datetime, timezone + +import numpy as np +from netCDF4 import Dataset +from pysolar.solar import get_altitude +from scipy.interpolate import griddata + +"""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 input files are initialized at 12 UTC only -- do not change) +hour = 0 # gfs forecast hour (0-24) in UTC +ref_lev = 10 # reference height above the canopy (m) +frp_src = 1 # frp data source for WAF (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.) + + +# ---------------------------------ATTENTION---------------------------------# +# If local FRP is used (frp_src=0,1), archived GBBEPx files since 2020 are # +# available for GMU HOPPER users. For users outside GMU, parameter "f_frp" # +# and function "read_frp_local" need to be modified accordingly. # +# Function read_frp_local is designed for reading GBBEPx by default. # +# # +# Recent GBBEPx v3 files (~ < 3 months) are available for download at: # +# https://www.ospo.noaa.gov/Products/land/gbbepx/ # +# GBBEPx_all01GRID.emissions_v003_'+YY+MM+DD+'.nc' # +# # +# Archived GBBEPx files since 2020 are available for GMU HOPPER users. # +# # +# 12 month climatological FRP will not be respective of actual conditions. # +# Only use it when users do not need the limited WAF application for fires. # +# ---------------------------------------------------------------------------# + + +starttime = datetime.now() +print("------------------------------------") +print("---- Global input pre-process start!", starttime.strftime("%Y/%m/%d %H:%M:%S")) +print("------------------------------------") + + +"""Settings""" +# date/time +YY = str(year) +MM = "%02d" % month +DD = "%02d" % day +HHI = "%02d" % houri +HH = "%02d" % hour + + +# domain +lat_lim = [-90, 90] +lon_lim = [0, 360] + + +# input/output files +f_met = ( + path + "/gfs.t" + HHI + "z." + YY + MM + DD + ".sfcf0" + HH + ".nc" +) # gfs met file +f_can = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + "01.sfcf000.nc" # canopy file +f_output = ( + path + "/gfs.canopy.t" + HHI + "z." + YY + MM + DD + ".global.f0" + HH + ".nc" +) # output file + +if (frp_src == 0) or (frp_src == 1): # local frp file + if int(YY + MM + DD) <= 20230510: # version 3 + f_frp = ( + "/groups/ESS/yli74/data/GBBEPx/ORI/GBBEPx_all01GRID.emissions_v003_" + + YY + + MM + + DD + + ".nc" + ) + else: # version 4 + f_frp = ( + "/groups/ESS/yli74/data/GBBEPx/ORI/GBBEPx-all01GRID_v4r0_blend_" + + YY + + MM + + DD + + ".nc " + ) +elif frp_src == 2: # climatological frp + f_frp = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + "01.sfcf000.nc" + + +# required variables +loclist = ["grid_xt", "lon", "grid_yt", "lat", "time"] +metlist = [ + "ugrd10m", + "vgrd10m", + "fricv", + "sfcr", + "vtype", + "sotyp", + "pressfc", + "dswrf", + "shtfl", + "tmpsfc", + "tmp2m", + "spfh2m", + "hpbl", + "prate_ave", +] +canlist = ["lai", "clu", "ffrac", "fh", "mol", "csz", "frp", "href"] + + +# constants +fill_value = 9.99e20 # fill value +den = 1.18 # air density (kg/m3) +Cp = 1004 # specific heat capacity of Air at 25C (J/kg/K) +K = 0.4 # Von Karman constant +g = 9.8 # gravitational acceleration (m/s2) + + +# functions +def read_varatt(var): + attname = var.ncattrs() + att = [var.getncattr(X) for X in attname] + return attname, att + + +def write_varatt(var, attname, att): + for X in np.arange(len(attname)): + if attname[X] == "_FillValue": + continue + elif attname[X] == "missing_value": + value = np.float32(att[X]) + value = np.round(value / 1e15) * 1e15 + var.setncattr(attname[X], value) + else: + var.setncattr(attname[X], att[X]) + + +def mapping(xgrid, ygrid, data, xdata, ydata, map_method, fvalue): + output = griddata( + (xdata, ydata), data, (xgrid, ygrid), method=map_method, fill_value=fvalue + ) + return output + + +def read_gfs_climatology(filename, lat, lon, varname): + readin = Dataset(filename) + + # map to met grids + yt = readin["lat"][:] + xt = readin["lon"][:] + data = np.squeeze(readin[varname][0, :, :]) + + DATA = mapping(lat, lon, data.flatten(), yt.flatten(), xt.flatten(), "linear", np.nan) + DATA[np.isnan(DATA)] = 0 + DATA[DATA < 0] = 0 + return DATA + + +def read_frp_local(filename, lat, lon, fill_value): + readin = Dataset(filename) + + # map to met grids + xt, yt = np.meshgrid(readin["Longitude"][:], readin["Latitude"][:]) + xt[xt < 0] = xt[xt < 0] + 360 + data = np.squeeze(readin["MeanFRP"][:]) + + DATA = mapping( + lat, lon, data.flatten(), yt.flatten(), xt.flatten(), "nearest", fill_value + ) + return DATA + + +"""Data Download""" +"""Download from servers if required files do not exist.""" +print("---- Checking required files...") +print("------------------------------------") + +# met file +if os.path.isfile(f_met) is True: + print("---- Met file found!") +else: + print("---- Cannot find met file. Downloading from AWS...") + subprocess.run( + [ + "wget", + "--no-check-certificate", + "--no-proxy", + "-O", + path + "/gfs.t" + HHI + "z." + YY + MM + DD + ".sfcf0" + HH + ".nc", + "https://nacc-in-the-cloud.s3.amazonaws.com/inputs/" + + YY + + MM + + DD + + "/gfs.t" + + HHI + + "z.sfcf0" + + HH + + ".nc", + ] + ) + if os.path.isfile(f_met) is True: + os.chmod(f_met, 0o0755) + print("---- Download complete!") + else: + print("---- No available met data. Terminated!") + exit() + +# can file +if os.path.isfile(f_can) is True: + print("---- Canopy file found!") +else: + print("---- Cannot find canopy file. Downloading from AWS...") + subprocess.run( + [ + "wget", + "--no-check-certificate", + "--no-proxy", + "-P", + path, + "https://nacc-in-the-cloud.s3.amazonaws.com/inputs/geo-files/gfs.canopy.t" + + HHI + + "z." + + YY + + MM + + "01.sfcf000.nc", + ] + ) + if os.path.isfile(f_can) is True: + os.chmod(f_can, 0o0755) + print("---- Download complete!") + else: + print("---- No available canopy data. Terminated!") + exit() + +# frp file +if frp_src == 0: # local source + if os.path.isfile(f_frp) is True: + os.system("cp " + f_frp + " " + path) + if int(YY + MM + DD) <= 20230510: + f_frp = path + "/GBBEPx_all01GRID.emissions_v003_" + YY + MM + DD + ".nc" + else: + f_frp = path + "/GBBEPx-all01GRID_v4r0_blend_" + YY + MM + DD + ".nc" + os.chmod(f_frp, 0o0755) + print("---- FRP file found!") + else: + print("---- No available FRP file. Terminated!") + exit() + +if frp_src == 1: # local source + if os.path.isfile(f_frp) is True: + os.system("cp " + f_frp + " " + path) + if int(YY + MM + DD) <= 20230510: + f_frp = path + "/GBBEPx_all01GRID.emissions_v003_" + YY + MM + DD + ".nc" + else: + f_frp = path + "/GBBEPx-all01GRID_v4r0_blend_" + YY + MM + DD + ".nc" + os.chmod(f_frp, 0o0755) + print("---- FRP file found!") + else: + print("---- No available FRP file. Switch to Climatology FRP...") + frp_src = 2 + f_frp = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + "01.sfcf000.nc" + +if frp_src == 2: # 12 month climatology frp + if os.path.isfile(f_frp) is True: + print("---- FRP file found!") + else: + print("---- Canot find FRP file. Downloading from AWS...") + subprocess.run( + [ + "wget", + "--no-check-certificate", + "--no-proxy", + "-P", + path, + "https://nacc-in-the-cloud.s3.amazonaws.com/inputs/geo-files/gfs.canopy.t" + + HHI + + "z." + + YY + + MM + + "01.sfcf000.nc", + ] + ) + if os.path.isfile(f_met) is True: + os.chmod(f_frp, 0o0755) + print("---- Download complete!") + else: + print("---- No available FRP data. Terminated!") + exit() + + print("-----------!!!WARNING!!!-------------") + print("---!!!Climatological FRP is used!!!--") + + +os.system("cp " + f_met + " " + f_output) # copy gfs met file for appending + + +"""Reading dimensions""" +print("------------------------------------") +print("---- Checking variable dimensions...") +print("------------------------------------") +readin = Dataset(f_met) +grid_yt = readin["grid_yt"][:] +grid_xt = readin["grid_xt"][:] +lat = readin["lat"][:] +lon = readin["lon"][:] +time = readin["time"][:] + + +# dimension sizes +ntime = len(time) +nlat = len(grid_yt) +nlon = len(grid_xt) + + +# var check +print("time", time.shape) +print("grid_yt", grid_yt.shape) +print("grid_xt", grid_xt.shape) +print("lat", lat.shape) +print("lon", lon.shape) + + +"""Adding canvar""" +print("------------------------------------") +print("---- Generating canopy variables...") +print("------------------------------------") + +for i in np.arange(len(canlist)): + varname = canlist[i] + + print("---- " + varname + " processing...") + + if varname == "lai": + ATTNAME = ["long_name", "units", "missing_value"] + ATT = ["Leaf area index", "m^2/m^2", fill_value] + DATA = read_gfs_climatology(f_can, lat, lon, "lai") + + elif varname == "clu": + ATTNAME = ["long_name", "units", "missing_value"] + ATT = ["Canopy clumping index", "none", fill_value] + DATA = read_gfs_climatology(f_can, lat, lon, "clu") + + elif varname == "ffrac": + ATTNAME = ["long_name", "units", "missing_value"] + ATT = ["Forest fraction of grid cell", "none", fill_value] + DATA = read_gfs_climatology(f_can, lat, lon, "ffrac") + + elif varname == "fh": + ATTNAME = ["long_name", "units", "missing_value"] + ATT = ["Canopy height above the surface", "m", fill_value] + DATA = read_gfs_climatology(f_can, lat, lon, "fh") + + elif varname == "mol": + # Reference: + # Essa 1999, ESTIMATION OF MONIN-OBUKHOV LENGTH USING RICHARDSON AND BULK RICHARDSON NUMBER + # https://inis.iaea.org/collection/NCLCollectionStore/_Public/37/118/37118528.pdf + ATTNAME = ["long_name", "units", "missing_value"] + ATT = ["Monin-Obukhov length", "m", fill_value] + + readin = Dataset(f_met) + t2m = np.squeeze(readin["tmp2m"][:]) + fricv = np.squeeze(readin["fricv"][:]) + shtfl = np.squeeze(readin["shtfl"][:]) + + DATA = (-1 * den * Cp * t2m * (fricv**3)) / (K * g * shtfl) + DATA[DATA > 500] = 500 + DATA[DATA < -500] = -500 + + del [readin, t2m, fricv, shtfl] + + elif varname == "csz": + ATTNAME = ["long_name", "units", "missing_value"] + ATT = ["Cosine of solar zenith angle", "none", fill_value] + + time_conv = datetime( + int(YY), int(MM), int(DD), int(HH), 0, 0, 0, tzinfo=timezone.utc + ) + sza = 90 - get_altitude(lat, lon, time_conv) + DATA = np.cos(sza * 0.0174532925) # degree to radian + + del [time_conv, sza] + + elif varname == "frp": + ATTNAME = ["long_name", "units", "missing_value"] + ATT = ["Mean fire radiative power", "MW", fill_value] + + if frp_src == 2: + DATA = read_gfs_climatology(f_can, "frp") + elif frp_src == 3: + DATA = np.empty(lat.shape) + DATA[:] = 1 + else: + DATA = read_frp_local(f_frp, lat, lon, fill_value) + + elif varname == "href": + ATTNAME = ["long_name", "units", "missing_value"] + ATT = ["Reference height above the canopy", "m", fill_value] + DATA = np.empty([nlat, nlon]) + DATA[:] = ref_lev + + # var check + print("Dimension/Attributes:") + print(DATA.shape) + print(ATTNAME) + print(ATT) + + # adding to output file + output = Dataset(f_output, "a") + + var = output.createVariable( + varname, "float", ("time", "grid_yt", "grid_xt"), fill_value=fill_value + ) + write_varatt(var, ATTNAME, ATT) + var[:] = DATA + + output.close() + + print("---- " + varname + " complete!") + + del [output, var] + del [varname, DATA, ATTNAME, ATT] + + +endtime = datetime.now() + + +print("------------------------------------") +print("---- Global input pre-process complete!", endtime.strftime("%Y/%m/%d %H:%M:%S")) +print("---- Process time:", str(endtime - starttime)) +print("------------------------------------") diff --git a/src/Makefile b/src/Makefile index eeca9da2..681d98d6 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,7 +1,7 @@ # # Use `FC= ... make` to select compiler # Use `DEBUG=1 ... make` for a debug build -# Use `NETCDF=1 ... make` to build with NetCDF using `nf-config` +# Use `NC=1 ... make` to build with NetCDF using `nf-config` # Compiler FC ?= gfortran @@ -14,33 +14,38 @@ $(info FC setting: '$(FC)') DEBUG ?= 0 # Default to NetCDF build -NETCDF ?= 1 +NC ?= 1 # Compile flags $(info DEBUG setting: '$(DEBUG)') ifeq ($(DEBUG), 1) - FCFLAGS := -g -Wall -Wextra -Wconversion -Og -pedantic -fcheck=bounds -fmax-errors=5 + ifeq ($(findstring gfortran,$(notdir $(FC))),gfortran) + FCFLAGS := -g -Wall -Wextra -Wconversion -Og -pedantic -fcheck=bounds -fmax-errors=5 -std=f2003 -fall-intrinsics + else ifeq ($(FC),ifort) + FCFLAGS := -g -warn all -check bounds -implicitnone -O0 -error-limit 5 + endif else ifeq ($(DEBUG), 0) FCFLAGS := -O3 else $(error invalid setting for DEBUG, should be 0 or 1 but is '$(DEBUG)') endif +$(info FCFLAGS: '$(FCFLAGS)') # NETCDF Settings here LIBS := INC := -$(info NETCDF setting: '$(NETCDF)') -ifeq ($(NETCDF), 1) +$(info NC setting: '$(NC)') +ifeq ($(NC), 1) NETCDF_FLIBS := $(shell nf-config --flibs) NETCDF_INC := -I$(shell nf-config --includedir) # LIBS += $(NETCDF_FLIBS) INC += $(NETCDF_INC) FCFLAGS += -DNETCDF -else ifeq ($(NETCDF), 0) +else ifeq ($(NC), 0) # else - $(error invalid setting for NETCDF, should be 0 or 1 but is '$(NETCDF)') + $(error invalid setting for NC, should be 0 or 1 but is '$(NC)') endif $(info LIBS: '$(LIBS)') $(info INC: '$(INC)') @@ -63,9 +68,11 @@ OBJS :=\ canopy_read_txt.o \ canopy_dxcalc_mod.o \ canopy_profile_mod.o \ + canopy_phot_mod.o \ + canopy_rad_mod.o \ + canopy_tleaf_mod.o \ canopy_wind_mod.o \ canopy_waf_mod.o \ - canopy_phot_mod.o \ canopy_eddy_mod.o \ canopy_bioparm_mod.o \ canopy_bioemi_mod.o \ @@ -74,7 +81,7 @@ OBJS :=\ canopy_dealloc.o \ canopy_app.o -ifeq ($(NETCDF), 0) +ifeq ($(NC), 0) _ncf_objs := canopy_check_input.o canopy_ncf_io_mod.o OBJS := $(filter-out $(_ncf_objs),$(OBJS)) endif diff --git a/src/canopy_alloc.F90 b/src/canopy_alloc.F90 index 46afb75a..517de2bf 100644 --- a/src/canopy_alloc.F90 +++ b/src/canopy_alloc.F90 @@ -25,8 +25,15 @@ SUBROUTINE canopy_alloc ! Allocate arrays for Internal Canopy Distribution Variables !------------------------------------------------------------------------------- - if(.not.allocated(zhc)) allocate(zhc(modlays)) - if(.not.allocated(fafraczInt)) allocate(fafraczInt(modlays)) + if(.not.allocated(zhc)) allocate(zhc(modlays)) + if(.not.allocated(fafraczInt)) allocate(fafraczInt(modlays)) + if(.not.allocated(fsun)) allocate(fsun(modlays)) + if(.not.allocated(tleaf_sun)) allocate(tleaf_sun(modlays)) + if(.not.allocated(tleaf_shade)) allocate(tleaf_shade(modlays)) + if(.not.allocated(tleaf_ave)) allocate(tleaf_ave(modlays)) + if(.not.allocated(ppfd_sun)) allocate(ppfd_sun(modlays)) + if(.not.allocated(ppfd_shade)) allocate(ppfd_shade(modlays)) + if(.not.allocated(ppfd_ave)) allocate(ppfd_ave(modlays)) !------------------------------------------------------------------------------- ! Allocate arrays for Canopy Wind Outputs diff --git a/src/canopy_bioemi_mod.F90 b/src/canopy_bioemi_mod.F90 index 85da7d7c..16d618ba 100644 --- a/src/canopy_bioemi_mod.F90 +++ b/src/canopy_bioemi_mod.F90 @@ -5,8 +5,10 @@ module canopy_bioemi_mod contains !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, CLU, COSZEN, SFCRAD, & - TEMP2, LU_OPT, VTYPE, MODRES, CCE, EMI_IND, EMI_OUT) + SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, PPFD_SUN, & + PPFD_SHADE, TLEAF_SUN, TLEAF_SHADE, TLEAF_AVE, LU_OPT, & + VTYPE, MODRES, CCE, VERT, CO2OPT, CO2SET, & + EMI_IND, EMI_OUT) !----------------------------------------------------------------------- @@ -14,7 +16,7 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, CLU, COSZEN, SFCRAD, & ! computes parameterized canopy biogenic emissions ! Preconditions: -! in-canopy FCLAI, model LAI, surface/skin temperature. +! in-canopy FCLAI, model LAI, etc. ! Subroutines and Functions Called: @@ -29,46 +31,38 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, CLU, COSZEN, SFCRAD, & !----------------------------------------------------------------------- ! Jan 2023 P.C. Campbell: Initial canopy isoprene only version ! Feb 2023 P.C. Campbell: Modified for multiple biogenic species +! Jul 2023 P.C. Campbell: Restructured to use FSUN, TLEAF, and PPFD +! as inputs !----------------------------------------------------------------------- !----------------------------------------------------------------------- use canopy_const_mod, ONLY: rk,rgasuniv !constants for canopy models - use canopy_utils_mod, ONLY: interp_linear1_internal - use canopy_phot_mod + use canopy_utils_mod, ONLY: interp_linear1_internal, GET_GAMMA_CO2 use canopy_bioparm_mod ! Arguments: ! IN/OUT - REAL(RK), INTENT( IN ) :: ZK(:) ! Input model heights (m) - REAL(RK), INTENT( IN ) :: FCLAI(:) ! Input Fractional (z) shapes of the + REAL(RK), INTENT( IN ) :: ZK(:) ! Model heights (m) + REAL(RK), INTENT( IN ) :: FCLAI(:) ! Fractional (z) shapes of the ! plant surface distribution (nondimensional), i.e., a Fractional Culmulative LAI - REAL(RK), INTENT( IN ) :: FCH ! Model input canopy height (m) - REAL(RK), INTENT( IN ) :: LAI ! Model input total Leaf Area Index - REAL(RK), INTENT( IN ) :: CLU ! Model input Clumping Index - REAL(RK), INTENT( IN ) :: COSZEN ! Model input Cosine Solar Zenith Angle - REAL(RK), INTENT( IN ) :: SFCRAD ! Model input Instantaneous surface downward shortwave flux (W/m2) - REAL(RK), INTENT( IN ) :: TEMP2 ! Model input 2-m Temperature (K) + REAL(RK), INTENT( IN ) :: FCH ! Canopy height (m) + REAL(RK), INTENT( IN ) :: LAI ! Total Leaf Area Index + REAL(RK), INTENT( IN ) :: FSUN(:) ! Sunlit/Shaded fraction from photolysis correction factor + REAL(RK), INTENT( IN ) :: PPFD_SUN(:) ! PPFD for sunlit leaves (umol phot/m2 s) + REAL(RK), INTENT( IN ) :: PPFD_SHADE(:) ! PPFD for shaded leaves (umol phot/m2 s) + REAL(RK), INTENT( IN ) :: TLEAF_SUN(:) ! Leaf temp for sunlit leaves (K) + REAL(RK), INTENT( IN ) :: TLEAF_SHADE(:) ! Leaf temp for shaded leaves (K) + REAL(RK), INTENT( IN ) :: TLEAF_AVE(:) ! Average Leaf temp (K) INTEGER, INTENT( IN ) :: LU_OPT ! integer for LU type from model mapped to Massman et al. (default = 0/VIIRS) INTEGER, INTENT( IN ) :: VTYPE ! Grid cell dominant vegetation type REAL(RK), INTENT( IN ) :: MODRES ! Canopy model input vertical resolution (m) REAL(RK), INTENT( IN ) :: CCE ! MEGAN Canopy environment coefficient. + INTEGER, INTENT( IN ) :: VERT ! MEGAN vertical integration option (default = 0/no integration) + INTEGER, INTENT( IN ) :: CO2OPT ! Option for co2 inhibition calculation + REAL(RK), INTENT( IN ) :: CO2SET ! User set atmospheric CO2 conc [ppmv] INTEGER, INTENT( IN ) :: EMI_IND ! Input biogenic emissions index REAL(RK), INTENT( OUT ) :: EMI_OUT(:) ! Output canopy layer volume emissions (kg m-3 s-1) ! Local Variables - REAL(RK) :: FSUN(SIZE(ZK)) ! Sunlit/Shaded fraction from photolysis correction factor - REAL(RK) :: PPFD_SUN(SIZE(ZK)) ! PPFD for sunlit leaves (umol phot/m2 s) - REAL(RK) :: PPFD_SHADE(SIZE(ZK)) ! PPFD for shaded leaves (umol phot/m2 s) - REAL(RK) :: ATEMP_SUN(SIZE(ZK)) ! Regression coefficient A for sun leaves (Silva et al., 2020) - REAL(RK) :: BTEMP_SUN(SIZE(ZK)) ! Regression coefficient B for sun leaves - REAL(RK) :: CTEMP_SUN(SIZE(ZK)) ! Regression coefficient C for sun leaves - REAL(RK) :: DTEMP_SUN(SIZE(ZK)) ! Regression coefficient D for sun leaves - REAL(RK) :: ATEMP_SHADE(SIZE(ZK)) ! Regression coefficient A for shade leaves - REAL(RK) :: BTEMP_SHADE(SIZE(ZK)) ! Regression coefficient B for shade leaves - REAL(RK) :: CTEMP_SHADE(SIZE(ZK)) ! Regression coefficient C for shade leaves - REAL(RK) :: DTEMP_SHADE(SIZE(ZK)) ! Regression coefficient D for shade leaves - REAL(RK) :: TLEAF_SUN(SIZE(ZK)) ! Leaf temp for sunlit leaves (K) - REAL(RK) :: TLEAF_SHADE(SIZE(ZK)) ! Leaf temp for shaded leaves (K) - REAL(RK) :: TLEAF_AVE(SIZE(ZK)) ! Average Leaf temp (K) REAL(RK) :: TLEAF24_AVE(SIZE(ZK)) ! Average Leaf temp over the past 24 hours (K) REAL(RK) :: TLEAF240_AVE(SIZE(ZK)) ! Average Leaf temp over the past 240 hours (K) REAL(RK) :: GammaTLEAF_SUN_NUM(SIZE(ZK)) ! Numerator in Tleaf sun activity factor @@ -92,176 +86,19 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, CLU, COSZEN, SFCRAD, & REAL(RK) :: E_OPT(SIZE(ZK)) ! maximum normalized emission capacity REAL(RK) :: TLEAF_OPT(SIZE(ZK)) ! Tleaf at which E_OPT occurs (K) REAL(RK) :: FLAI(SIZE(ZK)) ! Fractional LAI in layer + REAL(RK) :: VPGWT(SIZE(ZK)) ! MEGANv3-like in-canopy weighting factor + REAL(RK) :: GAUSS(SIZE(ZK)) ! MEGANv3-like in-canopy gaussian REAL(RK) :: CT1 ! Activation energy (kJ/mol) REAL(RK) :: CEO ! Empirical coefficient REAL(RK) :: EF ! Final Mapped Emission factor (EF) (ug/m2 hr) - integer i + REAL(RK) :: GAMMACO2 ! CO2 inhibition factor (isoprene only) + integer i, LAYERS ! Constant Canopy Parameters - REAL(RK), PARAMETER :: FRAC_PAR = 0.5_rk !Fraction of incoming solar irradiance that is PAR REAL(RK), PARAMETER :: PPFD0_SUN = 200.0 !Constant PPFDo sunlit (umol/m2 s) (Guenther et al.,2012) REAL(RK), PARAMETER :: PPFD0_SHADE = 50.0 !Constant PPFDo shaded (umol/m2 s) (Guenther et al.,2012) - REAL(RK), PARAMETER :: ATEMP_1_SUN = -13.891_rk !Linearized 2-m temp --> leaf temp parameters (Level 1 = - !top of canopy - REAL(RK), PARAMETER :: ATEMP_2_SUN = -12.322_rk !Based on Table 1 in Silva et al. (2022) - REAL(RK), PARAMETER :: ATEMP_3_SUN = -1.032_rk ! - REAL(RK), PARAMETER :: ATEMP_4_SUN = -5.172_rk ! - REAL(RK), PARAMETER :: ATEMP_5_SUN = -5.589_rk !... - REAL(RK), PARAMETER :: BTEMP_1_SUN = 1.064_rk !... - REAL(RK), PARAMETER :: BTEMP_2_SUN = 1.057_rk !... - REAL(RK), PARAMETER :: BTEMP_3_SUN = 1.031_rk !... - REAL(RK), PARAMETER :: BTEMP_4_SUN = 1.050_rk !... - REAL(RK), PARAMETER :: BTEMP_5_SUN = 1.051_rk !... - REAL(RK), PARAMETER :: CTEMP_1_SUN = 1.083_rk !Exponential 2-m PPFD --> PPFD parameters (Level 1 = - !top of canopy - REAL(RK), PARAMETER :: CTEMP_2_SUN = 1.096_rk !Based on Table 1 in Silva et al. (2022) - REAL(RK), PARAMETER :: CTEMP_3_SUN = 1.104_rk ! - REAL(RK), PARAMETER :: CTEMP_4_SUN = 1.098_rk ! - REAL(RK), PARAMETER :: CTEMP_5_SUN = 1.090_rk !... - REAL(RK), PARAMETER :: DTEMP_1_SUN = 0.002_rk !... - REAL(RK), PARAMETER :: DTEMP_2_SUN = -0.128_rk !... - REAL(RK), PARAMETER :: DTEMP_3_SUN = -0.298_rk !... - REAL(RK), PARAMETER :: DTEMP_4_SUN = -0.445_rk !... - REAL(RK), PARAMETER :: DTEMP_5_SUN = -0.535_rk !... - REAL(RK), PARAMETER :: ATEMP_1_SHADE = -12.846_rk !... - REAL(RK), PARAMETER :: ATEMP_2_SHADE = -11.343_rk !... - REAL(RK), PARAMETER :: ATEMP_3_SHADE = -1.068_rk !... - REAL(RK), PARAMETER :: ATEMP_4_SHADE = -5.551_rk !... - REAL(RK), PARAMETER :: ATEMP_5_SHADE = -5.955_rk !... - REAL(RK), PARAMETER :: BTEMP_1_SHADE = 1.060_rk !... - REAL(RK), PARAMETER :: BTEMP_2_SHADE = 1.053_rk !... - REAL(RK), PARAMETER :: BTEMP_3_SHADE = 1.031_rk !... - REAL(RK), PARAMETER :: BTEMP_4_SHADE = 1.051_rk !... - REAL(RK), PARAMETER :: BTEMP_5_SHADE = 1.053_rk !... - REAL(RK), PARAMETER :: CTEMP_1_SHADE = 0.871_rk !... - REAL(RK), PARAMETER :: CTEMP_2_SHADE = 0.890_rk !... - REAL(RK), PARAMETER :: CTEMP_3_SHADE = 0.916_rk !... - REAL(RK), PARAMETER :: CTEMP_4_SHADE = 0.941_rk !... - REAL(RK), PARAMETER :: CTEMP_5_SHADE = 0.956_rk !... - REAL(RK), PARAMETER :: DTEMP_1_SHADE = 0.015_rk !... - REAL(RK), PARAMETER :: DTEMP_2_SHADE = -0.141_rk !... - REAL(RK), PARAMETER :: DTEMP_3_SHADE = -0.368_rk !... - REAL(RK), PARAMETER :: DTEMP_4_SHADE = -0.592_rk !... - REAL(RK), PARAMETER :: DTEMP_5_SHADE = -0.743_rk !... - REAL(RK), PARAMETER :: CT2 = 230.0_rk !Deactivation energy (kJ/mol) (Guenther et al., 2012) -!Calculate photolyis shading/correction factor through canopy, i.e., the fraction of sunlit leaves downward through canopy -! `canopy_phot` gives relative direct beam irradiance, -! which, multiplied by clumping index, gives sunlit fraction (e.g., Bonan 2019, eq. 14.18) - - call canopy_phot(FCLAI, LAI, CLU, COSZEN, FSUN) - FSUN = FSUN * CLU - -! Use linear canopy temperature model based on Silva et al. (2020) to get approx. sun/shade leaf temperatures -! through canopy (ignores effect of wind speed on leaf boundary layer ~ 1 % error/bias) -!Citation: -!Silva, S. J., Heald, C. L., and Guenther, A. B.: Development of a reduced-complexity plant canopy -!physics surrogate model for use in chemical transport models: a case study with GEOS-Chem v12.3.0, -!Geosci. Model Dev., 13, 2569–2585, https://doi.org/10.5194/gmd-13-2569-2020, 2020. - do i=1, SIZE(ZK) !calculate linear change in parameters interpolated to Silva et al. 5 layer canopy regions - if (ZK(i) .gt. FCH) then ! above canopy, Tleaf = Tair - ATEMP_SUN(i) = 0.0 - BTEMP_SUN(i) = 1.0 - ATEMP_SHADE(i) = 0.0 - BTEMP_SHADE(i) = 1.0 - else if (ZK(i) .le. FCH .and. ZK(i) .gt. FCH*(4.0_rk/5.0_rk)) then !Level 1 - 2 - ATEMP_SUN(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & - (/ ATEMP_2_SUN,ATEMP_1_SUN /),ZK(i)) - BTEMP_SUN(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & - (/ BTEMP_2_SUN,BTEMP_1_SUN /),ZK(i)) - ATEMP_SHADE(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & - (/ ATEMP_2_SHADE,ATEMP_1_SHADE /),ZK(i)) - BTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & - (/ BTEMP_2_SHADE,BTEMP_1_SHADE /),ZK(i)) - else if (ZK(i) .le. FCH*(4.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(3.0_rk/5.0_rk)) then !Level 2 - 3 - ATEMP_SUN(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & - (/ ATEMP_3_SUN,ATEMP_2_SUN /),ZK(i)) - BTEMP_SUN(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & - (/ BTEMP_3_SUN,BTEMP_2_SUN /),ZK(i)) - ATEMP_SHADE(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & - (/ ATEMP_3_SHADE,ATEMP_2_SHADE /),ZK(i)) - BTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & - (/ BTEMP_3_SHADE,BTEMP_2_SHADE /),ZK(i)) - else if (ZK(i) .le. FCH*(3.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(2.0_rk/5.0_rk)) then !Level 3 - 4 - ATEMP_SUN(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & - (/ ATEMP_4_SUN,ATEMP_3_SUN /),ZK(i)) - BTEMP_SUN(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & - (/ BTEMP_4_SUN,BTEMP_3_SUN /),ZK(i)) - ATEMP_SHADE(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & - (/ ATEMP_4_SHADE,ATEMP_3_SHADE /),ZK(i)) - BTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & - (/ BTEMP_4_SHADE,BTEMP_3_SHADE /),ZK(i)) - else if (ZK(i) .le. FCH*(2.0_rk/5.0_rk) ) then !Level 4 - Bottom - ATEMP_SUN(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & - (/ ATEMP_5_SUN,ATEMP_4_SUN /),ZK(i)) - BTEMP_SUN(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & - (/ BTEMP_5_SUN,BTEMP_4_SUN /),ZK(i)) - ATEMP_SHADE(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & - (/ ATEMP_5_SHADE,ATEMP_4_SHADE /),ZK(i)) - BTEMP_SHADE(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & - (/ BTEMP_5_SHADE,BTEMP_4_SHADE /),ZK(i)) - end if - end do - - TLEAF_SUN = ATEMP_SUN + (BTEMP_SUN*TEMP2) - TLEAF_SHADE = ATEMP_SHADE + (BTEMP_SHADE*TEMP2) - TLEAF_AVE = (TLEAF_SUN*FSUN) + (TLEAF_SHADE*(1.0-FSUN)) ! average = sum sun and shade weighted by sunlit fraction - -! Use exponential PPFD model based on Silva et al. (2020) to get approx. sun/shade PPFD -! through canopy -!Citation: -!Silva, S. J., Heald, C. L., and Guenther, A. B.: Development of a reduced-complexity plant canopy -!physics surrogate model for use in chemical transport models: a case study with GEOS-Chem v12.3.0, -!Geosci. Model Dev., 13, 2569–2585, https://doi.org/10.5194/gmd-13-2569-2020, 2020. - do i=1, SIZE(ZK) !calculate linear change in parameters interpolated to Silva et al. 5 layer canopy regions - if (ZK(i) .gt. FCH) then ! above canopy, PPFD_leaf = PPFD_toc (toc=top of canopy) - CTEMP_SUN(i) = 0.0 - DTEMP_SUN(i) = 0.0 - CTEMP_SHADE(i) = 0.0 - DTEMP_SHADE(i) = 0.0 - else if (ZK(i) .le. FCH .and. ZK(i) .gt. FCH*(4.0_rk/5.0_rk)) then !Level 1 - 2 - CTEMP_SUN(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & - (/ CTEMP_2_SUN,CTEMP_1_SUN /),ZK(i)) - DTEMP_SUN(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & - (/ DTEMP_2_SUN,DTEMP_1_SUN /),ZK(i)) - CTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & - (/ CTEMP_2_SHADE,CTEMP_1_SHADE /),ZK(i)) - DTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & - (/ DTEMP_2_SHADE,DTEMP_1_SHADE /),ZK(i)) - else if (ZK(i) .le. FCH*(4.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(3.0_rk/5.0_rk)) then !Level 2 - 3 - CTEMP_SUN(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & - (/ CTEMP_3_SUN,CTEMP_2_SUN /),ZK(i)) - DTEMP_SUN(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & - (/ DTEMP_3_SUN,DTEMP_2_SUN /),ZK(i)) - CTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & - (/ CTEMP_3_SHADE,CTEMP_2_SHADE /),ZK(i)) - DTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & - (/ DTEMP_3_SHADE,DTEMP_2_SHADE /),ZK(i)) - else if (ZK(i) .le. FCH*(3.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(2.0_rk/5.0_rk)) then !Level 3 - 4 - CTEMP_SUN(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & - (/ CTEMP_4_SUN,CTEMP_3_SUN /),ZK(i)) - DTEMP_SUN(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & - (/ DTEMP_4_SUN,DTEMP_3_SUN /),ZK(i)) - CTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & - (/ CTEMP_4_SHADE,CTEMP_3_SHADE /),ZK(i)) - DTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & - (/ DTEMP_4_SHADE,DTEMP_3_SHADE /),ZK(i)) - else if (ZK(i) .le. FCH*(2.0_rk/5.0_rk) ) then !Level 4 - Bottom - CTEMP_SUN(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & - (/ CTEMP_5_SUN,CTEMP_4_SUN /),ZK(i)) - DTEMP_SUN(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & - (/ DTEMP_5_SUN,DTEMP_4_SUN /),ZK(i)) - CTEMP_SHADE(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & - (/ CTEMP_5_SHADE,CTEMP_4_SHADE /),ZK(i)) - DTEMP_SHADE(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & - (/ DTEMP_5_SHADE,DTEMP_4_SHADE /),ZK(i)) - end if - end do - - PPFD_SUN = FRAC_PAR * SFCRAD * EXP(CTEMP_SUN + DTEMP_SUN * LAI) !Silva et al. W/m2 --> umol m-2 s-1 - PPFD_SHADE = FRAC_PAR * SFCRAD * EXP(CTEMP_SHADE + DTEMP_SHADE * LAI) - ! Calculate maximum normalized emission capacity (E_OPT) and Tleaf at E_OPT TLEAF240_AVE = TLEAF_AVE !Assume instantaneous TLEAF estimate for TLEAF240 and TLEAF24 (could improve...) TLEAF24_AVE = TLEAF_AVE @@ -298,15 +135,89 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, CLU, COSZEN, SFCRAD, & GammaPPFD_AVE = (GammaPPFD_SUN*FSUN) + (GammaPPFD_SHADE*(1.0-FSUN)) ! average = sum sun and shade weighted by sunlit fraction +! Get CO2 inhibition factor for isoprene only + + if (EMI_IND .eq. 1) then !Isoprene + GAMMACO2 = GET_GAMMA_CO2(CO2OPT,CO2SET) + else + GAMMACO2 = 1.0_rk + end if + ! Calculate emissions profile in the canopy EMI_OUT = 0.0_rk ! set initial emissions profile to zero - do i=1, SIZE(ZK) - if (ZK(i) .gt. 0.0 .and. ZK(i) .le. FCH) then ! above ground level and at/below canopy top - FLAI(i) = ((FCLAI(i+1) - FCLAI(i)) * LAI)/MODRES !fractional LAI in each layer converted to LAD (m2 m-3) - EMI_OUT(i) = FLAI(i) * EF * GammaTLEAF_AVE(i) * GammaPPFD_AVE(i) * CCE ! (ug m-3 hr-1) - EMI_OUT(i) = EMI_OUT(i) * 2.7777777777778E-13_rk !TBD: convert emissions output to (kg m-3 s-1) - end if - end do + FLAI = 0.0_rk ! set initial fractional FLAI (LAD) profile to zero + + if (VERT .eq. 0) then !Full 3D leaf-level biogenic emissions (no averaging, summing, or integration) + do i=1, SIZE(ZK) + if (ZK(i) .gt. 0.0 .and. ZK(i) .le. FCH) then ! above ground level and at/below canopy top + FLAI(i) = ((FCLAI(i+1) - FCLAI(i)) * LAI)/MODRES !fractional LAI in each layer converted to LAD (m2 m-3) + EMI_OUT(i) = FLAI(i) * EF * GammaTLEAF_AVE(i) * GammaPPFD_AVE(i) * GAMMACO2 * CCE ! (ug m-3 hr-1) + EMI_OUT(i) = EMI_OUT(i) * 2.7777777777778E-13_rk !convert emissions output to (kg m-3 s-1) + end if + end do + else if (VERT .eq. 1) then !"MEGANv3-like": Use weighting factors normalized to plant distribution shape (FCLAI) + !across canopy layers + LAYERS = floor(FCH/MODRES) + 1 + do i=1, SIZE(ZK) + if (ZK(i) .gt. 0.0 .and. ZK(i) .le. FCH) then + FLAI(i) = ((FCLAI(i+1) - FCLAI(i)) * LAI)/MODRES !fractional LAI in each layer converted to LAD (m2 m-3) + end if + end do + do i=1, SIZE(ZK) + if (ZK(i) .gt. 0.0 .and. ZK(i) .le. FCH) then + VPGWT(i) = (FLAI(i))/sum(FLAI(1:LAYERS)) + end if + end do + EMI_OUT(SIZE(ZK)) = LAI * EF * SUM(GammaTLEAF_AVE(1:LAYERS) * GammaPPFD_AVE(1:LAYERS) * & + VPGWT(1:LAYERS)) * GAMMACO2 * CCE !put into top model layer (ug m-2 hr-1) + EMI_OUT = EMI_OUT * 2.7777777777778E-13_rk !convert emissions output to (kg m-2 s-1) + else if (VERT .eq. 2) then !"MEGANv3-like": Add weighted sum of activity coefficients using normal distribution + !across canopy layers using 5 layer numbers directly from MEGANv3 + !--warning: weights are not consistent with FCLAI distribution + !used for biomass distribution used for sunlit/shaded in Gamma TLEAF and GammaPPFD. + LAYERS = floor(FCH/MODRES) + 1 + do i=1, SIZE(ZK) + if (ZK(i) .gt. FCH) then + GAUSS(i) = 0.0 + else if (ZK(i) .le. FCH .and. ZK(i) .gt. FCH*(4.0_rk/5.0_rk)) then !Level 1 - 2 + GAUSS(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & + (/ 0.118464_rk,0.0_rk /),ZK(i)) + else if (ZK(i) .le. FCH*(4.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(3.0_rk/5.0_rk)) then !Level 2 - 3 + GAUSS(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & + (/ 0.239314_rk,0.118464_rk /),ZK(i)) + else if (ZK(i) .le. FCH*(3.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(2.0_rk/5.0_rk)) then !Level 3 - 4 + GAUSS(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & + (/ 0.284444_rk,0.239314_rk /),ZK(i)) + else if (ZK(i) .le. FCH*(2.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(1.0_rk/5.0_rk) ) then !Level 4 - Bottom + GAUSS(i) = interp_linear1_internal((/ FCH*(1.0_rk/5.0_rk),FCH*(2.0_rk/5.0_rk) /), & + (/ 0.239314_rk,0.284444_rk /),ZK(i)) + else if (ZK(i) .le. FCH*(1.0_rk/5.0_rk) ) then !Level 4 - Bottom + GAUSS(i) = interp_linear1_internal((/ ZK(1),FCH*(1.0_rk/5.0_rk) /), & + (/ 0.118464_rk,0.239314_rk /),ZK(i)) + end if + end do + + do i=1, SIZE(ZK) + VPGWT(i) = GAUSS(i)/sum(GAUSS(1:LAYERS)) + end do + EMI_OUT(SIZE(ZK)) = LAI * EF * SUM(GammaTLEAF_AVE(1:LAYERS) * GammaPPFD_AVE(1:LAYERS) * & + VPGWT(1:LAYERS)) * GAMMACO2 * CCE !put into top model layer (ug m-2 hr-1) + EMI_OUT = EMI_OUT * 2.7777777777778E-13_rk !convert emissions output to (kg m-2 s-1) + else if (VERT .eq. 3) then !"MEGANv3-like": Add weighted sum of activity coefficients equally + !across canopy layers + !--warning: weights are not consistent with FCLAI distribution + !used for biomass distribution used for sunlit/shaded in Gamma TLEAF and GammaPPFD. + LAYERS = floor(FCH/MODRES) + 1 + do i=1, SIZE(ZK) + VPGWT(i) = 1.0_rk/LAYERS + end do + EMI_OUT(SIZE(ZK)) = LAI * EF * SUM(GammaTLEAF_AVE(1:LAYERS) * GammaPPFD_AVE(1:LAYERS) * & + VPGWT(1:LAYERS)) * GAMMACO2 * CCE !put into top model layer (ug m-2 hr-1) + EMI_OUT = EMI_OUT * 2.7777777777778E-13_rk !convert emissions output to (kg m-2 s-1) + else + write(*,*) 'Wrong BIOVERT_OPT choice of ', VERT, ' in namelist...exiting' + call exit(2) + end if END SUBROUTINE CANOPY_BIO diff --git a/src/canopy_calcs.F90 b/src/canopy_calcs.F90 index 3f47929b..5ced80cd 100644 --- a/src/canopy_calcs.F90 +++ b/src/canopy_calcs.F90 @@ -15,7 +15,9 @@ SUBROUTINE canopy_calcs use canopy_utils_mod !main canopy utilities use canopy_dxcalc_mod !main canopy dx calculation use canopy_profile_mod !main canopy foliage profile routines - use canopy_wind_mod !main canopy wind components + use canopy_rad_mod !main canopy radiation sunlit/shaded routines + use canopy_tleaf_mod !main canopy leaf temperature sunlit/shaded routines + use canopy_wind_mod !main canopy components use canopy_waf_mod use canopy_phot_mod use canopy_eddy_mod @@ -84,6 +86,38 @@ SUBROUTINE canopy_calcs if (lu_opt .eq. 0 .or. lu_opt .eq. 1 ) then !VIIRS or MODIS if (vtyperef .gt. 0 .and. vtyperef .le. 10 .or. vtyperef .eq. 12) then !VIIRS or MODIS types +! ... check for ssg_opt from user namelist + if (vtyperef .ge. 6 .and. vtyperef .le. 10) then !VIIRS/MODIS shrubs/savannas/grasses (SSG) type + if (ssg_opt .eq. 0) then !use GEDI inputs for SSG heights (not likely captured...) + hcmref = hcmref + else if (ssg_opt .eq. 1) then !user set constant shrubs/savannas/grasslands height + hcmref = ssg_set + !recalculate + zhc = zk/hcmref + cansublays = floor(hcmref/modres) + else + write(*,*) 'Wrong SSG_OPT choice of ', ssg_opt, & + ' in namelist...exiting' + call exit(2) + end if + end if + +! ... check for crop_opt from user namelist + if (vtyperef .eq. 12) then !VIIRS/MODIS crop type + if (crop_opt .eq. 0) then !use GEDI inputs for crop height (not likely captured...) + hcmref = hcmref + else if (crop_opt .eq. 1) then !user set constant crop height + hcmref = crop_set + !recalculate + zhc = zk/hcmref + cansublays = floor(hcmref/modres) + else + write(*,*) 'Wrong CROP_OPT choice of ', crop_opt, & + ' in namelist...exiting' + call exit(2) + end if + end if + ! ... check for contiguous canopy conditions at each model grid cell if (hcmref .gt. fch_thresh .and. ffracref .gt. frt_thresh & .and. lairef .gt. lai_thresh) then @@ -105,6 +139,21 @@ SUBROUTINE canopy_calcs ubzref, z0ghc, lambdars, cdrag, pai, hcmref, hgtref, & z0ref, vtyperef, lu_opt, z0_opt, d_h, zo_h) +! ... calculate canopy radiation (sunlit and shade) profile + + call canopy_fsun_clu( fafraczInt, lairef, cluref, cszref, fsun) + +! ... calculate canopy leaf temperature (sun/shade) profile + + call canopy_tleaf_lin(zk, hcmref, tmp2mref, fsun, & + tleaf_sun, tleaf_shade, tleaf_ave) + +! ... calculate canopy Photosynthetic Photon Flux Density (PPFD) (sun/shade) profile + + call canopy_ppfd_exp(zk, hcmref, dswrfref, lairef, fsun, & + ppfd_sun, ppfd_shade, ppfd_ave) + +! ... **** User Canopy-App Options *** ! ... user option to calculate in-canopy wind speeds at height z and midflame WAF if (ifcanwind .or. ifcanwaf) then @@ -122,12 +171,26 @@ SUBROUTINE canopy_calcs ! ... determine midflamepoint and flame height from user or FRP calculation call canopy_flameh(flameh_opt, flameh_set, dx_2d(i,j), modres, & frpref, frp_fac, hcmref, midflamepoint, flameh_2d(i,j)) - - if (flameh_2d(i,j) .gt. 0.0 .and. flameh_2d(i,j) .le. hcmref) then - !only calculate WAF when flameh > 0 and <= FH - call canopy_waf(hcmref, lambdars, hgtref, flameh_2d(i,j), & - firetype, d_h, zo_h, canBOT(midflamepoint), & - canTOP(midflamepoint), waf_2d(i,j)) + if (firetype .eq. 0) then !forest/sub-canopy firetype + if (flameh_2d(i,j) .gt. 0.0) then !flameh must be > 0 + if (flameh_2d(i,j) .le. hcmref) then !only calculate when flameh <= FCH + call canopy_waf(hcmref, lambdars, hgtref, flameh_2d(i,j), & + firetype, d_h, zo_h, canBOT(midflamepoint), & + canTOP(midflamepoint), waf_2d(i,j)) + else + write(*,*) 'warning...sub-canopy type fire, but flameh > FCH, setting WAF=1' + waf_2d(i,j) = 1.0_rk + end if + end if + else !grass/crops, above-canopy firetype + if (flameh_2d(i,j) .gt. 0.0) then !flameh still must be > 0 + call canopy_waf(hcmref, lambdars, hgtref, flameh_2d(i,j), & + firetype, d_h, zo_h, canBOT(midflamepoint), & + canTOP(midflamepoint), waf_2d(i,j)) + end if + end if + if (waf_2d(i,j) .gt. 1.0_rk) then !Final check of WAF > 1, must be <=1 + waf_2d(i,j) = 1.0_rk end if end if @@ -152,80 +215,118 @@ SUBROUTINE canopy_calcs .and. cluref .gt. 0.0_rk) then !ISOP call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 1, emi_isop_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 1, emi_isop_3d(i,j,:)) !MYRC call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 2, emi_myrc_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 2, emi_myrc_3d(i,j,:)) !SABI call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 3, emi_sabi_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 3, emi_sabi_3d(i,j,:)) !LIMO call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 4, emi_limo_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 4, emi_limo_3d(i,j,:)) !CARE call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 5, emi_care_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 5, emi_care_3d(i,j,:)) !OCIM call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 6, emi_ocim_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 6, emi_ocim_3d(i,j,:)) !BPIN call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 7, emi_bpin_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 7, emi_bpin_3d(i,j,:)) !APIN call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 8, emi_apin_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 8, emi_apin_3d(i,j,:)) !MONO call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 9, emi_mono_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 9, emi_mono_3d(i,j,:)) !FARN call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 10, emi_farn_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 10, emi_farn_3d(i,j,:)) !CARY call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 11, emi_cary_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 11, emi_cary_3d(i,j,:)) !SESQ call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 12, emi_sesq_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 12, emi_sesq_3d(i,j,:)) !MBOL call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 13, emi_mbol_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 13, emi_mbol_3d(i,j,:)) !METH call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 14, emi_meth_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 14, emi_meth_3d(i,j,:)) !ACET call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 15, emi_acet_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 15, emi_acet_3d(i,j,:)) !CO call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 16, emi_co_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 16, emi_co_3d(i,j,:)) !BIDI VOC call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 17, emi_bvoc_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 17, emi_bvoc_3d(i,j,:)) !Stress VOC call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 18, emi_svoc_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 18, emi_svoc_3d(i,j,:)) !Other VOC call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 19, emi_ovoc_3d(i,j,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 19, emi_ovoc_3d(i,j,:)) end if end if @@ -296,6 +397,38 @@ SUBROUTINE canopy_calcs if (lu_opt .eq. 0 .or. lu_opt .eq. 1 ) then !VIIRS or MODIS if (vtyperef .gt. 0 .and. vtyperef .le. 10 .or. vtyperef .eq. 12) then !VIIRS or MODIS types +! ... check for ssg_opt from user namelist + if (vtyperef .ge. 6 .and. vtyperef .le. 10) then !VIIRS/MODIS shrubs/savannas/grasses (SSG) type + if (ssg_opt .eq. 0) then !use GEDI inputs for SSG heights (not likely captured...) + hcmref = hcmref + else if (ssg_opt .eq. 1) then !user set constant shrubs/savannas/grasslands height + hcmref = ssg_set + !recalculate + zhc = zk/hcmref + cansublays = floor(hcmref/modres) + else + write(*,*) 'Wrong SSG_OPT choice of ', ssg_opt, & + ' in namelist...exiting' + call exit(2) + end if + end if + +! ... check for crop_opt from user namelist + if (vtyperef .eq. 12) then !VIIRS/MODIS crop type + if (crop_opt .eq. 0) then !use GEDI inputs for crop height + hcmref = hcmref + else if (crop_opt .eq. 1) then !user set constant crop height + hcmref = crop_set + !recalculate + zhc = zk/hcmref + cansublays = floor(hcmref/modres) + else + write(*,*) 'Wrong CROP_OPT choice of ', crop_opt, & + ' in namelist...exiting' + call exit(2) + end if + end if + ! ... check for contiguous canopy conditions at each model grid cell if (hcmref .gt. fch_thresh .and. ffracref .gt. frt_thresh & .and. lairef .gt. lai_thresh) then @@ -317,6 +450,22 @@ SUBROUTINE canopy_calcs ubzref, z0ghc, lambdars, cdrag, pai, hcmref, hgtref, & z0ref, vtyperef, lu_opt, z0_opt, d_h, zo_h) +! ... calculate canopy radiation (sunlit and shade) profile + + call canopy_fsun_clu( fafraczInt, lairef, cluref, cszref, fsun) + +! ... calculate canopy leaf temperature (sun/shade) profile + + call canopy_tleaf_lin(zk, hcmref, tmp2mref, fsun, & + tleaf_sun, tleaf_shade, tleaf_ave) + +! ... calculate canopy Photosynthetic Photon Flux Density (PPFD) (sun/shade) profile + + call canopy_ppfd_exp(zk, hcmref, dswrfref, lairef, fsun, & + ppfd_sun, ppfd_shade, ppfd_ave) + +! ... **** User Canopy-App Options *** + ! ... user option to calculate in-canopy wind speeds at height z and midflame WAF if (ifcanwind .or. ifcanwaf) then @@ -334,12 +483,26 @@ SUBROUTINE canopy_calcs ! ... determine midflamepoint and flame height from user or FRP calculation call canopy_flameh(flameh_opt, flameh_set, dx(loc), modres, & frpref, frp_fac, hcmref, midflamepoint, flameh(loc)) - - if (flameh(loc) .gt. 0.0 .and. flameh(loc) .le. hcmref) then - !only calculate WAF when flameh > 0 - call canopy_waf(hcmref, lambdars, hgtref, flameh(loc), & - firetype, d_h, zo_h, canBOT(midflamepoint), & - canTOP(midflamepoint), waf(loc)) + if (firetype .eq. 0) then !forest/sub-canopy firetype + if (flameh(loc) .gt. 0.0) then !flameh must be > 0 + if (flameh(loc) .le. hcmref) then !only calculate when flameh <= FCH + call canopy_waf(hcmref, lambdars, hgtref, flameh(loc), & + firetype, d_h, zo_h, canBOT(midflamepoint), & + canTOP(midflamepoint), waf(loc)) + else + write(*,*) 'warning...sub-canopy type fire, but flameh > FCH, setting WAF=1' + waf(loc) = 1.0_rk + end if + end if + else !grass/crops, above-canopy firetype + if (flameh(loc) .gt. 0.0) then !flameh still must be > 0 + call canopy_waf(hcmref, lambdars, hgtref, flameh(loc), & + firetype, d_h, zo_h, canBOT(midflamepoint), & + canTOP(midflamepoint), waf(loc)) + end if + end if + if (waf(loc) .gt. 1.0_rk) then !Final check of WAF > 1, must be <=1 + waf(loc) = 1.0_rk end if end if @@ -364,80 +527,118 @@ SUBROUTINE canopy_calcs .and. cluref .gt. 0.0_rk) then !ISOP call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 1, emi_isop(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 1, emi_isop(loc,:)) !MYRC call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 2, emi_myrc(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 2, emi_myrc(loc,:)) !SABI call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 3, emi_sabi(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 3, emi_sabi(loc,:)) !LIMO call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 4, emi_limo(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 4, emi_limo(loc,:)) !CARE call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 5, emi_care(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 5, emi_care(loc,:)) !OCIM call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 6, emi_ocim(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 6, emi_ocim(loc,:)) !BPIN call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 7, emi_bpin(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 7, emi_bpin(loc,:)) !APIN call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 8, emi_apin(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 8, emi_apin(loc,:)) !MONO call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 9, emi_mono(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 9, emi_mono(loc,:)) !FARN call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 10, emi_farn(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 10, emi_farn(loc,:)) !CARY call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 11, emi_cary(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 11, emi_cary(loc,:)) !SESQ call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 12, emi_sesq(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 12, emi_sesq(loc,:)) !MBOL call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 13, emi_mbol(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 13, emi_mbol(loc,:)) !METH call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 14, emi_meth(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 14, emi_meth(loc,:)) !ACET call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 15, emi_acet(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 15, emi_acet(loc,:)) !CO call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 16, emi_co(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 16, emi_co(loc,:)) !BIDI VOC call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 17, emi_bvoc(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 17, emi_bvoc(loc,:)) !Stress VOC call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 18, emi_svoc(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 18, emi_svoc(loc,:)) !Other VOC call canopy_bio(zk, fafraczInt, hcmref, & - lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 19, emi_ovoc(loc,:)) + lairef, fsun, ppfd_sun, ppfd_shade, tleaf_sun, tleaf_shade, & + tleaf_ave, & + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 19, emi_ovoc(loc,:)) end if end if diff --git a/src/canopy_canopts_mod.F90 b/src/canopy_canopts_mod.F90 index 2fd44171..b3f1c29c 100644 --- a/src/canopy_canopts_mod.F90 +++ b/src/canopy_canopts_mod.F90 @@ -33,5 +33,12 @@ MODULE canopy_canopts_mod real(rk) :: z0ghc !ratio of ground roughness length to canopy top height real(rk) :: lambdars !Value representing influence of roughness sublayer (nondimensional) real(rk) :: bio_cce !MEGAN biogenic emission canopy environment coefficient. + integer :: biovert_opt !MEGAN vertical integration of emissions option (default = 0, off) + integer :: ssg_opt !Set default integer for shrubs/savanna/grassland vegtype option from GEDI or user (default = 0) + real(rk) :: ssg_set !Set default value for shrubs/savanna/grassland vegtype heights used in model (m) (Default = 1 m) + integer :: crop_opt !Set default integer for crop vegtype option from GEDI or user (default = 0) + real(rk) :: crop_set !Set default value for crop vegtype heights used in model (m) (Default = 3 m) + integer :: co2_opt !Set default integer for co2 inhibition option for biogenic isoprene emissions (default=0; Possell & Hewitt (2011)) + real(rk) :: co2_set !Set default value for atmospheric co2 concentration for co2_opt (m) (Default = 400.0 ppmv) END MODULE canopy_canopts_mod diff --git a/src/canopy_canvars_mod.F90 b/src/canopy_canvars_mod.F90 index 9365262c..5519e1f2 100644 --- a/src/canopy_canvars_mod.F90 +++ b/src/canopy_canvars_mod.F90 @@ -33,6 +33,13 @@ MODULE canopy_canvars_mod real(rk), allocatable :: fainc ( : ) ! incremental foliage shape function real(rk), allocatable :: fafracz ( : ) ! incremental fractional foliage shape function real(rk), allocatable :: fafraczInt ( : ) ! integral of incremental fractional foliage shape function + real(rk), allocatable :: fsun ( : ) ! Sunlit/Shaded fraction from photolysis correction factor + real(rk), allocatable :: tleaf_sun ( : ) ! Leaf temp for sunlit leaves (K) + real(rk), allocatable :: tleaf_shade ( : ) ! Leaf temp for shaded leaves (K) + real(rk), allocatable :: tleaf_ave ( : ) ! Average Leaf temp for sunlit and shaded leaves (K) + real(rk), allocatable :: ppfd_sun ( : ) ! PPFD for sunlit leaves (umol phot/m2 s) + real(rk), allocatable :: ppfd_shade ( : ) ! PPFD for shaded leaves (umol phot/m2 s) + real(rk), allocatable :: ppfd_ave ( : ) ! Average PPFD for sunlit and shaded leaves (umol phot/m2 s) real(rk), allocatable :: canBOT ( : ) ! Canopy bottom wind reduction factors real(rk), allocatable :: canTOP ( : ) ! Canopy top wind reduction factors real(rk), allocatable :: canWIND ( : , : ) ! canopy wind speeds (m/s) diff --git a/src/canopy_dealloc.F90 b/src/canopy_dealloc.F90 index 2094a096..29f015bc 100644 --- a/src/canopy_dealloc.F90 +++ b/src/canopy_dealloc.F90 @@ -25,11 +25,18 @@ SUBROUTINE canopy_dealloc ! Dellocate arrays for Canopy Distribution !------------------------------------------------------------------------------- - if(allocated(zk)) deallocate(zk) !allocated in canopy_readnml - if(allocated(zhc)) deallocate(zhc) - if(allocated(fainc)) deallocate(fainc) !allocated in canopy_profile - if(allocated(fafracz)) deallocate(fafracz) !allocated in canopy_profile - if(allocated(fafraczInt)) deallocate(fafraczInt) + if(allocated(zk)) deallocate(zk) !allocated in canopy_readnml + if(allocated(zhc)) deallocate(zhc) + if(allocated(fainc)) deallocate(fainc) !allocated in canopy_profile + if(allocated(fafracz)) deallocate(fafracz) !allocated in canopy_profile + if(allocated(fafraczInt)) deallocate(fafraczInt) + if(allocated(fsun)) deallocate(fsun) + if(allocated(tleaf_sun)) deallocate(tleaf_sun) + if(allocated(tleaf_shade)) deallocate(tleaf_shade) + if(allocated(tleaf_ave)) deallocate(tleaf_ave) + if(allocated(ppfd_sun)) deallocate(ppfd_sun) + if(allocated(ppfd_shade)) deallocate(ppfd_shade) + if(allocated(ppfd_ave)) deallocate(ppfd_ave) !------------------------------------------------------------------------------- ! Deallocate arrays for Canopy Wind diff --git a/src/canopy_init.F90 b/src/canopy_init.F90 index d7af1183..9a0bfa34 100644 --- a/src/canopy_init.F90 +++ b/src/canopy_init.F90 @@ -19,8 +19,15 @@ SUBROUTINE canopy_init ! Initialize arrays for Canopy Distribution !------------------------------------------------------------------------------- - if(allocated(zhc)) zhc(:) = fillreal - if(allocated(fafraczInt)) fafraczInt(:) = fillreal + if(allocated(zhc)) zhc(:) = fillreal + if(allocated(fafraczInt)) fafraczInt(:) = fillreal + if(allocated(fsun)) fsun(:) = fillreal + if(allocated(tleaf_sun)) tleaf_sun(:) = fillreal + if(allocated(tleaf_shade)) tleaf_shade(:) = fillreal + if(allocated(tleaf_ave)) tleaf_ave(:) = fillreal + if(allocated(ppfd_sun)) ppfd_sun(:) = fillreal + if(allocated(ppfd_shade)) ppfd_shade(:) = fillreal + if(allocated(ppfd_ave)) ppfd_ave(:) = fillreal !------------------------------------------------------------------------------- ! Initialize arrays for Canopy Wind diff --git a/src/canopy_profile_mod.F90 b/src/canopy_profile_mod.F90 index 5882a5a6..ee5506bd 100644 --- a/src/canopy_profile_mod.F90 +++ b/src/canopy_profile_mod.F90 @@ -245,13 +245,13 @@ SUBROUTINE CANOPY_ZPD( ZHC, FCLAI, UBZREF, Z0GHC, & REAL(RK), INTENT( IN ) :: ZHC ( : ) ! SUB-CANOPY Total z/h layers (nondimensional) REAL(RK), INTENT( IN ) :: FCLAI ( : ) ! SUB-CANOPY Fractional (z) shapes of the ! plant surface distribution, i.e., a Fractional Culmulative LAI - REAL(RK), INTENT( IN ) :: UBZREF ! Mean wind speed at zref-height of canopy top (m/s) + REAL(RK), INTENT( IN ) :: UBZREF ! Mean wind speed at reference height (ZREF) (m/s) REAL(RK), INTENT( IN ) :: Z0GHC ! Ratio of ground roughness length to canopy top height (nondimensional) REAL(RK), INTENT( IN ) :: LAMBDARS ! Value representing influence of roughness sublayer (nondimensional) REAL(RK), INTENT( IN ) :: CDRAG ! Drag coefficient (nondimensional) REAL(RK), INTENT( IN ) :: PAI ! Total plant/foliage area index (nondimensional) REAL(RK), INTENT( IN ) :: FCH ! Grid cell canopy height (m) - REAL(RK), INTENT( IN ) :: HREF ! Reference Height (m) above the canopy + REAL(RK), INTENT( IN ) :: HREF ! Reference Height above the canopy (ZREF = FCH + HREF; m) REAL(RK), INTENT( IN ) :: Z0_MOD ! Input model value of surface roughness length, z0 (m) INTEGER, INTENT( IN ) :: VTYPE ! Grid cell dominant vegetation type INTEGER, INTENT( IN ) :: LU_OPT ! integer for LU type from model mapped to Massman et al. (default = 0/VIIRS) diff --git a/src/canopy_rad_mod.F90 b/src/canopy_rad_mod.F90 new file mode 100644 index 00000000..163f7a5c --- /dev/null +++ b/src/canopy_rad_mod.F90 @@ -0,0 +1,170 @@ +module canopy_rad_mod + + implicit none + +contains + +!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + SUBROUTINE CANOPY_FSUN_CLU( FCLAI, LAI, CLU, COSZEN, FSUN) + +!----------------------------------------------------------------------- + +! Description: +! computes linear interpolation method for PPFD sun/shade in canopy. + +! Preconditions: +! in-canopy height, and model LAI, clumping index, and solar zenith angle + +! Subroutines and Functions Called: + +! Revision History: +! Prototype 06/23 by PCC +! Jun 2023 P.C. Campbell: Initial standalone PPFD linear subroutine based on +! Silva et al. (2020) exponential curve algorithms +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- + use canopy_const_mod, ONLY: RK !constants for canopy models + use canopy_phot_mod + +! Arguments: +! IN/OUT + REAL(RK), INTENT( IN ) :: FCLAI(:) ! Input Fractional (z) shapes of the + ! plant surface distribution (nondimensional), i.e., a Fractional Culmulative LAI + REAL(RK), INTENT( IN ) :: LAI ! Model input total Leaf Area Index + REAL(RK), INTENT( IN ) :: CLU ! Model input Clumping Index + REAL(RK), INTENT( IN ) :: COSZEN ! Model input Cosine Solar Zenith Angle + REAL(RK), INTENT( OUT ) :: FSUN(SIZE(FCLAI)) ! Sunlit/Shaded fraction from photolysis correction factor + +!Calculate photolyis shading/correction factor through canopy, i.e., the fraction of sunlit leaves downward through canopy +! `canopy_phot` gives relative direct beam irradiance, +! which, multiplied by clumping index, gives sunlit fraction (e.g., Bonan 2019, eq. 14.18) + + call canopy_phot(FCLAI, LAI, CLU, COSZEN, FSUN) + FSUN = FSUN * CLU + + END SUBROUTINE CANOPY_FSUN_CLU + +!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + SUBROUTINE CANOPY_PPFD_EXP( ZK, FCH, SFCRAD, LAI, FSUN, & + PPFD_SUN, PPFD_SHADE, PPFD_AVE) + +!----------------------------------------------------------------------- + +! Description: +! computes linear interpolation method for PPFD sun/shade in canopy. + +! Preconditions: +! in-canopy height, and model LAI, clumping index, and solar zenith angle + +! Subroutines and Functions Called: + +! Revision History: +! Prototype 06/23 by PCC +! Jun 2023 P.C. Campbell: Initial standalone PPFD linear subroutine based on +! Silva et al. (2020) exponential curve algorithms +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- + use canopy_const_mod, ONLY: RK !constants for canopy models + use canopy_utils_mod, ONLY: interp_linear1_internal + +! Arguments: +! IN/OUT + REAL(RK), INTENT( IN ) :: ZK(:) ! Input model heights (m) + REAL(RK), INTENT( IN ) :: FCH ! Model input canopy height (m) + REAL(RK), INTENT( IN ) :: SFCRAD ! Model input Instantaneous surface downward shortwave flux (W/m2) + REAL(RK), INTENT( IN ) :: LAI ! Model input total Leaf Area Index + REAL(RK), INTENT( IN ) :: FSUN(:) ! Sunlit/Shaded fraction from photolysis correction factor + REAL(RK), INTENT( OUT ) :: PPFD_SUN(SIZE(ZK)) ! PPFD for sunlit leaves (umol phot/m2 s) + REAL(RK), INTENT( OUT ) :: PPFD_SHADE(SIZE(ZK)) ! PPFD for shaded leaves (umol phot/m2 s) + REAL(RK), INTENT( OUT ) :: PPFD_AVE(SIZE(ZK)) ! Average PPFD for sunlit and shaded leaves (umol phot/m2 s) + +! LOCAL + REAL(RK), PARAMETER :: CTEMP_1_SUN = 1.083_rk !Exponential 2-m PPFD --> PPFD parameters (Level 1 = + !top of canopy + REAL(RK), PARAMETER :: CTEMP_2_SUN = 1.096_rk !Based on Table 1 in Silva et al. (2022) + REAL(RK), PARAMETER :: CTEMP_3_SUN = 1.104_rk ! + REAL(RK), PARAMETER :: CTEMP_4_SUN = 1.098_rk ! + REAL(RK), PARAMETER :: CTEMP_5_SUN = 1.090_rk !... + REAL(RK), PARAMETER :: DTEMP_1_SUN = 0.002_rk !... + REAL(RK), PARAMETER :: DTEMP_2_SUN = -0.128_rk !... + REAL(RK), PARAMETER :: DTEMP_3_SUN = -0.298_rk !... + REAL(RK), PARAMETER :: DTEMP_4_SUN = -0.445_rk !... + REAL(RK), PARAMETER :: DTEMP_5_SUN = -0.535_rk !... + REAL(RK), PARAMETER :: CTEMP_1_SHADE = 0.871_rk !... + REAL(RK), PARAMETER :: CTEMP_2_SHADE = 0.890_rk !... + REAL(RK), PARAMETER :: CTEMP_3_SHADE = 0.916_rk !... + REAL(RK), PARAMETER :: CTEMP_4_SHADE = 0.941_rk !... + REAL(RK), PARAMETER :: CTEMP_5_SHADE = 0.956_rk !... + REAL(RK), PARAMETER :: DTEMP_1_SHADE = 0.015_rk !... + REAL(RK), PARAMETER :: DTEMP_2_SHADE = -0.141_rk !... + REAL(RK), PARAMETER :: DTEMP_3_SHADE = -0.368_rk !... + REAL(RK), PARAMETER :: DTEMP_4_SHADE = -0.592_rk !... + REAL(RK), PARAMETER :: DTEMP_5_SHADE = -0.743_rk !... + + REAL(RK), PARAMETER :: FRAC_PAR = 0.5_rk !Fraction of incoming solar irradiance that is PAR + + REAL(RK) :: CTEMP_SUN(SIZE(ZK)) ! Regression coefficient C for sun leaves + REAL(RK) :: DTEMP_SUN(SIZE(ZK)) ! Regression coefficient D for sun leaves + REAL(RK) :: CTEMP_SHADE(SIZE(ZK)) ! Regression coefficient C for shade leaves + REAL(RK) :: DTEMP_SHADE(SIZE(ZK)) ! Regression coefficient D for shade leaves + + integer i + +! Use exponential PPFD model based on Silva et al. (2020) to get approx. sun/shade PPFD +! through canopy +!Citation: +!Silva, S. J., Heald, C. L., and Guenther, A. B.: Development of a reduced-complexity plant canopy +!physics surrogate model for use in chemical transport models: a case study with GEOS-Chem v12.3.0, +!Geosci. Model Dev., 13, 2569–2585, https://doi.org/10.5194/gmd-13-2569-2020, 2020. + do i=1, SIZE(ZK) !calculate linear change in parameters interpolated to Silva et al. 5 layer canopy regions + if (ZK(i) .gt. FCH) then ! above canopy, PPFD_leaf = PPFD_toc (toc=top of canopy) + CTEMP_SUN(i) = 0.0 + DTEMP_SUN(i) = 0.0 + CTEMP_SHADE(i) = 0.0 + DTEMP_SHADE(i) = 0.0 + else if (ZK(i) .le. FCH .and. ZK(i) .gt. FCH*(4.0_rk/5.0_rk)) then !Level 1 - 2 + CTEMP_SUN(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & + (/ CTEMP_2_SUN,CTEMP_1_SUN /),ZK(i)) + DTEMP_SUN(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & + (/ DTEMP_2_SUN,DTEMP_1_SUN /),ZK(i)) + CTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & + (/ CTEMP_2_SHADE,CTEMP_1_SHADE /),ZK(i)) + DTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & + (/ DTEMP_2_SHADE,DTEMP_1_SHADE /),ZK(i)) + else if (ZK(i) .le. FCH*(4.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(3.0_rk/5.0_rk)) then !Level 2 - 3 + CTEMP_SUN(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & + (/ CTEMP_3_SUN,CTEMP_2_SUN /),ZK(i)) + DTEMP_SUN(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & + (/ DTEMP_3_SUN,DTEMP_2_SUN /),ZK(i)) + CTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & + (/ CTEMP_3_SHADE,CTEMP_2_SHADE /),ZK(i)) + DTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & + (/ DTEMP_3_SHADE,DTEMP_2_SHADE /),ZK(i)) + else if (ZK(i) .le. FCH*(3.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(2.0_rk/5.0_rk)) then !Level 3 - 4 + CTEMP_SUN(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & + (/ CTEMP_4_SUN,CTEMP_3_SUN /),ZK(i)) + DTEMP_SUN(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & + (/ DTEMP_4_SUN,DTEMP_3_SUN /),ZK(i)) + CTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & + (/ CTEMP_4_SHADE,CTEMP_3_SHADE /),ZK(i)) + DTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & + (/ DTEMP_4_SHADE,DTEMP_3_SHADE /),ZK(i)) + else if (ZK(i) .le. FCH*(2.0_rk/5.0_rk) ) then !Level 4 - Bottom + CTEMP_SUN(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & + (/ CTEMP_5_SUN,CTEMP_4_SUN /),ZK(i)) + DTEMP_SUN(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & + (/ DTEMP_5_SUN,DTEMP_4_SUN /),ZK(i)) + CTEMP_SHADE(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & + (/ CTEMP_5_SHADE,CTEMP_4_SHADE /),ZK(i)) + DTEMP_SHADE(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & + (/ DTEMP_5_SHADE,DTEMP_4_SHADE /),ZK(i)) + end if + end do + + PPFD_SUN = FRAC_PAR * SFCRAD * EXP(CTEMP_SUN + DTEMP_SUN * LAI) !Silva et al. W/m2 --> umol m-2 s-1 + PPFD_SHADE = FRAC_PAR * SFCRAD * EXP(CTEMP_SHADE + DTEMP_SHADE * LAI) + PPFD_AVE = (PPFD_SUN*FSUN) + (PPFD_SHADE*(1.0-FSUN)) ! average = sum sun and shade weighted by sunlit fraction + + END SUBROUTINE CANOPY_PPFD_EXP + +end module canopy_rad_mod diff --git a/src/canopy_readnml.F90 b/src/canopy_readnml.F90 index a9eb4f4a..52e63989 100644 --- a/src/canopy_readnml.F90 +++ b/src/canopy_readnml.F90 @@ -23,7 +23,8 @@ SUBROUTINE canopy_readnml NAMELIST /userdefs/ infmt_opt, nlat, nlon, modlays, modres, href_opt, href_set, & z0ghc, lambdars, flameh_opt, flameh_set, frp_fac, ifcanwind, ifcanwaf, & ifcaneddy, ifcanphot, ifcanbio, pai_opt, pai_set, lu_opt, z0_opt, dx_opt, & - dx_set, lai_thresh, frt_thresh, fch_thresh, rsl_opt, bio_cce + dx_set, lai_thresh, frt_thresh, fch_thresh, rsl_opt, bio_cce, biovert_opt, & + ssg_opt, ssg_set, crop_opt, crop_set, co2_opt, co2_set !------------------------------------------------------------------------------- ! Error, warning, and informational messages. @@ -201,6 +202,41 @@ SUBROUTINE canopy_readnml bio_cce = 0.21_rk !------------------------------------------------------------------------------- +!------------------------------------------------------------------------------- +! Set default integer value for MEGAN vertical integration of emissions (0, off full leaf-level emissions) + biovert_opt = 0 +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default integer for shrubs/savanaa/grasslands vegtype option from GEDI or user (default = 0) + ssg_opt = 0 +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default value for shrubs/savanaa/grasslands vegtype heights used in model (m) (Default = 1 m) + ssg_set = 1.0_rk +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default integer for crop vegtype option from GEDI or user (default = 0) + crop_opt = 0 +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default value for crop vegtype heights used in model (m) (Default = 3 m) + crop_set = 3.0_rk +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default integer for co2 inhibition option for biogenic isoprene emissions (default = 0; Possell & Hewitt (2011)) + co2_opt = 0 +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default value for atmospheric co2 concentration for co2_opt (m) (Default = 400.0 ppmv) + co2_set = 400.0_rk +!------------------------------------------------------------------------------- + !------------------------------------------------------------------------------- ! Read namelist to get user definitions. Rewind namelist file after each ! read in case namelists are not in the correct order in the namelist. diff --git a/src/canopy_tleaf_mod.F90 b/src/canopy_tleaf_mod.F90 new file mode 100644 index 00000000..aac93106 --- /dev/null +++ b/src/canopy_tleaf_mod.F90 @@ -0,0 +1,126 @@ +module canopy_tleaf_mod + + implicit none + +contains + +!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + SUBROUTINE CANOPY_TLEAF_LIN( ZK, FCH, TEMP2, FSUN, & + TLEAF_SUN, TLEAF_SHADE, TLEAF_AVE) + +!----------------------------------------------------------------------- + +! Description: +! computes linear interpolation method for tleaf sun/shade in canopy. + +! Preconditions: +! canopy height,in-canopy model height, 2-m temp + +! Subroutines and Functions Called: + +! Revision History: +! Prototype 06/23 by PCC +! Jun 2023 P.C. Campbell: Initial standalone tleaf linear subroutine based on +! Silva et al. (2020) algorithms +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- + use canopy_const_mod, ONLY: RK !constants for canopy models + use canopy_utils_mod, ONLY: interp_linear1_internal + +! Arguments: +! IN/OUT + REAL(RK), INTENT( IN ) :: ZK(:) ! Input model heights (m) + REAL(RK), INTENT( IN ) :: FCH ! Model input canopy height (m) + REAL(RK), INTENT( IN ) :: TEMP2 ! Model input 2-m Temperature (K + REAL(RK), INTENT( IN ) :: FSUN(:) ! Sunlit/Shaded fraction from photolysis correction factor + REAL(RK), INTENT( OUT ) :: TLEAF_SUN(SIZE(ZK)) ! Leaf temp for sunlit leaves (K) + REAL(RK), INTENT( OUT ) :: TLEAF_SHADE(SIZE(ZK)) ! Leaf temp for shaded leaves (K) + REAL(RK), INTENT( OUT ) :: TLEAF_AVE(SIZE(ZK)) ! Ave Leaf temp for sun/shaded leaves (K) + +! LOCAL + !Linearized 2-m temp --> leaf temp parameters Based on Table 1 in Silva et al. (2020) + REAL(RK), PARAMETER :: ATEMP_1_SUN = -13.891_rk !Level 1 = top of canopy + REAL(RK), PARAMETER :: ATEMP_2_SUN = -12.322_rk !... + REAL(RK), PARAMETER :: ATEMP_3_SUN = -1.032_rk !... + REAL(RK), PARAMETER :: ATEMP_4_SUN = -5.172_rk !... + REAL(RK), PARAMETER :: ATEMP_5_SUN = -5.589_rk !... + REAL(RK), PARAMETER :: BTEMP_1_SUN = 1.064_rk !... + REAL(RK), PARAMETER :: BTEMP_2_SUN = 1.057_rk !... + REAL(RK), PARAMETER :: BTEMP_3_SUN = 1.031_rk !... + REAL(RK), PARAMETER :: BTEMP_4_SUN = 1.050_rk !... + REAL(RK), PARAMETER :: BTEMP_5_SUN = 1.051_rk !... + REAL(RK), PARAMETER :: ATEMP_1_SHADE = -12.846_rk !... + REAL(RK), PARAMETER :: ATEMP_2_SHADE = -11.343_rk !... + REAL(RK), PARAMETER :: ATEMP_3_SHADE = -1.068_rk !... + REAL(RK), PARAMETER :: ATEMP_4_SHADE = -5.551_rk !... + REAL(RK), PARAMETER :: ATEMP_5_SHADE = -5.955_rk !... + REAL(RK), PARAMETER :: BTEMP_1_SHADE = 1.060_rk !... + REAL(RK), PARAMETER :: BTEMP_2_SHADE = 1.053_rk !... + REAL(RK), PARAMETER :: BTEMP_3_SHADE = 1.031_rk !... + REAL(RK), PARAMETER :: BTEMP_4_SHADE = 1.051_rk !... + REAL(RK), PARAMETER :: BTEMP_5_SHADE = 1.053_rk !... + REAL(RK) :: ATEMP_SUN(SIZE(ZK)) ! Regression coefficient A for sun leaves (Silva et al., 2020) + REAL(RK) :: BTEMP_SUN(SIZE(ZK)) ! Regression coefficient B for sun leaves + REAL(RK) :: ATEMP_SHADE(SIZE(ZK)) ! Regression coefficient A for shade leaves + REAL(RK) :: BTEMP_SHADE(SIZE(ZK)) ! Regression coefficient B for shade leaves + + integer i + +! Use linear canopy temperature model based on Silva et al. (2020) to get approx. sun/shade leaf temperatures +! through canopy (ignores effect of wind speed on leaf boundary layer ~ 1 % error/bias) +!Citation: +!Silva, S. J., Heald, C. L., and Guenther, A. B.: Development of a reduced-complexity plant canopy +!physics surrogate model for use in chemical transport models: a case study with GEOS-Chem v12.3.0, +!Geosci. Model Dev., 13, 2569–2585, https://doi.org/10.5194/gmd-13-2569-2020, 2020. + do i=1, SIZE(ZK) !calculate linear change in parameters interpolated to Silva et al. 5 layer canopy regions + if (ZK(i) .gt. FCH) then ! above canopy, Tleaf = Tair + ATEMP_SUN(i) = 0.0 + BTEMP_SUN(i) = 1.0 + ATEMP_SHADE(i) = 0.0 + BTEMP_SHADE(i) = 1.0 + else if (ZK(i) .le. FCH .and. ZK(i) .gt. FCH*(4.0_rk/5.0_rk)) then !Level 1 - 2 + ATEMP_SUN(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & + (/ ATEMP_2_SUN,ATEMP_1_SUN /),ZK(i)) + BTEMP_SUN(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & + (/ BTEMP_2_SUN,BTEMP_1_SUN /),ZK(i)) + ATEMP_SHADE(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & + (/ ATEMP_2_SHADE,ATEMP_1_SHADE /),ZK(i)) + BTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), & + (/ BTEMP_2_SHADE,BTEMP_1_SHADE /),ZK(i)) + else if (ZK(i) .le. FCH*(4.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(3.0_rk/5.0_rk)) then !Level 2 - 3 + ATEMP_SUN(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & + (/ ATEMP_3_SUN,ATEMP_2_SUN /),ZK(i)) + BTEMP_SUN(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & + (/ BTEMP_3_SUN,BTEMP_2_SUN /),ZK(i)) + ATEMP_SHADE(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & + (/ ATEMP_3_SHADE,ATEMP_2_SHADE /),ZK(i)) + BTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), & + (/ BTEMP_3_SHADE,BTEMP_2_SHADE /),ZK(i)) + else if (ZK(i) .le. FCH*(3.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(2.0_rk/5.0_rk)) then !Level 3 - 4 + ATEMP_SUN(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & + (/ ATEMP_4_SUN,ATEMP_3_SUN /),ZK(i)) + BTEMP_SUN(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & + (/ BTEMP_4_SUN,BTEMP_3_SUN /),ZK(i)) + ATEMP_SHADE(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & + (/ ATEMP_4_SHADE,ATEMP_3_SHADE /),ZK(i)) + BTEMP_SHADE(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), & + (/ BTEMP_4_SHADE,BTEMP_3_SHADE /),ZK(i)) + else if (ZK(i) .le. FCH*(2.0_rk/5.0_rk) ) then !Level 4 - Bottom + ATEMP_SUN(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & + (/ ATEMP_5_SUN,ATEMP_4_SUN /),ZK(i)) + BTEMP_SUN(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & + (/ BTEMP_5_SUN,BTEMP_4_SUN /),ZK(i)) + ATEMP_SHADE(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & + (/ ATEMP_5_SHADE,ATEMP_4_SHADE /),ZK(i)) + BTEMP_SHADE(i) = interp_linear1_internal((/ ZK(1),FCH*(2.0_rk/5.0_rk) /), & + (/ BTEMP_5_SHADE,BTEMP_4_SHADE /),ZK(i)) + end if + end do + + TLEAF_SUN = ATEMP_SUN + (BTEMP_SUN*TEMP2) + TLEAF_SHADE = ATEMP_SHADE + (BTEMP_SHADE*TEMP2) + TLEAF_AVE = (TLEAF_SUN*FSUN) + (TLEAF_SHADE*(1.0-FSUN)) ! average = sum sun and shade weighted by sunlit fraction + + END SUBROUTINE CANOPY_TLEAF_LIN + +end module canopy_tleaf_mod diff --git a/src/canopy_utils_mod.F90 b/src/canopy_utils_mod.F90 index 312f07a3..ae1117bb 100644 --- a/src/canopy_utils_mod.F90 +++ b/src/canopy_utils_mod.F90 @@ -6,7 +6,7 @@ module canopy_utils_mod private public IntegrateTrapezoid,interp_linear1_internal,CalcPAI, & - CalcDX,CalcFlameH + CalcDX,CalcFlameH,GET_GAMMA_CO2 contains @@ -125,4 +125,130 @@ function CalcFlameH(frp, dx) end function + real(rk) function GET_GAMMA_CO2(co2_opt, co2_set) result( GAMMA_CO2 ) + ! !IROUTINE: get_gamma_co2 + ! + ! !DESCRIPTION: Function GET\_GAMMA\_CO2 computes the CO2 activity factor + ! associated with CO2 inhibition of isoprene emission. Called from + ! GET\_MEGAN\_EMISSIONS only. + !\\ + !\\ + ! !INTERFACE: + + ! !INPUT PARAMETERS: + integer, INTENT(IN) :: co2_opt ! Option for co2 inhibition calculation + ! 0=Possell & Hewitt (2011); + ! 1=Wilkinson et al. (2009) + ! >1=off + real(rk), INTENT(IN) :: co2_set ! User set atmospheric CO2 conc [ppmv] + ! + ! !RETURN VALUE: +! REAL(rk) :: GAMMA_CO2 ! CO2 activity factor [unitless] + ! + ! !LOCAL VARIABLES: + REAL(rk) :: CO2i ! Intercellular CO2 conc [ppmv] + REAL(rk) :: ISMAXi ! Asymptote for intercellular CO2 + REAL(rk) :: HEXPi ! Exponent for intercellular CO2 + REAL(rk) :: CSTARi ! Scaling coef for intercellular CO2 + REAL(rk) :: ISMAXa ! Asymptote for atmospheric CO2 + REAL(rk) :: HEXPa ! Exponent for atmospheric CO2 + REAL(rk) :: CSTARa ! Scaling coef for atmospheric CO2 + ! + ! !REMARKS: + ! References: + ! ============================================================================ + ! (1 ) Heald, C. L., Wilkinson, M. J., Monson, R. K., Alo, C. A., + ! Wang, G. L., and Guenther, A.: Response of isoprene emission + ! to ambient co(2) changes and implications for global budgets, + ! Global Change Biology, 15, 1127-1140, 2009. + ! (2 ) Wilkinson, M. J., Monson, R. K., Trahan, N., Lee, S., Brown, E., + ! Jackson, R. B., Polley, H. W., Fay, P. A., and Fall, R.: Leaf + ! isoprene emission rate as a function of atmospheric CO2 + ! concentration, Global Change Biology, 15, 1189-1200, 2009. + ! (3 ) Possell, M., and Hewitt, C. N.: Isoprene emissions from plants + ! are mediated by atmospheric co2 concentrations, Global Change + ! Biology, 17, 1595-1610, 2011. + + ! !REVISION HISTORY: + ! (1 ) Implemented in the standard code by A. Tai (Jun 2012). + ! See https://github.com/geoschem/hemco for complete history + !EOP + !------------------------------------------------------------------------------ + !BOC + + !----------------------- + ! Compute GAMMA_CO2 + !----------------------- + + !---------------------------------------------------------- + ! Choose between two alternative CO2 inhibition schemes + !---------------------------------------------------------- + + IF ( co2_opt .eq. 0 ) THEN + + ! Use empirical relationship of Possell & Hewitt (2011): + ! Empirical relationship of Possell & Hewitt (2011) based on nine + ! experimental studies including Wilkinson et al. (2009). + + GAMMA_CO2 = 8.9406_rk / ( 1.0_rk + 8.9406_rk * 0.0024_rk * co2_set ) + + + ELSEIF ( co2_opt .eq. 1 ) THEN + + ! Use parameterization of Wilkinson et al. (2009): + ! Semi-process-based parameterization of Wilkinson et al. (2009), + ! taking into account of sensitivity to intercellular CO2 + ! fluctuation, which is here set as a constant fraction of + ! atmospheric CO2. This is especially recommended for sub-ambient + ! CO2 concentrations:: + + + ! Parameters for intercellular CO2 using linear interpolation: + IF ( co2_set <= 600.0_rk ) THEN + ISMAXi = 1.036_rk - (1.036_rk - 1.072_rk) / & + (600.0_rk - 400.0_rk) * (600.0_rk - co2_set) + HEXPi = 2.0125_rk - (2.0125_rk - 1.7000_rk) / & + (600.0_rk - 400.0_rk) * (600.0_rk - co2_set) + CSTARi = 1150.0_rk - (1150.0_rk - 1218.0_rk) / & + (600.0_rk - 400.0_rk) * (600.0_rk - co2_set) + ELSEIF ( co2_set > 600.0_rk .AND. co2_set < 800.0_rk ) THEN + ISMAXi = 1.046_rk - (1.046_rk - 1.036_rk) / & + (800.0_rk - 600.0_rk) * (800.0_rk - co2_set) + HEXPi = 1.5380_rk - (1.5380_rk - 2.0125_rk) / & + (800.0_rk - 600.0_rk) * (800.0_rk - co2_set) + CSTARi = 2025.0_rk - (2025.0_rk - 1150.0_rk) / & + (800.0_rk - 600.0_rk) * (800.0_rk - co2_set) + ELSE + ISMAXi = 1.014_rk - (1.014_rk - 1.046_rk) / & + (1200.0_rk - 800.0_rk) * (1200.0_rk - co2_set) + HEXPi = 2.8610_rk - (2.8610_rk - 1.5380_rk) / & + (1200.0_rk - 800.0_rk) * (1200.0_rk - co2_set) + CSTARi = 1525.0_rk - (1525.0_rk - 2025.0_rk) / & + (1200.0_rk - 800.0_rk) * (1200.0_rk - co2_set) + ENDIF + + ! Parameters for atmospheric CO2: + ISMAXa = 1.344_rk + HEXPa = 1.4614_rk + CSTARa = 585.0_rk + + ! For now, set CO2_Ci = 0.7d0 * CO2_Ca as recommended by Heald + ! et al. (2009): + CO2i = 0.7_rk * co2_set + + ! Compute GAMMA_CO2: + GAMMA_CO2 = ( ISMAXi - ISMAXi * CO2i**HEXPi / & + ( CSTARi**HEXPi + CO2i**HEXPi ) ) & + * ( ISMAXa - ISMAXa * ( 0.7_rk * co2_set )**HEXPa / & + ( CSTARa**HEXPa + ( 0.7_rk * co2_set )**HEXPa ) ) + + ELSE + + ! No CO2 inhibition scheme is used; GAMMA_CO2 set to unity: + GAMMA_CO2 = 1.0_rk + + ENDIF + + end function GET_GAMMA_CO2 + end module canopy_utils_mod diff --git a/src/canopy_wind_mod.F90 b/src/canopy_wind_mod.F90 index 53f87d51..6b6e846e 100644 --- a/src/canopy_wind_mod.F90 +++ b/src/canopy_wind_mod.F90 @@ -33,12 +33,12 @@ SUBROUTINE CANOPY_WIND_MOST( HCM, ZK, FAFRACK, UBZREF, Z0GHC, & REAL(RK), INTENT( IN ) :: ZK ! Above/Below canopy height, z (m) REAL(RK), INTENT( IN ) :: FAFRACK ! Fractional (z) shapes of the ! plant surface distribution (nondimensional) - REAL(RK), INTENT( IN ) :: UBZREF ! Mean wind speed at reference height (m/s) + REAL(RK), INTENT( IN ) :: UBZREF ! Mean wind speed at reference height (ZREF) (m/s) REAL(RK), INTENT( IN ) :: Z0GHC ! Ratio of ground roughness length to canopy top height (nondimensional) REAL(RK), INTENT( IN ) :: LAMBDARS ! Value representing influence of roughness sublayer (nondimensional) REAL(RK), INTENT( IN ) :: CDRAG ! Drag coefficient (nondimensional) REAL(RK), INTENT( IN ) :: PAI ! Total plant/foliage area index (nondimensional) - REAL(RK), INTENT( IN ) :: HREF ! Reference Height above canopy associated with ref wind speed (m) + REAL(RK), INTENT( IN ) :: HREF ! Reference Height above the canopy (ZREF = HCM + HREF; m) REAL(RK), INTENT( IN ) :: D_H ! Zero plane displacement heights (nondimensional) REAL(RK), INTENT( IN ) :: ZO_H ! Surface (soil+veg) roughness lengths (nondimensional) REAL(RK), INTENT( OUT ) :: CANBOT_OUT ! Canopy bottom wind reduction factor = canbot (nondimensional)