RRFS-SD Reader

[bbakernoaa]

It came to my attention that the RRFS-SD verification has been using the _rrfs_cmaq_mm.py reader. While this works it is probably better to create a separate one that is a simpler use case than the complex chemistry of cmaq itself.

Specially there are three variables for air composition within the files.

smoke (pm25)
dust (pm25)
coarse_pm (PM10 coarse mode only)

We need to aggregate these into PM25 and PM10

PM25 = smoke + dust

PM10 = smoke + dust + PM10

@jordanschnell can you confirm the variable names here

[zmoon]

Comparing current develop to what @jordanschnell was using:

--- monetio/models/_rrfs_cmaq_mm.py     2023-12-13 19:10:27.241244093 +0000
+++ /scratch1/BMC/wrf-chem/Jordan/miniconda3/envs/melodies-monet/lib/python3.9/site-packages/monetio/models/_rrfs_
cmaq_mm.py      2023-08-05 01:17:24.000000000 +0000
@@ -19,7 +19,8 @@
     var_list=None,
     fname_pm25=None,
     surf_only=False,
-    **kwargs,
+    convert_pm25=True,
+    **kwargs
 ):
     # Like WRF-chem add var list that just determines whether to calculate sums or not to speed this up.
     """Method to open RFFS-CMAQ dyn* netcdf files.
@@ -107,7 +108,7 @@
         var_list.append("tmp")
         var_list.append("pressfc")
         var_list.append("dpres")
-        var_list.append("hgtsfc")
+#        var_list.append("hgtsfc")
         var_list.append("delz")

         # Remove duplicates just in case:
@@ -119,7 +120,7 @@
         # If variables in pm25 files are included remove these as these are not in the main file
         # And will be added later.
         for pm25_var in [
-            "PM25_TOT",
+#            "PM25_TOT",
             "PM25_TOT_NSOM",
             "PM25_EC",
             "PM25_NH4",
@@ -153,8 +154,6 @@
         ]

     if fname_pm25 is not None:
-        from ..util import _try_merge_exact
-
         # Add the processed pm2.5 species.
         dset_pm25 = xr.open_mfdataset(fname_pm25, concat_dim="time", combine="nested", **kwargs)
         dset_pm25 = dset_pm25.drop(
@@ -164,7 +163,7 @@
         # same pressure levels from the model dynf* files.
         # Attributes are formatted differently in pm25 file so remove attributes and use those from dynf* files.
         dset_pm25.attrs = {}
-        dset = _try_merge_exact(dset, dset_pm25, right_name="PM2.5")
+        dset = dset.merge(dset_pm25)

     # Standardize some variable names
     dset = dset.rename(
@@ -178,26 +177,23 @@
             "tmp": "temperature_k",  # standard temperature (kelvin)
             "pressfc": "surfpres_pa",
             "dpres": "dp_pa",  # Change names so standard surfpres_pa and dp_pa
-            "hgtsfc": "surfalt_m",
+#            "hgtsfc": "surfalt_m",
             "delz": "dz_m",
         }
     )  # Optional, but when available include altitude info

     # Calculate pressure. This has to go before sorting because ak and bk
     # are not sorted as they are in attributes
-    dset["pres_pa_mid"] = _calc_pressure(dset)
+    dset["pres_pa_mid"] = _calc_pressure(dset,surf_only)

     # Adjust pressure levels for all models such that the surface is first.
-    if np.all(np.diff(dset.z.values) > 0):  # increasing pressure
-        dset = dset.isel(z=slice(None, None, -1))  # -> decreasing
-    if np.all(np.diff(dset.z_i.values) > 0):  # increasing pressure
-        dset = dset.isel(z_i=slice(None, None, -1))  # -> decreasing
-    dset["dz_m"] = dset["dz_m"] * -1.0  # Change to positive values.
+    dset = dset.sortby("z", ascending=False)
+    dset = dset.sortby("z_i", ascending=False)

     # Note this altitude calcs needs to always go after resorting.
     # Altitude calculations are all optional, but for each model add values that are easy to calculate.
-    if not surf_only:
-        dset["alt_msl_m_full"] = _calc_hgt(dset)
+#    dset["alt_msl_m_full"] = _calc_hgt(dset)
+#    dset["dz_m"] = dset["dz_m"] * -1.0  # Change to positive values.

     # Set coordinates
     dset = dset.reset_index(
@@ -226,9 +222,10 @@
     for i in dset.variables:
         if "units" in dset[i].attrs:
             if "ug/kg" in dset[i].attrs["units"]:
-                # ug/kg -> ug/m3 using dry air density
-                dset[i] = dset[i] * dset["pres_pa_mid"] / dset["temperature_k"] / 287.05535
-                dset[i].attrs["units"] = r"$\mu g m^{-3}$"
+                if convert_pm25:
+                    # ug/kg -> ug/m3 using dry air density
+                    dset[i] = dset[i] * dset["pres_pa_mid"] / dset["temperature_k"] / 287.05535
+                    dset[i].attrs["units"] = r"$\mu g m^{-3}$"

     # add lazy diagnostic variables
     # Note that because there are so many species to sum. Summing the aerosols is slowing down the code.
@@ -259,7 +256,7 @@
     if "pm25_om" in list_calc_sum:
         dset = add_lazy_om_pm25(dset, dict_sum)
     # Change the times to pandas format
-    dset["time"] = dset.indexes["time"].to_datetimeindex(unsafe=True)
+# JLS   dset["time"] = dset.indexes["time"].to_datetimeindex(unsafe=True)
     # Turn off warning for now. This is just because the model is in julian time

     # Drop extra variables that were part of sum, but are not in original var_list
@@ -1074,7 +1071,7 @@
     return z


-def _calc_pressure(dset):
+def _calc_pressure(dset,surf_only):
     """Calculate the mid-layer pressure in Pa from surface pressure
     and ak and bk constants.

@@ -1094,14 +1091,19 @@
     xarray.DataArray
         Mid-layer pressure with attributes.
     """
-    p = dset.dp_pa.copy().load()  # Have to load into memory here so can assign levels.
-    psfc = dset.surfpres_pa.copy().load()
+    pres = dset.dp_pa.copy().load()  # Have to load into memory here so can assign levels.
+    srfpres = dset.surfpres_pa.copy().load()
     for k in range(len(dset.z)):
-        pres_2 = dset.ak[k + 1] + psfc * dset.bk[k + 1]
-        pres_1 = dset.ak[k] + psfc * dset.bk[k]
-        p[:, k, :, :] = (pres_2 - pres_1) / np.log(pres_2 / pres_1)
-
-    p.name = "pres_pa_mid"
-    p.attrs["units"] = "pa"
-    p.attrs["long_name"] = "Pressure Mid Layer in Pa"
-    return p
+        if surf_only:
+           pres_2 = 0.0 + srfpres * 0.9978736
+           pres_1 = 0.0 + srfpres * 1.0
+        else:
+           pres_2 = dset.ak[k + 1] + srfpres * dset.bk[k + 1]
+           pres_1 = dset.ak[k] + srfpres * dset.bk[k]
+
+        pres[:, k, :, :] = (pres_2 - pres_1) / np.log(pres_2 / pres_1)
+
+    pres.name = "pres_pa_mid"
+    pres.attrs["units"] = "pa"
+    pres.attrs["long_name"] = "Pressure Mid Layer in Pa"
+    return pres