{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e13b606f-3d73-4d58-8216-8909a6abd059",
   "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": 1,
   "id": "5153111d-056c-4f70-9fa6-dbbf232ae8e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#set to ignore warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3ea59f3e-32d7-4453-ac2c-1b3037b3288d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#observations directory\n",
    "obsdir = '/global/cscratch1/sd/avarble/eagles/observations/ena/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "33b3f606-df7b-45ec-b753-556a1ef2fa6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#constants\n",
    "G = 9.8      #gravity\n",
    "Cp = 1005.7  #specific heat of air at constant pressure\n",
    "Rd = 287.    #dry air constant\n",
    "Rv = 461.    #moist air constant\n",
    "lv = 2.477e6 #latent heat of vaporization at 10 C\n",
    "epsilon = Rd/Rv\n",
    "pres_const = 85000. #Pa (used by Bennartz (2007, JGR)), could use cloud top pressure, but shouldn't alter estimates much\n",
    "Q = 2.   #scattering efficiency\n",
    "k = 0.74 #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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "fd11c216-19bb-441d-8391-656e58be269f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/global/homes/a/avarble/.conda/envs/myenv/lib/python3.7/site-packages/xarray/core/nanops.py:142: RuntimeWarning: Mean of empty slice\n",
      "  return np.nanmean(a, axis=axis, dtype=dtype)\n"
     ]
    }
   ],
   "source": [
    "#loop through 0.5-deg VISST satellite retrieval files, coarsen to 1 degree, and write to files for analyses\n",
    "files = glob.glob(os.path.join(obsdir, 'visst/grid/enavisstgridm*minnisX1.c1.*.cdf'))\n",
    "files2 = sorted(files)\n",
    "\n",
    "for file in files2:\n",
    "    data = xr.open_dataset(file)\n",
    "    tmp = file.split('.')[0]\n",
    "    strdataset = tmp.split('/')[-1]\n",
    "    strdate = file.split('.')[2]\n",
    "    strtime = file.split('.')[3]\n",
    "    base_time = data['base_time']\n",
    "    time_offset = data['time_offset']\n",
    "    time = data['time']\n",
    "\n",
    "    lon_1deg = data['longitude'].coarsen(lon=2, boundary=\"exact\").mean()\n",
    "    lat_1deg = data['latitude'].coarsen(lat=2, boundary=\"exact\").mean()\n",
    "\n",
    "    bb_lw_flux_1deg = data['broadband_longwave_flux'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    bb_sw_albedo_1deg = data['broadband_shortwave_albedo'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    cf_cldtype_1deg = data['cloud_percentage'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    cf_level_1deg = data['cloud_percentage_level'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    irtemp_1deg = data['ir_temperature'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    sfc_lw_down_1deg = data['surface_down_longwave_flux'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    sfc_sw_down_1deg = data['surface_down_shortwave_flux'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    sfc_lw_net_1deg = data['surface_net_longwave_flux'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    sfc_sw_net_1deg = data['surface_net_shortwave_flux'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    vis_reflectance_1deg = data['visible_reflectance'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    irtemp_clear_1deg = data['clearsky_ir_temperature'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    vis_reflectance_clear_1deg = data['clearsky_vis_reflectance'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    solar_zenith_1deg = data['solar_zenith_angle'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    viewing_zenith_1deg = data['viewing_zenith_angle'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    azimuth_1deg = data['azimuth_angle'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    scan_time_1deg = data['scan_time'].coarsen(lat=2, lon=2, boundary=\"exact\").mean()\n",
    "    \n",
    "    #cloud type and phase variables (weight by cloud fraction)\n",
    "    cf = data['cloud_percentage']\n",
    "    cf_sum_1deg = data['cloud_percentage'].coarsen(lat=2, lon=2, boundary=\"exact\").sum() # units of %\n",
    "\n",
    "    cbh = data['cloud_height_base']\n",
    "    cbh_sum = cbh*cf\n",
    "    cbh_sum_1deg = cbh_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cbh_1deg = cbh_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    cbh_sd = data['cloud_height_base_sd']\n",
    "    cbh_sd_sum = cbh_sd*cf\n",
    "    cbh_sd_sum_1deg = cbh_sd_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cbh_sd_1deg = cbh_sd_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    cch = data['cloud_height_center']\n",
    "    cch_sum = cch*cf\n",
    "    cch_sum_1deg = cch_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cch_1deg = cch_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    cch_sd = data['cloud_height_center_sd']\n",
    "    cch_sd_sum = cch_sd*cf\n",
    "    cch_sd_sum_1deg = cch_sd_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cch_sd_1deg = cch_sd_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    cth = data['cloud_height_top']\n",
    "    cth_sum = cth*cf\n",
    "    cth_sum_1deg = cth_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cth_1deg = cth_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    cth_sd = data['cloud_height_top_sd']\n",
    "    cth_sd_sum = cth_sd*cf\n",
    "    cth_sd_sum_1deg = cth_sd_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cth_sd_1deg = cth_sd_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    cbp = data['cloud_pressure_base']\n",
    "    cbp_sum = cbp*cf\n",
    "    cbp_sum_1deg = cbp_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cbp_1deg = cbp_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    cbp_sd = data['cloud_pressure_base_sd']\n",
    "    cbp_sd_sum = cbp_sd*cf\n",
    "    cbp_sd_sum_1deg = cbp_sd_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cbp_sd_1deg = cbp_sd_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    ccp = data['cloud_pressure_center']\n",
    "    ccp_sum = ccp*cf\n",
    "    ccp_sum_1deg = ccp_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    ccp_1deg = ccp_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    ccp_sd = data['cloud_pressure_center_sd']\n",
    "    ccp_sd_sum = ccp_sd*cf\n",
    "    ccp_sd_sum_1deg = ccp_sd_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    ccp_sd_1deg = ccp_sd_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    ctp = data['cloud_pressure_top']\n",
    "    ctp_sum = ctp*cf\n",
    "    ctp_sum_1deg = ctp_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    ctp_1deg = ctp_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    ctp_sd = data['cloud_pressure_top_sd']\n",
    "    ctp_sd_sum = ctp_sd*cf\n",
    "    ctp_sd_sum_1deg = ctp_sd_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    ctp_sd_1deg = ctp_sd_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    ct = data['cloud_temperature']\n",
    "    ct_sum = ct*cf\n",
    "    ct_sum_1deg = ct_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    ct_1deg = ct_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    ct_sd = data['cloud_temperature_sd']\n",
    "    ct_sd_sum = ct_sd*cf\n",
    "    ct_sd_sum_1deg = ct_sd_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    ct_sd_1deg = ct_sd_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    ir_emit = data['ir_emit']\n",
    "    ir_emit_sum = ir_emit*cf\n",
    "    ir_emit_sum_1deg = ir_emit_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    ir_emit_1deg = ir_emit_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    ir_emit_sd = data['ir_emit_sd']\n",
    "    ir_emit_sd_sum = ir_emit_sd*cf\n",
    "    ir_emit_sd_sum_1deg = ir_emit_sd_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    ir_emit_sd_1deg = ir_emit_sd_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    cod_lin = data['optical_depth_linear']\n",
    "    cod_lin_sum = cod_lin*cf\n",
    "    cod_lin_sum_1deg = cod_lin_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cod_lin_1deg = cod_lin_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    cod_lin_sd = data['optical_depth_linear_sd']\n",
    "    cod_lin_sd_sum = cod_lin_sd*cf\n",
    "    cod_lin_sd_sum_1deg = cod_lin_sd_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cod_lin_sd_1deg = cod_lin_sd_sum_1deg/cf_sum_1deg\n",
    "\n",
    "    cod_log = data['optical_depth_log']\n",
    "    cod_log_sum = cod_log*cf\n",
    "    cod_log_sum_1deg = cod_log_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cod_log_1deg = cod_log_sum_1deg/cf_sum_1deg \n",
    "\n",
    "    particle_size = data['particle_size']\n",
    "    particle_size_sum = particle_size*np.array(cf[:,:,:,1:])\n",
    "    particle_size_sum_1deg = particle_size_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    particle_size_1deg = particle_size_sum_1deg/np.array(cf_sum_1deg[:,:,:,1:])\n",
    "\n",
    "    particle_size_sd = data['particle_size_sd']\n",
    "    particle_size_sd_sum = particle_size_sd*np.array(cf[:,:,:,1:])\n",
    "    particle_size_sd_sum_1deg = particle_size_sd_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    particle_size_sd_1deg = particle_size_sd_sum_1deg/np.array(cf_sum_1deg[:,:,:,1:])\n",
    "\n",
    "    water_path = data['water_path']\n",
    "    water_path_sum = water_path*np.array(cf[:,:,:,1:])\n",
    "    water_path_sum_1deg = water_path_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    water_path_1deg = water_path_sum_1deg/np.array(cf_sum_1deg[:,:,:,1:])\n",
    "\n",
    "    water_path_sd = data['water_path_sd']\n",
    "    water_path_sd_sum = water_path_sd*np.array(cf[:,:,:,1:])\n",
    "    water_path_sd_sum_1deg = water_path_sd_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    water_path_sd_1deg = water_path_sd_sum_1deg/np.array(cf_sum_1deg[:,:,:,1:])\n",
    "    \n",
    "    #level variables (weight by cloud fraction)\n",
    "    cf_lev = data['cloud_percentage_level']\n",
    "    cf_lev_sum_1deg = data['cloud_percentage_level'].coarsen(lat=2, lon=2, boundary=\"exact\").sum() # units of %\n",
    "\n",
    "    cbh_level = data['cloud_height_base_level']\n",
    "    cbh_level_sum = cbh_level*cf_lev\n",
    "    cbh_level_sum_1deg = cbh_level_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cbh_level_1deg = cbh_level_sum_1deg/cf_lev_sum_1deg\n",
    "\n",
    "    cch_level = data['cloud_height_center_level']\n",
    "    cch_level_sum = cch_level*cf_lev\n",
    "    cch_level_sum_1deg = cch_level_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cch_level_1deg = cch_level_sum_1deg/cf_lev_sum_1deg\n",
    "\n",
    "    cth_level = data['cloud_height_top_level']\n",
    "    cth_level_sum = cth_level*cf_lev\n",
    "    cth_level_sum_1deg = cth_level_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cth_level_1deg = cth_level_sum_1deg/cf_lev_sum_1deg\n",
    "\n",
    "    cbp_level = data['cloud_pressure_base_level']\n",
    "    cbp_level_sum = cbp_level*cf_lev\n",
    "    cbp_level_sum_1deg = cbp_level_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cbp_level_1deg = cbp_level_sum_1deg/cf_lev_sum_1deg\n",
    "\n",
    "    ccp_level = data['cloud_pressure_center_level']\n",
    "    ccp_level_sum = ccp_level*cf_lev\n",
    "    ccp_level_sum_1deg = ccp_level_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    ccp_level_1deg = ccp_level_sum_1deg/cf_lev_sum_1deg\n",
    "\n",
    "    ctp_level = data['cloud_pressure_top_level']\n",
    "    ctp_level_sum = ctp_level*cf_lev\n",
    "    ctp_level_sum_1deg = ctp_level_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    ctp_level_1deg = ctp_level_sum_1deg/cf_lev_sum_1deg\n",
    "\n",
    "    cbt_level = data['cloud_temperature_base_level']\n",
    "    cbt_level_sum = cbt_level*cf_lev\n",
    "    cbt_level_sum_1deg = cbt_level_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cbt_level_1deg = cbt_level_sum_1deg/cf_lev_sum_1deg\n",
    "\n",
    "    cct_level = data['cloud_temperature_center_level']\n",
    "    cct_level_sum = cct_level*cf_lev\n",
    "    cct_level_sum_1deg = cct_level_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cct_level_1deg = cct_level_sum_1deg/cf_lev_sum_1deg\n",
    "\n",
    "    ctt_level = data['cloud_temperature_top_level']\n",
    "    ctt_level_sum = ctt_level*cf_lev\n",
    "    ctt_level_sum_1deg = ctt_level_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    ctt_level_1deg = ctt_level_sum_1deg/cf_lev_sum_1deg\n",
    "\n",
    "    cod_lin_level = data['optical_depth_linear_level']\n",
    "    cod_lin_level_sum = cod_lin_level*cf_lev\n",
    "    cod_lin_level_sum_1deg = cod_lin_level_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cod_lin_level_1deg = cod_lin_level_sum_1deg/cf_lev_sum_1deg\n",
    "\n",
    "    cod_log_level = data['optical_depth_log_level']\n",
    "    cod_log_level_sum = cod_log_level*cf_lev\n",
    "    cod_log_level_sum_1deg = cod_log_level_sum.coarsen(lat=2, lon=2, boundary=\"exact\").sum()\n",
    "    cod_log_level_1deg = cod_log_level_sum_1deg/cf_lev_sum_1deg\n",
    "    \n",
    "    #cloud droplet number concentration retrieval following Bennartz (2007, JGR)\n",
    "    iwp = water_path_1deg[:,:,:,0] #IWP\n",
    "    cod = cod_lin_1deg[:,:,:,2]    #cloud optical depth\n",
    "    lwp = water_path_1deg[:,:,:,1] #LWP (g/m2)\n",
    "    ctt = ct_1deg[:,:,:,2]         #cloud top temperature (K)\n",
    "    \n",
    "    rho_air = pres_const/(Rd*ctt) #dry air density\n",
    "    es = 611.2*np.exp(17.62*(ctt-273.15)/(243.12 + ctt - 273.15))\n",
    "    ws = epsilon*es/(pres_const - es)\n",
    "    gamma_w = G*((1 + lv*ws/(Rd*ctt))/(Cp + lv**2*ws*epsilon/(Rd*ctt**2))) #moist adiabatic lapse rate in cloud\n",
    "    gamma_ad = (((epsilon + ws)*ws*lv*gamma_w)/(Rd*ctt**2) - (G*ws*pres_const/(Rd*ctt*(pres_const - es))))*rho_air #adiabatic LWC lapse rate in cloud\n",
    "    H_1deg = (2.*1e-3*lwp/(0.8*gamma_ad))**0.5 #cloud depth with 80% adiabaticity\n",
    "    Nd_1deg = 1e-6*(cod**3/k)*((2*(1e-3*lwp))**(-2.5))*((0.6*np.pi*Q)**(-3))*((3./(4.*np.pi*rho_liq))**(-2))*((0.8*gamma_ad)**0.5) #cloud drop concentration (80% adiabatic)\n",
    "    H_ad_1deg = (2.*1e-3*lwp/(gamma_ad))**0.5 #cloud depth with 100% adiabaticity\n",
    "    Nd_ad_1deg = 1e-6*(cod**3/k)*((2*(1e-3*lwp))**(-2.5))*((0.6*np.pi*Q)**(-3))*((3./(4.*np.pi*rho_liq))**(-2))*((gamma_ad)**0.5) #adiabatic cloud drop concentration\n",
    "\n",
    "    #remove columns with ice or bad values\n",
    "    H_array_1deg = np.array(H_1deg)\n",
    "    H_ad_array_1deg = np.array(H_ad_1deg)\n",
    "    Nd_array_1deg = np.array(Nd_1deg)\n",
    "    Nd_ad_array_1deg = np.array(Nd_ad_1deg)    \n",
    "    ind = np.array(iwp > 0)\n",
    "    H_array_1deg[ind] = np.nan\n",
    "    H_ad_array_1deg[ind] = np.nan\n",
    "    Nd_array_1deg[ind] = np.nan\n",
    "    Nd_ad_array_1deg[ind] = np.nan    \n",
    "    ind = np.isinf(Nd_array_1deg)\n",
    "    H_array_1deg[ind] = np.nan\n",
    "    H_ad_array_1deg[ind] = np.nan\n",
    "    Nd_array_1deg[ind] = np.nan\n",
    "    Nd_ad_array_1deg[ind] = np.nan\n",
    "    \n",
    "    ds = xr.Dataset({'base_time': ('time', np.arange(1)), 'base_time': base_time,\n",
    "                     'time_offset': ('time', time_offset),\n",
    "                     'broadband_longwave_flux': (('time','lat','lon','scn_type'), np.float32(bb_lw_flux_1deg)),\n",
    "                     'broadband_shortwave_albedo': (('time','lat','lon','scn_type'), np.float32(bb_sw_albedo_1deg)),\n",
    "                     'cloud_percentage': (('time','lat','lon','cld_type'), np.float32(cf_cldtype_1deg)),\n",
    "                     'cloud_percentage_level': (('time','lat','lon','level'), np.float32(cf_level_1deg)),\n",
    "                     'ir_temperature': (('time','lat','lon','scn_type'), np.float32(irtemp_1deg)),\n",
    "                     'surface_down_longwave_flux': (('time','lat','lon'), np.float32(sfc_lw_down_1deg)),\n",
    "                     'surface_down_shortwave_flux': (('time','lat','lon'), np.float32(sfc_sw_down_1deg)),\n",
    "                     'surface_net_longwave_flux': (('time','lat','lon'), np.float32(sfc_lw_net_1deg)),\n",
    "                     'surface_net_shortwave_flux': (('time','lat','lon'), np.float32(sfc_sw_net_1deg)),\n",
    "                     'visible_reflectance': (('time','lat','lon','scn_type'), np.float32(vis_reflectance_1deg)),\n",
    "                     'clearsky_ir_temperature': (('time','lat','lon'), np.float32(irtemp_clear_1deg)),\n",
    "                     'clearsky_visible_reflectance': (('time','lat','lon'), np.float32(vis_reflectance_clear_1deg)),\n",
    "                     'solar_zenith_angle': (('time','lat','lon'), np.float32(solar_zenith_1deg)),\n",
    "                     'viewing_zenith_angle': (('time','lat','lon'), np.float32(viewing_zenith_1deg)),\n",
    "                     'azimuth_angle': (('time','lat','lon'), np.float32(azimuth_1deg)),\n",
    "                     'scan_time': (('time','lat','lon'), np.float32(scan_time_1deg)),\n",
    "                     'cloud_height_base': (('time','lat','lon','cld_type'), np.float32(cbh_1deg)),\n",
    "                     'cloud_height_base_sd': (('time','lat','lon','cld_type'), np.float32(cbh_sd_1deg)),\n",
    "                     'cloud_height_center': (('time','lat','lon','cld_type'), np.float32(cch_1deg)),\n",
    "                     'cloud_height_center_sd': (('time','lat','lon','cld_type'), np.float32(cch_sd_1deg)),\n",
    "                     'cloud_height_top': (('time','lat','lon','cld_type'), np.float32(cth_1deg)),\n",
    "                     'cloud_height_top_sd': (('time','lat','lon','cld_type'), np.float32(cth_sd_1deg)),\n",
    "                     'cloud_pressure_base': (('time','lat','lon','cld_type'), np.float32(cbp_1deg)),\n",
    "                     'cloud_pressure_base_sd': (('time','lat','lon','cld_type'), np.float32(cbp_sd_1deg)),\n",
    "                     'cloud_pressure_center': (('time','lat','lon','cld_type'), np.float32(ccp_1deg)),\n",
    "                     'cloud_pressure_center_sd': (('time','lat','lon','cld_type'), np.float32(ccp_sd_1deg)),\n",
    "                     'cloud_pressure_top': (('time','lat','lon','cld_type'), np.float32(ctp_1deg)),\n",
    "                     'cloud_pressure_top_sd': (('time','lat','lon','cld_type'), np.float32(ctp_sd_1deg)),\n",
    "                     'cloud_temperature': (('time','lat','lon','cld_type'), np.float32(ct_1deg)),\n",
    "                     'cloud_temperature_sd': (('time','lat','lon','cld_type'), np.float32(ct_sd_1deg)),\n",
    "                     'ir_emit': (('time','lat','lon','cld_type'), np.float32(ir_emit_1deg)),\n",
    "                     'ir_emit_sd': (('time','lat','lon','cld_type'), np.float32(ir_emit_sd_1deg)),\n",
    "                     'optical_depth_linear': (('time','lat','lon','cld_type'), np.float32(cod_lin_1deg)),\n",
    "                     'optical_depth_linear_sd': (('time','lat','lon','cld_type'), np.float32(cod_lin_sd_1deg)),\n",
    "                     'optical_depth_log': (('time','lat','lon','cld_type'), np.float32(cod_log_1deg)),\n",
    "                     'particle_size': (('time','lat','lon','cld_phase'), np.float32(particle_size_1deg)),\n",
    "                     'particle_size_sd': (('time','lat','lon','cld_phase'), np.float32(particle_size_sd_1deg)),\n",
    "                     'water_path': (('time','lat','lon','cld_phase'), np.float32(water_path_1deg)),\n",
    "                     'water_path_sd': (('time','lat','lon','cld_phase'), np.float32(water_path_sd_1deg)),\n",
    "                     'cloud_height_base_level': (('time','lat','lon','level'), np.float32(cbh_level_1deg)),\n",
    "                     'cloud_height_center_level': (('time','lat','lon','level'), np.float32(cch_level_1deg)),\n",
    "                     'cloud_height_top_level': (('time','lat','lon','level'), np.float32(cth_level_1deg)),\n",
    "                     'cloud_pressure_base_level': (('time','lat','lon','level'), np.float32(cbp_level_1deg)),\n",
    "                     'cloud_pressure_center_level': (('time','lat','lon','level'), np.float32(ccp_level_1deg)),\n",
    "                     'cloud_pressure_top_level': (('time','lat','lon','level'), np.float32(ctp_level_1deg)),\n",
    "                     'cloud_temperature_base_level': (('time','lat','lon','level'), np.float32(cbt_level_1deg)),\n",
    "                     'cloud_temperature_center_level': (('time','lat','lon','level'), np.float32(cct_level_1deg)),\n",
    "                     'cloud_temperature_top_level': (('time','lat','lon','level'), np.float32(ctt_level_1deg)),\n",
    "                     'optical_depth_linear_level': (('time','lat','lon','level'), np.float32(cod_lin_level_1deg)),\n",
    "                     'optical_depth_log_level': (('time','lat','lon','level'), np.float32(cod_log_level_1deg)),\n",
    "                     'CDNC': (('time','lat','lon'), np.float32(Nd_array_1deg)),\n",
    "                     'CDNC_adiabatic': (('time','lat','lon'), np.float32(Nd_ad_array_1deg)),\n",
    "                     'H': (('time','lat','lon'), np.float32(H_array_1deg)),\n",
    "                     'H_adiabatic': (('time','lat','lon'), np.float32(H_ad_array_1deg))},\n",
    "                     coords={'time': ('time', time),\n",
    "                             'latitude': ('lat', np.float32(lat_1deg)),\n",
    "                             'longitude': ('lon', np.float32(lon_1deg)),\n",
    "                             'scan_type': ('scn_type', ['all sky', 'clear (cloudy = total - clear)']),\n",
    "                             'cloud_type': ('cld_type', ['all','ice','liquid','supercooled liquid']),\n",
    "                             'cloud_phase': ('cld_phase', ['ice','liquid','supercooled liquid']),\n",
    "                             'altitude_level': ('level', ['total','low (0-2 km)','mid (2-6 km)','high (> 6 km)'])})\n",
    "\n",
    "    ds['base_time'].attrs[\"ancillary_variables\"] = \"time_offset\"\n",
    "    ds['time_offset'].attrs[\"long_name\"] = \"Time offset from base_time\"\n",
    "    ds['time_offset'].attrs[\"ancillary_variables\"] = \"base_time\"\n",
    "    ds['longitude'].attrs[\"long_name\"] = \"Longitude\"\n",
    "    ds['longitude'].attrs[\"units\"] = \"degrees\"\n",
    "    ds['latitude'].attrs[\"long_name\"] = \"Latitude\"\n",
    "    ds['latitude'].attrs[\"units\"] = \"degrees\"\n",
    "    ds['broadband_longwave_flux'].attrs[\"long_name\"] = \"Broadband TOA Upwelling Longwave Flux\"\n",
    "    ds['broadband_longwave_flux'].attrs[\"units\"] = \"W m-2\"\n",
    "    ds['broadband_longwave_flux'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['broadband_shortwave_albedo'].attrs[\"long_name\"] = \"Broadband TOA Shortwave Albedo\"\n",
    "    ds['broadband_shortwave_albedo'].attrs[\"units\"] = \"%\"\n",
    "    ds['broadband_shortwave_albedo'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_percentage'].attrs[\"long_name\"] = \"Cloud Fractional Coverage\"\n",
    "    ds['cloud_percentage'].attrs[\"units\"] = \"%\"\n",
    "    ds['cloud_percentage'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_percentage_level'].attrs[\"long_name\"] = \"Cloud Fractional Coverage\"\n",
    "    ds['cloud_percentage_level'].attrs[\"units\"] = \"%\"\n",
    "    ds['cloud_percentage_level'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['ir_temperature'].attrs[\"long_name\"] = \"Infrared Temperature\"\n",
    "    ds['ir_temperature'].attrs[\"units\"] = \"K\"\n",
    "    ds['ir_temperature'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['surface_down_longwave_flux'].attrs[\"long_name\"] = \"Estimated Surface Downwelling Longwave Flux\"\n",
    "    ds['surface_down_longwave_flux'].attrs[\"units\"] = \"W m-2\"\n",
    "    ds['surface_down_longwave_flux'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['surface_down_shortwave_flux'].attrs[\"long_name\"] = \"Estimated Surface Downwelling Shortwave Flux\"\n",
    "    ds['surface_down_shortwave_flux'].attrs[\"units\"] = \"W m-2\"\n",
    "    ds['surface_down_shortwave_flux'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['surface_net_longwave_flux'].attrs[\"long_name\"] = \"Estimated Surface Net Longwave Flux\"\n",
    "    ds['surface_net_longwave_flux'].attrs[\"units\"] = \"W m-2\"\n",
    "    ds['surface_net_longwave_flux'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['surface_net_shortwave_flux'].attrs[\"long_name\"] = \"Estimated Surface Net Shortwave Flux\"\n",
    "    ds['surface_net_shortwave_flux'].attrs[\"units\"] = \"W m-2\"\n",
    "    ds['surface_net_shortwave_flux'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['visible_reflectance'].attrs[\"long_name\"] = \"Visible Reflectance\"\n",
    "    ds['visible_reflectance'].attrs[\"units\"] = \"Fraction\"\n",
    "    ds['visible_reflectance'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['clearsky_ir_temperature'].attrs[\"long_name\"] = \"Clear Sky Infrared Temperature\"\n",
    "    ds['clearsky_ir_temperature'].attrs[\"units\"] = \"K\"\n",
    "    ds['clearsky_ir_temperature'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['clearsky_visible_reflectance'].attrs[\"long_name\"] = \"Clear Sky Visible Reflectance\"\n",
    "    ds['clearsky_visible_reflectance'].attrs[\"units\"] = \"Fraction\"\n",
    "    ds['clearsky_visible_reflectance'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['solar_zenith_angle'].attrs[\"long_name\"] = \"Solar Zenith Angle\"\n",
    "    ds['solar_zenith_angle'].attrs[\"units\"] = \"Degrees\"\n",
    "    ds['solar_zenith_angle'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['viewing_zenith_angle'].attrs[\"long_name\"] = \"Viewing Zenith Angle\"\n",
    "    ds['viewing_zenith_angle'].attrs[\"units\"] = \"Degrees\"\n",
    "    ds['viewing_zenith_angle'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['azimuth_angle'].attrs[\"long_name\"] = \"Azimuth Angle\"\n",
    "    ds['azimuth_angle'].attrs[\"units\"] = \"Degrees\"\n",
    "    ds['azimuth_angle'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['scan_time'].attrs[\"long_name\"] = \"Average Scan Time\"\n",
    "    ds['scan_time'].attrs[\"units\"] = \"Hour\"\n",
    "    ds['scan_time'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_height_base'].attrs[\"long_name\"] = \"Cloud Base Height\"\n",
    "    ds['cloud_height_base'].attrs[\"units\"] = \"km\"\n",
    "    ds['cloud_height_base'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_height_base_sd'].attrs[\"long_name\"] = \"Standard Deviation of Cloud Base Height\"\n",
    "    ds['cloud_height_base_sd'].attrs[\"units\"] = \"km\"\n",
    "    ds['cloud_height_base_sd'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_height_center'].attrs[\"long_name\"] = \"Cloud Center Height\"\n",
    "    ds['cloud_height_center'].attrs[\"units\"] = \"km\"\n",
    "    ds['cloud_height_center'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_height_center_sd'].attrs[\"long_name\"] = \"Standard Deviation of Cloud Center Height\"\n",
    "    ds['cloud_height_center_sd'].attrs[\"units\"] = \"km\"\n",
    "    ds['cloud_height_center_sd'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_height_top'].attrs[\"long_name\"] = \"Cloud Top Height\"\n",
    "    ds['cloud_height_top'].attrs[\"units\"] = \"km\"\n",
    "    ds['cloud_height_top'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_height_top_sd'].attrs[\"long_name\"] = \"Standard Deviation of Cloud Top Height\"\n",
    "    ds['cloud_height_top_sd'].attrs[\"units\"] = \"km\"\n",
    "    ds['cloud_height_top_sd'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_pressure_base'].attrs[\"long_name\"] = \"Cloud Base Pressure\"\n",
    "    ds['cloud_pressure_base'].attrs[\"units\"] = \"hPa\"\n",
    "    ds['cloud_pressure_base'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_pressure_base_sd'].attrs[\"long_name\"] = \"Standard Deviation of Cloud Base Pressure\"\n",
    "    ds['cloud_pressure_base_sd'].attrs[\"units\"] = \"hPa\"\n",
    "    ds['cloud_pressure_base_sd'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_pressure_center'].attrs[\"long_name\"] = \"Cloud Center Pressure\"\n",
    "    ds['cloud_pressure_center'].attrs[\"units\"] = \"hPa\"\n",
    "    ds['cloud_pressure_center'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_pressure_center_sd'].attrs[\"long_name\"] = \"Standard Deviation of Cloud Center Pressure\"\n",
    "    ds['cloud_pressure_center_sd'].attrs[\"units\"] = \"hPa\"\n",
    "    ds['cloud_pressure_center_sd'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_pressure_top'].attrs[\"long_name\"] = \"Cloud Top Pressure\"\n",
    "    ds['cloud_pressure_top'].attrs[\"units\"] = \"hPa\"\n",
    "    ds['cloud_pressure_top'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_pressure_top_sd'].attrs[\"long_name\"] = \"Standard Deviation of Cloud Top Pressure\"\n",
    "    ds['cloud_pressure_top_sd'].attrs[\"units\"] = \"hPa\"\n",
    "    ds['cloud_pressure_top_sd'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_temperature'].attrs[\"long_name\"] = \"Cloud Temperature\"\n",
    "    ds['cloud_temperature'].attrs[\"units\"] = \"K\"\n",
    "    ds['cloud_temperature'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_temperature_sd'].attrs[\"long_name\"] = \"Standard Deviation of Cloud Temperature\"\n",
    "    ds['cloud_temperature_sd'].attrs[\"units\"] = \"K\"\n",
    "    ds['cloud_temperature_sd'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['ir_emit'].attrs[\"long_name\"] = \"Average Infrared Emissivity\"\n",
    "    ds['ir_emit'].attrs[\"units\"] = \"unitless\"\n",
    "    ds['ir_emit'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['ir_emit_sd'].attrs[\"long_name\"] = \"Standard Deviation of Infrared Emissivity\"\n",
    "    ds['ir_emit_sd'].attrs[\"units\"] = \"unitless\"\n",
    "    ds['ir_emit_sd'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['optical_depth_linear'].attrs[\"long_name\"] = \"Linear Average of Cloud Optical Depth\"\n",
    "    ds['optical_depth_linear'].attrs[\"units\"] = \"unitless\"\n",
    "    ds['optical_depth_linear'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['optical_depth_linear_sd'].attrs[\"long_name\"] = \"Standard Deviation of Cloud Optical Depth\"\n",
    "    ds['optical_depth_linear_sd'].attrs[\"units\"] = \"unitless\"\n",
    "    ds['optical_depth_linear_sd'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['optical_depth_log'].attrs[\"long_name\"] = \"Log Average of Cloud Optical Depth\"\n",
    "    ds['optical_depth_log'].attrs[\"units\"] = \"unitless\"\n",
    "    ds['optical_depth_log'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['particle_size'].attrs[\"long_name\"] = \"Average Effective Particle Size\"\n",
    "    ds['particle_size'].attrs[\"units\"] = \"um\"\n",
    "    ds['particle_size'].attrs[\"description\"] = \"From VISST gridded product. Radius for liquid and diameter for ice.\"\n",
    "    ds['particle_size_sd'].attrs[\"long_name\"] = \"Standard Deviation of Effective Particle Size\"\n",
    "    ds['particle_size_sd'].attrs[\"units\"] = \"um\"\n",
    "    ds['particle_size_sd'].attrs[\"description\"] = \"From VISST gridded product. Radius for liquid and diameter for ice.\"\n",
    "    ds['water_path'].attrs[\"long_name\"] = \"Cloud Water Path\"\n",
    "    ds['water_path'].attrs[\"units\"] = \"g m-2\"\n",
    "    ds['water_path'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['water_path_sd'].attrs[\"long_name\"] = \"Standard Deviation of Cloud Water Path\"\n",
    "    ds['water_path_sd'].attrs[\"units\"] = \"g m-2\"\n",
    "    ds['water_path_sd'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_height_base_level'].attrs[\"long_name\"] = \"Cloud Base Height for Select Layers\"\n",
    "    ds['cloud_height_base_level'].attrs[\"units\"] = \"km\"\n",
    "    ds['cloud_height_base_level'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_height_center_level'].attrs[\"long_name\"] = \"Cloud Center Height for Select Layers\"\n",
    "    ds['cloud_height_center_level'].attrs[\"units\"] = \"km\"\n",
    "    ds['cloud_height_center_level'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_height_top_level'].attrs[\"long_name\"] = \"Cloud Top Height for Select Layers\"\n",
    "    ds['cloud_height_top_level'].attrs[\"units\"] = \"km\"\n",
    "    ds['cloud_height_top_level'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_pressure_base_level'].attrs[\"long_name\"] = \"Cloud Base Pressure for Select Layers\"\n",
    "    ds['cloud_pressure_base_level'].attrs[\"units\"] = \"hPa\"\n",
    "    ds['cloud_pressure_base_level'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_pressure_center_level'].attrs[\"long_name\"] = \"Cloud Center Pressure for Select Layers\"\n",
    "    ds['cloud_pressure_center_level'].attrs[\"units\"] = \"hPa\"\n",
    "    ds['cloud_pressure_center_level'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_pressure_top_level'].attrs[\"long_name\"] = \"Cloud Top Pressure for Select Layers\"\n",
    "    ds['cloud_pressure_top_level'].attrs[\"units\"] = \"hPa\"\n",
    "    ds['cloud_pressure_top_level'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_temperature_base_level'].attrs[\"long_name\"] = \"Cloud Base Temperature for Select Layers\"\n",
    "    ds['cloud_temperature_base_level'].attrs[\"units\"] = \"K\"\n",
    "    ds['cloud_temperature_base_level'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_temperature_center_level'].attrs[\"long_name\"] = \"Cloud Center Temperature for Select Layers\"\n",
    "    ds['cloud_temperature_center_level'].attrs[\"units\"] = \"K\"\n",
    "    ds['cloud_temperature_center_level'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['cloud_temperature_top_level'].attrs[\"long_name\"] = \"Cloud Top Temperature for Select Layers\"\n",
    "    ds['cloud_temperature_top_level'].attrs[\"units\"] = \"K\"\n",
    "    ds['cloud_temperature_top_level'].attrs[\"description\"] = \"From VISST gridded product\"\n",
    "    ds['optical_depth_linear_level'].attrs[\"long_name\"] = \"Cloud Optical Depth Linear Average for Select Layers\"\n",
    "    ds['optical_depth_linear_level'].attrs[\"units\"] = \"unitless\"\n",
    "    ds['optical_depth_linear_level'].attrs[\"description\"] = \"From VISST gridded product\"    \n",
    "    ds['optical_depth_log_level'].attrs[\"long_name\"] = \"Cloud Optical Depth Logarithmic Average for Select Layers\"\n",
    "    ds['optical_depth_log_level'].attrs[\"units\"] = \"unitless\"\n",
    "    ds['optical_depth_log_level'].attrs[\"description\"] = \"From VISST gridded product\"    \n",
    "    ds['CDNC'].attrs[\"long_name\"] = \"Cloud Droplet Concentration\"\n",
    "    ds['CDNC'].attrs[\"units\"] = \"cm-3\"\n",
    "    ds['CDNC'].attrs[\"description\"] = \"Retrieved following Bennartz using VISST product Meteosat data assuming adiabaticity = 0.8.\"\n",
    "    ds['CDNC_adiabatic'].attrs[\"long_name\"] = \"Adiabatic Cloud Droplet Concentration\"\n",
    "    ds['CDNC_adiabatic'].attrs[\"units\"] = \"cm-3\"\n",
    "    ds['CDNC_adiabatic'].attrs[\"description\"] = \"Retrieved following Bennartz using VISST product Meteosat data assuming adiabaticity = 1.\"\n",
    "    ds['H'].attrs[\"long_name\"] = \"Cloud Depth\"\n",
    "    ds['H'].attrs[\"units\"] = \"m\"\n",
    "    ds['H'].attrs[\"description\"] = \"Estimated using satellite-retrieved LWP and moist adiabatic lapse rate with adiabaticity = 0.8.\"\n",
    "    ds['H_adiabatic'].attrs[\"long_name\"] = \"Adiabatic Cloud Depth\"\n",
    "    ds['H_adiabatic'].attrs[\"units\"] = \"m\"\n",
    "    ds['H_adiabatic'].attrs[\"description\"] = \"Estimated using satellite-retrieved LWP and moist adiabatic lapse rate with adiabaticity = 1.\"\n",
    "\n",
    "    ds.attrs[\"description\"] = \"These data are averaged to 1 degree from the pre-existing 0.5-degree VISST grid files with added CDNC and cloud depth retrievals. Cloud droplet concentration retrievals are only valid in overcast single layer liquid cloud situations with sufficiently high solar angles. The retrieval follows Bennartz with an adiabaticity of 0.8, except k is set to 0.74 rather than 0.8. The adiabatic retrievals assume adiabaticity = 1.\"\n",
    "    ds.attrs[\"contact\"] = \"Adam Varble, adam.varble@pnnl.gov\"\n",
    "    ds.attrs[\"date\"] = \"3 January 2022\"\n",
    "\n",
    "    outfile = strdataset+'.c1.cdnc.1deg.'+strdate+'.cdf'\n",
    "    ds.to_netcdf(obsdir+'visst/1deg/'+outfile, mode='w')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3af8d9d-dd3e-43ba-a6e0-0df820331e1f",
   "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": 5
}
