Merge pull request #92 from quaz115/develop

Leaf Age Response function for all Biogenic VOCs consistent with MEGAN

README.md

@@ -224,6 +224,8 @@ You can also [generate global inputs using Python (see python/global_data_proces
 | `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`) |
+| `leafage_opt`   | user-set options for applying leaf-age response to biogenic VOC emissions based on [Guenther et al. 2006](https://doi.org/10.5194/acp-6-3181-2006) (default is off i.e., `leafage_opt=1`, the corresponding $\gamma$ is set to 1). If turned on (`leafage_opt=0`), leafage $\gamma$ is calculated and the lai_tstep option needs to be set to ensure correct interpolation in this leafage_opt calculation. |
+| `lai_tstep`     | user-defined options for the number of seconds in the interval at which LAI (Leaf Area Index) data is provided to the model. For instance, if LAI data is given on a daily basis, lai_tstep would be set to the number of seconds in a day (86,400 seconds). If LAI data is provided monthly, then lai_tstep would represent the total number of seconds in that month (e.g., 2,592,000 seconds for a 30-day month). This parameter helps in determining the frequency of LAI input and is crucial for interpolating LAI values to the model's hourly timesteps when the model's timestep (time_intvl) is smaller than the LAI input interval.  |
 | `lai_thresh`    | user-set real value of LAI threshold for contiguous canopy (m2/m2)                 |
 | `frt_thresh`    | user-set real value of forest fraction threshold for contiguous canopy             |
 | `fch_thresh`    | user-set real value of canopy height  threshold for contiguous canopy (m)          |

input/namelist.canopy

@@ -49,6 +49,8 @@
   crop_set    =  3.0
   co2_opt     =  0
   co2_set     =  400.0
+  lai_tstep   =  86400
+  leafage_opt =  1
   lai_thresh  =  0.1
   frt_thresh  =  0.1
   fch_thresh  =  0.5

python/examples.ipynb

@@ -882,6 +882,63 @@
    "source": [
     "ds.emi_apin.mean(\"z\", keep_attrs=True).plot(size=4, marker=\"o\")"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0a445f3c-52a2-4a3f-89be-c06011b91a76",
+   "metadata": {},
+   "source": [
+    "### Leaf age parameterization"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a0eee654-2333-408e-a4f6-0bb9ed574600",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "config = {\n",
+    "    \"filenames\": {\n",
+    "        \"file_vars\": [\n",
+    "            \"../input/point_file_20220630.sfcf023.txt\",\n",
+    "            \"../input/point_file_20220701.sfcf000.txt\",\n",
+    "            \"../input/point_file_20220701.sfcf001.txt\",\n",
+    "        ],\n",
+    "    },\n",
+    "    \"userdefs\": {\n",
+    "        \"infmt_opt\": 1,\n",
+    "        \"ntime\": 3,\n",
+    "        \"nlat\": 1,\n",
+    "        \"nlon\": 1,\n",
+    "        \"leafage_opt\": 0,  # on\n",
+    "    },\n",
+    "}\n",
+    "ds0 = run(config=config)\n",
+    "\n",
+    "config[\"userdefs\"][\"leafage_opt\"] = 1  # off\n",
+    "ds1 = run(config=config)\n",
+    "\n",
+    "ds = xr.concat([ds0, ds1], dim=\"leafage_opt\")\n",
+    "ds"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2439516f-bd60-495f-bbc5-5dfe70555c03",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(5, 6), sharex=True)\n",
+    "\n",
+    "ds.emi_apin.mean(\"z\", keep_attrs=True).plot(hue=\"leafage_opt\", marker=\"o\", ax=ax1)\n",
+    "ds.emi_isop.mean(\"z\", keep_attrs=True).plot(hue=\"leafage_opt\", marker=\"o\", ax=ax2)"
+   ]
   }
  ],
  "metadata": {
@@ -900,7 +957,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.10"
+   "version": "3.11.5"
   }
  },
  "nbformat": 4,

python/global_data_process.py

@@ -1,5 +1,6 @@
 """
-Created on Sun Jun 4 15:45:29 2023
+Created on Sun Jun 4 2023
+Updated on Mon Oct 16 2023: Use daily gfs.canopy files
 
 Author: Wei-Ting Hung
 """
@@ -14,7 +15,7 @@
 from scipy.interpolate import griddata
 
 """User Options"""
-path = "/scratch/pcampbe8/canopy-app/input"  # work directory
+path = "/scratch1/BMC/rcm2/qrasool/canopy-app-LeafAgeTimeFeature/input"  # "/scratch/pcampbe8/canopy-app/input"  # work directory
 year = 2022  # year
 month = 7  # month
 day = 1  # day
@@ -65,7 +66,7 @@
 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_can = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + DD + ".sfcf000.nc"  # canopy file
 f_output = (
     path + "/gfs.canopy.t" + HHI + "z." + YY + MM + DD + ".global.f0" + HH + ".nc"
 )  # output file
@@ -88,7 +89,7 @@
             + ".nc "
         )
 elif frp_src == 2:  # climatological frp
-    f_frp = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + "01.sfcf000.nc"
+    f_frp = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + DD + ".sfcf000.nc"
 
 
 # required variables
@@ -109,7 +110,7 @@
     "hpbl",
     "prate_ave",
 ]
-canlist = ["lai", "clu", "ffrac", "fh", "mol", "csz", "frp", "href"]
+canlist = ["lai", "clu", "ffrac", "fh", "pavd", "mol", "csz", "frp", "href"]
 
 
 # constants
@@ -149,12 +150,35 @@ def mapping(xgrid, ygrid, data, xdata, ydata, map_method, fvalue):
 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, :, :])
+    if varname == "pavd":
+        # map to met grids
+        yt = readin["lat"][:]
+        xt = readin["lon"][:]
+        data = np.squeeze(readin[varname][:])
+
+        DATA = np.empty([data.shape[0], lat.shape[0], lat.shape[1]])
+
+        for ll in np.arange(data.shape[0]):
+            DATA[ll, :, :] = mapping(
+                lat,
+                lon,
+                data[ll, :, :].flatten(),
+                yt.flatten(),
+                xt.flatten(),
+                "linear",
+                np.nan,
+            )
+
+    else:
+        # 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 = mapping(lat, lon, data.flatten(), yt.flatten(), xt.flatten(), "linear", np.nan)
     DATA[np.isnan(DATA)] = 0
     DATA[DATA < 0] = 0
     return DATA
@@ -226,7 +250,8 @@ def read_frp_local(filename, lat, lon, fill_value):
             + "z."
             + YY
             + MM
-            + "01.sfcf000.nc",
+            + DD
+            + ".sfcf000.nc",
         ]
     )
     if os.path.isfile(f_can) is True:
@@ -262,7 +287,7 @@ def read_frp_local(filename, lat, lon, fill_value):
     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"
+        f_frp = path + "/gfs.canopy.t" + HHI + "z." + YY + MM + DD + ".sfcf000.nc"
 
 if frp_src == 2:  # 12 month climatology frp
     if os.path.isfile(f_frp) is True:
@@ -281,7 +306,8 @@ def read_frp_local(filename, lat, lon, fill_value):
                 + "z."
                 + YY
                 + MM
-                + "01.sfcf000.nc",
+                + DD
+                + ".sfcf000.nc",
             ]
         )
         if os.path.isfile(f_met) is True:
@@ -354,6 +380,11 @@ def read_frp_local(filename, lat, lon, fill_value):
         ATT = ["Canopy height above the surface", "m", fill_value]
         DATA = read_gfs_climatology(f_can, lat, lon, "fh")
 
+    elif varname == "pavd":
+        ATTNAME = ["long_name", "units", "missing_value"]
+        ATT = ["Plant area volume density profile", "m2/m3", fill_value]
+        DATA = read_gfs_climatology(f_can, lat, lon, "pavd")
+
     elif varname == "mol":
         # Reference:
         # Essa 1999, ESTIMATION OF MONIN-OBUKHOV LENGTH USING RICHARDSON AND BULK RICHARDSON NUMBER
@@ -409,15 +440,48 @@ def read_frp_local(filename, lat, lon, fill_value):
     print(ATT)
 
     # adding to output file
-    output = Dataset(f_output, "a")
+    if varname == "pavd":
+        output = Dataset(f_output, "a")
+        output.createDimension("level", DATA.shape[0])
+
+        var_bot = output.createVariable("layer_bottom", "i4", ("level",))
+        var_top = output.createVariable("layer_top", "i4", ("level",))
+        var = output.createVariable(
+            varname,
+            "float",
+            ("time", "level", "grid_yt", "grid_xt"),
+            fill_value=fill_value,
+        )
 
-    var = output.createVariable(
-        varname, "float", ("time", "grid_yt", "grid_xt"), fill_value=fill_value
-    )
-    write_varatt(var, ATTNAME, ATT)
-    var[:] = DATA
+        write_varatt(var, ATTNAME, ATT)
+        write_varatt(
+            var_bot,
+            ["long_name", "units"],
+            ["height of the layer bottom above the ground", "m"],
+        )
+        write_varatt(
+            var_top,
+            ["long_name", "units"],
+            ["height of the layer top above the ground", "m"],
+        )
+
+        var_bot[:] = np.arange(0, 65 + 1, 5)
+        var_top[:] = np.arange(5, 70 + 1, 5)
+        var[:] = DATA
+
+        output.close()
+        del [var_bot, var_top]
+
+    else:
+        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()
+        output.close()
 
     print("---- " + varname + " complete!")
 

src/canopy_app.F90

@@ -85,7 +85,7 @@ program canopy_app
 ! Main canopy model calculations.
 !-------------------------------------------------------------------------------
 
-        call canopy_calcs
+        call canopy_calcs(nn)
 
 !-------------------------------------------------------------------------------
 ! Write model output of canopy model calculations.