From dd09a7c31d80c5056e58bf1715d07124590ab36a Mon Sep 17 00:00:00 2001
From: Daniel Kennedy <djk2120@columbia.edu>
Date: Mon, 27 Sep 2021 15:14:00 -0600
Subject: [PATCH] bfb script

---
 analysis/check_bfb.ipynb | 304 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 304 insertions(+)
 create mode 100644 analysis/check_bfb.ipynb

diff --git a/analysis/check_bfb.ipynb b/analysis/check_bfb.ipynb
new file mode 100644
index 0000000..54b13a9
--- /dev/null
+++ b/analysis/check_bfb.ipynb
@@ -0,0 +1,304 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "1a318faf-f427-47da-889f-d577fea42163",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import xarray as xr\n",
+    "import glob\n",
+    "import matplotlib\n",
+    "import matplotlib.pyplot as plt\n",
+    "import cftime\n",
+    "import dask\n",
+    "from dask_jobqueue import PBSCluster\n",
+    "from dask.distributed import Client\n",
+    "import statsmodels.api as sm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "id": "7187399b-ea09-40d5-b5ea-2df757ac1c37",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_ensemble(files,data_vars,p=True):\n",
+    "\n",
+    "    def preprocess(ds):\n",
+    "        return ds[data_vars]\n",
+    "\n",
+    "    #read in the dataset\n",
+    "    ds = xr.open_mfdataset(files,combine='nested',concat_dim='ens',\n",
+    "                           parallel=p,preprocess=preprocess)\n",
+    "\n",
+    "    #fix up time dimension\n",
+    "    htape='h0'\n",
+    "    #if htape=='h0' or htape=='h1':\n",
+    "    #    ds['time'] = xr.cftime_range(str(2005),periods=len(ds.time),freq='MS') #fix time bug\n",
+    "\n",
+    "    #specify extra variables    \n",
+    "    if htape=='h0':\n",
+    "        extras     = ['grid1d_lat','grid1d_lon']\n",
+    "    elif htape=='h1':\n",
+    "        extras     = ['pfts1d_lat','pfts1d_lon','pfts1d_wtgcell','pfts1d_itype_veg']\n",
+    "    \n",
+    "    #add in some extra variables\n",
+    "    ds0 = xr.open_dataset(files[0])\n",
+    "    for extra in extras:\n",
+    "        ds[extra]=ds0[extra]\n",
+    "\n",
+    "    return ds"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "id": "53f265a9-c0a4-4052-a9d9-280261af0449",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Setup your PBSCluster\n",
+    "ncores=1\n",
+    "nmem='25GB'\n",
+    "cluster = PBSCluster(\n",
+    "    cores=ncores, # The number of cores you want\n",
+    "    memory=nmem, # Amount of memory\n",
+    "    processes=1, # How many processes\n",
+    "    queue='casper', # The type of queue to utilize (/glade/u/apps/dav/opt/usr/bin/execcasper)\n",
+    "    local_directory='$TMPDIR', # Use your local directory\n",
+    "    resource_spec='select=1:ncpus='+str(ncores)+':mem='+nmem, # Specify resources\n",
+    "    project='P93300641', # Input your project ID here\n",
+    "    walltime='03:00:00', # Amount of wall time\n",
+    "    interface='ib0', # Interface to use\n",
+    ")\n",
+    "\n",
+    "# Scale up\n",
+    "cluster.scale(20)\n",
+    "\n",
+    "# Setup your client\n",
+    "client = Client(cluster)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "f3be8bac-5847-4e99-a0db-f4f12d3c066a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#fetch the paraminfo\n",
+    "csv = '/glade/scratch/djk2120/PPEn11/SP_bfb_test.csv' \n",
+    "paramkey = pd.read_csv(csv)\n",
+    "\n",
+    "#fetch the sparsegrid landarea\n",
+    "la_file = '/glade/scratch/djk2120/PPEn08/sparsegrid_landarea.nc'\n",
+    "la = xr.open_dataset(la_file).landarea  #km2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 79,
+   "id": "02baf223-a08c-4cba-906a-330571789134",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "kdir = '/glade/scratch/oleson/'\n",
+    "keys  = []; params = []\n",
+    "files = []\n",
+    "for key,param in zip(paramkey.key,paramkey.param):\n",
+    "    thisdir = kdir+'PPEn11_CTL2010SP_'+key+'/run/'\n",
+    "    rfile   = glob.glob(thisdir+'*.clm2.r.*.nc')\n",
+    "    if len(rfile)>0:\n",
+    "        keys.append(key)\n",
+    "        params.append(param)\n",
+    "        h0 = glob.glob(thisdir+'*.clm2.h0.*.nc')\n",
+    "        files.append(h0[0])\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "id": "82e3b26a-f2ad-4bb1-bbd9-8f1511c151dc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "datavars = ['FPSN','EFLX_LH_TOT','FSA','TV','TSOI_10CM','SOILWATER_10CM']\n",
+    "ds =get_ensemble(files,datavars)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 74,
+   "id": "cbfaf107-9924-45a5-b77c-abd25122bc95",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "FPSN 75\n",
+      "EFLX_LH_TOT 75\n",
+      "FSA 75\n",
+      "TV 75\n",
+      "TSOI_10CM 75\n",
+      "SOILWATER_10CM 75\n"
+     ]
+    }
+   ],
+   "source": [
+    "nens = len(ds.ens)\n",
+    "bfb_all = np.zeros(nens)+1\n",
+    "for f in datavars:\n",
+    "\n",
+    "    x0  = ds[v].sel(ens=0)\n",
+    "    isnan = np.tile(np.isnan(x0),[nens,1,1])\n",
+    "    bfb_grid = (ds[v]==x0).values\n",
+    "    bfb_grid[isnan]=1                                 #ignore nans\n",
+    "    bfb = bfb_grid.sum(axis=(1,2))==24*400            #all gridcells / all times must be BFB\n",
+    "    print(f,bfb.sum())\n",
+    "    \n",
+    "    bfb_all = bfb_all*bfb"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 80,
+   "id": "fd01fa9a-2a9b-4950-b785-c8591f13638b",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array(['default', 'grperc', 'br_mr', 'lmr_intercept_atkin',\n",
+       "       'FUN_fracfixers', 'fun_cn_flex_a', 'fun_cn_flex_b',\n",
+       "       'fun_cn_flex_c', 'kc_nonmyc', 'kn_nonmyc', 'akc_active',\n",
+       "       'akn_active', 'ekc_active', 'ekn_active', 'stem_leaf',\n",
+       "       'croot_stem', 'flivewd', 'frootcn', 'leaf_long', 'lwtop_ann',\n",
+       "       'ndays_off', 'ndays_on', 'tau_cwd', 'tau_l1', 'tau_l2_l3',\n",
+       "       'tau_s1', 'tau_s2', 'tau_s3', 'q10_mr', 'minpsi_hr', 'maxpsi_hr',\n",
+       "       'rf_l1s1_bgc', 'rf_l2s1_bgc', 'rf_l3s2_bgc', 'rf_s2s1_bgc',\n",
+       "       'rf_s2s3_bgc', 'rf_s3s1_bgc', 'cn_s3_bgc', 'decomp_depth_efolding',\n",
+       "       'max_altdepth_cryoturbation', 'max_altmultiplier_cryoturb',\n",
+       "       'cryoturb_diffusion_k', 'som_diffus', 'k_nitr_max_perday',\n",
+       "       'denitrif_respiration_coefficient',\n",
+       "       'denitrif_respiration_exponent',\n",
+       "       'denitrif_nitrateconc_coefficient',\n",
+       "       'denitrif_nitrateconc_exponent', 'r_mort', 'fsr_pft', 'fd_pft',\n",
+       "       'prh30', 'ignition_efficiency', 'cc_dstem', 'cc_leaf', 'cc_lstem',\n",
+       "       'cc_other', 'fm_droot', 'fm_leaf', 'fm_lroot', 'fm_lstem',\n",
+       "       'fm_other', 'fm_root', 'KCN', 'LF', 'FR', 'Q10', 'CWD',\n",
+       "       'perched_baseflow_scalar', 'xdrdt', 'frootcn_max', 'frootcn_min',\n",
+       "       'leafcn_max', 'leafcn_min', 'fm_dstem'], dtype='<U32')"
+      ]
+     },
+     "execution_count": 80,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.array(params)[bfb]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 81,
+   "id": "3553bb44-82fd-497a-8361-709f524c71fb",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array(['taulnir', 'taulvis', 'tausnir', 'tausvis', 'rholnir', 'rholvis',\n",
+       "       'rhosnir', 'rhosvis', 'xl', 'displar', 'dleaf', 'z0mr', 'csoilc',\n",
+       "       'cv', 'a_coef', 'a_exp', 'zlnd', 'zsno', 'd_max',\n",
+       "       'frac_sat_soil_dsl_init', 'lai_dl', 'z_dl', 'zetamaxstable',\n",
+       "       'wind_min', 'tkd_sand', 'tkd_clay', 'tkd_om', 'tkm_om', 'pd',\n",
+       "       'csol_om', 'csol_sand', 'csol_clay', 'bsw_sf', 'hksat_sf',\n",
+       "       'sucsat_sf', 'watsat_sf', 'baseflow_scalar',\n",
+       "       'maximum_leaf_wetted_fraction', 'interception_fraction',\n",
+       "       'aq_sp_yield_min', 'fff', 'liq_canopy_storage_scalar',\n",
+       "       'snow_canopy_storage_scalar', 'e_ice', 'n_baseflow', 'n_melt_coef',\n",
+       "       'accum_factor', 'eta0_vionnet', 'drift_gs', 'ssi', 'wimp',\n",
+       "       'upplim_destruct_metamorph', 'wind_snowcompact_fact', 'rho_max',\n",
+       "       'tau_ref', 'snowcan_unload_wind_fact', 'snowcan_unload_temp_fact',\n",
+       "       'snw_rds_refrz', 'scvng_fct_mlt_sf', 'ceta', 'medlynslope',\n",
+       "       'medlynintercept', 'fnps', 'theta_psii', 'theta_ip', 'theta_cj',\n",
+       "       'kc25_coef', 'ko25_coef', 'cp25_yr2000', 'tpu25ratio', 'kp25ratio',\n",
+       "       'lmrse', 'slatop', 'jmaxb0', 'jmaxb1', 'wc2wjb0',\n",
+       "       'enzyme_turnover_daily', 'relhExp', 'minrelh', 'luna_theta_cj',\n",
+       "       'kmax', 'krmax', 'psi50', 'ck', 'rootprof_beta', 'fbw', 'nstem',\n",
+       "       'rstem', 'wood_density', 'froot_leaf', 'leafcn', 'vcmaxha',\n",
+       "       'jmaxha', 'tpuha', 'lmrha', 'kcha', 'koha', 'cpha', 'vcmaxhd',\n",
+       "       'jmaxhd', 'tpuhd', 'lmrhd', 'vcmaxse_sf', 'jmaxse_sf', 'tpuse_sf',\n",
+       "       'jmax25top_sf', 'om_frac_sf', 'slopebeta', 'slopemax', 'pc', 'mu',\n",
+       "       'C2_liq_Brun89', 'fnr', 'act25'], dtype='<U32')"
+      ]
+     },
+     "execution_count": 81,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.array(params)[~bfb]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 82,
+   "id": "3a9ce84a-35f7-44d1-a561-b0c81874617c",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "114"
+      ]
+     },
+     "execution_count": 82,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "(~bfb).sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ce87c02f-477a-4037-b1c6-b7dec42a56ad",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python [conda env:miniconda3-ppe-py]",
+   "language": "python",
+   "name": "conda-env-miniconda3-ppe-py-py"
+  },
+  "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.7.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}