{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import dask\n",
    "import pandas as pd\n",
    "import datetime\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Directory containing model output\n",
    "e3smdir = '/global/cscratch1/sd/plma/e3sm_scratch/EAMv1_ENA/latlon/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Years and months to process in loops\n",
    "stryear = ['2016','2017','2018','2019','2020']\n",
    "intyear = [2016, 2017, 2018, 2019, 2020]\n",
    "strmonth = ['01','02','03','04','05','06','07','08','09','10','11','12']\n",
    "intmonth = [1,2,3,4,5,6,7,8,9,10,11,12]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ENA is located at 39° 5′ 27.96″ N, 28° 1′ 46.92″ W\n",
    "#In the E3SM output, the ENA site is encapsulated by longitude ID = 37 and latitude ID = 22\n",
    "#These are the IDs of longitude and latitude extracted from the E3SM output\n",
    "#They are spaced about 7-8 points apart because the output was interpolated onto a higher resolution grid even though it was run at ~1 deg\n",
    "lonid=[29,37,45]\n",
    "latid=[14,22,29]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Constants\n",
    "G = 9.8\n",
    "Cp = 1005.7\n",
    "Rd = 287.\n",
    "Rv = 461.\n",
    "lv = 2.477e6 #at 10 C\n",
    "epsilon = Rd/Rv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set this to ignore warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Year: 2016\n",
      "Month: 01\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 02\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 03\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 04\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 05\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 06\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 07\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 08\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 09\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 10\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 11\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 12\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Year: 2017\n",
      "Month: 01\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 02\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 03\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 04\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 05\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 06\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 07\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 08\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 09\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 10\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 11\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 12\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Year: 2018\n",
      "Month: 01\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 02\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 03\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 04\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 05\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 06\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 07\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 08\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 09\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 10\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 11\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 12\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Year: 2019\n",
      "Month: 01\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 02\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 03\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 04\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 05\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 06\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 07\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 08\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 09\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 10\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 11\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 12\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Year: 2020\n",
      "Month: 01\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 02\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 03\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 04\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 06\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Month: 07\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 08\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 09\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 10\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 11\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n",
      "Month: 12\n",
      "Opening E3SM files...\n",
      "Reading E3SM output...\n",
      "Computing additional variables...\n",
      "Interpolating atmospheric state variables...\n",
      "Computing cloud base, top, and depth...\n",
      "Saving cloud macrophysical properties...\n"
     ]
    }
   ],
   "source": [
    "#Loop through each year and month to compute relevant variables and save to files to be used in analyses\n",
    "for ii in np.arange(5):\n",
    "    print('Year: '+stryear[ii])\n",
    "    for jj in np.arange(12):\n",
    "        print('Month: '+strmonth[jj])\n",
    "        print('Opening E3SM files...')\n",
    "        e3smfiles = glob.glob(os.path.join(e3smdir, 'EAMv1_ENA.cam.h1.'+stryear[ii]+'-'+strmonth[jj]+'*.nc'))\n",
    "        e3smdata = xr.open_mfdataset(e3smfiles, decode_times=False, combine='by_coords', parallel=True)\n",
    "        \n",
    "        #Define a time array for the model output that xarray can use\n",
    "        e3sm_ds = pd.date_range(strmonth[jj]+'/01/'+stryear[ii], periods = np.size(e3smdata['time']), freq ='H')\n",
    "        e3sm_times = xr.DataArray(e3sm_ds, dims=[\"time\"])\n",
    "        \n",
    "        print('Reading E3SM output...')\n",
    "        ptop=45 #level includes through 6 km to encompass liquid clouds in 3D variables\n",
    "\n",
    "        #time and space\n",
    "        e3sm_time = e3smdata['time']\n",
    "        e3sm_lon = e3smdata['lon'][lonid]\n",
    "        e3sm_lat = e3smdata['lat'][latid]\n",
    "        e3sm_lev = e3smdata['lev'][ptop:72]\n",
    "        e3sm_ilev = e3smdata['ilev'][ptop:73]\n",
    "        e3sm_elev = e3smdata['PHIS'][0,latid,lonid].values/9.80665\n",
    "        e3sm_zw = e3smdata['ZM_CLUBB'][:,ptop:,latid,lonid] + e3sm_elev #dynamics vertical levels\n",
    "        e3sm_z = e3smdata['ZT_CLUBB'][:,ptop:72,latid,lonid] + e3sm_elev #mass vertical levels\n",
    "        \n",
    "        #variables for computing pressure\n",
    "        hyai = e3smdata['hyai'][:,ptop:73]\n",
    "        hybi = e3smdata['hybi'][:,ptop:73]\n",
    "        hyam = e3smdata['hyam'][:,ptop:72]\n",
    "        hybm = e3smdata['hybm'][:,ptop:72]\n",
    "        P0 = e3smdata['P0']\n",
    "        PS = e3smdata['PS'][:,latid,lonid]\n",
    "        \n",
    "        #CCN\n",
    "        e3sm_ccn1 = e3smdata['CCN3'][:,71,latid,lonid]     #surface 0.1% CCN\n",
    "        e3sm_ccn2 = e3smdata['CCN4'][:,71,latid,lonid]     #surface 0.2% CCN\n",
    "        e3sm_ccn5 = e3smdata['CCN5'][:,71,latid,lonid]     #surface 0.5% CCN\n",
    "\n",
    "        #COSP outputted MODIS retrievals\n",
    "        e3sm_modis_cf_liq = e3smdata['CLWMODIS'][:,latid,lonid]     #liquid cloud fraction\n",
    "        e3sm_modis_cf_ice = e3smdata['CLIMODIS'][:,latid,lonid]     #ice cloud fraction\n",
    "        # e3sm_modis_cf = e3smdata['CLTMODIS'][:,latid,lonid]         #total cloud fraction\n",
    "        e3sm_modis_iwp = e3smdata['IWPMODIS'][:,latid,lonid]        #ice water path\n",
    "        e3sm_modis_lwp = e3smdata['LWPMODIS'][:,latid,lonid]        #liquid water path\n",
    "        e3sm_modis_liq_re = e3smdata['REFFCLWMODIS'][:,latid,lonid] #liquid cloud effective radius\n",
    "        e3sm_modis_cod = e3smdata['TAUWMODIS'][:,latid,lonid]       #cloud optical depth\n",
    "        # e3sm_modis_ctp = e3smdata['PCTMODIS'][:,latid,lonid]        #cloud top pressure\n",
    "\n",
    "        #grid-scale cloud properties\n",
    "        e3sm_grid_lwp = e3smdata['TGCLDLWP'][:,latid,lonid]             #liquid water path (includes convection)\n",
    "        # e3sm_grid_twp = e3smdata['TGCLDCWP'][:,latid,lonid]             #total water path (liquid + ice)\n",
    "        # e3sm_cdnc_column = e3smdata['CDNUMC'][:,latid,lonid]            #column integrated cloud droplet number 1/m^2\n",
    "        e3sm_grid_cf = e3smdata['CLOUD'][:,ptop:72,latid,lonid]         #cloud fraction at each height\n",
    "        e3sm_grid_cf_strat = e3smdata['CLOUDFRAC_CLUBB'][:,ptop:72,latid,lonid] #cloud fraction at each height from CLUBB\n",
    "        e3sm_cld_low = e3smdata['CLDLOW'][:,latid,lonid]                #low level cloud fraction\n",
    "        # e3sm_mp_strat_re_liq = e3smdata['REL'][:,ptop:72,latid,lonid]   #stratiform effective radius from MP scheme\n",
    "        e3sm_cld = e3smdata['CLDLIQ'][:,ptop:72,latid,lonid].load()     #grid point cloud liquid mass (kg/kg)\n",
    "        e3sm_cld_num = e3smdata['NUMLIQ'][:,ptop:72,latid,lonid].load() #grid point cloud liquid number (1/kg)\n",
    "        e3sm_grid_rain = e3smdata['RAINQM'][:,ptop:72,latid,lonid]      #grid point rain mass (kg/kg)\n",
    "        e3sm_precip_conv = 3.6e6*e3smdata['PRECC'][:,latid,lonid]       #surface convective rain rate (mm/h)\n",
    "        e3sm_precip_strat = 3.6e6*e3smdata['PRECL'][:,latid,lonid]      #surface stratiform rain rate\n",
    "        e3sm_rain = e3sm_precip_conv + e3sm_precip_strat                #surface total rain rate\n",
    "        \n",
    "        #grid-scale thermodynamics\n",
    "        e3sm_t = e3smdata['T'][:,ptop:72,latid,lonid] #grid point temperature (K)\n",
    "        e3sm_q = e3smdata['Q'][:,ptop:72,latid,lonid] #grid point specific humidity (kg/kg)\n",
    "    \n",
    "        #Radiation\n",
    "        #TOA\n",
    "        e3sm_swup_clear = e3smdata['FSUTOAC'][:,latid,lonid]  #TOA clear sky upwelling broadband SW\n",
    "        e3sm_swup = e3smdata['FSUTOA'][:,latid,lonid]         #TOA upwelling broadband SW\n",
    "        e3sm_swnet_clear = e3smdata['FSNTOAC'][:,latid,lonid] #TOA clear sky net broadband SW\n",
    "        # e3sm_swnet = e3smdata['FSNTOA'][:,latid,lonid]       #TOA net broadband SW\n",
    "        e3sm_swdown = e3sm_swnet_clear + e3sm_swup_clear      #TOA downwelling broadband SW \n",
    "        e3sm_albedo_tot = 1e2*e3sm_swup/e3sm_swdown           #total albedo\n",
    "        e3sm_albedo_clr = 1e2*e3sm_swup_clear/e3sm_swdown     #clear sky albedo\n",
    "        e3sm_albedo_cld = e3sm_albedo_tot - e3sm_albedo_clr   #cloud albedo\n",
    "        #surface\n",
    "        e3sm_sw_down_sfc = e3smdata['FSDS'][:,latid,lonid]   #downwelling solar flux at surface\n",
    "        e3sm_sw_down_clr = e3smdata['FSDSC'][:,latid,lonid]  #downwelling clear sky solar flux at surface\n",
    "        e3sm_cre = e3sm_sw_down_clr - e3sm_sw_down_sfc       #cloud radiative effect\n",
    "        e3sm_cld_eff_albedo = 1e2*(e3sm_sw_down_clr - e3sm_sw_down_sfc)/e3sm_sw_down_clr #surface cloud effective albedo\n",
    "        \n",
    "        e3smdata.close()\n",
    "        \n",
    "        print('Computing additional variables...')\n",
    "        e3sm_modis_iwp_rel = e3sm_modis_iwp/e3sm_modis_cf_ice #MODIS cloud IWP not including clear sky\n",
    "        e3sm_modis_lwp_rel = e3sm_modis_lwp/e3sm_modis_cf_liq #MODIS cloud LWP not including clear sky\n",
    "        e3sm_modis_cod_rel = e3sm_modis_cod/e3sm_modis_cf_liq #MODIS cloud optical depth not including clear sky\n",
    "        e3sm_modis_liq_re_rel = e3sm_modis_liq_re/e3sm_modis_cf_liq #MODIS cloud effective radius not including clear sky\n",
    "        e3sm_grid_lwp_rel = e3sm_grid_lwp/e3sm_cld_low        #Grid cloud LWP not including clear sky\n",
    "        \n",
    "        #compute and expand pressure array for faster calculations\n",
    "        PS_lev = PS.expand_dims(dim={\"lev\": 27}, axis=1)\n",
    "        PS_ilev = PS.expand_dims(dim={\"ilev\": 28}, axis=1)\n",
    "        hyam_latlon = hyam.expand_dims(dim={\"lat\": 3, \"lon\": 3}, axis=[2, 3])\n",
    "        hybm_latlon = hybm.expand_dims(dim={\"lat\": 3, \"lon\": 3}, axis=[2, 3])\n",
    "        hyai_latlon = hyai.expand_dims(dim={\"lat\": 3, \"lon\": 3}, axis=[2, 3])\n",
    "        hybi_latlon = hybi.expand_dims(dim={\"lat\": 3, \"lon\": 3}, axis=[2, 3])\n",
    "        e3sm_pm = hyam_latlon*100000. + hybm_latlon*PS_lev #P0 (reference pressure) is always 100000 Pa\n",
    "        e3sm_pi = hyai_latlon*100000. + hybi_latlon*PS_ilev\n",
    "        e3sm_pi.rename('Pi')\n",
    "        e3sm_pm.rename('Pm')\n",
    "\n",
    "        #Compute relative humidity, potential temperature, and dry air density\n",
    "        e3sm_rh = e3sm_q*100./(0.622*(611.*np.exp((2.5e6/462)*(1/273.15 - 1/e3sm_t)))/\n",
    "                       (e3sm_pm - (611.*np.exp((2.5e6/462)*(1/273.15 - 1/e3sm_t)))))\n",
    "        e3sm_theta = e3sm_t*(1000./(1e-2*e3sm_pm))**0.286\n",
    "        e3sm_rho_air = e3sm_pm/(287.04*e3sm_t)\n",
    "        \n",
    "        #Compute cloud and rain water paths from integrating grid point values\n",
    "        delpi = e3sm_pi[:,1:,:,:].values - e3sm_pi[:,:np.size(e3sm_ilev)-1,:,:].values\n",
    "        \n",
    "        tmprain = e3sm_grid_rain*delpi/9.80665\n",
    "        e3sm_rain_water_path = np.sum(tmprain, axis=1)   #rain water path from column integrated grid point values\n",
    "        tmpcloud = e3sm_cld*delpi/9.80665\n",
    "        e3sm_cloud_water_path = np.sum(tmpcloud, axis=1) #cloud water path from column integrated grid point values\n",
    "        e3sm_total_water_path = e3sm_cloud_water_path + e3sm_rain_water_path #column integrated total water path\n",
    "        \n",
    "        #convert water paths to cloud only values (removing clear sky)\n",
    "        e3sm_rain_water_path_rel = e3sm_rain_water_path/e3sm_cld_low\n",
    "        e3sm_cloud_water_path_rel = e3sm_cloud_water_path/e3sm_cld_low\n",
    "        e3sm_total_water_path_rel = e3sm_total_water_path/e3sm_cld_low\n",
    "        \n",
    "        #grid scale LWP does not include rain so add that here for comparison with observations affected by rain\n",
    "        e3sm_grid_twp = e3sm_grid_lwp + e3sm_rain_water_path\n",
    "        e3sm_grid_twp_rel = e3sm_grid_lwp_rel + e3sm_rain_water_path_rel\n",
    "                                       \n",
    "        print('Interpolating atmospheric state variables...')\n",
    "        #load xarrays for faster numpy computations\n",
    "        e3sm_t.load()\n",
    "        e3sm_q.load()\n",
    "        e3sm_theta.load()\n",
    "        e3sm_z.load()\n",
    "        e3sm_rh.load()\n",
    "        e3sm_pm.load()\n",
    "        e3sm_pi.load()    \n",
    "        e3sm_tw = xr.zeros_like(e3sm_zw)\n",
    "        e3sm_tw = e3sm_tw.rename('Ti')\n",
    "        e3sm_qw = xr.zeros_like(e3sm_zw)\n",
    "        e3sm_qw = e3sm_qw.rename('Qi')\n",
    "        e3sm_tw.load()\n",
    "        e3sm_qw.load()\n",
    "        e3sm_theta_700 = np.zeros((np.size(e3sm_time), np.size(e3sm_lat), np.size(e3sm_lon)))\n",
    "        e3sm_t_700 = np.zeros((np.size(e3sm_time), np.size(e3sm_lat), np.size(e3sm_lon)))\n",
    "        e3sm_z_700 = np.zeros((np.size(e3sm_time), np.size(e3sm_lat), np.size(e3sm_lon)))\n",
    "        e3sm_rh_700 = np.zeros((np.size(e3sm_time), np.size(e3sm_lat), np.size(e3sm_lon)))\n",
    "        e3sm_rh_750 = np.zeros((np.size(e3sm_time), np.size(e3sm_lat), np.size(e3sm_lon)))\n",
    "        e3sm_rh_800 = np.zeros((np.size(e3sm_time), np.size(e3sm_lat), np.size(e3sm_lon)))\n",
    "        for kk in np.arange(np.size(e3sm_time)):\n",
    "            for i in np.arange(np.size(e3sm_lon)):\n",
    "                for j in np.arange(np.size(e3sm_lat)):\n",
    "                    e3sm_tw[kk,:,j,i] = np.interp(e3sm_pi[kk,:,j,i], e3sm_pm[kk,:,j,i], e3sm_t[kk,:,j,i]) #temperature on dynamics vertical levels\n",
    "                    e3sm_qw[kk,:,j,i] = np.interp(e3sm_pi[kk,:,j,i], e3sm_pm[kk,:,j,i], e3sm_q[kk,:,j,i]) #specific humidity on dynamics vertical levels\n",
    "                    e3sm_theta_700[kk,j,i] = np.interp(70000, e3sm_pm[kk,:,j,i], e3sm_theta[kk,:,j,i])    #700-mb potential temperature\n",
    "                    e3sm_t_700[kk,j,i] = np.interp(70000, e3sm_pm[kk,:,j,i], e3sm_t[kk,:,j,i])            #700-mb temperature\n",
    "                    e3sm_z_700[kk,j,i] = np.interp(70000, e3sm_pm[kk,:,j,i], e3sm_z[kk,:,j,i])            #700-mb height\n",
    "                    e3sm_rh_700[kk,j,i] = np.interp(70000, e3sm_pm[kk,:,j,i], e3sm_rh[kk,:,j,i])          #700-mb relative humidity\n",
    "                    e3sm_rh_750[kk,j,i] = np.interp(75000, e3sm_pm[kk,:,j,i], e3sm_rh[kk,:,j,i])          #750-mb relative humidity\n",
    "                    e3sm_rh_800[kk,j,i] = np.interp(80000, e3sm_pm[kk,:,j,i], e3sm_rh[kk,:,j,i])          #800-mb relative humidity\n",
    "            \n",
    "        #compute cloud base and top temperature and height, and cloud depth using dynamics vertical levels (bounding mass points)\n",
    "        #Get cloud base by grabbing temperature of lowest cloud fraction (dynamics level interpolated temperature), multiplying it by the cloud fraction,\n",
    "        #and then stepping up to the next level and doing the same thing until the cloud fraction normalized to the total low level cloud fraction is 1.\n",
    "        print('Computing cloud base, top, and depth...')\n",
    "        e3sm_cld_low_repeat = e3sm_cld_low.expand_dims(dim={\"lev\": 27}, axis=1) #expand low level cloud fraction array for faster vector computations\n",
    "        e3sm_grid_cf_rel = e3sm_grid_cf/e3sm_cld_low_repeat #cloud fraction at each height normalized by the total low cloud fraction\n",
    "        e3sm_grid_cf_flip = e3sm_grid_cf_rel[:,::-1,:,:] #flip array to go from the surface upward\n",
    "        e3sm_grid_cf_sum = e3sm_grid_cf_flip.cumsum(dim='lev') #cumulative sum of normalized cloud fraction moving upward\n",
    "        e3sm_grid_cf_sum_filt = e3sm_grid_cf_sum.where(e3sm_grid_cf_sum <= 1, other=1) #Set all cumulative values > 1 to 1\n",
    "        e3sm_grid_cf_sum_diff = e3sm_grid_cf_sum_filt.diff(\"lev\", n=1) #Difference each level to find cloud contributions to cloud base from each level\n",
    "        \n",
    "        #cloud base height and temperature calculations\n",
    "        e3sm_zw = e3sm_zw.rename({'ilev': 'lev'})\n",
    "        e3sm_tw = e3sm_tw.rename({'ilev': 'lev'})\n",
    "        e3sm_pi = e3sm_pi.rename({'ilev': 'lev'})\n",
    "        e3sm_qw = e3sm_qw.rename({'ilev': 'lev'})\n",
    "        e3sm_zw_flip = e3sm_zw[:,::-1,:,:]\n",
    "        e3sm_tw_flip = e3sm_tw[:,::-1,:,:]\n",
    "        e3sm_pi_flip = e3sm_pi[:,::-1,:,:]\n",
    "        e3sm_qw_flip = e3sm_qw[:,::-1,:,:]\n",
    "        e3sm_cloud_base_z = np.sum(e3sm_zw_flip[:,1:np.size(e3sm_lev),:,:].values*e3sm_grid_cf_sum_diff, axis=1)\n",
    "        e3sm_cloud_base_t = np.sum(e3sm_tw_flip[:,1:np.size(e3sm_lev),:,:].values*e3sm_grid_cf_sum_diff, axis=1)\n",
    "        e3sm_cloud_base_p = np.sum(e3sm_pi_flip[:,1:np.size(e3sm_lev),:,:].values*e3sm_grid_cf_sum_diff, axis=1)\n",
    "        e3sm_cloud_base_q = np.sum(e3sm_qw_flip[:,1:np.size(e3sm_lev),:,:].values*e3sm_grid_cf_sum_diff, axis=1)\n",
    "        \n",
    "        #following cloud base but for cloud top height and temperature weighting by normalized cloud fraction (top down instead of bottom up)\n",
    "        e3sm_grid_cf_sum_top = e3sm_grid_cf_rel.cumsum(dim='lev')\n",
    "        e3sm_grid_cf_sum_top_filt = e3sm_grid_cf_sum_top.where(e3sm_grid_cf_sum_top <= 1, other=1)\n",
    "        e3sm_grid_cf_sum_top_diff = e3sm_grid_cf_sum_top_filt.diff(\"lev\",n=1).load()\n",
    "        e3sm_grid_cf_sum_top_diff[:,0,:,:] = e3sm_grid_cf_sum_top_diff[:,0,:,:] + e3sm_grid_cf_sum_top_filt[:,0,:,:]\n",
    "        e3sm_cloud_top_z = np.sum(e3sm_zw[:,1:np.size(e3sm_lev),:,:].values*e3sm_grid_cf_sum_top_diff, axis=1)\n",
    "        e3sm_cloud_top_t = np.sum(e3sm_tw[:,1:np.size(e3sm_lev),:,:].values*e3sm_grid_cf_sum_top_diff, axis=1)\n",
    "        \n",
    "        #cloud depth (this is a maximum since it does not allow cloud depth to be less than the vertical grid spacing at the cloud level)\n",
    "        e3sm_cloud_depth = e3sm_cloud_top_z - e3sm_cloud_base_z\n",
    "        \n",
    "        #compute cloud base, top, depth using mass levels now\n",
    "        #this can be viewed as an estimate of the minimum cloud depth assuming if vertical resolution were higher that cloud depths would decrease because\n",
    "        #they would not need to always step between discrete levels\n",
    "        e3sm_z2 = e3sm_z[:,:np.size(e3sm_lev)-1,:,:]\n",
    "        e3sm_z3 = e3sm_z[:,1:,:,:]\n",
    "        e3sm_t2 = e3sm_t[:,:np.size(e3sm_lev)-1,:,:]\n",
    "        e3sm_t3 = e3sm_t[:,1:,:,:]\n",
    "        e3sm_p2 = e3sm_pm[:,:np.size(e3sm_lev)-1,:,:]\n",
    "        e3sm_q2 = e3sm_q[:,:np.size(e3sm_lev)-1,:,:]        \n",
    "        e3sm_cloud_base_z2 = np.sum(e3sm_z2[:,::-1,:,:].values*e3sm_grid_cf_sum_diff, axis=1)\n",
    "        e3sm_cloud_base_t2 = np.sum(e3sm_t2[:,::-1,:,:].values*e3sm_grid_cf_sum_diff, axis=1)\n",
    "        e3sm_cloud_base_p2 = np.sum(e3sm_p2[:,::-1,:,:].values*e3sm_grid_cf_sum_diff, axis=1)\n",
    "        e3sm_cloud_base_q2 = np.sum(e3sm_q2[:,::-1,:,:].values*e3sm_grid_cf_sum_diff, axis=1)\n",
    "        e3sm_cloud_top_z2 = np.sum(e3sm_z3[:,:,:,:].values*e3sm_grid_cf_sum_top_diff, axis=1)\n",
    "        e3sm_cloud_top_t2 = np.sum(e3sm_t3[:,:,:,:].values*e3sm_grid_cf_sum_top_diff, axis=1)\n",
    "        e3sm_cloud_depth2 = e3sm_cloud_top_z2 - e3sm_cloud_base_z2\n",
    "        \n",
    "        #average the minimum and maximum cloud depth estimates from dynamics and mass vertical levels to get a best estimate depth, base, and top\n",
    "        #again, this is done to decrease the biases that come from the discrete model vertiacl grid not fully resolving the cloud base and top\n",
    "        e3sm_cloud_depth_avg = 0.5*(e3sm_cloud_depth + e3sm_cloud_depth2)\n",
    "        e3sm_cloud_base_z_avg = 0.5*(e3sm_cloud_base_z + e3sm_cloud_base_z2)\n",
    "        e3sm_cloud_base_t_avg = 0.5*(e3sm_cloud_base_t + e3sm_cloud_base_t2)\n",
    "        e3sm_cloud_base_p_avg = 0.5*(e3sm_cloud_base_p + e3sm_cloud_base_p2)\n",
    "        e3sm_cloud_base_q_avg = 0.5*(e3sm_cloud_base_q + e3sm_cloud_base_q2)\n",
    "        e3sm_cloud_top_z_avg = 0.5*(e3sm_cloud_top_z + e3sm_cloud_top_z2)\n",
    "        e3sm_cloud_top_t_avg = 0.5*(e3sm_cloud_top_t + e3sm_cloud_top_t2)\n",
    "        \n",
    "        print('Computing thermodynamic indices...')\n",
    "        #variable for cloud-surface coupling indices\n",
    "        e3sm_theta_sfc = e3sm_theta[:,-1,:,:] #surface potential temperature\n",
    "        e3sm_q_sfc = e3sm_q[:,-1,:,:] #surface water vapor mixing ratio\n",
    "        e3sm_q_cb = e3sm_cloud_base_q_avg #cloud base water vapor mixing ratio\n",
    "        e3sm_theta_cb = e3sm_cloud_base_t_avg*(1000./(1e-2*e3sm_cloud_base_p_avg))**0.286 #cloud base potential temperature\n",
    "        \n",
    "        #Compute surface dew point to get LCL height for EIS\n",
    "        b = 18.678\n",
    "        c = 257.14\n",
    "        d = 234.5\n",
    "        gamma_m = np.log(1e-2*e3sm_rh[:,-1,:,:]*np.exp((b - ((e3sm_t[:,-1,:,:] - 273.15)/d))*((e3sm_t[:,-1,:,:] - 273.15)/(c + (e3sm_t[:,-1,:,:] - 273.15)))))\n",
    "        e3sm_dp_sfc = c*gamma_m/(b - gamma_m) + 273.15 #surface dew point\n",
    "        \n",
    "        #Compute LCL height\n",
    "        T0 = 0.5*(e3sm_t[:,-1,:,:]-273.15 + e3sm_t_700-273.15)\n",
    "        p0 = 85000 #Pa\n",
    "        es = 611.2*np.exp(17.62*T0/(243.12 + T0))\n",
    "        qs = epsilon*es/(p0 - es)\n",
    "        gamma = (G/Cp)*(1 - (1 + lv*qs/(Rd*(T0 + 273.15)))/(1 + lv**2*qs/(Cp*Rv*(T0 + 273.15)**2))) #moist adiabatic potential temperature lapse rate\n",
    "        e3sm_z_lcl = 125.*(e3sm_t[:,-1,:,:] - e3sm_dp_sfc) + e3sm_z[:,-1,:,:] #LCL height\n",
    "\n",
    "        e3sm_lts = e3sm_theta_700 - e3sm_theta_sfc #lower tropospheric stability\n",
    "        e3sm_eis = e3sm_lts - gamma*(e3sm_z_700 - e3sm_z_lcl) #estimated inversion strength (EIS)\n",
    "        \n",
    "        print('Computing cloud drop number concentrations...')\n",
    "        #cloud droplet numbers here are for cloudy area only (note that this is stratiform only, while column integrated LWP and cloud fraction include convective clouds)\n",
    "        e3sm_cdnc_rel = e3sm_cld_num*e3sm_rho_air/e3sm_grid_cf_strat/1e6 #converting units to #/cc and in cloud only\n",
    "        e3sm_cdnc_rel = e3sm_cdnc_rel.where(e3sm_grid_cf_strat > 0, other = 0) #set points with no cloud fraction to 0\n",
    "\n",
    "        #compute the cloud layer mean drop number by weighting it at each level by the cloud fraction at that level\n",
    "        e3sm_cf_column = e3sm_grid_cf_strat.sum(dim='lev')\n",
    "        e3sm_cf_column_wt = e3sm_grid_cf_strat\n",
    "        e3sm_cdnc_rel_avg_tmp = e3sm_cdnc_rel.dot(e3sm_grid_cf_strat, dims='lev')\n",
    "        e3sm_cdnc_rel_avg = np.divide(e3sm_cdnc_rel_avg_tmp, e3sm_cf_column) #cloud layer average drop concentration (note: not weighted to mass)\n",
    "                \n",
    "        #compute cloud drop concentration from cloud optical depth, LWP, and cloud effective radius like MODIS\n",
    "        pres_const = 85000. #Pa (used by Bennartz, 2007, JGR); could use cloud top pressure instead, but shouldn't alter estimates much\n",
    "        Q = 2. #typical assumed scattering efficiency\n",
    "        k = 0.74 #related to the drop size dispersion; Bennartz (2007, JGR) uses 0.8 +/- 0.1 but can be 0.5-0.9 depending on cloud type\n",
    "        rho_liq = 1000. #liquid water density\n",
    "        C1 = 0.05789\n",
    "        rho_air = pres_const/(Rd*e3sm_cloud_top_t_avg) #air density at cloud top\n",
    "        es = 611.2*np.exp(17.62*(e3sm_cloud_top_t_avg-273.15)/(243.12 + e3sm_cloud_top_t_avg - 273.15))\n",
    "        ws = epsilon*es/(pres_const - es)\n",
    "        gamma_w = G*((1 + lv*ws/(Rd*e3sm_cloud_top_t_avg))/(Cp + lv**2*ws*epsilon/(Rd*e3sm_cloud_top_t_avg**2))) #cloud moist adiabatic lapse rate\n",
    "        gamma_ad = (((epsilon + ws)*ws*lv*gamma_w)/(Rd*e3sm_cloud_top_t_avg**2) - (G*ws*pres_const/(Rd*e3sm_cloud_top_t_avg*(pres_const - es))))*rho_air #cloud adiabatic LWC lapse rate\n",
    "        \n",
    "        #Assume 80% adiabatic fraction to compute cloud depth (H) and cloud droplet number concentration\n",
    "        H = (2.*e3sm_modis_lwp_rel/(0.8*gamma_ad))**0.5\n",
    "        e3sm_cdnc_remote_toa = 1e-6*(((e3sm_modis_cod_rel)**3)/k)*((2*e3sm_modis_lwp_rel)**(-2.5))*((0.6*np.pi*Q)**(-3))*((3./(4.*np.pi*rho_liq))**(-2))*((0.8*gamma_ad)**0.5)\n",
    "        \n",
    "        #Assume fully adiabatic for cloud depth and cloud droplet number concentration\n",
    "        H_ad = (2.*e3sm_modis_lwp_rel/(gamma_ad))**0.5\n",
    "        e3sm_cdnc_remote_toa_ad = 1e-6*((e3sm_modis_cod_rel**3)/k)*((2*e3sm_modis_lwp_rel)**(-2.5))*((0.6*np.pi*Q)**(-3))*((3./(4.*np.pi*rho_liq))**(-2))*((gamma_ad)**0.5)\n",
    "        \n",
    "        #compute cloud droplet number following the ARM surface-based retrieval for several different LWP estimates (including rain and/or convection or not)\n",
    "        e3sm_cdnc_remote_sfc_gridlwp = 1e-6*(C1/k)*(2.**0.5)*(rho_liq**2)*((e3sm_modis_cod_rel)**3)/(e3sm_grid_lwp_rel**2)/e3sm_cloud_depth_avg\n",
    "        e3sm_cdnc_remote_sfc_gridtwp = 1e-6*(C1/k)*(2.**0.5)*(rho_liq**2)*((e3sm_modis_cod_rel)**3)/(e3sm_grid_twp_rel**2)/e3sm_cloud_depth_avg\n",
    "        e3sm_cdnc_remote_sfc_twp = 1e-6*(C1/k)*(2.**0.5)*(rho_liq**2)*((e3sm_modis_cod_rel)**3)/(e3sm_total_water_path_rel**2)/e3sm_cloud_depth_avg\n",
    "        e3sm_cdnc_remote_sfc_modislwp = 1e-6*(C1/k)*(2.**0.5)*(rho_liq**2)*((e3sm_modis_cod_rel)**3)/(e3sm_modis_lwp_rel**2)/e3sm_cloud_depth_avg\n",
    "        \n",
    "        print('Computing adiabatic LWP...')  \n",
    "        e3sm_cloud_base_tC = e3sm_cloud_base_t_avg - 273.15\n",
    "        rho_air = e3sm_cloud_base_p_avg/(Rd*e3sm_cloud_base_t_avg)\n",
    "        es = 611.2*np.exp(17.62*e3sm_cloud_base_tC/(243.12 + e3sm_cloud_base_tC))\n",
    "        ws = epsilon*es/(e3sm_cloud_base_p_avg - es)\n",
    "        gamma_w = G*((1 + lv*ws/(Rd*e3sm_cloud_base_t_avg))/(Cp + lv**2*ws*epsilon/(Rd*e3sm_cloud_base_t_avg**2)))\n",
    "        gamma_ad = (((epsilon + ws)*ws*lv*gamma_w)/(Rd*e3sm_cloud_base_t_avg**2) - (G*ws*e3sm_cloud_base_p_avg/(Rd*e3sm_cloud_base_t_avg*(e3sm_cloud_base_p_avg - es))))*rho_air\n",
    "        e3sm_lwp_ad = 0.5*gamma_ad*e3sm_cloud_depth_avg**2 #Cloud adiabatic LWP using cloud base temperature (following method used for observations)\n",
    "        \n",
    "        #Grab solar zenith angles from 1-deg CERES retrievals to E3SM model output locations and times\n",
    "        print('CERES solar angle...')\n",
    "        ceresdata = xr.open_dataset(os.path.join('ceres/CERES_SYN1deg-1H_Terra-Aqua-MODIS_Ed4.1_Subset_'+stryear[ii]+'0101-'+stryear[ii]+'1231.nc'))\n",
    "#        ceres_lonid = 5 #longitude closest to ENA in CERES file\n",
    "#        ceres_latid = 6 #latitude closest to ENA in CERES file\n",
    "        ceres_lonid = [3,5,6]\n",
    "        ceres_latid = [4,6,7]\n",
    "        e3sm_solar_angle = ceresdata['solar_zen_angle_1h'][np.where(ceresdata['time'].dt.month == intmonth[jj])[0],ceres_latid,ceres_lonid]\n",
    "        if intmonth[jj] == 2:\n",
    "            e3sm_solar_angle = e3sm_solar_angle[np.where(e3sm_solar_angle['time'].dt.day < 29)[0],:,:] #solar zenith angle\n",
    "\n",
    "        #Write variables to files to be used in analyses\n",
    "        print('Saving cloud macrophysical properties...')\n",
    "        ds = xr.Dataset({#'base_time': ('time', np.arange(1)), 'base_time': base_time,\n",
    "                         'time_offset': ('time', e3sm_times),\n",
    "                         'cloud_fraction_modis': (('time','lat','lon'), np.float32(e3sm_modis_cf_liq)),\n",
    "                         'cloud_fraction_grid': (('time','lat','lon'), np.float32(1e2*e3sm_cld_low)),\n",
    "                         'rain_rate': (('time','lat','lon'), np.float32(e3sm_rain)),\n",
    "                         'rain_rate_conv': (('time','lat','lon'), np.float32(e3sm_precip_conv)),\n",
    "                         'rain_rate_strat': (('time','lat','lon'), np.float32(e3sm_precip_strat)),\n",
    "                         'cloud_depth': (('time','lat','lon'), np.float32(e3sm_cloud_depth_avg)),\n",
    "                         'cloud_base_height': (('time','lat','lon'), np.float32(e3sm_cloud_base_z_avg)),\n",
    "                         'cloud_base_temperature': (('time','lat','lon'), np.float32(e3sm_cloud_base_t_avg)),\n",
    "                         'cloud_base_pressure': (('time','lat','lon'), np.float32(e3sm_cloud_base_p_avg)),\n",
    "                         'cloud_top_height': (('time','lat','lon'), np.float32(e3sm_cloud_top_z_avg)),\n",
    "                         'cloud_top_temperature': (('time','lat','lon'), np.float32(e3sm_cloud_top_t_avg))},\n",
    "                         coords={'time': ('time', e3sm_times),\n",
    "                                 'lon': ('lon', e3sm_lon),\n",
    "                                 'lat': ('lat', e3sm_lat)})\n",
    "\n",
    "        #assign attributes\n",
    "        ds['time'].attrs[\"long_name\"] = \"UTC Datetime\"\n",
    "        ds['time'].attrs[\"standard_name\"] = \"time\"\n",
    "        ds['time_offset'].attrs[\"long_name\"] = \"UTC Datetime\"\n",
    "        ds['time_offset'].attrs[\"ancillary_variables\"] = \"base_time\"\n",
    "        ds['cloud_fraction_modis'].attrs[\"long_name\"] = \"Cloud Fraction\"\n",
    "        ds['cloud_fraction_modis'].attrs[\"units\"] = \"%\"\n",
    "        ds['cloud_fraction_modis'].attrs[\"description\"] = \"Using MODIS simulator\"\n",
    "        ds['cloud_fraction_grid'].attrs[\"long_name\"] = \"Cloud Fraction\"\n",
    "        ds['cloud_fraction_grid'].attrs[\"units\"] = \"%\"\n",
    "        ds['cloud_fraction_grid'].attrs[\"description\"] = \"Grid scale\"\n",
    "        ds['rain_rate'].attrs[\"long_name\"] = \"Total surface rain rate\"\n",
    "        ds['rain_rate'].attrs[\"units\"] = \"mm/h\"\n",
    "        ds['rain_rate'].attrs[\"description\"] = \"Instantaneous\"\n",
    "        ds['rain_rate_conv'].attrs[\"long_name\"] = \"Convective surface rain rate\"\n",
    "        ds['rain_rate_conv'].attrs[\"units\"] = \"mm/h\"\n",
    "        ds['rain_rate_conv'].attrs[\"description\"] = \"Instantaneous\"\n",
    "        ds['rain_rate_strat'].attrs[\"long_name\"] = \"Stratiform surface rain rate\"\n",
    "        ds['rain_rate_strat'].attrs[\"units\"] = \"mm/h\"\n",
    "        ds['rain_rate_strat'].attrs[\"description\"] = \"Instantaneous\"\n",
    "        ds['cloud_depth'].attrs[\"long_name\"] = \"Cloud Depth\"\n",
    "        ds['cloud_depth'].attrs[\"units\"] = \"m\"\n",
    "        ds['cloud_depth'].attrs[\"description\"] = \"Using estimate cloud base and top heights\"\n",
    "        ds['cloud_base_height'].attrs[\"long_name\"] = \"Cloud Base Height\"\n",
    "        ds['cloud_base_height'].attrs[\"units\"] = \"m\"\n",
    "        ds['cloud_base_height'].attrs[\"description\"] = \"Weighting height by cloud fraction from bottom up until low cloud fraction is reached, assuming maximum cloud overlap by height\"\n",
    "        ds['cloud_base_pressure'].attrs[\"long_name\"] = \"Cloud Base Pressure\"\n",
    "        ds['cloud_base_pressure'].attrs[\"units\"] = \"Pa\"\n",
    "        ds['cloud_base_pressure'].attrs[\"description\"] = \"Matched to Cloud Base Height Estimate\"\n",
    "        ds['cloud_base_temperature'].attrs[\"long_name\"] = \"Cloud Base Temperature\"\n",
    "        ds['cloud_base_temperature'].attrs[\"units\"] = \"C\"\n",
    "        ds['cloud_base_temperature'].attrs[\"description\"] = \"Matched to Cloud Base Height Estimate\"\n",
    "        ds['cloud_top_height'].attrs[\"long_name\"] = \"Cloud Top Height\"\n",
    "        ds['cloud_top_height'].attrs[\"units\"] = \"m\"\n",
    "        ds['cloud_top_height'].attrs[\"description\"] = \"Weighting height by cloud fraction from top down until low cloud fraction is reached, assuming maximum cloud overlap by height\"\n",
    "        ds['cloud_top_temperature'].attrs[\"long_name\"] = \"Cloud Top Temperature\"\n",
    "        ds['cloud_top_temperature'].attrs[\"units\"] = \"C\"\n",
    "        ds['cloud_top_temperature'].attrs[\"description\"] = \"Matched to Cloud Base Height Estimate\"\n",
    "\n",
    "        ds.attrs[\"description\"] = \"Cloud depth retrievals from grid average cloud fraction profiles are a work in progress.\"\n",
    "        ds.attrs[\"contact\"] = \"Adam Varble, adam.varble@pnnl.gov\"\n",
    "        ds.attrs[\"date\"] = \"11 March 2023\"\n",
    "\n",
    "#        outfile = 'E3SM_1Deg_ENA_Cloud_Macro_Properties_'+stryear[ii]+'-'+strmonth[jj]+'.nc'\n",
    "        outfile = 'E3SM_1Deg_ENA_3x3_Cloud_Macro_Properties_'+stryear[ii]+'-'+strmonth[jj]+'.nc'\n",
    "        ds.to_netcdf(outfile, mode='w')\n",
    "\n",
    "        print('Saving cloud drop number concentration properties...')\n",
    "        ds = xr.Dataset({#'base_time': ('time', np.arange(1)), 'base_time': base_time,\n",
    "                         'time_offset': ('time', e3sm_times),\n",
    "                         'cloud_drop_concentration': (('time','lat','lon'), np.float32(e3sm_cdnc_rel_avg)),\n",
    "                         'cloud_drop_concentration_remote_toa': (('time','lat','lon'), np.float32(e3sm_cdnc_remote_toa)),\n",
    "                         'cloud_drop_concentration_remote_toa_ad': (('time','lat','lon'), np.float32(e3sm_cdnc_remote_toa_ad)),\n",
    "                         'cloud_drop_concentration_remote_sfc_gridlwp': (('time','lat','lon'), np.float32(e3sm_cdnc_remote_sfc_gridlwp)),\n",
    "                         'cloud_drop_concentration_remote_sfc_gridtwp': (('time','lat','lon'), np.float32(e3sm_cdnc_remote_sfc_gridtwp)),\n",
    "                         'cloud_drop_concentration_remote_sfc_twp': (('time','lat','lon'), np.float32(e3sm_cdnc_remote_sfc_twp)),\n",
    "                         'cloud_drop_concentration_remote_sfc_modislwp': (('time','lat','lon'), np.float32(e3sm_cdnc_remote_sfc_modislwp))},\n",
    "                         coords={'time': ('time', e3sm_times),\n",
    "                                 'lon': ('lon', e3sm_lon),\n",
    "                                 'lat': ('lat', e3sm_lat)})\n",
    "\n",
    "        #assign attributes\n",
    "        ds['time'].attrs[\"long_name\"] = \"UTC Datetime\"\n",
    "        ds['time'].attrs[\"standard_name\"] = \"time\"\n",
    "        ds['time_offset'].attrs[\"long_name\"] = \"UTC Datetime\"\n",
    "        ds['time_offset'].attrs[\"ancillary_variables\"] = \"base_time\"\n",
    "        ds['cloud_drop_concentration'].attrs[\"long_name\"] = \"Cloud Droplet Concentration\"\n",
    "        ds['cloud_drop_concentration'].attrs[\"units\"] = \" cm-3\"\n",
    "        ds['cloud_drop_concentration'].attrs[\"description\"] = \"Retrieved by weighting each cloud layer droplet number by cloud fraction. Does not include convective scheme contributions.\"\n",
    "        ds['cloud_drop_concentration_remote_toa'].attrs[\"long_name\"] = \"Cloud Droplet Concentration\"\n",
    "        ds['cloud_drop_concentration_remote_toa'].attrs[\"units\"] = \" cm-3\"\n",
    "        ds['cloud_drop_concentration_remote_toa'].attrs[\"description\"] = \"Retrieved following the Bennartz method using optical depth, cloud lwp, and cloud effective radius with adiabaticity of 0.8\"\n",
    "        ds['cloud_drop_concentration_remote_toa_ad'].attrs[\"long_name\"] = \"Cloud Droplet Concentration\"\n",
    "        ds['cloud_drop_concentration_remote_toa_ad'].attrs[\"units\"] = \" cm-3\"\n",
    "        ds['cloud_drop_concentration_remote_toa_ad'].attrs[\"description\"] = \"Retrieved following the Bennartz method using optical depth, cloud lwp, and cloud effective radius with adiabaticity of 1\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_gridlwp'].attrs[\"long_name\"] = \"Cloud Droplet Concentration\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_gridlwp'].attrs[\"units\"] = \" cm-3\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_gridlwp'].attrs[\"description\"] = \"Retrieved following the ARM NDROP method using MODIS simulator cloud optical depth, E3SM lwp, and cloud depth\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_gridtwp'].attrs[\"long_name\"] = \"Cloud Droplet Concentration\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_gridtwp'].attrs[\"units\"] = \" cm-3\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_gridtwp'].attrs[\"description\"] = \"Retrieved following the ARM NDROP method using MODIS simulator cloud optical depth, E3SM lwp plus rain water path (twp), and cloud depth\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_twp'].attrs[\"long_name\"] = \"Cloud Droplet Concentration\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_twp'].attrs[\"units\"] = \" cm-3\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_twp'].attrs[\"description\"] = \"Retrieved following the ARM NDROP method using MODIS simulator cloud optical depth, E3SM vertically integrated cloud and rain water contents, and cloud depth\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_modislwp'].attrs[\"long_name\"] = \"Cloud Droplet Concentration\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_modislwp'].attrs[\"units\"] = \" cm-3\"\n",
    "        ds['cloud_drop_concentration_remote_sfc_modislwp'].attrs[\"description\"] = \"Retrieved following the ARM NDROP method using MODIS simulator cloud optical depth, MODIS simulator lwp, and cloud depth\"\n",
    "\n",
    "        ds.attrs[\"description\"] = \"Surface remote rerieval uses MODIS simulator COD, cloud depth estimated from E3SM predicted cloud fraction, and either E3SM or MODIS estimated LWP. TOA retrievals use MODIS variables.\"\n",
    "        ds.attrs[\"contact\"] = \"Adam Varble, adam.varble@pnnl.gov\"\n",
    "        ds.attrs[\"date\"] = \"11 November 2022\"\n",
    "\n",
    "#        outfile = 'E3SM_1Deg_ENA_CDNC_Properties_'+stryear[ii]+'-'+strmonth[jj]+'.nc'\n",
    "        outfile = 'E3SM_1Deg_ENA_3x3_CDNC_Properties_'+stryear[ii]+'-'+strmonth[jj]+'.nc'\n",
    "        ds.to_netcdf(outfile, mode='w')\n",
    "\n",
    "        print('Saving LWP properties...')\n",
    "        ds = xr.Dataset({#'base_time': ('time', np.arange(1)), 'base_time': base_time,\n",
    "                         'time_offset': ('time', e3sm_times),\n",
    "                         'liquid_water_path_grid': (('time','lat','lon'), np.float32(e3sm_grid_lwp)),\n",
    "                         'liquid_water_path_cloud': (('time','lat','lon'), np.float32(e3sm_grid_lwp_rel)),\n",
    "                         'liquid_water_path_cloud_adiabatic': (('time','lat','lon'), np.float32(e3sm_lwp_ad)),\n",
    "                         'liquid_water_path_modis': (('time','lat','lon'), np.float32(e3sm_modis_lwp_rel)),\n",
    "                         'liquid_water_path_modis_cloud': (('time','lat','lon'), np.float32(e3sm_modis_lwp_rel)),\n",
    "                         'ice_water_path_modis': (('time','lat','lon'), np.float32(e3sm_modis_iwp/1e2)),\n",
    "                         'rain_water_path_grid': (('time','lat','lon'), np.float32(e3sm_rain_water_path)),\n",
    "                         'cloud_water_path_grid': (('time','lat','lon'), np.float32(e3sm_cloud_water_path)),\n",
    "                         'rain_water_path_cloud': (('time','lat','lon'), np.float32(e3sm_rain_water_path_rel)),\n",
    "                         'cloud_water_path_cloud': (('time','lat','lon'), np.float32(e3sm_cloud_water_path_rel))},\n",
    "                         coords={'time': ('time', e3sm_times),\n",
    "                                 'lon': ('lon', e3sm_lon),\n",
    "                                 'lat': ('lat', e3sm_lat)})\n",
    "\n",
    "        #assign attributes\n",
    "        ds['time'].attrs[\"long_name\"] = \"UTC Datetime\"\n",
    "        ds['time'].attrs[\"standard_name\"] = \"time\"\n",
    "        ds['time_offset'].attrs[\"long_name\"] = \"UTC Datetime\"\n",
    "        ds['time_offset'].attrs[\"ancillary_variables\"] = \"base_time\"\n",
    "        ds['liquid_water_path_grid'].attrs[\"long_name\"] = \"Grid Cloud Liquid Water Path\"\n",
    "        ds['liquid_water_path_grid'].attrs[\"units\"] = \"kg m-2\"\n",
    "        ds['liquid_water_path_grid'].attrs[\"description\"] = \"Grid box mean including convection (averages in non-cloudy area)\"\n",
    "        ds['liquid_water_path_cloud'].attrs[\"long_name\"] = \"Cloud Liquid Water Path\"\n",
    "        ds['liquid_water_path_cloud'].attrs[\"units\"] = \"kg m-2\"\n",
    "        ds['liquid_water_path_cloud'].attrs[\"description\"] = \"Cloud mean including convection (excludes non-cloudy area)\"\n",
    "        ds['liquid_water_path_cloud_adiabatic'].attrs[\"long_name\"] = \"Adiabatic Cloud Liquid Water Path\"\n",
    "        ds['liquid_water_path_cloud_adiabatic'].attrs[\"units\"] = \"kg m-2\"\n",
    "        ds['liquid_water_path_cloud_adiabatic'].attrs[\"description\"] = \"Computed following Braun et al. (2018, JGR) using cloud depth and lwp lapse rate computed from cloud base properties\"\n",
    "        ds['liquid_water_path_modis'].attrs[\"long_name\"] = \"MODIS Grid Liquid Water Path\"\n",
    "        ds['liquid_water_path_modis'].attrs[\"units\"] = \"kg m-2\"\n",
    "        ds['liquid_water_path_modis'].attrs[\"description\"] = \"Using MODIS simulator\"\n",
    "        ds['liquid_water_path_modis_cloud'].attrs[\"long_name\"] = \"MODIS Cloud Liquid Water Path\"\n",
    "        ds['liquid_water_path_modis_cloud'].attrs[\"units\"] = \"kg m-2\"\n",
    "        ds['liquid_water_path_modis_cloud'].attrs[\"description\"] = \"Using MODIS simulator\"\n",
    "        ds['ice_water_path_modis'].attrs[\"long_name\"] = \"MODIS Grid Ice Water Path\"\n",
    "        ds['ice_water_path_modis'].attrs[\"units\"] = \"kg m-2\"\n",
    "        ds['ice_water_path_modis'].attrs[\"description\"] = \"Using MODIS simulator\"\n",
    "        ds['rain_water_path_grid'].attrs[\"long_name\"] = \"Grid Rain Water Path\"\n",
    "        ds['rain_water_path_grid'].attrs[\"units\"] = \"kg m-2\"\n",
    "        ds['rain_water_path_grid'].attrs[\"description\"] = \"Grid box mean excluding convection computed by vertically integrating RAINQM\"\n",
    "        ds['cloud_water_path_grid'].attrs[\"long_name\"] = \"Grid Cloud Liquid Water Path\"\n",
    "        ds['cloud_water_path_grid'].attrs[\"units\"] = \"kg m-2\"\n",
    "        ds['cloud_water_path_grid'].attrs[\"description\"] = \"Grid box mean excluding convection computed by vertically integrating CLOUD\"\n",
    "        ds['rain_water_path_cloud'].attrs[\"long_name\"] = \"Rain Water Path\"\n",
    "        ds['rain_water_path_cloud'].attrs[\"units\"] = \"kg m-2\"\n",
    "        ds['rain_water_path_cloud'].attrs[\"description\"] = \"Cloud mean excluding convection computed by vertically integrating RAINQM\"\n",
    "        ds['cloud_water_path_cloud'].attrs[\"long_name\"] = \"Cloud Liquid Water Path\"\n",
    "        ds['cloud_water_path_cloud'].attrs[\"units\"] = \"kg m-2\"\n",
    "        ds['cloud_water_path_cloud'].attrs[\"description\"] = \"Cloud mean excluding convection computed by vertically integrating CLOUD\"\n",
    "        \n",
    "        ds.attrs[\"description\"] = \"Adiabatic lwp is still a work in progress since it is heavily influenced by cloud depth, which is difficult to diagnose.\"\n",
    "        ds.attrs[\"contact\"] = \"Adam Varble, adam.varble@pnnl.gov\"\n",
    "        ds.attrs[\"date\"] = \"9 November 2022\"\n",
    "\n",
    "#        outfile = 'E3SM_1Deg_ENA_LWP_Properties_'+stryear[ii]+'-'+strmonth[jj]+'.nc'\n",
    "        outfile = 'E3SM_1Deg_ENA_3x3_LWP_Properties_'+stryear[ii]+'-'+strmonth[jj]+'.nc'\n",
    "        ds.to_netcdf(outfile, mode='w')\n",
    "\n",
    "        print('Saving radiation properties...')\n",
    "        ds = xr.Dataset({#'base_time': ('time', np.arange(1)), 'base_time': base_time,\n",
    "                         'time_offset': ('time', e3sm_times),\n",
    "                         'cloud_optical_depth_modis': (('time','lat','lon'), np.float32(e3sm_modis_cod_rel)),\n",
    "                         'cloud_effective_radius_modis': (('time','lat','lon'), np.float32(e3sm_modis_liq_re_rel)),\n",
    "                         'solar_zenith': (('time','lat','lon'), np.float32(e3sm_solar_angle)),\n",
    "                         'cre': (('time','lat','lon'), np.float32(e3sm_cre)),\n",
    "                         'cloud_effective_albedo': (('time','lat','lon'), np.float32(e3sm_cld_eff_albedo)),\n",
    "                         'sw_down': (('time','lat','lon'), np.float32(e3sm_sw_down_sfc)),\n",
    "                         'sw_down_clear': (('time','lat','lon'), np.float32(e3sm_sw_down_clr)),\n",
    "                         'sw_up': (('time','lat','lon'), np.float32(e3sm_swup)),\n",
    "                         'sw_up_clear': (('time','lat','lon'), np.float32(e3sm_swup_clear)), \n",
    "                         'albedo_total': (('time','lat','lon'), np.float32(e3sm_albedo_tot)),\n",
    "                         'albedo_clear': (('time','lat','lon'), np.float32(e3sm_albedo_clr)),\n",
    "                         'albedo_cloud': (('time','lat','lon'), np.float32(e3sm_albedo_cld))},\n",
    "                         coords={'time': ('time', e3sm_times),\n",
    "                                 'lon': ('lon', e3sm_lon),\n",
    "                                 'lat': ('lat', e3sm_lat)})\n",
    "\n",
    "        #assign attributes\n",
    "        ds['time'].attrs[\"long_name\"] = \"UTC Datetime\"\n",
    "        ds['time'].attrs[\"standard_name\"] = \"time\"\n",
    "        ds['time_offset'].attrs[\"long_name\"] = \"UTC Datetime\"\n",
    "        ds['time_offset'].attrs[\"ancillary_variables\"] = \"base_time\"\n",
    "        ds['cloud_optical_depth_modis'].attrs[\"long_name\"] = \"Cloud Optical Depth\"\n",
    "        ds['cloud_optical_depth_modis'].attrs[\"units\"] = \" \"\n",
    "        ds['cloud_optical_depth_modis'].attrs[\"description\"] = \"Retrieved by MODIS Simulator\"\n",
    "        ds['cloud_effective_radius_modis'].attrs[\"long_name\"] = \"Cloud Effective Radius\"\n",
    "        ds['cloud_effective_radius_modis'].attrs[\"units\"] = \"m\"\n",
    "        ds['cloud_effective_radius_modis'].attrs[\"description\"] = \"Using MODIS simulator\"\n",
    "        ds['solar_zenith'].attrs[\"long_name\"] = \"Solar Zenith Angle\"\n",
    "        ds['solar_zenith'].attrs[\"units\"] = \" \"\n",
    "        ds['solar_zenith'].attrs[\"description\"] = \"From CERES\"\n",
    "        ds['cre'].attrs[\"long_name\"] = \"Cloud Radiative Effect\"\n",
    "        ds['cre'].attrs[\"units\"] = \"W m-2\"\n",
    "        ds['cre'].attrs[\"description\"] = \"Using downwelling surface SW radiation\"\n",
    "        ds['cloud_effective_albedo'].attrs[\"long_name\"] = \"Cloud Effective Albedo\"\n",
    "        ds['cloud_effective_albedo'].attrs[\"units\"] = \"%\"\n",
    "        ds['cloud_effective_albedo'].attrs[\"description\"] = \"Using downwelling surface SW radiation\"\n",
    "        ds['sw_down'].attrs[\"long_name\"] = \"Downwelling Broadband Shortwave Radiation\"\n",
    "        ds['sw_down'].attrs[\"units\"] = \"W m-2\"\n",
    "        ds['sw_down'].attrs[\"description\"] = \"At the surface\"\n",
    "        ds['sw_down_clear'].attrs[\"long_name\"] = \"Clear Sky Downwelling Broadband Shortwave Irradiance\"\n",
    "        ds['sw_down_clear'].attrs[\"units\"] = \"W m-2\"\n",
    "        ds['sw_down_clear'].attrs[\"description\"] = \"At the surface\"\n",
    "        ds['sw_up'].attrs[\"long_name\"] = \"Upwelling Broadband Shortwave Radiation\"\n",
    "        ds['sw_up'].attrs[\"units\"] = \"W m-2\"\n",
    "        ds['sw_up'].attrs[\"description\"] = \"At TOA\"\n",
    "        ds['sw_up_clear'].attrs[\"long_name\"] = \"Clear Skey Upwelling Broadband Shortwave Radiation\"\n",
    "        ds['sw_up_clear'].attrs[\"units\"] = \"W m-2\"\n",
    "        ds['sw_up_clear'].attrs[\"description\"] = \"At TOA\"\n",
    "        ds['albedo_total'].attrs[\"long_name\"] = \"Grid Albedo\"\n",
    "        ds['albedo_total'].attrs[\"units\"] = \"%\"\n",
    "        ds['albedo_total'].attrs[\"description\"] = \"At TOA, incoporating the entire grid box inclusive of both cloudy and clear area\"\n",
    "        ds['albedo_clear'].attrs[\"long_name\"] = \"Clear Sky Albedo\"\n",
    "        ds['albedo_clear'].attrs[\"units\"] = \"%\"\n",
    "        ds['albedo_clear'].attrs[\"description\"] = \"At TOA, if no clouds\"\n",
    "        ds['albedo_cloud'].attrs[\"long_name\"] = \"Grid Albedo\"\n",
    "        ds['albedo_cloud'].attrs[\"units\"] = \"%\"\n",
    "        ds['albedo_cloud'].attrs[\"description\"] = \"At TOA, only cloud component\"\n",
    "\n",
    "        ds.attrs[\"description\"] = \" \"\n",
    "        ds.attrs[\"contact\"] = \"Adam Varble, adam.varble@pnnl.gov\"\n",
    "        ds.attrs[\"date\"] = \"11 May 2022\"\n",
    "\n",
    "#        outfile = 'E3SM_1Deg_ENA_Radiation_Properties_'+stryear[ii]+'-'+strmonth[jj]+'.nc'\n",
    "        outfile = 'E3SM_1Deg_ENA_3x3_Radiation_Properties_'+stryear[ii]+'-'+strmonth[jj]+'.nc'\n",
    "        ds.to_netcdf(outfile, mode='w')\n",
    "\n",
    "        print('Saving atmospheric state properties...')\n",
    "        ds = xr.Dataset({#'base_time': ('time', np.arange(1)), 'base_time': base_time,\n",
    "                         'time_offset': ('time', e3sm_times),\n",
    "                         'qv_sfc': (('time','lat','lon'), np.float32(e3sm_q_sfc)),\n",
    "                         'qv_cb': (('time','lat','lon'), np.float32(e3sm_q_cb)),\n",
    "                         'theta_sfc': (('time','lat','lon'), np.float32(e3sm_theta_sfc)),\n",
    "                         'theta_cb': (('time','lat','lon'), np.float32(e3sm_theta_cb)),\n",
    "                         'z_lcl': (('time','lat','lon'), np.float32(e3sm_z_lcl)),\n",
    "                         'LTS': (('time','lat','lon'), np.float32(e3sm_lts)),\n",
    "                         'EIS': (('time','lat','lon'), np.float32(e3sm_eis)),\n",
    "                         'rh_700': (('time','lat','lon'), np.float32(e3sm_rh_700)),\n",
    "                         'rh_750': (('time','lat','lon'), np.float32(e3sm_rh_750)),\n",
    "                         'rh_800': (('time','lat','lon'), np.float32(e3sm_rh_800)),\n",
    "                         'ccn1': (('time','lat','lon'), np.float32(e3sm_ccn1)),\n",
    "                         'ccn2': (('time','lat','lon'), np.float32(e3sm_ccn2)),\n",
    "                         'ccn5': (('time','lat','lon'), np.float32(e3sm_ccn5))},\n",
    "                         coords={'time': ('time', e3sm_times),\n",
    "                                 'lon': ('lon', e3sm_lon),\n",
    "                                 'lat': ('lat', e3sm_lat)})\n",
    "\n",
    "        #assign attributes\n",
    "        ds['time'].attrs[\"long_name\"] = \"UTC Datetime\"\n",
    "        ds['time'].attrs[\"standard_name\"] = \"time\"\n",
    "        ds['time_offset'].attrs[\"long_name\"] = \"UTC Datetime\"\n",
    "        ds['time_offset'].attrs[\"ancillary_variables\"] = \"base_time\"\n",
    "        ds['qv_sfc'].attrs[\"long_name\"] = \"Near Surface Water Vapor Mixing Ratio\"\n",
    "        ds['qv_sfc'].attrs[\"units\"] = \"kg/kg\"\n",
    "        ds['qv_sfc'].attrs[\"description\"] = \" \"\n",
    "        ds['qv_cb'].attrs[\"long_name\"] = \"Cloud Base Water Vapor Mixing Ratio\"\n",
    "        ds['qv_cb'].attrs[\"units\"] = \"kg/kg\"\n",
    "        ds['qv_cb'].attrs[\"description\"] = \"Using Environmental Water Vapor Estimate Matched to Cloud Base Height Estimate\"\n",
    "        ds['theta_sfc'].attrs[\"long_name\"] = \"Near Surface Potential Temperature\"\n",
    "        ds['theta_sfc'].attrs[\"units\"] = \"K\"\n",
    "        ds['theta_sfc'].attrs[\"description\"] = \" \"\n",
    "        ds['theta_cb'].attrs[\"long_name\"] = \"Cloud Base Potential Temperature\"\n",
    "        ds['theta_cb'].attrs[\"units\"] = \"K\"\n",
    "        ds['theta_cb'].attrs[\"description\"] = \"Matching temperature to cloud base height estimate\"\n",
    "        ds['z_lcl'].attrs[\"long_name\"] = \"Lifted Condensation Level\"\n",
    "        ds['z_lcl'].attrs[\"units\"] = \"m\"\n",
    "        ds['z_lcl'].attrs[\"description\"] = \"For well-mixed layer. Estimated Using Near Surface Temperature and Dew Point Difference\"\n",
    "        ds['LTS'].attrs[\"long_name\"] = \"Lower Tropospheric Stability\"\n",
    "        ds['LTS'].attrs[\"units\"] = \"K\"\n",
    "        ds['LTS'].attrs[\"description\"] = \" \"\n",
    "        ds['EIS'].attrs[\"long_name\"] = \"Estimated Inversion Strength\"\n",
    "        ds['EIS'].attrs[\"units\"] = \"K\"\n",
    "        ds['EIS'].attrs[\"description\"] = \"Estimated Using Method From Wood and Bretherton (2006, J. Climate)\"\n",
    "        ds['rh_700'].attrs[\"long_name\"] = \"700-hPa Relative Humidity\"\n",
    "        ds['rh_700'].attrs[\"units\"] = \"%\"\n",
    "        ds['rh_700'].attrs[\"description\"] = \" \"\n",
    "        ds['rh_750'].attrs[\"long_name\"] = \"750-hPa Relative Humidity\"\n",
    "        ds['rh_750'].attrs[\"units\"] = \"%\"\n",
    "        ds['rh_750'].attrs[\"description\"] = \" \"\n",
    "        ds['rh_800'].attrs[\"long_name\"] = \"800-hPa Relative Humidity\"\n",
    "        ds['rh_800'].attrs[\"units\"] = \"%\"\n",
    "        ds['rh_800'].attrs[\"description\"] = \" \"\n",
    "        ds['ccn1'].attrs[\"long_name\"] = \"0.1% Cloud Condensation Nuclei\"\n",
    "        ds['ccn1'].attrs[\"units\"] = \"cm-3\"\n",
    "        ds['ccn1'].attrs[\"description\"] = \"At the surface\"\n",
    "        ds['ccn2'].attrs[\"long_name\"] = \"0.2% Cloud Condensation Nuclei\"\n",
    "        ds['ccn2'].attrs[\"units\"] = \"cm-3\"\n",
    "        ds['ccn2'].attrs[\"description\"] = \"At the surface\"\n",
    "        ds['ccn5'].attrs[\"long_name\"] = \"0.5% Cloud Condensation Nuclei\"\n",
    "        ds['ccn5'].attrs[\"units\"] = \"cm-3\"\n",
    "        ds['ccn5'].attrs[\"description\"] = \"At the surface\"\n",
    "\n",
    "        ds.attrs[\"description\"] = \"LCL height is estimated from surface theta and dew point assuming a convectively mixed layer\"\n",
    "        ds.attrs[\"contact\"] = \"Adam Varble, adam.varble@pnnl.gov\"\n",
    "        ds.attrs[\"date\"] = \"9 November 2022\"\n",
    "        \n",
    "#        outfile = 'E3SM_1Deg_ENA_Environment_Properties_'+stryear[ii]+'-'+strmonth[jj]+'.nc'\n",
    "        outfile = 'E3SM_1Deg_ENA_3x3_Environment_Properties_'+stryear[ii]+'-'+strmonth[jj]+'.nc'\n",
    "        ds.to_netcdf(outfile, mode='w')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "MyEnv",
   "language": "python",
   "name": "myenv"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
