ncar_veg_params.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import xarray as xr\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import collections\n",
    "import warnings \n",
    "from netCDF4 import default_fillvals\n",
    "from scipy.stats import hmean\n",
    "\n",
    "# import veg functions\n",
    "from parameter_functions import (calculate_cv_pft, map_pft_to_nldas_class, is_overstory, \n",
    "                                 calc_root_fract, calc_root_depth_rz1, calc_root_depth_rz2)\n",
    "\n",
    "from parameter_functions_v2 import (create_empty_arrays, \n",
    "                                    create_parameter_dataset, calculate_nveg_pfts)\n",
    "\n",
    "from veg_parameter_functions import create_veg_parameter_dataset\n",
    "\n",
    "# define fillvals\n",
    "fillval_f = default_fillvals['f8']\n",
    "fillval_i = default_fillvals['i4']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a notebook for deriving VIC 5 vegetation parameters for the NCAR RAL group, to be used for high-resolution simulations over Alaska. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Set resolution for filenames__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating parameters at 12km\n"
     ]
    }
   ],
   "source": [
    "res = '12km'\n",
    "\n",
    "# set file extensions \n",
    "if res == \"50km\":\n",
    "    grid = 'wr50a_ar9v4'\n",
    "elif res == \"25km\":\n",
    "    grid = 'wr25b_ar9v4'\n",
    "elif res == \"12km\":\n",
    "    grid = 'ncar_ral'\n",
    "    \n",
    "if res == \"50km\":\n",
    "    nj = 205\n",
    "    ni = 275\n",
    "elif res == \"25km\":\n",
    "    nj = 413\n",
    "    ni = 551\n",
    "elif res == \"12km\":\n",
    "    nj = 209\n",
    "    ni = 299\n",
    "    \n",
    "num_veg = 17\n",
    "print(\"calculating parameters at %s\" %res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Set domain file__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "if res == '50km':\n",
    "    domain = xr.open_dataset(os.path.join('/u/home/gergel/data/parameters', \n",
    "                                      'domain.lnd.wr50a_ar9v4.100920.nc'))\n",
    "elif res == \"25km\":\n",
    "    domain = xr.open_dataset(os.path.join('/u/home/gergel/data/parameters', \n",
    "                                      'domain.lnd.wr25b_ar9v4.170413.nc'))\n",
    "elif res == \"12km\":\n",
    "    domain = xr.open_dataset(os.path.join('/u/home/gergel/data/parameters', \n",
    "                                      'alaska_vic_domain_ncar.nc'))\n",
    "masknan_vals = domain['mask'].where(domain['mask'] == 1).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "old_params = xr.open_dataset(os.path.join('/u/home/gergel/data', \n",
    "                                          'vic_params_wr50a_vic5.0.dev_20160328.nc'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Options__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if set to True, add additional organic_fract options to parameter file, including soil organic fraction, \n",
    "# soil particle density of OM, and bulk density of OM \n",
    "organic_fract = True\n",
    "max_snow_albedo = True\n",
    "bulk_density_comb = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "load CLM PFTs to use for vegetation parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load regridded PFTs \n",
    "pfts_data_dir = '/u/home/gergel/data/parameters/pfts/regridded_pfts'\n",
    "pfts_filename = 'mksrf_landuse_rc2000_c110913_ncar.nc'\n",
    "veg_data = xr.open_dataset(os.path.join(pfts_data_dir, pfts_filename))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Calculate Cv from PFTs__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv = xr.apply_ufunc(calculate_cv_pft, \n",
    "                    veg_data['PCT_PFT'].where(domain.mask == 1),\n",
    "                    dask='allowed',\n",
    "                    vectorize=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Calculate number of active PFTs, `Nveg`__ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nveg = xr.apply_ufunc(calculate_nveg_pfts,\n",
    "                      veg_data['PCT_PFT'].where(domain.mask == 1),\n",
    "                      dask='allowed',\n",
    "                      input_core_dims=[['pft']],\n",
    "                      vectorize=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Load LAI and vegetation height, `MONTHLY_LAI` and `MONTHLY_HEIGHT_TOP`__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "lai_file = xr.open_dataset(os.path.join('/u/home/gergel/data/parameters/lai/regridded_lai', \n",
    "                                   'mksrf_lai_78pfts_simyr2005.c170413_ncar_lai.nc'))\n",
    "veg_height_file = xr.open_dataset(os.path.join('/u/home/gergel/data/parameters/lai/regridded_lai', \n",
    "                                'mksrf_lai_78pfts_simyr2005.c170413_ncar_veg_height.nc'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LAI and veg_height from CLM and `PCT_PFT` from CLM have a different number of PFTs (`PCT_PFT` has one more PFT, 17 vs 16). The extra PFT in `PCT_PFT` has `PCT_PFT` = 0 over the entire RASM domain, so I just slice the LAI and veg_height from the 0th PFT (water/bare soil) and concatenate it for the 16th PFT. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "lai_slice = lai_file['MONTHLY_LAI'].isel(pft = 0)\n",
    "vegheight_slice = veg_height_file['MONTHLY_HEIGHT_TOP'].isel(pft=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "lai = xr.concat([lai_file['MONTHLY_LAI'], lai_slice], dim='pft')\n",
    "veg_height = xr.concat([veg_height_file['MONTHLY_HEIGHT_TOP'], vegheight_slice], dim='pft')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`veg_rough` = 0.123 * `veg_height`\n",
    "\n",
    "`displacement` = 0.67 * `veg_height`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "veg_rough = 0.123 * veg_height\n",
    "displacement = 0.67 * veg_height\n",
    "\n",
    "displacement.values[displacement.values == 0] = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Change dims and order of dims of LAI array__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "lai = lai.rename({'time': 'month', 'pft': 'veg_class'})\n",
    "lai = lai.transpose('veg_class', 'month', 'nj', 'ni')\n",
    "\n",
    "veg_rough = veg_rough.rename({'time': 'month', 'pft': 'veg_class'})\n",
    "veg_rough = veg_rough.transpose('veg_class', 'month', 'nj', 'ni')\n",
    "\n",
    "displacement = displacement.rename({'time': 'month', 'pft': 'veg_class'})\n",
    "displacement = displacement.transpose('veg_class', 'month', 'nj', 'ni')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Note__: map albedo, root zone fraction and root zone depth based on vegetation type. see `deriving_new_parameters_v2.xlsx` sheet titled `PFT-NLDAS Mapping` for mapping between NLDAS vegetation classes (used in old VIC 5 parameters) and CLM PFTs. This mapping is based on obvious relationships and some approximations (used for PFTs 8-11)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Create Dataset for variables and define data_vars__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr_months, arr_nlayer, \\\n",
    "arr_rootzone, arr_veg_classes, \\\n",
    "arr_veg_classes_rootzone, arr_veg_classes_month = create_empty_arrays(domain, nj, ni, num_veg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create DataSet\n",
    "params = create_veg_parameter_dataset(domain, old_params, nj, ni, num_veg, max_snow_albedo)\n",
    "\n",
    "# fill in values\n",
    "params['Cv'].values = cv.values\n",
    "params['Nveg'].values = Nveg.values\n",
    "params['LAI'].values = lai.values.reshape(num_veg, 12, nj, ni)\n",
    "params['displacement'].values = displacement.values.reshape(num_veg, 12, nj, ni)\n",
    "params['veg_rough'].values = veg_rough.values.reshape(num_veg, 12, nj, ni)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "define `trunk_ratio`, `rarc`, `rmin`, `wind_h`, `RGL`, `rad_atten`, `wind_atten`, `overstory`, `max_snow_albedo`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# trunk ratio, rarc, rad_atten\n",
    "trunk_ratio = np.copy(arr_veg_classes)\n",
    "params['trunk_ratio'].values = trunk_ratio * 0.2\n",
    "# adjust for bare soil \n",
    "params['trunk_ratio'].values[0, :, :] = 0.0\n",
    "\n",
    "rarc = np.copy(arr_veg_classes)\n",
    "params['rarc'].values = rarc * 60\n",
    "# adjust for bare soil\n",
    "params['rarc'].values[0, :, :] = 100\n",
    "\n",
    "rad_atten = np.copy(arr_veg_classes)\n",
    "params['rad_atten'].values = rad_atten * 0.5\n",
    "# adjust for bare soil \n",
    "params['rad_atten'].values[0, :, :] = 0.0\n",
    "\n",
    "wind_atten = np.copy(arr_veg_classes)\n",
    "params['wind_atten'].values = wind_atten * 0.5\n",
    "# adjust for bare soil \n",
    "params['wind_atten'].values[0, :, :] = 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "if max_snow_albedo == True:\n",
    "    # max_albedo\n",
    "    for pft in veg_data.pft.values:\n",
    "        # get nldas mapping from pft\n",
    "        nldas = map_pft_to_nldas_class(pft)\n",
    "        if nldas == 0:\n",
    "            max_alb = 0.34\n",
    "        elif nldas == 1:\n",
    "            max_alb = 0.37\n",
    "        elif nldas == 2:\n",
    "            max_alb = 0.35\n",
    "        elif nldas == 3: \n",
    "            max_alb = 0.35\n",
    "        elif nldas == 4: \n",
    "            max_alb = 0.44\n",
    "        elif nldas == 5:\n",
    "            max_alb = 0.69\n",
    "        elif nldas == 6:\n",
    "            max_alb = 0.43\n",
    "        elif nldas == 7:\n",
    "            max_alb = 0.56\n",
    "        elif nldas == 8:\n",
    "            max_alb = 0.70\n",
    "        elif nldas == 9:\n",
    "            max_alb = 0.65\n",
    "        elif nldas == 10:\n",
    "            max_alb = 0.46\n",
    "        elif nldas == 11:\n",
    "            max_alb = 0.84\n",
    "        params['max_snow_albedo'].values[pft, :, :] = np.ones((1, nj, ni)) * max_alb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rmin, wind_h\n",
    "for pft in veg_data.pft.values:\n",
    "    # get nldas mapping from pft\n",
    "    nldas = map_pft_to_nldas_class(pft)\n",
    "    if nldas >= 0 and nldas <= 3:\n",
    "        rmin = np.asscalar(old_params.rmin.isel(veg_class=0).mean())\n",
    "        wind_h = np.asscalar(old_params.wind_h.isel(veg_class=0).mean())\n",
    "    elif nldas == 4:\n",
    "        rmin = np.asscalar(old_params.rmin.isel(veg_class=4).mean())\n",
    "        wind_h = np.asscalar(old_params.wind_h.isel(veg_class=4).mean())\n",
    "    elif nldas >= 5 and nldas <= 6:\n",
    "        rmin = np.asscalar(old_params.rmin.isel(veg_class=5).mean())\n",
    "        wind_h = np.asscalar(old_params.wind_h.isel(veg_class=5).mean())\n",
    "    elif nldas >= 7 and nldas <= 8:\n",
    "        rmin = np.asscalar(old_params.rmin.isel(veg_class=7).mean())\n",
    "        wind_h = np.asscalar(old_params.wind_h.isel(veg_class=7).mean())\n",
    "    elif nldas == 9:\n",
    "        rmin = np.asscalar(old_params.rmin.isel(veg_class=9).mean())\n",
    "        wind_h = np.asscalar(old_params.wind_h.isel(veg_class=9).mean())\n",
    "    elif nldas == 10:\n",
    "        rmin = np.asscalar(old_params.rmin.isel(veg_class=10).mean())\n",
    "        wind_h = np.asscalar(old_params.wind_h.isel(veg_class=10).mean())\n",
    "    elif nldas == 11:\n",
    "        rmin = np.asscalar(old_params.rmin.isel(veg_class=11).mean())\n",
    "        wind_h = np.asscalar(old_params.wind_h.isel(veg_class=11).mean())\n",
    "    params['rmin'].values[pft, :, :] = np.ones((1, nj, ni)) * rmin\n",
    "    params['wind_h'].values[pft, :, :] = np.ones((1, nj, ni)) * wind_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# RGL\n",
    "for pft in veg_data.pft.values:\n",
    "    # get nldas mapping from pft\n",
    "    nldas = map_pft_to_nldas_class(pft)\n",
    "    if nldas >= 0 and nldas <= 3:\n",
    "        rgl = np.asscalar(old_params.wind_h.isel(veg_class=0).mean())\n",
    "    elif nldas >= 4 and nldas <= 5:\n",
    "        rgl = np.asscalar(old_params.wind_h.isel(veg_class=4).mean())\n",
    "    elif nldas >= 6 and nldas <= 8:\n",
    "        rgl = np.asscalar(old_params.wind_h.isel(veg_class=6).mean())\n",
    "    elif nldas >= 9 and nldas <= 10:\n",
    "        rgl = np.asscalar(old_params.wind_h.isel(veg_class=9).mean())\n",
    "    elif nldas == 11:\n",
    "        rgl = np.asscalar(old_params.wind_h.isel(veg_class=11).mean())\n",
    "    params['RGL'].values[pft, :, :] = np.ones((1, nj, ni)) * rgl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# overstory\n",
    "overstory = np.copy(arr_veg_classes)\n",
    "for pft in veg_data.pft.values:\n",
    "    nldas = map_pft_to_nldas_class(pft)\n",
    "    if nldas > 6:\n",
    "        # no overstory\n",
    "        overstory = 0.0\n",
    "    else: \n",
    "        overstory = 1.0\n",
    "    params['overstory'].values[pft, :, :] = overstory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "root fraction and root depth "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "root_depth_rz1 = xr.apply_ufunc(calc_root_depth_rz1,\n",
    "                           params['Cv'].where(domain.mask == 1), \n",
    "                           dask='allowed',\n",
    "                           vectorize=True)\n",
    "root_depth_rz2 = xr.apply_ufunc(calc_root_depth_rz2,\n",
    "                           params['Cv'].where(domain.mask == 1), \n",
    "                           dask='allowed',\n",
    "                           vectorize=True)\n",
    "root_depth = xr.concat([root_depth_rz1, root_depth_rz2],\n",
    "                      dim='root_zone').transpose('veg_class', 'root_zone', 'nj', 'ni')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# root fract \n",
    "\n",
    "rz = 0\n",
    "for pft in veg_data.pft.values:\n",
    "    if pft == 0:\n",
    "        root_fract_rz1 = xr.apply_ufunc(calc_root_fract,\n",
    "                                        params['Cv'].isel(veg_class=pft),\n",
    "                                        str(pft),\n",
    "                                        str(rz),\n",
    "                                        dask='allowed',\n",
    "                                        vectorize=True)\n",
    "    else: \n",
    "        root_fract_rz1 = xr.concat([root_fract_rz1, xr.apply_ufunc(calc_root_fract,\n",
    "                                                                   params['Cv'].isel(veg_class=pft),\n",
    "                                                                   str(pft),\n",
    "                                                                   str(rz),\n",
    "                                                                   dask='allowed',\n",
    "                                                                   vectorize=True)],\n",
    "                                  dim='veg_class')\n",
    "rz = 1\n",
    "for pft in veg_data.pft.values:\n",
    "    if pft == 0:\n",
    "        root_fract_rz2 = xr.apply_ufunc(calc_root_fract,\n",
    "                                        params['Cv'].isel(veg_class=pft), \n",
    "                                        str(pft),\n",
    "                                        str(rz),\n",
    "                                        dask='allowed',\n",
    "                                        vectorize=True)\n",
    "    else: \n",
    "        root_fract_rz2 = xr.concat([root_fract_rz2, xr.apply_ufunc(calc_root_fract,\n",
    "                                                                   params['Cv'].isel(veg_class=pft), \n",
    "                                                                   str(pft),\n",
    "                                                                   str(rz),\n",
    "                                                                   dask='allowed',\n",
    "                                                                   vectorize=True)],\n",
    "                                  dim='veg_class')\n",
    "        \n",
    "root_fract = xr.concat([root_fract_rz1, root_fract_rz2], dim='root_zone').transpose('veg_class', 'root_zone', 'nj', 'ni')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "params['root_depth'].values = root_depth\n",
    "params['root_fract'].values = root_fract"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__albedo__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# loop over pft classes and months \n",
    "for pft in veg_data.pft.values:\n",
    "    for month in old_params.month.values:\n",
    "        nldas = map_pft_to_nldas_class(pft)\n",
    "        if nldas == 0 or nldas == 1:\n",
    "            albedo = 0.12\n",
    "        elif nldas >= 2 and nldas <= 5:\n",
    "            albedo = 0.18\n",
    "        elif nldas >= 6 and nldas <= 8:\n",
    "            albedo = 0.19\n",
    "        elif nldas == 9:\n",
    "            albedo = 0.2\n",
    "        elif nldas == 10: \n",
    "            albedo = 0.12\n",
    "        elif nldas == 11: \n",
    "            albedo = 0.2\n",
    "        params['albedo'].values[pft, month-1, :, :] = np.ones((1, 1, nj, ni)) * albedo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "if max_snow_albedo == True:\n",
    "    veg_class_vars = ['Cv', 'trunk_ratio', 'rarc', 'rmin', 'wind_h', 'RGL', 'rad_atten',\n",
    "                  'wind_atten', 'albedo', 'LAI', 'overstory',\n",
    "                  'root_depth', 'root_fract', 'displacement', 'veg_rough', 'max_snow_albedo']\n",
    "else: \n",
    "    veg_class_vars = ['Cv', 'trunk_ratio', 'rarc', 'rmin', 'wind_h', 'RGL', 'rad_atten',\n",
    "                  'wind_atten', 'albedo', 'LAI', 'overstory',\n",
    "                  'root_depth', 'root_fract', 'displacement', 'veg_rough']\n",
    "\n",
    "for veg_class_var in veg_class_vars: \n",
    "    if ((veg_class_var == \"LAI\") or (veg_class_var == \"albedo\") or (veg_class_var == \"root_depth\") \n",
    "    or (veg_class_var == \"root_fract\") or (veg_class_var == \"veg_rough\") or \n",
    "        (veg_class_var == \"displacement\")):\n",
    "        bare = np.copy(params[veg_class_var].isel(veg_class=0))\n",
    "        last = np.copy(params[veg_class_var].isel(veg_class=16))\n",
    "        params[veg_class_var].values[0, :, :, :] = last\n",
    "        params[veg_class_var].values[16, :, :, :] = bare\n",
    "    else:\n",
    "        bare = np.copy(params[veg_class_var].isel(veg_class=0))\n",
    "        last = np.copy(params[veg_class_var].isel(veg_class=16))\n",
    "        params[veg_class_var].values[0, :, :] = last\n",
    "        params[veg_class_var].values[16, :, :] = bare\n",
    "params['root_fract'].values[16, :, :, :] = 0\n",
    "params['root_depth'].values[16, :, :, :] = 0\n",
    "params['displacement'].values[16, :, :, :] = 0\n",
    "params['veg_rough'].values[16, :, :, :] = 0\n",
    "params['overstory'].values[16, :, :] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adjust data vars that need adjusting \n",
    "params['Nveg'].values = params['Nveg'].where(domain.mask == 1)\n",
    "params['trunk_ratio'] = params['trunk_ratio'].where(domain.mask == 1)\n",
    "params['rarc'] = params['rarc'].where(domain.mask == 1)\n",
    "params['rmin'] = params['rmin'].where(domain.mask == 1)\n",
    "params['wind_h'] = params['wind_h'].where(domain.mask == 1)\n",
    "params['RGL'] = params['RGL'].where(domain.mask == 1)\n",
    "params['rad_atten'] = params['rad_atten'].where(domain.mask == 1)\n",
    "params['wind_atten'] = params['wind_atten'].where(domain.mask == 1)\n",
    "params['root_depth'] = params['root_depth'].where(domain.mask == 1)\n",
    "params['root_fract'] = params['root_fract'].where(domain.mask == 1)\n",
    "params['albedo'] = params['albedo'].where(domain.mask == 1)\n",
    "params['LAI'] = params['LAI'].where(domain.mask == 1)\n",
    "params['overstory'] = params['overstory'].where(domain.mask == 1)\n",
    "params['displacement'] = params['displacement'].where(domain.mask == 1)\n",
    "params['veg_rough'] = params['veg_rough'].where(domain.mask == 1)\n",
    "if max_snow_albedo == True:\n",
    "    params['max_snow_albedo'] = params['max_snow_albedo'].where(domain.mask == 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Organic Fract is True, including soil_dens_org and organic_fract in parameters file\n",
      "Max Snow Albedo is True, including veg-dep snow albedo in parameters file\n",
      "Bulk Density Comb is True, including combined bulk density in parameters file\n",
      "saved new parameters to /u/home/gergel/data/parameters/new_vic5_params_ncar_ral_all_options.nc\n"
     ]
    }
   ],
   "source": [
    "encoding_params = {'Nveg': {'dtype': 'int32', \"_FillValue\": fillval_i},\n",
    "                   'overstory': {'dtype': 'int32', \"_FillValue\": fillval_i}}\n",
    "\n",
    "direc = '/u/home/gergel/data/parameters'\n",
    "if organic_fract == True and bulk_density_comb == True and max_snow_albedo == True:\n",
    "    filename = 'new_vic5_params_%s_%s.nc' %(grid, \"all_options\")\n",
    "elif organic_fract == False and bulk_density_comb == True and max_snow_albedo == False:\n",
    "    filename = 'new_vic5_params_%s_%s.nc' %(grid, \"bulk_density\")\n",
    "else: \n",
    "    filename = 'new_vic5_params_%s_%s.nc' %(grid, \"no_options\")\n",
    "new_params_file = os.path.join(direc, filename)\n",
    "params.to_netcdf(new_params_file, format='NETCDF4_CLASSIC', encoding=encoding_params)\n",
    "\n",
    "if organic_fract:\n",
    "    print(\"Organic Fract is True, including soil_dens_org and organic_fract in parameters file\")\n",
    "if max_snow_albedo:\n",
    "    print(\"Max Snow Albedo is True, including veg-dep snow albedo in parameters file\")\n",
    "if bulk_density_comb:\n",
    "    print(\"Bulk Density Comb is True, including combined bulk density in parameters file\")\n",
    "    \n",
    "print(\"saved new parameters to %s\" %new_params_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}