#!/usr/bin/env python
# coding: utf-8

# In[1]:


import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma
from scipy.optimize import curve_fit
import pandas as pd
import csv
import os
import matplotlib.dates as mdates
import matplotlib.cm as cm
from matplotlib.colors import Normalize
from colorspacious import cspace_converter
import matplotlib as mpl
import pandas as pd
import netCDF4 as nc
from sklearn import linear_model
import statistics
from numpy import ma
import pandas as pd
import scipy.stats as stats
import seaborn as sns
import statsmodels.api as sm
from statsmodels.formula.api import ols
from sklearn.linear_model import LinearRegression

# Standardize the data (important for PCA)
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()

from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression


# In[2]:


fold = 'newrad2/'
desc = 'pafog_figs/newrad2/'
df_other =  pd.read_csv(fold + 'other.csv')
df_sat = pd.read_csv(fold + 'sat.csv')
df_sat = df_sat.reset_index(drop=True)
df_tb_ten_xm2 = pd.read_csv(fold + 'tb_ten_xm2.csv')
df_tb_xm23 = pd.read_csv(fold + 'tb_xm23.csv')
df_ave_met = pd.read_csv(fold + 'ave_met.csv')
df_ave_rad = pd.read_csv(fold + 'ave_rad.csv')
df_ave_ten_xm2 = pd.read_csv(fold + 'ave_ten_xm2.csv')
df_ave_ten_xm3 = pd.read_csv(fold + 'ave_ten_xm3.csv')
df_belave_xm23 = pd.read_csv(fold + 'belave_xm23.csv')
df_incldave_xm23 = pd.read_csv(fold + 'incldave_xm23.csv')
df_sum_met = pd.read_csv(fold + 'sum_met.csv')
df_sum_rad = pd.read_csv(fold + 'sum_rad.csv')
df_sum_ten_xm2 = pd.read_csv(fold + 'sum_ten_xm2.csv')
df_sum_ten_xm3 = pd.read_csv(fold + 'sum_ten_xm3.csv')
df_sum_xm23 = pd.read_csv(fold + 'sum_xm23.csv')
df_tb_ten_xm3 = pd.read_csv(fold + 'tb_ten_xm3.csv')
df_tb_met = pd.read_csv(fold + 'tb_met.csv')
df_tb_rad = pd.read_csv(fold + 'tb_rad.csv')
df_col = pd.read_csv(fold + 'col.csv')
df_col = df_col.reset_index(drop=True)
df_sed = pd.read_csv(fold + 'sed.csv')
df_sed = df_sed.reset_index(drop=True)

df_inputs = pd.read_csv(fold + 'bin_scheme_inputs.csv')

df_aall = pd.concat([df_inputs, df_other, df_sat, df_tb_ten_xm2, df_tb_xm23, df_ave_rad, df_ave_met, df_ave_ten_xm2, df_ave_ten_xm3, df_belave_xm23, df_incldave_xm23, df_sum_xm23, df_sum_met, df_sum_rad, df_sum_ten_xm3, df_sum_ten_xm2, df_tb_rad, df_tb_met, df_tb_ten_xm3, df_col], axis = 1).reindex(df_inputs.index)
df_aall['cap_f'] = df_aall['xm_col_f'] / df_aall['CRH_f']
df_aall['bl_cap_f'] = df_aall['xm_bl_f'] / df_aall['CRH_bl_f']
df_aall['xm_FT_f'] = df_aall['xm_col_f'] - df_aall['xm_bl_f']
df_aall['CRH_FT_f'] = (df_aall['xm_col_f'] - df_aall['xm_bl_f']) / ((df_aall['xm_col_f'] / df_aall['CRH_f']) - (df_aall['xm_bl_f'] / df_aall['CRH_bl_f']))
df_aall['FT_cap_f'] = df_aall['xm_FT_f'] / df_aall['CRH_FT_f']
df_aall['entrainment_z'] = df_aall['delta_z_blt'] / 86.4 - df_aall['w_blt'] * 1000
df_aall['entrainment_zh'] = df_aall['dz_blt_h'] / 86.4 - df_aall['w_blt_h'] * 1000
df_aall['duration'] = df_aall['diss_hr'] - df_aall['onset_hr']
df_aall['t_dur'] = (df_aall['t_diss_i'] - df_aall['t_onset_i']) / 12
df_aall['nl_dur'] = (df_aall['noliq_diss_i'] - df_aall['noliq_onset_i']) / 12
df_aall['ent_xm_ave2'] = df_aall['ent_delt_ave'] * df_aall['entrainment_z'] / 1000
df_aall['ent_xm_ave_h2'] = df_aall['ent_delt_ave_h'] * df_aall['entrainment_zh'] / 1000
df_aall['ent_sed'] = df_aall['ent_xm_ave2'] * (86400) - df_aall['sedt_f']
df_aall['ent_sed_h'] = df_aall['ent_xm_ave_h2'] * (86400 / 2) - (df_aall['sedt_h'])
df_all = df_aall.drop(columns = ['Unnamed: 0', 'bindir', 'inpdir', 'outdir'])
df_out = pd.concat([df_other, df_sat, df_tb_ten_xm2, df_tb_xm23, df_ave_rad, df_ave_met, df_ave_ten_xm2, df_ave_ten_xm3, df_belave_xm23, df_incldave_xm23, df_sum_xm23, df_sum_met, df_sum_rad, df_sum_ten_xm3, df_sum_ten_xm2, df_tb_rad, df_tb_met, df_tb_ten_xm3, df_col], axis = 1).reindex(df_inputs.index)

fog = df_all['onset_hr'] > 0.0
cond_outlier = df_all['diss_hr'] < 24.0
df_fog = df_all[fog]
df_nout = df_all[cond_outlier]
print(np.shape(df_fog), np.shape(df_nout))
# cond_outlier = df_all['entrainment_z'] > 1.5
cond_adzblt = df_all['blt_delta_abs'] > -1.
cond_constant = df_all['blt_delta_abs'] == 0.
df_constant = df_inputs[cond_constant]
cond_rev = df_all['blt_delta_abs'] > np.abs(df_all['delta_z_blt'])

cldN_list = sorted(df_inputs['cldN'].unique())
dts_list = sorted(df_inputs['dtsurf'].unique())
wmax_list = sorted(df_inputs['wmax'].unique())
Ug_list = sorted(df_inputs['Ug'].unique())

dtsurf__2  = df_all['dtsurf'] == dts_list[0]
dtsurf__1  = df_all['dtsurf'] == dts_list[1]
dtsurf_1  = df_all['dtsurf'] == dts_list[3]
dtsurf_2  = df_all['dtsurf'] == dts_list[4]
base_dtsurf = df_all['dtsurf'] == dts_list[2]
high_dtsurf = df_all['dtsurf'] > dts_list[2]
low_dtsurf = df_all['dtsurf'] < dts_list[2]

print(np.mean(df_all[low_dtsurf]['dtsurf']))



base_Ug = df_all['Ug'] == Ug_list[2]
high_Ug = df_all['Ug'] > Ug_list[2]
low_Ug = df_all['Ug'] < Ug_list[2]
Ug_5 = df_all['Ug'] == Ug_list[0]
Ug_7 = df_all['Ug'] == Ug_list[1]
Ug_12 = df_all['Ug'] == Ug_list[3]
Ug_15 = df_all['Ug'] == Ug_list[4]

base_wmax = df_all['wmax'] == wmax_list[2]
high_wmax = df_all['wmax'] > wmax_list[2]
low_wmax = df_all['wmax'] < wmax_list[2]
wmax_35 = df_all['wmax'] == wmax_list[0]
wmax_325 = df_all['wmax'] == wmax_list[1]
wmax_275 = df_all['wmax'] == wmax_list[3]
wmax_25 = df_all['wmax'] == wmax_list[4]


# In[3]:


fogprop_cldN = np.zeros((289, 10))
fog_cldN = np.zeros(np.shape(cldN_list))
fogdur_cldN = np.zeros(np.shape(cldN_list))
fogprop_a_cldN = np.zeros((289, 10))
fog_a_cldN = np.zeros(np.shape(cldN_list))
fogdur_a_cldN = np.zeros(np.shape(cldN_list))

for j in range(np.shape(fogprop_cldN)[1]):
    cldN_cond = df_all['cldN'] == cldN_list[j]
    df_cldN_fog = df_all[cldN_cond & fog & cond_outlier]
    df_cond = df_all[cldN_cond & cond_outlier]
    df_cldN_a_fog = df_all[cldN_cond & fog]
    fog_a_cldN[j] = np.shape(df_cldN_a_fog)[0] / 125.
    fog_cldN[j] = np.shape(df_cldN_fog)[0] / np.shape(df_cond)[0]
    fogdur_a_cldN[j] = np.mean(df_cldN_a_fog['diss_hr']) - np.mean(df_cldN_a_fog['onset_hr'])
    fogdur_cldN[j] = np.mean(df_cldN_fog['diss_hr']) - np.mean(df_cldN_fog['onset_hr'])
    for i in range(np.shape(fogprop_cldN)[0]):
        on = df_all['onset_hr'] <= i / 12
        diss = df_all['diss_hr'] <= i / 12
        fognum = df_all[on & ~diss & cldN_cond & cond_outlier].shape[0]
        totnum = df_all[cldN_cond & cond_outlier].shape[0]
        fognum_a = df_all[on & ~diss & cldN_cond].shape[0]
        totnum_a = df_all[cldN_cond].shape[0]
        fogprop_a_cldN[i,j] = fognum_a/totnum_a
        fogprop_cldN[i,j] = fognum/totnum
        

fogprop_a_cldN_low_Ug = np.zeros((289, 10))
fog_a_cldN_low_Ug = np.zeros(np.shape(cldN_list))
fogdur_a_cldN_low_Ug = np.zeros(np.shape(cldN_list))
fogprop_cldN_low_Ug = np.zeros((289, 10))
fog_cldN_low_Ug = np.zeros(np.shape(cldN_list))
fogdur_cldN_low_Ug = np.zeros(np.shape(cldN_list))

for j in range(np.shape(fogprop_cldN)[1]):
    cldN_cond = df_all['cldN'] == cldN_list[j]
    df_cldN_fog = df_all[cldN_cond & fog & low_Ug & cond_outlier]
    df_cond = df_all[cldN_cond & low_Ug & cond_outlier]
    df_cldN_a_fog = df_all[cldN_cond & fog & low_Ug]
    fog_a_cldN_low_Ug[j] = np.shape(df_cldN_a_fog)[0] / 50.
    fog_cldN_low_Ug[j] = np.shape(df_cldN_fog)[0] / np.shape(df_cond)[0]
    fogdur_a_cldN_low_Ug[j] = np.mean(df_cldN_a_fog['diss_hr']) - np.mean(df_cldN_a_fog['onset_hr'])
    fogdur_cldN_low_Ug[j] = np.mean(df_cldN_fog['diss_hr']) - np.mean(df_cldN_fog['onset_hr'])
    for i in range(np.shape(fogprop_cldN)[0]):
        on = df_all['onset_hr'] <= i / 12
        diss = df_all['diss_hr'] <= i / 12
        fognum = df_all[on & ~diss & cldN_cond & low_Ug & cond_outlier].shape[0]
        totnum = df_all[cldN_cond & low_Ug & cond_outlier].shape[0]
        fognum_a = df_all[on & ~diss & cldN_cond & low_Ug].shape[0]
        totnum_a = df_all[cldN_cond & low_Ug].shape[0]
        if totnum_a == 0:
            fogprop_a_cldN_low_Ug[i,j] = float("NaN")
        else:
            fogprop_a_cldN_low_Ug[i,j] = fognum_a/totnum_a
        if totnum == 0:
            fogprop_cldN_low_Ug[i,j] = float("NaN")
        else:
            fogprop_cldN_low_Ug[i,j] = fognum/totnum
        

fogprop_a_cldN_high_Ug = np.zeros((289, 10))
fog_a_cldN_high_Ug = np.zeros(np.shape(cldN_list))
fogdur_a_cldN_high_Ug = np.zeros(np.shape(cldN_list))
fogprop_cldN_high_Ug = np.zeros((289, 10))
fog_cldN_high_Ug = np.zeros(np.shape(cldN_list))
fogdur_cldN_high_Ug = np.zeros(np.shape(cldN_list))

for j in range(np.shape(fogprop_cldN)[1]):
    cldN_cond = df_all['cldN'] == cldN_list[j]
    df_cldN_fog = df_all[cldN_cond & fog & high_Ug & cond_outlier]
    df_cond = df_all[cldN_cond & high_Ug & cond_outlier]
    df_cldN_a_fog = df_all[cldN_cond & fog & high_Ug]
    fog_a_cldN_high_Ug[j] = np.shape(df_cldN_a_fog)[0] / 50.
    fog_cldN_high_Ug[j] = np.shape(df_cldN_fog)[0] / np.shape(df_cond)[0]
    fogdur_a_cldN_high_Ug[j] = np.mean(df_cldN_a_fog['diss_hr']) - np.mean(df_cldN_a_fog['onset_hr'])
    fogdur_cldN_high_Ug[j] = np.mean(df_cldN_fog['diss_hr']) - np.mean(df_cldN_fog['onset_hr'])
    for i in range(np.shape(fogprop_cldN)[0]):
        on = df_all['onset_hr'] <= i / 12
        diss = df_all['diss_hr'] <= i / 12
        fognum = df_all[on & ~diss & cldN_cond & high_Ug & cond_outlier].shape[0]
        totnum = df_all[cldN_cond & high_Ug & cond_outlier].shape[0]
        fognum_a = df_all[on & ~diss & cldN_cond & high_Ug].shape[0]
        totnum_a = df_all[cldN_cond & high_Ug].shape[0]
        if totnum_a == 0:
            fogprop_a_cldN_high_Ug[i,j] = float("NaN")
        else:
            fogprop_a_cldN_high_Ug[i,j] = fognum_a/totnum_a
        if totnum == 0:
            fogprop_cldN_high_Ug[i,j] = float("NaN")
        else:
            fogprop_cldN_high_Ug[i,j] = fognum/totnum
        

fogprop_a_cldN_low_wmax = np.zeros((289, 10))
fog_a_cldN_low_wmax = np.zeros(np.shape(cldN_list))
fogdur_a_cldN_low_wmax = np.zeros(np.shape(cldN_list))
fogprop_cldN_low_wmax = np.zeros((289, 10))
fog_cldN_low_wmax = np.zeros(np.shape(cldN_list))
fogdur_cldN_low_wmax = np.zeros(np.shape(cldN_list))

for j in range(np.shape(fogprop_cldN)[1]):
    cldN_cond = df_all['cldN'] == cldN_list[j]
    df_cldN_fog = df_all[cldN_cond & fog & low_wmax & cond_outlier]
    df_cond = df_all[cldN_cond & low_wmax & cond_outlier]
    df_cldN_a_fog = df_all[cldN_cond & fog & low_wmax]
    fog_a_cldN_low_wmax[j] = np.shape(df_cldN_a_fog)[0] / 50.
    fog_cldN_low_wmax[j] = np.shape(df_cldN_fog)[0] / np.shape(df_cond)[0]
    fogdur_a_cldN_low_wmax[j] = np.mean(df_cldN_a_fog['diss_hr']) - np.mean(df_cldN_a_fog['onset_hr'])
    fogdur_cldN_low_wmax[j] = np.mean(df_cldN_fog['diss_hr']) - np.mean(df_cldN_fog['onset_hr'])
    for i in range(np.shape(fogprop_cldN)[0]):
        on = df_all['onset_hr'] <= i / 12
        diss = df_all['diss_hr'] <= i / 12
        fognum = df_all[on & ~diss & cldN_cond & low_wmax & cond_outlier].shape[0]
        totnum = df_all[cldN_cond & low_wmax & cond_outlier].shape[0]
        fognum_a = df_all[on & ~diss & cldN_cond & low_wmax].shape[0]
        totnum_a = df_all[cldN_cond & low_wmax].shape[0]
        if totnum_a == 0:
            fogprop_a_cldN_low_wmax[i,j] = float("NaN")
        else:
            fogprop_a_cldN_low_wmax[i,j] = fognum_a/totnum_a
        if totnum == 0:
            fogprop_cldN_low_wmax[i,j] = float("NaN")
        else:
            fogprop_cldN_low_wmax[i,j] = fognum/totnum


fogprop_a_cldN_high_wmax = np.zeros((289, 10))
fog_a_cldN_high_wmax = np.zeros(np.shape(cldN_list))
fogdur_a_cldN_high_wmax = np.zeros(np.shape(cldN_list))
fogprop_cldN_high_wmax = np.zeros((289, 10))
fog_cldN_high_wmax = np.zeros(np.shape(cldN_list))
fogdur_cldN_high_wmax = np.zeros(np.shape(cldN_list))

for j in range(np.shape(fogprop_cldN)[1]):
    cldN_cond = df_all['cldN'] == cldN_list[j]
    df_cldN_fog = df_all[cldN_cond & fog & high_wmax & cond_outlier]
    df_cond = df_all[cldN_cond & high_wmax & cond_outlier]
    df_cldN_a_fog = df_all[cldN_cond & fog & high_wmax]
    fog_a_cldN_high_wmax[j] = np.shape(df_cldN_a_fog)[0] / 50.
    fog_cldN_high_wmax[j] = np.shape(df_cldN_fog)[0] / np.shape(df_cond)[0]
    fogdur_a_cldN_high_wmax[j] = np.mean(df_cldN_a_fog['diss_hr']) - np.mean(df_cldN_a_fog['onset_hr'])
    fogdur_cldN_high_wmax[j] = np.mean(df_cldN_fog['diss_hr']) - np.mean(df_cldN_fog['onset_hr'])
    for i in range(np.shape(fogprop_cldN)[0]):
        on = df_all['onset_hr'] <= i / 12
        diss = df_all['diss_hr'] <= i / 12
        fognum = df_all[on & ~diss & cldN_cond & high_wmax & cond_outlier].shape[0]
        totnum = df_all[cldN_cond & high_wmax & cond_outlier].shape[0]
        fognum_a = df_all[on & ~diss & cldN_cond & high_wmax].shape[0]
        totnum_a = df_all[cldN_cond & high_wmax].shape[0]
        if totnum_a == 0:
            fogprop_a_cldN_high_wmax[i,j] = float("NaN")
        else:
            fogprop_a_cldN_high_wmax[i,j] = fognum_a/totnum_a
        if totnum == 0:
            fogprop_cldN_high_wmax[i,j] = float("NaN")
        else:
            fogprop_cldN_high_wmax[i,j] = fognum/totnum
        

fogprop_a_cldN_low_dtsurf = np.zeros((289, 10))
fog_a_cldN_low_dtsurf = np.zeros(np.shape(cldN_list))
fogdur_a_cldN_low_dtsurf = np.zeros(np.shape(cldN_list))
fogprop_cldN_low_dtsurf = np.zeros((289, 10))
fog_cldN_low_dtsurf = np.zeros(np.shape(cldN_list))
fogdur_cldN_low_dtsurf = np.zeros(np.shape(cldN_list))

for j in range(np.shape(fogprop_cldN)[1]):
    cldN_cond = df_all['cldN'] == cldN_list[j]
    df_cldN_fog = df_all[cldN_cond & fog & low_dtsurf & cond_outlier]
    df_cond = df_all[cldN_cond & low_dtsurf & cond_outlier]
    df_cldN_a_fog = df_all[cldN_cond & fog & low_dtsurf]
    fog_a_cldN_low_dtsurf[j] = np.shape(df_cldN_a_fog)[0] / 50.
    fog_cldN_low_dtsurf[j] = np.shape(df_cldN_fog)[0] / np.shape(df_cond)[0]
    fogdur_a_cldN_low_dtsurf[j] = np.mean(df_cldN_a_fog['diss_hr']) - np.mean(df_cldN_a_fog['onset_hr'])
    fogdur_cldN_low_dtsurf[j] = np.mean(df_cldN_fog['diss_hr']) - np.mean(df_cldN_fog['onset_hr'])
    for i in range(np.shape(fogprop_cldN)[0]):
        on = df_all['onset_hr'] <= i / 12
        diss = df_all['diss_hr'] <= i / 12
        fognum = df_all[on & ~diss & cldN_cond & low_dtsurf & cond_outlier].shape[0]
        totnum = df_all[cldN_cond & low_dtsurf & cond_outlier].shape[0]
        fognum_a = df_all[on & ~diss & cldN_cond & low_dtsurf].shape[0]
        totnum_a = df_all[cldN_cond & low_dtsurf].shape[0]
        if totnum_a == 0:
            fogprop_a_cldN_low_dtsurf[i,j] = float("NaN")
        else:
            fogprop_a_cldN_low_dtsurf[i,j] = fognum_a/totnum_a
        if totnum == 0:
            fogprop_cldN_low_dtsurf[i,j] = float("NaN")
        else:
            fogprop_cldN_low_dtsurf[i,j] = fognum/totnum


fogprop_a_cldN_high_dtsurf = np.zeros((289, 10))
fog_a_cldN_high_dtsurf = np.zeros(np.shape(cldN_list))
fogdur_a_cldN_high_dtsurf = np.zeros(np.shape(cldN_list))
fogprop_cldN_high_dtsurf = np.zeros((289, 10))
fog_cldN_high_dtsurf = np.zeros(np.shape(cldN_list))
fogdur_cldN_high_dtsurf = np.zeros(np.shape(cldN_list))

for j in range(np.shape(fogprop_cldN)[1]):
    cldN_cond = df_all['cldN'] == cldN_list[j]
    df_cldN_fog = df_all[cldN_cond & fog & high_dtsurf & cond_outlier]
    df_cond = df_all[cldN_cond & high_dtsurf & cond_outlier]
    df_cldN_a_fog = df_all[cldN_cond & fog & high_dtsurf]
    fog_a_cldN_high_dtsurf[j] = np.shape(df_cldN_a_fog)[0] / 50.
    fog_cldN_high_dtsurf[j] = np.shape(df_cldN_fog)[0] / np.shape(df_cond)[0]
    fogdur_a_cldN_high_dtsurf[j] = np.mean(df_cldN_a_fog['diss_hr']) - np.mean(df_cldN_a_fog['onset_hr'])
    fogdur_cldN_high_dtsurf[j] = np.mean(df_cldN_fog['diss_hr']) - np.mean(df_cldN_fog['onset_hr'])
    for i in range(np.shape(fogprop_cldN)[0]):
        on = df_all['onset_hr'] <= i / 12
        diss = df_all['diss_hr'] <= i / 12
        fognum = df_all[on & ~diss & cldN_cond & high_dtsurf & cond_outlier].shape[0]
        totnum = df_all[cldN_cond & high_dtsurf & cond_outlier].shape[0]
        fognum_a = df_all[on & ~diss & cldN_cond & high_dtsurf].shape[0]
        totnum_a = df_all[cldN_cond & high_dtsurf].shape[0]
        if totnum_a == 0:
            fogprop_a_cldN_high_dtsurf[i,j] = float("NaN")
        else:
            fogprop_a_cldN_high_dtsurf[i,j] = fognum_a/totnum_a
        if totnum == 0:
            fogprop_cldN_high_dtsurf[i,j] = float("NaN")
        else:
            fogprop_cldN_high_dtsurf[i,j] = fognum/totnum

        
fogamount_a_cldN = np.nanmean(fogprop_a_cldN, axis = 0)
fogamount_a_cldN_low_wmax = np.nanmean(fogprop_a_cldN_low_wmax, axis = 0)
fogamount_a_cldN_high_wmax = np.nanmean(fogprop_a_cldN_high_wmax, axis = 0)
fogamount_a_cldN_low_Ug = np.nanmean(fogprop_a_cldN_low_Ug, axis = 0)
fogamount_a_cldN_high_Ug = np.nanmean(fogprop_a_cldN_high_Ug, axis = 0) 
fogamount_a_cldN_low_dtsurf = np.nanmean(fogprop_a_cldN_low_dtsurf, axis = 0)
fogamount_a_cldN_high_dtsurf = np.nanmean(fogprop_a_cldN_high_dtsurf, axis = 0) 

fogtime_a_cldN = np.nansum(fogprop_a_cldN * np.arange(np.shape(fogprop_a_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_a_cldN, axis = 0) /12
fogtime_a_cldN_low_wmax = np.nansum(fogprop_a_cldN_low_wmax * np.arange(np.shape(fogprop_a_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_a_cldN_low_wmax, axis = 0) / 12
fogtime_a_cldN_high_wmax = np.nansum(fogprop_a_cldN_high_wmax * np.arange(np.shape(fogprop_a_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_a_cldN_high_wmax, axis = 0) / 12
fogtime_a_cldN_low_Ug = np.nansum(fogprop_a_cldN_low_Ug * np.arange(np.shape(fogprop_a_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_a_cldN_low_Ug, axis = 0) / 12
fogtime_a_cldN_high_Ug = np.nansum(fogprop_a_cldN_high_Ug * np.arange(np.shape(fogprop_a_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_a_cldN_high_Ug, axis = 0) / 12
fogtime_a_cldN_low_dtsurf = np.nansum(fogprop_a_cldN_low_dtsurf * np.arange(np.shape(fogprop_a_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_a_cldN_low_dtsurf, axis = 0) / 12
fogtime_a_cldN_high_dtsurf = np.nansum(fogprop_a_cldN_high_dtsurf * np.arange(np.shape(fogprop_a_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_a_cldN_high_dtsurf, axis = 0) / 12

fogamount_cldN = np.nanmean(fogprop_cldN, axis = 0)
fogamount_cldN_low_wmax = np.nanmean(fogprop_cldN_low_wmax, axis = 0)
fogamount_cldN_high_wmax = np.nanmean(fogprop_cldN_high_wmax, axis = 0)
fogamount_cldN_low_Ug = np.nanmean(fogprop_cldN_low_Ug, axis = 0)
fogamount_cldN_high_Ug = np.nanmean(fogprop_cldN_high_Ug, axis = 0) 
fogamount_cldN_low_dtsurf = np.nanmean(fogprop_cldN_low_dtsurf, axis = 0)
fogamount_cldN_high_dtsurf = np.nanmean(fogprop_cldN_high_dtsurf, axis = 0) 

fogtime_cldN = np.nansum(fogprop_cldN * np.arange(np.shape(fogprop_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_cldN, axis = 0) /12
fogtime_cldN_low_wmax = np.nansum(fogprop_cldN_low_wmax * np.arange(np.shape(fogprop_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_cldN_low_wmax, axis = 0) / 12
fogtime_cldN_high_wmax = np.nansum(fogprop_cldN_high_wmax * np.arange(np.shape(fogprop_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_cldN_high_wmax, axis = 0) / 12
fogtime_cldN_low_Ug = np.nansum(fogprop_cldN_low_Ug * np.arange(np.shape(fogprop_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_cldN_low_Ug, axis = 0) / 12
fogtime_cldN_high_Ug = np.nansum(fogprop_cldN_high_Ug * np.arange(np.shape(fogprop_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_cldN_high_Ug, axis = 0) / 12
fogtime_cldN_low_dtsurf = np.nansum(fogprop_cldN_low_dtsurf * np.arange(np.shape(fogprop_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_cldN_low_dtsurf, axis = 0) / 12
fogtime_cldN_high_dtsurf = np.nansum(fogprop_cldN_high_dtsurf * np.arange(np.shape(fogprop_cldN)[0])[:,np.newaxis], axis = 0) / np.nansum(fogprop_cldN_high_dtsurf, axis = 0) / 12


# In[4]:


import warnings; warnings.simplefilter('ignore')

df_ExpDur = df_all[['Ug', 'wmax', 'dtsurf', 'cldN', 'duration', 't_dur', 'nl_dur']]

df_ED_vars = df_ExpDur[cond_outlier]
df_ED_h_Ug = df_ExpDur[cond_outlier & high_Ug]
df_ED_l_Ug = df_ExpDur[cond_outlier & low_Ug]
df_ED_h_wmax = df_ExpDur[cond_outlier & high_wmax]
df_ED_l_wmax = df_ExpDur[cond_outlier & low_wmax]
df_ED_h_dtsurf = df_ExpDur[cond_outlier & high_dtsurf]
df_ED_l_dtsurf = df_ExpDur[cond_outlier & low_dtsurf]

df_ED_bg = df_ED_vars.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_h_Ug = df_ED_h_Ug.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_l_Ug = df_ED_l_Ug.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_h_wmax = df_ED_h_wmax.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_l_wmax = df_ED_l_wmax.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_h_dtsurf = df_ED_h_dtsurf.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_l_dtsurf = df_ED_l_dtsurf.groupby(['Ug', 'wmax', 'dtsurf']).mean()

f_durs = np.empty(np.shape(df_ED_bg)[0])
f_durweights = np.empty(np.shape(df_ED_bg)[0])
f_durs_h_Ug = np.empty(np.shape(df_ED_bg_h_Ug)[0])
f_durweights_h_Ug = np.empty(np.shape(df_ED_bg_h_Ug)[0])
f_durs_l_Ug = np.empty(np.shape(df_ED_bg_l_Ug)[0])
f_durweights_l_Ug = np.empty(np.shape(df_ED_bg_l_Ug)[0])
f_durs_h_wmax = np.empty(np.shape(df_ED_bg_h_wmax)[0])
f_durweights_h_wmax = np.empty(np.shape(df_ED_bg_h_wmax)[0])
f_durs_l_wmax = np.empty(np.shape(df_ED_bg_l_wmax)[0])
f_durweights_l_wmax = np.empty(np.shape(df_ED_bg_l_wmax)[0])
f_durs_h_dtsurf = np.empty(np.shape(df_ED_bg_h_dtsurf)[0])
f_durweights_h_dtsurf = np.empty(np.shape(df_ED_bg_h_dtsurf)[0])
f_durs_l_dtsurf = np.empty(np.shape(df_ED_bg_l_dtsurf)[0])
f_durweights_l_dtsurf = np.empty(np.shape(df_ED_bg_l_dtsurf)[0])

f_durs_cldN = np.empty((np.shape(df_ED_bg)[0], 10))
f_durs_cldN[:] = np.nan
f_durs_cldN_h_Ug = np.empty((np.shape(df_ED_bg_h_Ug)[0], 10))
f_durs_cldN_h_Ug[:] = np.nan
f_durs_cldN_l_Ug = np.empty((np.shape(df_ED_bg_l_Ug)[0], 10))
f_durs_cldN_l_Ug[:] = np.nan
f_durs_cldN_h_wmax = np.empty((np.shape(df_ED_bg_h_wmax)[0], 10))
f_durs_cldN_h_wmax[:] = np.nan
f_durs_cldN_l_wmax = np.empty((np.shape(df_ED_bg_l_wmax)[0], 10))
f_durs_cldN_l_wmax[:] = np.nan
f_durs_cldN_h_dtsurf = np.empty((np.shape(df_ED_bg_h_dtsurf)[0], 10))
f_durs_cldN_h_dtsurf[:] = np.nan
f_durs_cldN_l_dtsurf = np.empty((np.shape(df_ED_bg_l_dtsurf)[0], 10))
f_durs_cldN_l_dtsurf[:] = np.nan

t_durs = np.empty(np.shape(df_ED_bg)[0])
t_durweights = np.empty(np.shape(df_ED_bg)[0])

t_durs_cldN = np.empty((np.shape(df_ED_bg)[0], 10))
t_durs_cldN[:] = np.nan

nl_durs = np.empty(np.shape(df_ED_bg)[0])
nl_durweights = np.empty(np.shape(df_ED_bg)[0])

nl_durs_cldN = np.empty((np.shape(df_ED_bg)[0], 10))
nl_durs_cldN[:] = np.nan

for i in range(np.shape(df_ED_bg)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg.index[i]
    cond_Ug = df_ED_vars['Ug'] == Ug
    cond_wmax = df_ED_vars['wmax'] == wmax
    cond_dtsurf = df_ED_vars['dtsurf'] == dtsurf
    cond_fog = df_ED_vars['duration'] > 0.0
    cond_t = df_ED_vars['t_dur'] > 0.0
    cond_nl = df_ED_vars['t_dur'] > 0.0

    df_this_bg = df_ED_vars[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_vars[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    df_this_bg_t_fog = df_ED_vars[cond_Ug & cond_wmax & cond_dtsurf & cond_t].reset_index()
    df_this_bg_nl_fog = df_ED_vars[cond_Ug & cond_wmax & cond_dtsurf & cond_nl].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    this_t_fog_cldN = df_this_bg_t_fog['cldN'].to_numpy()
    this_nl_fog_cldN = df_this_bg_nl_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights[i] = np.shape(f_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN[i,j] = f_df_cldN['duration'][index] - f_durs[i]
            else:
                f_durs_cldN[i,j] = np.nan
                
    if np.shape(this_t_fog_cldN)[0] > 0:
        t_maxcldN = np.nanmax(this_t_fog_cldN)
        t_mincldN = np.nanmin(this_t_fog_cldN)
        t_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == t_maxcldN)[0] + 1])
        t_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == t_mincldN)[0] - 1])
        t_maxNcond = df_this_bg['cldN'] <= cldN_list[t_maxNi]
        t_minNcond = df_this_bg['cldN'] >= cldN_list[t_minNi]
        t_df_cldN = df_this_bg[t_maxNcond & t_minNcond].reset_index()
        
        t_durs[i] = np.nanmean(t_df_cldN['t_dur'])
        t_durweights[i] = np.shape(t_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in t_df_cldN['cldN'].to_numpy():
                index = np.where(t_df_cldN['cldN'] == cldN_list[j])[0]
                t_durs_cldN[i,j] = t_df_cldN['t_dur'][index] - t_durs[i]
            else:
                t_durs_cldN[i,j] = np.nan
                
    if np.shape(this_nl_fog_cldN)[0] > 0:
        nl_maxcldN = np.nanmax(this_nl_fog_cldN)
        nl_mincldN = np.nanmin(this_nl_fog_cldN)
        nl_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == nl_maxcldN)[0] + 1])
        nl_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == nl_mincldN)[0] - 1])
        nl_maxNcond = df_this_bg['cldN'] <= cldN_list[nl_maxNi]
        nl_minNcond = df_this_bg['cldN'] >= cldN_list[nl_minNi]
        nl_df_cldN = df_this_bg[nl_maxNcond & nl_minNcond].reset_index()
        
        nl_durs[i] = np.nanmean(nl_df_cldN['nl_dur'])
        nl_durweights[i] = np.shape(nl_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in nl_df_cldN['cldN'].to_numpy():
                index = np.where(nl_df_cldN['cldN'] == cldN_list[j])[0]
                nl_durs_cldN[i,j] = nl_df_cldN['nl_dur'][index] - nl_durs[i]
            else:
                nl_durs_cldN[i,j] = np.nan

for i in range(np.shape(df_ED_bg_h_Ug)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_h_Ug.index[i]
    cond_Ug = df_ED_h_Ug['Ug'] == Ug
    cond_wmax = df_ED_h_Ug['wmax'] == wmax
    cond_dtsurf = df_ED_h_Ug['dtsurf'] == dtsurf
    cond_fog = df_ED_h_Ug['duration'] > 0.0

    df_this_bg = df_ED_h_Ug[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_h_Ug[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_h_Ug[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_h_Ug[i] = np.shape(f_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_h_Ug[i,j] = f_df_cldN['duration'][index] - f_durs_h_Ug[i]
            else:
                f_durs_cldN_h_Ug[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_l_Ug)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_l_Ug.index[i]
    cond_Ug = df_ED_l_Ug['Ug'] == Ug
    cond_wmax = df_ED_l_Ug['wmax'] == wmax
    cond_dtsurf = df_ED_l_Ug['dtsurf'] == dtsurf
    cond_fog = df_ED_l_Ug['duration'] > 0.0

    df_this_bg = df_ED_l_Ug[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_l_Ug[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_l_Ug[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_l_Ug[i] = np.shape(f_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_l_Ug[i,j] = f_df_cldN['duration'][index] - f_durs_l_Ug[i]
            else:
                f_durs_cldN_l_Ug[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_h_wmax)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_h_wmax.index[i]
    cond_Ug = df_ED_h_wmax['Ug'] == Ug
    cond_wmax = df_ED_h_wmax['wmax'] == wmax
    cond_dtsurf = df_ED_h_wmax['dtsurf'] == dtsurf
    cond_fog = df_ED_h_wmax['duration'] > 0.0

    df_this_bg = df_ED_h_wmax[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_h_wmax[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_h_wmax[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_h_wmax[i] = np.shape(f_df_cldN)[0]
        # print(f_durs_h_wmax[i], f_durweights_h_wmax[i])
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_h_wmax[i,j] = f_df_cldN['duration'][index] - f_durs_h_wmax[i]
            else:
                f_durs_cldN_h_wmax[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_l_wmax)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_l_wmax.index[i]
    cond_Ug = df_ED_l_wmax['Ug'] == Ug
    cond_wmax = df_ED_l_wmax['wmax'] == wmax
    cond_dtsurf = df_ED_l_wmax['dtsurf'] == dtsurf
    cond_fog = df_ED_l_wmax['duration'] > 0.0

    df_this_bg = df_ED_l_wmax[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_l_wmax[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_l_wmax[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_l_wmax[i] = np.shape(f_df_cldN)[0]
        print(f_durs_l_wmax[i], f_durweights_l_wmax[i])
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_l_wmax[i,j] = f_df_cldN['duration'][index] - f_durs_l_wmax[i]
            else:
                f_durs_cldN_l_wmax[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_h_dtsurf)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_h_dtsurf.index[i]
    cond_Ug = df_ED_h_dtsurf['Ug'] == Ug
    cond_wmax = df_ED_h_dtsurf['wmax'] == wmax
    cond_dtsurf = df_ED_h_dtsurf['dtsurf'] == dtsurf
    cond_fog = df_ED_h_dtsurf['duration'] > 0.0

    df_this_bg = df_ED_h_dtsurf[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_h_dtsurf[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_h_dtsurf[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_h_dtsurf[i] = np.shape(f_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_h_dtsurf[i,j] = f_df_cldN['duration'][index] - f_durs_h_dtsurf[i]
            else:
                f_durs_cldN_h_dtsurf[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_l_dtsurf)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_l_dtsurf.index[i]
    cond_Ug = df_ED_l_dtsurf['Ug'] == Ug
    cond_wmax = df_ED_l_dtsurf['wmax'] == wmax
    cond_dtsurf = df_ED_l_dtsurf['dtsurf'] == dtsurf
    cond_fog = df_ED_l_dtsurf['duration'] > 0.0

    df_this_bg = df_ED_l_dtsurf[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_l_dtsurf[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_l_dtsurf[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_l_dtsurf[i] = np.shape(f_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_l_dtsurf[i,j] = f_df_cldN['duration'][index] - f_durs_l_dtsurf[i]
            else:
                f_durs_cldN_l_dtsurf[i,j] = np.nan


                    
f_durchange = np.nanmean(f_durs_cldN, axis = 0)
t_durchange = np.nanmean(t_durs_cldN, axis = 0)
nl_durchange = np.nanmean(nl_durs_cldN, axis = 0)
f_durstd = np.nanstd(f_durs_cldN, axis = 0)
t_durstd = np.nanstd(t_durs_cldN, axis = 0)
nl_durstd = np.nanstd(nl_durs_cldN, axis = 0)
f_dcount = np.count_nonzero(~np.isnan(f_durs_cldN), axis=0)
t_dcount = np.count_nonzero(~np.isnan(t_durs_cldN), axis=0)
nl_dcount = np.count_nonzero(~np.isnan(nl_durs_cldN), axis=0)

f_durchange_h_Ug = np.nanmean(f_durs_cldN_h_Ug, axis = 0)
f_durchange_l_Ug = np.nanmean(f_durs_cldN_l_Ug, axis = 0)
f_durchange_h_wmax = np.nanmean(f_durs_cldN_h_wmax, axis = 0)
f_durchange_l_wmax = np.nanmean(f_durs_cldN_l_wmax, axis = 0)
f_durchange_h_dtsurf = np.nanmean(f_durs_cldN_h_dtsurf, axis = 0)
f_durchange_l_dtsurf = np.nanmean(f_durs_cldN_l_dtsurf, axis = 0)

f_durstd_h_Ug = np.nanstd(f_durs_cldN_h_Ug, axis = 0)
f_durstd_l_Ug = np.nanstd(f_durs_cldN_l_Ug, axis = 0)
f_durstd_h_wmax = np.nanstd(f_durs_cldN_h_wmax, axis = 0)
f_durstd_l_wmax = np.nanstd(f_durs_cldN_l_wmax, axis = 0)
f_durstd_h_dtsurf = np.nanstd(f_durs_cldN_h_dtsurf, axis = 0)
f_durstd_l_dtsurf = np.nanstd(f_durs_cldN_l_dtsurf, axis = 0)

f_dcount_h_Ug = np.count_nonzero(~np.isnan(f_durs_cldN_h_Ug), axis=0)
f_dcount_l_Ug = np.count_nonzero(~np.isnan(f_durs_cldN_l_Ug), axis=0)
f_dcount_h_wmax = np.count_nonzero(~np.isnan(f_durs_cldN_h_wmax), axis=0)
f_dcount_l_wmax = np.count_nonzero(~np.isnan(f_durs_cldN_l_wmax), axis=0)
f_dcount_h_dtsurf = np.count_nonzero(~np.isnan(f_durs_cldN_h_dtsurf), axis=0)
f_dcount_l_dtsurf = np.count_nonzero(~np.isnan(f_durs_cldN_l_dtsurf), axis=0)

wmean_f_dur = np.nansum(f_durs * f_durweights) / np.nansum(f_durweights)
wmean_t_dur = np.nansum(t_durs * t_durweights) / np.nansum(t_durweights)
wmean_nl_dur = np.nansum(nl_durs * nl_durweights) / np.nansum(nl_durweights)

wmean_f_dur_h_Ug = np.nansum(f_durs_h_Ug * f_durweights_h_Ug) / np.nansum(f_durweights_h_Ug)
wmean_f_dur_l_Ug = np.nansum(f_durs_l_Ug * f_durweights_l_Ug) / np.nansum(f_durweights_l_Ug)
wmean_f_dur_h_wmax = np.nansum(f_durs_h_wmax * f_durweights_h_wmax) / np.nansum(f_durweights_h_wmax)
wmean_f_dur_l_wmax = np.nansum(f_durs_l_wmax * f_durweights_l_wmax) / np.nansum(f_durweights_l_wmax)
wmean_f_dur_h_dtsurf = np.nansum(f_durs_h_dtsurf * f_durweights_h_dtsurf) / np.nansum(f_durweights_h_dtsurf)
wmean_f_dur_l_dtsurf = np.nansum(f_durs_l_dtsurf * f_durweights_l_dtsurf) / np.nansum(f_durweights_l_dtsurf)


# In[5]:


df_ED_vars_a = df_ExpDur.copy()
df_ED_h_Ug_a = df_ExpDur[high_Ug]
df_ED_l_Ug_a = df_ExpDur[low_Ug]
df_ED_h_wmax_a = df_ExpDur[high_wmax]
df_ED_l_wmax_a = df_ExpDur[low_wmax]
df_ED_h_dtsurf_a = df_ExpDur[high_dtsurf]
df_ED_l_dtsurf_a = df_ExpDur[low_dtsurf]

df_ED_bg_a = df_ED_vars_a.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_h_Ug_a = df_ED_h_Ug_a.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_l_Ug_a = df_ED_l_Ug_a.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_h_wmax_a = df_ED_h_wmax_a.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_l_wmax_a = df_ED_l_wmax_a.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_h_dtsurf_a = df_ED_h_dtsurf_a.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_l_dtsurf_a = df_ED_l_dtsurf_a.groupby(['Ug', 'wmax', 'dtsurf']).mean()

f_durs_a = np.empty(np.shape(df_ED_bg_a)[0])
f_durweights_a = np.empty(np.shape(df_ED_bg_a)[0])
f_durs_h_Ug_a = np.empty(np.shape(df_ED_bg_h_Ug_a)[0])
f_durweights_h_Ug_a = np.empty(np.shape(df_ED_bg_h_Ug_a)[0])
f_durs_l_Ug_a = np.empty(np.shape(df_ED_bg_l_Ug_a)[0])
f_durweights_l_Ug_a = np.empty(np.shape(df_ED_bg_l_Ug_a)[0])
f_durs_h_wmax_a = np.empty(np.shape(df_ED_bg_h_wmax_a)[0])
f_durweights_h_wmax_a = np.empty(np.shape(df_ED_bg_h_wmax_a)[0])
f_durs_l_wmax_a = np.empty(np.shape(df_ED_bg_l_wmax_a)[0])
f_durweights_l_wmax_a = np.empty(np.shape(df_ED_bg_l_wmax_a)[0])
f_durs_h_dtsurf_a = np.empty(np.shape(df_ED_bg_h_dtsurf_a)[0])
f_durweights_h_dtsurf_a = np.empty(np.shape(df_ED_bg_h_dtsurf_a)[0])
f_durs_l_dtsurf_a = np.empty(np.shape(df_ED_bg_l_dtsurf_a)[0])
f_durweights_l_dtsurf_a = np.empty(np.shape(df_ED_bg_l_dtsurf_a)[0])

f_durs_cldN_a = np.empty((np.shape(df_ED_bg_a)[0], 10))
f_durs_cldN_a[:] = np.nan
f_durs_cldN_h_Ug_a = np.empty((np.shape(df_ED_bg_h_Ug_a)[0], 10))
f_durs_cldN_h_Ug_a[:] = np.nan
f_durs_cldN_l_Ug_a = np.empty((np.shape(df_ED_bg_l_Ug_a)[0], 10))
f_durs_cldN_l_Ug_a[:] = np.nan
f_durs_cldN_h_wmax_a = np.empty((np.shape(df_ED_bg_h_wmax_a)[0], 10))
f_durs_cldN_h_wmax_a[:] = np.nan
f_durs_cldN_l_wmax_a = np.empty((np.shape(df_ED_bg_l_wmax_a)[0], 10))
f_durs_cldN_l_wmax_a[:] = np.nan
f_durs_cldN_h_dtsurf_a = np.empty((np.shape(df_ED_bg_h_dtsurf_a)[0], 10))
f_durs_cldN_h_dtsurf_a[:] = np.nan
f_durs_cldN_l_dtsurf_a = np.empty((np.shape(df_ED_bg_l_dtsurf_a)[0], 10))
f_durs_cldN_l_dtsurf_a[:] = np.nan

t_durs_a = np.empty(np.shape(df_ED_bg_a)[0])
t_durweights_a = np.empty(np.shape(df_ED_bg_a)[0])

t_durs_cldN_a = np.empty((np.shape(df_ED_bg_a)[0], 10))
t_durs_cldN_a[:] = np.nan

nl_durs_a = np.empty(np.shape(df_ED_bg_a)[0])
nl_durweights_a = np.empty(np.shape(df_ED_bg_a)[0])

nl_durs_cldN_a = np.empty((np.shape(df_ED_bg_a)[0], 10))
nl_durs_cldN_a[:] = np.nan

for i in range(np.shape(df_ED_bg_a)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_a.index[i]
    cond_Ug = df_ED_vars_a['Ug'] == Ug
    cond_wmax = df_ED_vars_a['wmax'] == wmax
    cond_dtsurf = df_ED_vars_a['dtsurf'] == dtsurf
    cond_fog = df_ED_vars_a['duration'] > 0.0
    cond_t = df_ED_vars_a['t_dur'] > 0.0
    cond_nl = df_ED_vars_a['t_dur'] > 0.0

    df_this_bg = df_ED_vars_a[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_vars_a[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    df_this_bg_t_fog = df_ED_vars_a[cond_Ug & cond_wmax & cond_dtsurf & cond_t].reset_index()
    df_this_bg_nl_fog = df_ED_vars_a[cond_Ug & cond_wmax & cond_dtsurf & cond_nl].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    this_t_fog_cldN = df_this_bg_t_fog['cldN'].to_numpy()
    this_nl_fog_cldN = df_this_bg_nl_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_a[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_a[i] = np.shape(f_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_a[i,j] = f_df_cldN['duration'][index] - f_durs_a[i]
            else:
                f_durs_cldN_a[i,j] = np.nan
                
    if np.shape(this_t_fog_cldN)[0] > 0:
        t_maxcldN = np.nanmax(this_t_fog_cldN)
        t_mincldN = np.nanmin(this_t_fog_cldN)
        t_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == t_maxcldN)[0] + 1])
        t_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == t_mincldN)[0] - 1])
        t_maxNcond = df_this_bg['cldN'] <= cldN_list[t_maxNi]
        t_minNcond = df_this_bg['cldN'] >= cldN_list[t_minNi]
        t_df_cldN = df_this_bg[t_maxNcond & t_minNcond].reset_index()
        
        t_durs_a[i] = np.nanmean(t_df_cldN['t_dur'])
        t_durweights_a[i] = np.shape(t_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in t_df_cldN['cldN'].to_numpy():
                index = np.where(t_df_cldN['cldN'] == cldN_list[j])[0]
                t_durs_cldN_a[i,j] = t_df_cldN['t_dur'][index] - t_durs_a[i]
            else:
                t_durs_cldN_a[i,j] = np.nan
                
    if np.shape(this_nl_fog_cldN)[0] > 0:
        nl_maxcldN = np.nanmax(this_nl_fog_cldN)
        nl_mincldN = np.nanmin(this_nl_fog_cldN)
        nl_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == nl_maxcldN)[0] + 1])
        nl_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == nl_mincldN)[0] - 1])
        nl_maxNcond = df_this_bg['cldN'] <= cldN_list[nl_maxNi]
        nl_minNcond = df_this_bg['cldN'] >= cldN_list[nl_minNi]
        nl_df_cldN = df_this_bg[nl_maxNcond & nl_minNcond].reset_index()
        
        nl_durs_a[i] = np.nanmean(nl_df_cldN['nl_dur'])
        nl_durweights_a[i] = np.shape(nl_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in nl_df_cldN['cldN'].to_numpy():
                index = np.where(nl_df_cldN['cldN'] == cldN_list[j])[0]
                nl_durs_cldN_a[i,j] = nl_df_cldN['nl_dur'][index] - nl_durs_a[i]
            else:
                nl_durs_cldN_a[i,j] = np.nan

for i in range(np.shape(df_ED_bg_h_Ug_a)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_h_Ug_a.index[i]
    cond_Ug = df_ED_h_Ug_a['Ug'] == Ug
    cond_wmax = df_ED_h_Ug_a['wmax'] == wmax
    cond_dtsurf = df_ED_h_Ug_a['dtsurf'] == dtsurf
    cond_fog = df_ED_h_Ug_a['duration'] > 0.0

    df_this_bg = df_ED_h_Ug_a[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_h_Ug_a[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_h_Ug_a[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_h_Ug_a[i] = np.shape(f_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_h_Ug_a[i,j] = f_df_cldN['duration'][index] - f_durs_h_Ug_a[i]
            else:
                f_durs_cldN_h_Ug_a[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_l_Ug_a)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_l_Ug_a.index[i]
    cond_Ug = df_ED_l_Ug_a['Ug'] == Ug
    cond_wmax = df_ED_l_Ug_a['wmax'] == wmax
    cond_dtsurf = df_ED_l_Ug_a['dtsurf'] == dtsurf
    cond_fog = df_ED_l_Ug_a['duration'] > 0.0

    df_this_bg = df_ED_l_Ug_a[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_l_Ug_a[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_l_Ug_a[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_l_Ug_a[i] = np.shape(f_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_l_Ug_a[i,j] = f_df_cldN['duration'][index] - f_durs_l_Ug_a[i]
            else:
                f_durs_cldN_l_Ug_a[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_h_wmax_a)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_h_wmax_a.index[i]
    cond_Ug = df_ED_h_wmax_a['Ug'] == Ug
    cond_wmax = df_ED_h_wmax_a['wmax'] == wmax
    cond_dtsurf = df_ED_h_wmax_a['dtsurf'] == dtsurf
    cond_fog = df_ED_h_wmax_a['duration'] > 0.0

    df_this_bg = df_ED_h_wmax_a[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_h_wmax_a[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_h_wmax_a[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_h_wmax_a[i] = np.shape(f_df_cldN)[0]
        # print(f_durs_h_wmax[i], f_durweights_h_wmax[i])
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_h_wmax_a[i,j] = f_df_cldN['duration'][index] - f_durs_h_wmax_a[i]
            else:
                f_durs_cldN_h_wmax_a[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_l_wmax_a)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_l_wmax_a.index[i]
    cond_Ug = df_ED_l_wmax_a['Ug'] == Ug
    cond_wmax = df_ED_l_wmax_a['wmax'] == wmax
    cond_dtsurf = df_ED_l_wmax_a['dtsurf'] == dtsurf
    cond_fog = df_ED_l_wmax_a['duration'] > 0.0

    df_this_bg = df_ED_l_wmax_a[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_l_wmax_a[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_l_wmax_a[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_l_wmax_a[i] = np.shape(f_df_cldN)[0]
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_l_wmax_a[i,j] = f_df_cldN['duration'][index] - f_durs_l_wmax_a[i]
            else:
                f_durs_cldN_l_wmax_a[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_h_dtsurf_a)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_h_dtsurf_a.index[i]
    cond_Ug = df_ED_h_dtsurf_a['Ug'] == Ug
    cond_wmax = df_ED_h_dtsurf_a['wmax'] == wmax
    cond_dtsurf = df_ED_h_dtsurf_a['dtsurf'] == dtsurf
    cond_fog = df_ED_h_dtsurf_a['duration'] > 0.0

    df_this_bg = df_ED_h_dtsurf_a[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_h_dtsurf_a[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_h_dtsurf_a[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_h_dtsurf_a[i] = np.shape(f_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_h_dtsurf_a[i,j] = f_df_cldN['duration'][index] - f_durs_h_dtsurf_a[i]
            else:
                f_durs_cldN_h_dtsurf_a[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_l_dtsurf_a)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_l_dtsurf_a.index[i]
    cond_Ug = df_ED_l_dtsurf_a['Ug'] == Ug
    cond_wmax = df_ED_l_dtsurf_a['wmax'] == wmax
    cond_dtsurf = df_ED_l_dtsurf_a['dtsurf'] == dtsurf
    cond_fog = df_ED_l_dtsurf_a['duration'] > 0.0

    df_this_bg = df_ED_l_dtsurf_a[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_l_dtsurf_a[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs_l_dtsurf_a[i] = np.nanmean(f_df_cldN['duration'])
        f_durweights_l_dtsurf_a[i] = np.shape(f_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN_l_dtsurf_a[i,j] = f_df_cldN['duration'][index] - f_durs_l_dtsurf_a[i]
            else:
                f_durs_cldN_l_dtsurf_a[i,j] = np.nan


                    
f_durchange_a = np.nanmean(f_durs_cldN_a, axis = 0)
t_durchange_a = np.nanmean(t_durs_cldN_a, axis = 0)
nl_durchange_a = np.nanmean(nl_durs_cldN_a, axis = 0)
f_durstd_a = np.nanstd(f_durs_cldN_a, axis = 0)
t_durstd_a = np.nanstd(t_durs_cldN_a, axis = 0)
nl_durstd_a = np.nanstd(nl_durs_cldN_a, axis = 0)
f_dcount_a = np.count_nonzero(~np.isnan(f_durs_cldN_a), axis=0)
t_dcount_a = np.count_nonzero(~np.isnan(t_durs_cldN_a), axis=0)
nl_dcount_a = np.count_nonzero(~np.isnan(nl_durs_cldN_a), axis=0)

f_durchange_h_Ug_a = np.nanmean(f_durs_cldN_h_Ug_a, axis = 0)
f_durchange_l_Ug_a = np.nanmean(f_durs_cldN_l_Ug_a, axis = 0)
f_durchange_h_wmax_a = np.nanmean(f_durs_cldN_h_wmax_a, axis = 0)
f_durchange_l_wmax_a = np.nanmean(f_durs_cldN_l_wmax_a, axis = 0)
f_durchange_h_dtsurf_a = np.nanmean(f_durs_cldN_h_dtsurf_a, axis = 0)
f_durchange_l_dtsurf_a = np.nanmean(f_durs_cldN_l_dtsurf_a, axis = 0)

f_durstd_h_Ug_a = np.nanstd(f_durs_cldN_h_Ug_a, axis = 0)
f_durstd_l_Ug_a = np.nanstd(f_durs_cldN_l_Ug_a, axis = 0)
f_durstd_h_wmax_a = np.nanstd(f_durs_cldN_h_wmax_a, axis = 0)
f_durstd_l_wmax_a = np.nanstd(f_durs_cldN_l_wmax_a, axis = 0)
f_durstd_h_dtsurf_a = np.nanstd(f_durs_cldN_h_dtsurf_a, axis = 0)
f_durstd_l_dtsurf_a = np.nanstd(f_durs_cldN_l_dtsurf_a, axis = 0)

f_dcount_h_Ug_a = np.count_nonzero(~np.isnan(f_durs_cldN_h_Ug_a), axis=0)
f_dcount_l_Ug_a = np.count_nonzero(~np.isnan(f_durs_cldN_l_Ug_a), axis=0)
f_dcount_h_wmax_a = np.count_nonzero(~np.isnan(f_durs_cldN_h_wmax_a), axis=0)
f_dcount_l_wmax_a = np.count_nonzero(~np.isnan(f_durs_cldN_l_wmax_a), axis=0)
f_dcount_h_dtsurf_a = np.count_nonzero(~np.isnan(f_durs_cldN_h_dtsurf_a), axis=0)
f_dcount_l_dtsurf_a = np.count_nonzero(~np.isnan(f_durs_cldN_l_dtsurf_a), axis=0)

wmean_f_dur_a = np.nansum(f_durs_a * f_durweights_a) / np.nansum(f_durweights_a)
wmean_t_dur_a = np.nansum(t_durs_a * t_durweights_a) / np.nansum(t_durweights_a)
wmean_nl_dur_a = np.nansum(nl_durs_a * nl_durweights_a) / np.nansum(nl_durweights_a)

wmean_f_dur_h_Ug_a = np.nansum(f_durs_h_Ug_a * f_durweights_h_Ug_a) / np.nansum(f_durweights_h_Ug_a)
wmean_f_dur_l_Ug_a = np.nansum(f_durs_l_Ug_a * f_durweights_l_Ug_a) / np.nansum(f_durweights_l_Ug_a)
wmean_f_dur_h_wmax_a = np.nansum(f_durs_h_wmax_a * f_durweights_h_wmax_a) / np.nansum(f_durweights_h_wmax_a)
wmean_f_dur_l_wmax_a = np.nansum(f_durs_l_wmax_a * f_durweights_l_wmax_a) / np.nansum(f_durweights_l_wmax_a)
wmean_f_dur_h_dtsurf_a = np.nansum(f_durs_h_dtsurf_a * f_durweights_h_dtsurf_a) / np.nansum(f_durweights_h_dtsurf_a)
wmean_f_dur_l_dtsurf_a = np.nansum(f_durs_l_dtsurf_a * f_durweights_l_dtsurf_a) / np.nansum(f_durweights_l_dtsurf_a)


# In[6]:


x = np.array([100, 150, 200, 300, 400, 500, 750, 1000, 1250, 1500])

fig, (ax2, ax3) = plt.subplots(1, 2)
fig.set_size_inches(14,7)
plt.subplots_adjust(left=0.1, right=0.95, bottom=0.1, top=0.925, wspace=0.25, hspace=0.2)

# ax1.plot(x, fogamount_a_cldN, color = 'black', label = 'Overall')
# ax1.plot(x, fogamount_a_cldN_high_wmax, color = 'blue')
# ax1.plot(x, fogamount_a_cldN_low_wmax, color = 'blue', linestyle = 'dashed')
# ax1.plot(x, fogamount_a_cldN_high_Ug, color = 'green')
# ax1.plot(x, fogamount_a_cldN_low_Ug, color = 'green', linestyle = 'dashed')
# ax1.plot(x, fogamount_a_cldN_high_dtsurf, color = 'red')
# ax1.plot(x, fogamount_a_cldN_low_dtsurf, color = 'red', linestyle = 'dashed')
# ax1.legend(fontsize = 14)

ax2.plot(x, fog_a_cldN, color = 'black', label = 'Overall')
ax2.plot(x, fog_a_cldN_high_wmax, color = 'blue', label = 'Weak Subsidence')
ax2.plot(x, fog_a_cldN_low_wmax, color = 'blue', linestyle = 'dashed', label = 'Strong Subsidence')
ax2.plot(x, fog_a_cldN_high_Ug, color = 'green', label = 'High Wind')
ax2.plot(x, fog_a_cldN_low_Ug, color = 'green', linestyle = 'dashed', label = 'Low Wind')
ax2.plot(x, fog_a_cldN_high_dtsurf, color = 'red', label = 'Warming Surface')
ax2.plot(x, fog_a_cldN_low_dtsurf, color = 'red', linestyle = 'dashed', label = 'Cooling Surface')
ax2.legend(fontsize = 14)

ax3.plot(x, wmean_f_dur_a + f_durchange_a, color = 'black', label = 'Overall')
ax3.plot(x, wmean_f_dur_h_wmax_a + f_durchange_h_wmax_a, color = 'blue', label = 'Weak Subsidence')
ax3.plot(x, wmean_f_dur_l_wmax_a + f_durchange_l_wmax_a, color = 'blue', linestyle = 'dashed', label = 'Strong Subsidence')
ax3.plot(x, wmean_f_dur_h_Ug_a + f_durchange_h_Ug_a, color = 'green', label = 'High Wind')
ax3.plot(x, wmean_f_dur_l_Ug_a + f_durchange_l_Ug_a, color = 'green', linestyle = 'dashed', label = 'Low Wind')
ax3.plot(x, wmean_f_dur_h_dtsurf_a + f_durchange_h_dtsurf_a, color = 'red', label = 'Warming Surface')
ax3.plot(x, wmean_f_dur_l_dtsurf_a + f_durchange_l_dtsurf_a, color = 'red', linestyle = 'dashed', label = 'Cooling Surface')
# ax3.legend(fontsize = 14)

fig.suptitle('Fog Formation Fraction and Expected Duration vs Na Under Different Conditions, no Outliers', fontsize = 20)
# ax1.set_title('Simulation Fogginess', fontsize = 20)
# ax1.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 18)
# ax1.set_ylabel('Average Fog Proportion', fontsize = 18)
# ax1.tick_params(axis='both', which='both', labelsize=14)


# ax2.set_title('Fog Formation Fraction', fontsize = 20)
ax2.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 20)
ax2.set_ylabel('Proportion of Simulations', fontsize = 20)
ax2.tick_params(axis='both', which='both', labelsize=14)

# ax3.set_title('Expected Fog Duration', fontsize = 20)
ax3.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 20)
ax3.set_ylabel('Fog Duration (hours)', fontsize = 20)
ax3.tick_params(axis='both', which='both', labelsize=14)

fname = desc + 'F_Regimes_a.png'
plt.savefig(fname)
plt.show()
plt.close()

fig, (ax2, ax3) = plt.subplots(1, 2)
fig.set_size_inches(14,7)
plt.subplots_adjust(left=0.1, right=0.95, bottom=0.1, top=0.925, wspace=0.25, hspace=0.2)

# ax1.plot(x, fogamount_cldN, color = 'black', label = 'Overall')
# ax1.plot(x, fogamount_cldN_high_wmax, color = 'blue')
# ax1.plot(x, fogamount_cldN_low_wmax, color = 'blue', linestyle = 'dashed')
# ax1.plot(x, fogamount_cldN_high_Ug, color = 'green')
# ax1.plot(x, fogamount_cldN_low_Ug, color = 'green', linestyle = 'dashed')
# ax1.plot(x, fogamount_cldN_high_dtsurf, color = 'red')
# ax1.plot(x, fogamount_cldN_low_dtsurf, color = 'red', linestyle = 'dashed')
# ax1.legend(fontsize = 14)

ax2.plot(x, fog_cldN, color = 'black', label = 'Overall')
ax2.plot(x, fog_cldN_high_wmax, color = 'blue', label = 'Weak Subsidence')
ax2.plot(x, fog_cldN_low_wmax, color = 'blue', linestyle = 'dashed', label = 'Strong Subsidence')
ax2.plot(x, fog_cldN_high_Ug, color = 'green', label = 'High Wind')
ax2.plot(x, fog_cldN_low_Ug, color = 'green', linestyle = 'dashed', label = 'Low Wind')
ax2.plot(x, fog_cldN_high_dtsurf, color = 'red', label = 'Warming Surface')
ax2.plot(x, fog_cldN_low_dtsurf, color = 'red', linestyle = 'dashed', label = 'Cooling Surface')
ax2.legend(fontsize = 14)

ax3.plot(x, wmean_f_dur + f_durchange, color = 'black', label = 'Overall')
ax3.plot(x, wmean_f_dur_h_wmax + f_durchange_h_wmax, color = 'blue', label = 'Weak Subsidence')
ax3.plot(x, wmean_f_dur_l_wmax + f_durchange_l_wmax, color = 'blue', linestyle = 'dashed', label = 'Strong Subsidence')
ax3.plot(x, wmean_f_dur_h_Ug + f_durchange_h_Ug, color = 'green', label = 'High Wind')
ax3.plot(x, wmean_f_dur_l_Ug + f_durchange_l_Ug, color = 'green', linestyle = 'dashed', label = 'Low Wind')
ax3.plot(x, wmean_f_dur_h_dtsurf + f_durchange_h_dtsurf, color = 'red', label = 'Warming Surface')
ax3.plot(x, wmean_f_dur_l_dtsurf + f_durchange_l_dtsurf, color = 'red', linestyle = 'dashed', label = 'Cooling Surface')
# ax3.legend(fontsize = 14)

fig.suptitle('Fog Formation Fraction and Expected Duration vs Na Under Different Conditions, no Outliers', fontsize = 20)

# ax1.set_title('Simulation Fogginess', fontsize = 20)
# ax1.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 18)
# ax1.set_ylabel('Average Fog Proportion', fontsize = 18)
# ax1.tick_params(axis='both', which='both', labelsize=14)


# ax2.set_title('Fog Formation Fraction', fontsize = 20)
ax2.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 20)
ax2.set_ylabel('Proportion of Simulations', fontsize = 20)
ax2.tick_params(axis='both', which='both', labelsize=14)

# ax3.set_title('Expected Fog Duration', fontsize = 20)
ax3.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 20)
ax3.set_ylabel('Fog Duration (hours)', fontsize = 20)
ax3.tick_params(axis='both', which='both', labelsize=14)

fname = desc + 'F_Regimes.png'
plt.savefig(fname)
plt.show()
plt.close()


# In[7]:


df_fig_fog = df_all[fog]
print(np.shape(df_fig_fog))
df_nout_fog = df_all[fog & cond_outlier]

th_max = 15
th_min = 5

df_gb_cldN_f = df_fig_fog.groupby(['cldN']).mean().reset_index()
df_gb_cldN_fn = df_nout_fog.groupby(['cldN']).mean().reset_index()
df_gb_ug_f = df_fig_fog.groupby(['Ug']).mean().reset_index()
df_gb_wmax_f = df_fig_fog.groupby(['wmax']).mean().reset_index()
df_gb_dtsurf_f = df_fig_fog.groupby(['dtsurf']).mean().reset_index()

gb_cldN_f = (df_gb_cldN_f.index.to_numpy() - min(df_gb_cldN_f.index.to_numpy())) / (max(df_gb_cldN_f.index.to_numpy()) - min(df_gb_cldN_f.index.to_numpy()))
gb_cldN_fn = (df_gb_cldN_fn.index.to_numpy() - min(df_gb_cldN_fn.index.to_numpy())) / (max(df_gb_cldN_fn.index.to_numpy()) - min(df_gb_cldN_fn.index.to_numpy()))
gb_ug_f = (df_gb_ug_f.index.to_numpy() - min(df_gb_ug_f.index.to_numpy())) / (max(df_gb_ug_f.index.to_numpy()) - min(df_gb_ug_f.index.to_numpy()))
gb_wmax_f = - (df_gb_wmax_f.index.to_numpy() - max(df_gb_wmax_f.index.to_numpy())) / (max(df_gb_wmax_f.index.to_numpy()) - min(df_gb_wmax_f.index.to_numpy()))
gb_dtsurf_f = (df_gb_dtsurf_f.index.to_numpy() - min(df_gb_dtsurf_f.index.to_numpy())) / (max(df_gb_dtsurf_f.index.to_numpy()) - min(df_gb_dtsurf_f.index.to_numpy()))
    
wmax_wmax_f = (gb_wmax_f * (th_max-th_min)) + th_min
ug_ug_f = (gb_ug_f * (th_max-th_min)) + th_min
dtsurf_dtsurf_f = (gb_dtsurf_f * (th_max-th_min)) + th_min
cldN_cldN_f = (gb_cldN_f * (th_max-th_min)) + th_min
cldN_cldN_fn = (gb_cldN_f * (th_max-th_min)) + th_min

factor_f = max([max(df_gb_ug_f['duration']), max(df_gb_wmax_f['duration']), max(df_gb_dtsurf_f['duration']), max(df_gb_cldN_f['duration'])])


# In[8]:


fig, ax2 = plt.subplots()
fig.set_size_inches(7,5.5)


c2 = ax2.scatter(df_fig_fog['max_sfc_xm2'] * 1000, (df_fig_fog['xm23_col_m'] * 1000), 
                 c = df_fig_fog['cldN'], cmap = 'rainbow', s = 20)

cbar2 = fig.colorbar(c2)

ax2.scatter(df_gb_cldN_fn['max_sfc_xm2'] * 1000, df_gb_cldN_fn['xm23_col_m'] * 1000, 
            marker = 's', c = cldN_cldN_fn, cmap = 'rainbow', edgecolor = 'black', linewidth = 3, s = 120)

# ax2.set_title('Fog Liquid Water Path vs Max Surface CWC', fontsize = 20)
ax2.set_xlabel('Max Surface CWC (g/kg)', fontsize = 18)
ax2.set_ylabel('Liquid Water Path (g/m$^2$)', fontsize = 18)
cbar2.set_label('Aerosol Number Concentration (#/cm$^3$)', fontsize = 16)
ax2.tick_params(axis = 'both', which = 'both', labelsize = 14)

fname = desc + 'fig02.png'
plt.savefig(fname)
plt.show()
plt.close()


# In[9]:


fig, ax2 = plt.subplots()
fig.set_size_inches(10,6)


c2 = ax2.scatter(df_fig_fog['sfc_xm2_h'] * 1000, (df_fig_fog['xm23_col_h']), 
                 c = df_fig_fog['cldN'], cmap = 'rainbow')

cbar2 = fig.colorbar(c2)

ax2.scatter(df_gb_cldN_f['sfc_xm2_h'] * 1000, df_gb_cldN_f['xm23_col_h'], 
            marker = 's', c = cldN_cldN_fn, cmap = 'rainbow', edgecolor = 'black', linewidth = 3, s = 160)

ax2.set_title('Fog Liquid Water Path vs Surface CLWC', fontsize = 20)
ax2.set_xlabel('Max Surface CLWC (g/kg)', fontsize = 18)
ax2.set_ylabel('Liquid Water Path (mm)', fontsize = 18)
cbar2.set_label('Aerosol Number Concentration (#/cm^3)', fontsize = 16)
ax2.tick_params(axis = 'both', which = 'both', labelsize = 14)

fname = desc + 'outliers_h.png'
plt.savefig(fname)
plt.show()
plt.close()


# In[10]:


df_Area_mCLWC_depc_all = df_all[['cldN', 'duration', 'max_sfc_xm2', 'dur_max_subsat_thickness', 'sedc_o', 'sedc_d', 'dur_bot_xm2_m3', 'onset_hr', 'diss_hr']]
df_A = df_Area_mCLWC_depc_all[fog & cond_outlier & cond_adzblt].reset_index()

df_A['dep_vel'] = (df_A['sedc_d'] - df_A['sedc_o']) / ((df_A['diss_hr'] - df_A['onset_hr']) * 3600 * (df_A['dur_bot_xm2_m3']))

cldN_100 = df_A['cldN'] == 100.
cldN_150 = df_A['cldN'] == 150.
cldN_200 = df_A['cldN'] == 200.
cldN_300 = df_A['cldN'] == 300.
cldN_400 = df_A['cldN'] == 400.
cldN_500 = df_A['cldN'] == 500.
cldN_750 = df_A['cldN'] == 750.
cldN_1000 = df_A['cldN'] == 1000.
cldN_1250 = df_A['cldN'] == 1250.
cldN_1500 = df_A['cldN'] == 1500.

df_A_100 = df_A[cldN_100].reset_index()
df_A_150 = df_A[cldN_150].reset_index()
df_A_200 = df_A[cldN_200].reset_index()
df_A_300 = df_A[cldN_300].reset_index()
df_A_400 = df_A[cldN_400].reset_index()
df_A_500 = df_A[cldN_500].reset_index()
df_A_750 = df_A[cldN_750].reset_index()
df_A_1000 = df_A[cldN_1000].reset_index()
df_A_1250 = df_A[cldN_1250].reset_index()
df_A_1500 = df_A[cldN_1500].reset_index()

xp = sm.add_constant(np.linspace(0, 200))

x_100 = df_A_100['duration'] ** 2
x_100_with_intercept = sm.add_constant(x_100)
y_100 = df_A_100['max_sfc_xm2']
x_150 = df_A_150['duration'] ** 2
x_150_with_intercept = sm.add_constant(x_150)
y_150 = df_A_150['max_sfc_xm2']
x_200 = df_A_200['duration'] ** 2
x_200_with_intercept = sm.add_constant(x_200)
y_200 = df_A_200['max_sfc_xm2']
x_300 = df_A_300['duration'] ** 2
x_300_with_intercept = sm.add_constant(x_300)
y_300 = df_A_300['max_sfc_xm2']
x_400 = df_A_400['duration'] ** 2
x_400_with_intercept = sm.add_constant(x_400)
y_400 = df_A_400['max_sfc_xm2']
x_500 = df_A_500['duration'] ** 2
x_500_with_intercept = sm.add_constant(x_500)
y_500 = df_A_500['max_sfc_xm2']
x_750 = df_A_750['duration'] ** 2
x_750_with_intercept = sm.add_constant(x_750)
y_750 = df_A_750['max_sfc_xm2']
x_1000 = df_A_1000['duration'] ** 2
x_1000_with_intercept = sm.add_constant(x_1000)
y_1000 = df_A_1000['max_sfc_xm2']
x_1250 = df_A_1250['duration'] ** 2
x_1250_with_intercept = sm.add_constant(x_1250)
y_1250 = df_A_1250['max_sfc_xm2']
x_1500 = df_A_1500['duration'] ** 2
x_1500_with_intercept = sm.add_constant(x_1500)
y_1500 = df_A_1500['max_sfc_xm2']

# Fit the linear regression model
model_100 = sm.OLS(y_100, x_100_with_intercept).fit()
model_150 = sm.OLS(y_150, x_150_with_intercept).fit()
model_200 = sm.OLS(y_200, x_200_with_intercept).fit()
model_300 = sm.OLS(y_300, x_300_with_intercept).fit()
model_400 = sm.OLS(y_400, x_400_with_intercept).fit()
model_500 = sm.OLS(y_500, x_500_with_intercept).fit()
model_750 = sm.OLS(y_750, x_750_with_intercept).fit()
model_1000 = sm.OLS(y_1000, x_1000_with_intercept).fit()
model_1250 = sm.OLS(y_1250, x_1250_with_intercept).fit()
model_1500 = sm.OLS(y_1500, x_1500_with_intercept).fit()

# Predict y values
p_100 = model_100.get_prediction(xp)
p_150 = model_150.get_prediction(xp)
p_200 = model_200.get_prediction(xp)
p_300 = model_300.get_prediction(xp)
p_400 = model_400.get_prediction(xp)
p_500 = model_500.get_prediction(xp)
p_750 = model_750.get_prediction(xp)
p_1000 = model_1000.get_prediction(xp)
p_1250 = model_1250.get_prediction(xp)
p_1500 = model_1500.get_prediction(xp)

yp_100 = model_100.predict(xp)
yp_150 = p_150.predicted_mean
yp_200 = p_200.predicted_mean
yp_300 = p_300.predicted_mean
yp_400 = p_400.predicted_mean
yp_500 = p_500.predicted_mean
yp_750 = p_750.predicted_mean
yp_1000 = p_1000.predicted_mean
yp_1250 = p_1250.predicted_mean
yp_1500 = p_1500.predicted_mean

# Get the confidence intervals for the regression coefficients
ci_100 = p_100.conf_int(alpha=0.05)  # 95% confidence intervals
ci_150 = p_150.conf_int(alpha=0.05)
ci_200 = p_200.conf_int(alpha=0.05)
ci_300 = p_300.conf_int(alpha=0.05)
ci_400 = p_400.conf_int(alpha=0.05)
ci_500 = p_500.conf_int(alpha=0.05)
ci_750 = p_750.conf_int(alpha=0.05)
ci_1000 = p_1000.conf_int(alpha=0.05)
ci_1250 = p_1250.conf_int(alpha=0.05)
ci_1500 = p_1500.conf_int(alpha=0.05)

YP = np.vstack((yp_100, yp_150, yp_200, yp_300, yp_400, yp_500, yp_750, yp_1000, yp_1250, yp_1500))
CI_0 = np.vstack((ci_100[:,0], ci_150[:,0], ci_200[:,0], ci_300[:,0], ci_400[:,0], ci_500[:,0], ci_750[:,0], 
                  ci_1000[:,0], ci_1250[:,0], ci_1500[:,0]))
CI_1 = np.vstack((ci_100[:,1], ci_150[:,1], ci_200[:,1], ci_300[:,1], ci_400[:,1], ci_500[:,1], ci_750[:,1], 
                  ci_1000[:,1], ci_1250[:,1], ci_1500[:,1]))

xp = xp[:,1]

yp100_2 = model_100.predict(x_100_with_intercept)
yp150_2 = model_150.predict(x_150_with_intercept)
yp200_2 = model_200.predict(x_200_with_intercept)
yp300_2 = model_300.predict(x_300_with_intercept)
yp400_2 = model_400.predict(x_400_with_intercept)
yp500_2 = model_500.predict(x_500_with_intercept)
yp750_2 = model_750.predict(x_750_with_intercept)
yp1000_2 = model_1000.predict(x_1000_with_intercept)
yp1250_2 = model_1250.predict(x_1250_with_intercept)
yp1500_2 = model_1500.predict(x_1500_with_intercept)

r2_100 = 1 - np.sum((y_100 - yp100_2) ** 2) / np.sum((y_100 - np.mean(y_100)) ** 2)
r2_150 = 1 - np.sum((y_150 - yp150_2) ** 2) / np.sum((y_150 - np.mean(y_150)) ** 2)
r2_200 = 1 - np.sum((y_200 - yp200_2) ** 2) / np.sum((y_200 - np.mean(y_200)) ** 2)
r2_300 = 1 - np.sum((y_300 - yp300_2) ** 2) / np.sum((y_300 - np.mean(y_300)) ** 2)
r2_400 = 1 - np.sum((y_400 - yp400_2) ** 2) / np.sum((y_400 - np.mean(y_400)) ** 2)
r2_500 = 1 - np.sum((y_500 - yp500_2) ** 2) / np.sum((y_500 - np.mean(y_500)) ** 2)
r2_750 = 1 - np.sum((y_750 - yp750_2) ** 2) / np.sum((y_750 - np.mean(y_750)) ** 2)
r2_1000 = 1 - np.sum((y_1000 - yp1000_2) ** 2) / np.sum((y_1000 - np.mean(y_1000)) ** 2)
r2_1250 = 1 - np.sum((y_1250 - yp1250_2) ** 2) / np.sum((y_1250 - np.mean(y_1250)) ** 2)
r2_1500 = 1 - np.sum((y_1500 - yp1500_2) ** 2) / np.sum((y_1500 - np.mean(y_1500)) ** 2)

std_100 = np.mean((y_100 - yp100_2) ** 2) ** (0.5)
std_150 = np.mean((y_150 - yp150_2) ** 2) ** (0.5)
std_200 = np.mean((y_200 - yp200_2) ** 2) ** (0.5)
std_300 = np.mean((y_300 - yp300_2) ** 2) ** (0.5)
std_400 = np.mean((y_400 - yp400_2) ** 2) ** (0.5)
std_500 = np.mean((y_500 - yp500_2) ** 2) ** (0.5)
std_750 = np.mean((y_750 - yp750_2) ** 2) ** (0.5)
std_1000 = np.mean((y_1000 - yp1000_2) ** 2) ** (0.5)
std_1250 = np.mean((y_1250 - yp1250_2) ** 2) ** (0.5)
std_1500 = np.mean((y_1500 - yp1500_2) ** 2) ** (0.5)

std_a = [std_100, std_150, std_200, std_300, std_400, std_500, std_750, std_1000, std_1250, std_1500]

print(r2_100, r2_150, r2_200, r2_300, r2_400, r2_500, r2_750, r2_1000, r2_1250, r2_1500)
print(np.sum((y_100 - yp100_2) ** 2), np.sum(((y_100-np.mean(y_100)) ** 2)))


# In[11]:


# lablist for newrad2
lablist = ["N$_a$ = 100, R$^2$ = 0.07", "N$_a$ = 100, R$^2$ = 0.35", "N$_a$ = 100, R$^2$ = 0.39", 
           "N$_a$ = 100, R$^2$ = 0.60", "N$_a$ = 100, R$^2$ = 0.72", "N$_a$ = 100, R$^2$ = 0.79", 
           "N$_a$ = 100, R$^2$ = 0.82", "N$_a$ = 100, R$^2$ = 0.79", "N$_a$ = 100, R$^2$ = 0.82", 
           "N$_a$ = 100, R$^2$ = 0.85"]

# lablist for constxm
# lablist = ["N$_a$ = 100, R$^2$ = 0.00", "N$_a$ = 150, R$^2$ = 0.10", "N$_a$ = 200, R$^2$ = 0.22", 
#            "N$_a$ = 300, R$^2$ = 0.48", "N$_a$ = 400, R$^2$ = 0.72", "N$_a$ = 500, R$^2$ = 0.80", 
#            "N$_a$ = 750, R$^2$ = 0.79", "N$_a$ = 1000, R$^2$ = 0.81", "N$_a$ = 1250, R$^2$ = 0.85", 
#            "N$_a$ = 1500, R$^2$ = 0.86"]


fig_gb_cldN = df_A.groupby(['cldN']).agg(['mean', 'std']).reset_index()

# Create a colormap object from the Viridis colormap
viridis_cmap = plt.get_cmap('rainbow')

# Define the number of levels (in this case, 10)
num_levels = 10

# Create an array of equally spaced values from 0 to 1 (representing the levels)
levels = [17, 26, 34, 51, 68, 85, 128, 170, 212, 255] 

# Sample colors from the Viridis colormap at the specified levels
colors = viridis_cmap(levels)

# Create a figure and plot the data with the specified colors
fig, (ax1, ax2) = plt.subplots(1,2)
fig.set_size_inches(16,7)
plt.subplots_adjust(left=0.1, right=0.95, bottom=0.1, top=0.92, wspace=0.2, hspace=0.2)

# ax1.scatter((fig_gb_cldN['diss_hr'] - fig_gb_cldN['onset_hr']) ** 2, fig_gb_cldN['max_sfc_xm2'] * 1000, c = fig_gb_cldN['cldN'], cmap = 'viridis', s = 80)

for i in range(num_levels):
    ax1.plot(xp, YP[i,:] * 1000, label=lablist[i], color=colors[i], linewidth = 3)
    ax1.fill_between(xp, CI_0[i,:] * 1000, CI_1[i,:] * 1000, color=colors[i], alpha = 0.25)

# Add a legend to distinguish the levels
ax1.legend(fontsize = 15)
ax1.set_ylabel('Max Surface CWC (g/kg)', fontsize = 20)
ax1.set_xlabel('Duration Squared (Hours$^2$)', fontsize = 20)
# ax1.set_title('Trends of Max Sfc CWC Vs. Duration Squared', fontsize = 20)
ax1.tick_params(axis='both', which='both', labelsize=15)
ax1.text(0.9, 0.95, '(A)', transform = ax1.transAxes, fontsize = 18, fontweight = 'bold')

scatter = ax2.scatter(fig_gb_cldN['dur_max_subsat_thickness']['mean'], fig_gb_cldN['dep_vel']['mean'], s = 100, c = colors, cmap = 'rainbow')



for i in range(len(fig_gb_cldN['dur_max_subsat_thickness']['mean'])):
    ax2.errorbar(fig_gb_cldN['dur_max_subsat_thickness']['mean'][i], 
                fig_gb_cldN['dep_vel']['mean'][i], 
                xerr=fig_gb_cldN['dur_max_subsat_thickness']['std'][i], 
                yerr=fig_gb_cldN['dep_vel']['std'][i], fmt='none', 
                ecolor=colors[i])

ax2.set_ylabel('Surface Deposition Velocity (m/s)', fontsize = 20)
ax2.set_xlabel('Max Fog Subsaturated Layer Thickness (m)', fontsize = 20)
# ax2.set_title('Sfc Dep. Velocity vs. Max Fog Subsat. Thickness', fontsize = 20)
ax2.tick_params(axis='both', which='both', labelsize=15)
ax2.text(0.9, 0.95, '(B)', transform = ax2.transAxes, fontsize = 18, fontweight = 'bold')
    
fname = desc + 'fig05.png'
plt.savefig(fname)
plt.show()
plt.close()



# Create a figure and plot the data with the specified colors
fig, (ax1, ax2) = plt.subplots(1,2)
fig.set_size_inches(16,7)
plt.subplots_adjust(left=0.1, right=0.95, bottom=0.1, top=0.92, wspace=0.2, hspace=0.2)

# ax1.scatter((fig_gb_cldN['diss_hr'] - fig_gb_cldN['onset_hr']) ** 2, fig_gb_cldN['max_sfc_xm2'] * 1000, c = fig_gb_cldN['cldN'], cmap = 'viridis', s = 80)

for i in range(num_levels):
    ax1.plot(xp, YP[i,:] * 1000, label=lablist[i], color=colors[i], linewidth = 3)
    ax1.fill_between(xp, (YP[i,:] - 1 * std_a[i]) * 1000, (YP[i,:] + 1 * std_a[i]) * 1000, color=colors[i], alpha = 0.25)

# Add a legend to distinguish the levels
ax1.legend(fontsize = 14)
ax1.set_ylabel('Max Surface CLWC (g/kg)', fontsize = 20)
ax1.set_xlabel('Duration Squared (Hours^2)', fontsize = 20)
ax1.set_title('Trends of Max Sfc CLWC Vs. Duration Squared', fontsize = 20)
ax1.tick_params(axis='both', which='both', labelsize=15)

scatter = ax2.scatter(fig_gb_cldN['dur_max_subsat_thickness']['mean'], fig_gb_cldN['dep_vel']['mean'], s = 100, c = colors, cmap = 'rainbow')



for i in range(len(fig_gb_cldN['dur_max_subsat_thickness']['mean'])):
    ax2.errorbar(fig_gb_cldN['dur_max_subsat_thickness']['mean'][i], 
                fig_gb_cldN['dep_vel']['mean'][i], 
                xerr=fig_gb_cldN['dur_max_subsat_thickness']['std'][i], 
                yerr=fig_gb_cldN['dep_vel']['std'][i], fmt='none', 
                ecolor=colors[i])

ax2.set_ylabel('Surface Deposition Velocity (m/s)', fontsize = 20)
ax2.set_xlabel('Max Fog Subsaturated Layer Thickness (m)', fontsize = 20)
ax2.set_title('Sfc Dep. Velocity vs. Max Fog Subsat. Thickness', fontsize = 20)
ax2.tick_params(axis='both', which='both', labelsize=14)
    
fname = desc + 'Area_mCLWC_depc_2.png'
# plt.savefig(fname)
plt.show()
plt.close()


# In[12]:


df_this_test = df_all[cond_outlier & cond_adzblt].reset_index()
df_nout = df_all[cond_outlier].reset_index()

fcond = df_this_test['duration'] > 0
df_this_test_f = df_this_test[fcond]

df_tt_ccldN = df_this_test.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_tt_ccldN_f = df_this_test_f.groupby(['Ug', 'wmax', 'dtsurf']).mean()

cldNcorr = np.zeros(np.shape(df_tt_ccldN)[0])
cldN_satke_corr = np.zeros(np.shape(df_tt_ccldN)[0])
cldN_entV_corr = np.zeros(np.shape(df_tt_ccldN)[0])
w_mean_cldN = np.zeros(np.shape(df_tt_ccldN)[0])
cldNcorr_b = np.zeros(np.shape(df_tt_ccldN)[0])
nf = np.zeros(np.shape(df_tt_ccldN)[0])
fogdur = np.zeros(np.shape(df_tt_ccldN)[0])
cldNimpact = np.zeros(np.shape(df_tt_ccldN)[0])
cldN_satke_impact = np.zeros(np.shape(df_tt_ccldN)[0])
cldN_entV_impact = np.zeros(np.shape(df_tt_ccldN)[0])
cldNimpact_b = np.zeros(np.shape(df_tt_ccldN)[0])
j = 0

for i in df_tt_ccldN.index:
    Ugcond = df_this_test['Ug'] == i[0]
    wmaxcond = df_this_test['wmax'] == i[1]
    dtscond = df_this_test['dtsurf'] == i[2]
    df_i = df_this_test[Ugcond & wmaxcond & dtscond]
    df_i_f = df_this_test[Ugcond & wmaxcond & dtscond & fcond]
    if np.nanmean(df_i['duration'] > 0.0):
        fogdur[j] = np.nanmean(df_i_f['duration'])
    cldcorr = np.corrcoef(df_i['cldN'], df_i['duration'])[0,1]
    satke_corr = np.corrcoef(df_i['cldN'], df_i['sat_tke'])[0,1]
    entV_corr = np.corrcoef(df_i['cldN'], df_i['entrainment_z'])[0,1]
    
    nt = np.shape(df_i)[0]
    if nt > 1.5:
        factor = 9. / (nt-1)
        cld_impact = np.std(df_i['duration']) * cldcorr * factor
        satke_impact = np.std(df_i['sat_tke']) * satke_corr * factor
        entV_impact = np.std(df_i['entrainment_z']) * entV_corr * factor
    else:
        cld_impact = np.nan
        satke_impact = np.nan
        entV_impact = np.nan
    cldNimpact[j] = cld_impact
    cldN_satke_impact[j] = satke_impact
    cldN_entV_impact[j] = entV_impact
    
    nf[j] = np.shape(df_i_f)[0]
    cldNcorr[j] = cldcorr
    cldN_satke_corr[j] = satke_corr
    cldN_entV_corr[j] = entV_corr
    
    
    if np.shape(df_i_f)[0] > 0:
        maxcldN = np.nanmax(df_i_f['cldN'])
        mincldN = np.nanmin(df_i_f['cldN'])
        maxNi = np.min([np.where(cldN_list == cldN_list[-1])[0], np.where(cldN_list == maxcldN)[0] + 1])
        minNi = np.max([np.where(cldN_list == cldN_list[0])[0], np.where(cldN_list == mincldN)[0] - 1])
        maxNcond = df_i['cldN'] <= cldN_list[maxNi]
        minNcond = df_i['cldN'] >= cldN_list[minNi]
        df_i_cldN = df_i[maxNcond & minNcond]
        cldcorrb = np.corrcoef(df_i_cldN['cldN'], df_i_cldN['duration'])[0,1]
        cldNcorr_b[j] = cldcorrb
        
        if np.shape(df_i_cldN)[0] > 1.5:
            factor_b = 9. / (np.shape(df_i_cldN)[0] - 1)
            cld_impact_b = np.std(df_i_cldN['duration']) * cldcorrb * factor_b
        else:
            cld_impact_b = np.nan    
        
        cldNimpact_b[j] = cld_impact_b
        
        wave_cldN = np.sum(df_i_f['cldN'] * df_i_f['duration']) / np.sum(df_i_f['duration'])
        w_mean_cldN[j] = wave_cldN
    else:
        cldNcorr_b[j] = float("NaN")
        cldNimpact_b[j] = float("NaN")
        w_mean_cldN[j] = float("NaN")
    
    j = j+1
    
df_tt_ccldN['cldN_corr'] = cldNcorr
df_tt_ccldN['cldN_satke_corr'] = cldN_satke_corr
df_tt_ccldN['cldN_entV_corr'] = cldN_entV_corr
df_tt_ccldN['cldN_corr_b'] = cldNcorr_b

df_tt_ccldN['cldN_impact'] = cldNimpact
df_tt_ccldN['cldN_satke_impact'] = cldN_satke_impact
df_tt_ccldN['cldN_entV_impact'] = cldN_entV_impact
df_tt_ccldN['cldN_impact_b'] = cldNimpact_b

df_tt_ccldN['fog_duration'] = fogdur
df_tt_ccldN['w_mean_cldN'] = w_mean_cldN
df_tt_ccldN['alpha'] = df_tt_ccldN['fog_duration'] / np.nanmax(df_tt_ccldN['fog_duration'])
df_tt_ccldN['nf'] = nf

nancond2 = df_tt_ccldN['cldN_corr'].notna()
df_tt_ccldN_ccr = df_tt_ccldN[nancond2]


# In[13]:


cldN_list = np.array([100, 150, 200, 300, 400, 500, 750, 1000, 1250, 1500])

df_vars = df_nout[['Ug', 'wmax', 'dtsurf', 'cldN', 'duration', 't_dur', 'nl_dur', 'sedt_f', 'ent_xm_ave', 
                     'ent_xm_ave2', 'ent_sed', 'entrainment_z', 'xm_bl_z_m', 'm_bl_zrat_m', 'xm_bl_zrat_m', 
                     'xm23_bl_z_m', 'xm23_bl_zrat_m', 'RH_maxz_m', 'xm23_maxz_m', 'xm_bl_z_h', 'm_bl_zrat_h', 
                     'xm_bl_zrat_h', 'xm23_bl_z_h', 'xm23_bl_zrat_h', 'RH_maxz_h', 'xm23_maxz_h', 'sedt_h', 
                     'ent_xm_ave_h', 'ent_xm_ave_h2', 'ent_sed_h', 'entrainment_zh', 'dz_blt_h']]

durs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
h_dzblt = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_h_dzblt = np.empty(np.shape(df_tt_ccldN_ccr)[0])

f_durs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_durweights = np.empty(np.shape(df_tt_ccldN_ccr)[0])
t_durs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
t_durweights = np.empty(np.shape(df_tt_ccldN_ccr)[0])
nl_durs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
nl_durweights = np.empty(np.shape(df_tt_ccldN_ccr)[0])

sedts = np.empty(np.shape(df_tt_ccldN_ccr)[0])
ent_seds = np.empty(np.shape(df_tt_ccldN_ccr)[0])
ent_xms = np.empty(np.shape(df_tt_ccldN_ccr)[0])
ent_vs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
sedts_h = np.empty(np.shape(df_tt_ccldN_ccr)[0])
ent_seds_h = np.empty(np.shape(df_tt_ccldN_ccr)[0])
ent_xms_h = np.empty(np.shape(df_tt_ccldN_ccr)[0])
ent_vs_h = np.empty(np.shape(df_tt_ccldN_ccr)[0])
xm_bl_zs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
xm23_bl_zs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
xm_bl_zrats = np.empty(np.shape(df_tt_ccldN_ccr)[0])
xm23_bl_zrats = np.empty(np.shape(df_tt_ccldN_ccr)[0])
RH_maxzs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
xm23_maxzs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
bg_weights = np.empty(np.shape(df_tt_ccldN_ccr)[0])
h_xm_bl_zs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
h_xm23_bl_zs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
h_xm_bl_zrats = np.empty(np.shape(df_tt_ccldN_ccr)[0])
h_xm23_bl_zrats = np.empty(np.shape(df_tt_ccldN_ccr)[0])
h_RH_maxzs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
h_xm23_maxzs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
h_bg_weights = np.empty(np.shape(df_tt_ccldN_ccr)[0])

f_sedts = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_ent_seds = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_ent_xms = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_ent_vs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_sedts_h = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_ent_seds_h = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_ent_xms_h = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_ent_vs_h = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_xm_bl_zs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_xm23_bl_zs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_xm_bl_zrats = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_xm23_bl_zrats = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_RH_maxzs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_xm23_maxzs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_h_xm_bl_zs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_h_xm23_bl_zs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_h_xm_bl_zrats = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_h_xm23_bl_zrats = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_h_RH_maxzs = np.empty(np.shape(df_tt_ccldN_ccr)[0])
f_h_xm23_maxzs = np.empty(np.shape(df_tt_ccldN_ccr)[0])

durs_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
durs_cldN[:] = np.nan

h_dzblt_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
h_dzblt_cldN[:] = np.nan
f_h_dzblt_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_h_dzblt_cldN[:] = np.nan

f_durs_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_durs_cldN[:] = np.nan
t_durs_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
t_durs_cldN[:] = np.nan
nl_durs_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
nl_durs_cldN[:] = np.nan

sedt_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
sedt_cldN[:] = np.nan
ent_xm_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
ent_xm_cldN[:] = np.nan
ent_sed_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
ent_sed_cldN[:] = np.nan
ent_v_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
ent_v_cldN[:] = np.nan
sedt_cldN_h = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
sedt_cldN_h[:] = np.nan
ent_xm_cldN_h = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
ent_xm_cldN_h[:] = np.nan
ent_sed_cldN_h = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
ent_sed_cldN_h[:] = np.nan
ent_v_cldN_h = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
ent_v_cldN_h[:] = np.nan
xm_bl_z_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
xm_bl_z_cldN[:] = np.nan
xm23_bl_z_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
xm23_bl_z_cldN[:] = np.nan
RH_maxz_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
RH_maxz_cldN[:] = np.nan
xm23_maxz_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
xm23_maxz_cldN[:] = np.nan
xm_bl_zrat_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
xm_bl_zrat_cldN[:] = np.nan
xm23_bl_zrat_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
xm23_bl_zrat_cldN[:] = np.nan
h_xm_bl_z_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
h_xm_bl_z_cldN[:] = np.nan
h_xm23_bl_z_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
h_xm23_bl_z_cldN[:] = np.nan
h_RH_maxz_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
h_RH_maxz_cldN[:] = np.nan
h_xm23_maxz_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
h_xm23_maxz_cldN[:] = np.nan
h_xm_bl_zrat_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
h_xm_bl_zrat_cldN[:] = np.nan
h_xm23_bl_zrat_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
h_xm23_bl_zrat_cldN[:] = np.nan

f_sedt_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_sedt_cldN[:] = np.nan
f_ent_xm_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_ent_xm_cldN[:] = np.nan
f_ent_sed_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_ent_sed_cldN[:] = np.nan
f_ent_v_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_ent_v_cldN[:] = np.nan
f_sedt_cldN_h = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_sedt_cldN_h[:] = np.nan
f_ent_xm_cldN_h = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_ent_xm_cldN_h[:] = np.nan
f_ent_sed_cldN_h = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_ent_sed_cldN_h[:] = np.nan
f_ent_v_cldN_h = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_ent_v_cldN_h[:] = np.nan
f_xm_bl_z_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_xm_bl_z_cldN[:] = np.nan
f_xm23_bl_z_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_xm23_bl_z_cldN[:] = np.nan
f_RH_maxz_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_RH_maxz_cldN[:] = np.nan
f_xm23_maxz_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_xm23_maxz_cldN[:] = np.nan
f_xm_bl_zrat_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_xm_bl_zrat_cldN[:] = np.nan
f_xm23_bl_zrat_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_xm23_bl_zrat_cldN[:] = np.nan
f_h_xm_bl_z_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_h_xm_bl_z_cldN[:] = np.nan
f_h_xm23_bl_z_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_h_xm23_bl_z_cldN[:] = np.nan
f_h_RH_maxz_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_h_RH_maxz_cldN[:] = np.nan
f_h_xm23_maxz_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_h_xm23_maxz_cldN[:] = np.nan
f_h_xm_bl_zrat_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_h_xm_bl_zrat_cldN[:] = np.nan
f_h_xm23_bl_zrat_cldN = np.empty((np.shape(df_tt_ccldN_ccr)[0], 10))
f_h_xm23_bl_zrat_cldN[:] = np.nan

for i in range(np.shape(df_tt_ccldN_ccr)[0]):
    (Ug, wmax, dtsurf) =  df_tt_ccldN_ccr.index[i]
    cond_Ug = df_vars['Ug'] == Ug
    cond_wmax = df_vars['wmax'] == wmax
    cond_dtsurf = df_vars['dtsurf'] == dtsurf
    cond_fog = df_vars['duration'] > 0.0
    cond_t_fog = df_vars['t_dur'] > 0.0
    cond_nl_fog = df_vars['nl_dur'] > 0.0

    
    df_this_bg = df_vars[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_vars[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    df_this_bg_t_fog = df_vars[cond_Ug & cond_wmax & cond_dtsurf & cond_t_fog].reset_index()
    df_this_bg_nl_fog = df_vars[cond_Ug & cond_wmax & cond_dtsurf & cond_nl_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    this_t_fog_cldN = df_this_bg_t_fog['cldN'].to_numpy()
    this_nl_fog_cldN = df_this_bg_nl_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    bg_maxNi = max(np.where(cldN_list == bg_maxcldN)[0])
    bg_minNi = min(np.where(cldN_list == bg_mincldN)[0])
    bg_maxNcond = df_this_bg['cldN'] <= cldN_list[bg_maxNi]
    bg_minNcond = df_this_bg['cldN'] >= cldN_list[bg_minNi]
    bg_df_cldN = df_this_bg[bg_maxNcond & bg_minNcond].reset_index()
    
    durs[i] = np.nanmean(bg_df_cldN['duration'])
    h_dzblt[i] = np.nanmean(bg_df_cldN['dz_blt_h'])
    bg_weights[i] = np.shape(bg_df_cldN)[0]
    sedts[i] = np.nanmean(bg_df_cldN['sedt_f'])
    ent_xms[i] = np.nanmean(bg_df_cldN['ent_xm_ave2'])
    ent_seds[i] = np.nanmean(bg_df_cldN['ent_sed'])
    ent_vs[i] = np.nanmean(bg_df_cldN['entrainment_z'])
    sedts_h[i] = np.nanmean(bg_df_cldN['sedt_h'])
    ent_xms_h[i] = np.nanmean(bg_df_cldN['ent_xm_ave_h2'])
    ent_seds_h[i] = np.nanmean(bg_df_cldN['ent_sed_h'])
    ent_vs_h[i] = np.nanmean(bg_df_cldN['entrainment_zh'])
    xm_bl_zs[i] = np.nanmean(bg_df_cldN['xm_bl_z_m'])
    xm23_bl_zs[i] = np.nanmean(bg_df_cldN['xm23_bl_z_m'])
    xm_bl_zrats[i] = np.nanmean(bg_df_cldN['xm_bl_zrat_m'])
    xm23_bl_zrats[i] = np.nanmean(bg_df_cldN['xm23_bl_zrat_m'])
    RH_maxzs[i] = np.nanmean(bg_df_cldN['RH_maxz_m'])
    xm23_maxzs[i] = np.nanmean(bg_df_cldN['xm23_maxz_m'])
    h_xm_bl_zs[i] = np.nanmean(bg_df_cldN['xm_bl_z_h'])
    h_xm23_bl_zs[i] = np.nanmean(bg_df_cldN['xm23_bl_z_h'])
    h_xm_bl_zrats[i] = np.nanmean(bg_df_cldN['xm_bl_zrat_h'])
    h_xm23_bl_zrats[i] = np.nanmean(bg_df_cldN['xm23_bl_zrat_h'])
    h_RH_maxzs[i] = np.nanmean(bg_df_cldN['RH_maxz_h'])
    h_xm23_maxzs[i] = np.nanmean(bg_df_cldN['xm23_maxz_h'])
    
    for j in range(10):
        if cldN_list[j] in bg_df_cldN['cldN'].to_numpy():
            index = np.where(bg_df_cldN['cldN'] == cldN_list[j])[0]
            durs_cldN[i,j] = bg_df_cldN['duration'][index] - durs[i]
            h_dzblt_cldN[i,j] = bg_df_cldN['dz_blt_h'][index] - h_dzblt[i]
            sedt_cldN[i,j] = bg_df_cldN['sedt_f'][index] - sedts[i]
            ent_xm_cldN[i,j] = bg_df_cldN['ent_xm_ave2'][index] - ent_xms[i]
            ent_sed_cldN[i,j] = bg_df_cldN['ent_sed'][index] - ent_seds[i]
            ent_v_cldN[i,j] = bg_df_cldN['entrainment_z'][index] - ent_vs[i]
            sedt_cldN_h[i,j] = bg_df_cldN['sedt_h'][index] - sedts_h[i]
            ent_xm_cldN_h[i,j] = bg_df_cldN['ent_xm_ave_h2'][index] - ent_xms_h[i]
            ent_sed_cldN_h[i,j] = bg_df_cldN['ent_sed_h'][index] - ent_seds_h[i]
            ent_v_cldN_h[i,j] = bg_df_cldN['entrainment_zh'][index] - ent_vs_h[i]
            xm_bl_z_cldN[i,j] = bg_df_cldN['xm_bl_z_m'][index] - xm_bl_zs[i]
            xm23_bl_z_cldN[i,j] = bg_df_cldN['xm23_bl_z_m'][index] - xm23_bl_zs[i]
            xm_bl_zrat_cldN[i,j] = bg_df_cldN['xm_bl_zrat_m'][index] - xm_bl_zrats[i]
            xm23_bl_zrat_cldN[i,j] = bg_df_cldN['xm23_bl_zrat_m'][index] - xm23_bl_zrats[i]
            RH_maxz_cldN[i,j] = bg_df_cldN['RH_maxz_m'][index] - RH_maxzs[i]
            xm23_maxz_cldN[i,j] = bg_df_cldN['xm23_maxz_m'][index] - xm23_maxzs[i]
            h_xm_bl_z_cldN[i,j] = bg_df_cldN['xm_bl_z_h'][index] - h_xm_bl_zs[i]
            h_xm23_bl_z_cldN[i,j] = bg_df_cldN['xm23_bl_z_h'][index] - h_xm23_bl_zs[i]
            h_xm_bl_zrat_cldN[i,j] = bg_df_cldN['xm_bl_zrat_h'][index] - h_xm_bl_zrats[i]
            h_xm23_bl_zrat_cldN[i,j] = bg_df_cldN['xm23_bl_zrat_h'][index] - h_xm23_bl_zrats[i]
            h_RH_maxz_cldN[i,j] = bg_df_cldN['RH_maxz_h'][index] - h_RH_maxzs[i]
            h_xm23_maxz_cldN[i,j] = bg_df_cldN['xm23_maxz_h'][index] - h_xm23_maxzs[i]
        else:
            durs_cldN[i,j] = np.nan
            h_dzblt_cldN[i,j] = np.nan
            sedt_cldN[i,j] = np.nan
            ent_xm_cldN[i,j] = np.nan
            ent_sed_cldN[i,j] = np.nan
            ent_v_cldN[i,j] = np.nan
            sedt_cldN_h[i,j] = np.nan
            ent_xm_cldN_h[i,j] = np.nan
            ent_sed_cldN_h[i,j] = np.nan
            ent_v_cldN_h[i,j] = np.nan
            xm_bl_z_cldN[i,j] = np.nan
            xm23_bl_z_cldN[i,j] = np.nan
            xm_bl_zrat_cldN[i,j] = np.nan
            xm23_bl_zrat_cldN[i,j] = np.nan
            RH_maxz_cldN[i,j] = np.nan
            xm23_maxz_cldN[i,j] = np.nan
            h_xm_bl_z_cldN[i,j] = np.nan
            h_xm23_bl_z_cldN[i,j] = np.nan
            h_xm_bl_zrat_cldN[i,j] = np.nan
            h_xm23_bl_zrat_cldN[i,j] = np.nan
            h_RH_maxz_cldN[i,j] = np.nan
            h_xm23_maxz_cldN[i,j] = np.nan
    
    if np.shape(this_fog_cldN)[0] > 0:
        f_maxcldN = np.nanmax(this_fog_cldN)
        f_mincldN = np.nanmin(this_fog_cldN)
        f_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == f_maxcldN)[0] + 1])
        f_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == f_mincldN)[0] - 1])
        f_maxNcond = df_this_bg['cldN'] <= cldN_list[f_maxNi]
        f_minNcond = df_this_bg['cldN'] >= cldN_list[f_minNi]
        f_df_cldN = df_this_bg[f_maxNcond & f_minNcond].reset_index()
        
        f_durs[i] = np.nanmean(f_df_cldN['duration'])
        f_h_dzblt[i] = np.nanmean(f_df_cldN['dz_blt_h'])
        f_durweights[i] = np.shape(f_df_cldN)[0]
        f_sedts[i] = np.nanmean(f_df_cldN['sedt_f'])
        f_ent_xms[i] = np.nanmean(f_df_cldN['ent_xm_ave2'])
        f_ent_seds[i] = np.nanmean(f_df_cldN['ent_sed'])
        f_ent_vs[i] = np.nanmean(f_df_cldN['entrainment_z'])
        f_sedts_h[i] = np.nanmean(f_df_cldN['sedt_h'])
        f_ent_xms_h[i] = np.nanmean(f_df_cldN['ent_xm_ave_h2'])
        f_ent_seds_h[i] = np.nanmean(f_df_cldN['ent_sed_h'])
        f_ent_vs_h[i] = np.nanmean(f_df_cldN['entrainment_zh'])
        f_xm_bl_zs[i] = np.nanmean(f_df_cldN['xm_bl_z_m'])
        f_xm23_bl_zs[i] = np.nanmean(f_df_cldN['xm23_bl_z_m'])
        f_xm_bl_zrats[i] = np.nanmean(f_df_cldN['xm_bl_zrat_m'])
        f_xm23_bl_zrats[i] = np.nanmean(f_df_cldN['xm23_bl_zrat_m'])
        f_RH_maxzs[i] = np.nanmean(f_df_cldN['RH_maxz_m'])
        f_xm23_maxzs[i] = np.nanmean(f_df_cldN['xm23_maxz_m'])
        f_h_xm_bl_zs[i] = np.nanmean(f_df_cldN['xm_bl_z_h'])
        f_h_xm23_bl_zs[i] = np.nanmean(f_df_cldN['xm23_bl_z_h'])
        f_h_xm_bl_zrats[i] = np.nanmean(f_df_cldN['xm_bl_zrat_h'])
        f_h_xm23_bl_zrats[i] = np.nanmean(f_df_cldN['xm23_bl_zrat_h'])
        f_h_RH_maxzs[i] = np.nanmean(f_df_cldN['RH_maxz_h'])
        f_h_xm23_maxzs[i] = np.nanmean(f_df_cldN['xm23_maxz_h'])
        
        for j in range(10):
            if cldN_list[j] in f_df_cldN['cldN'].to_numpy():
                index = np.where(f_df_cldN['cldN'] == cldN_list[j])[0]
                f_durs_cldN[i,j] = f_df_cldN['duration'][index] - f_durs[i]
                f_h_dzblt_cldN[i,j] = f_df_cldN['dz_blt_h'][index] - f_h_dzblt[i]
                f_sedt_cldN[i,j] = f_df_cldN['sedt_f'][index] - f_sedts[i]
                f_ent_xm_cldN[i,j] = f_df_cldN['ent_xm_ave2'][index] - f_ent_xms[i]
                f_ent_sed_cldN[i,j] = f_df_cldN['ent_sed'][index] - f_ent_seds[i]
                f_ent_v_cldN[i,j] = f_df_cldN['entrainment_z'][index] - f_ent_vs[i]
                f_sedt_cldN_h[i,j] = f_df_cldN['sedt_h'][index] - f_sedts_h[i]
                f_ent_xm_cldN_h[i,j] = f_df_cldN['ent_xm_ave_h2'][index] - f_ent_xms_h[i]
                f_ent_sed_cldN_h[i,j] = f_df_cldN['ent_sed_h'][index] - f_ent_seds_h[i]
                f_ent_v_cldN_h[i,j] = f_df_cldN['entrainment_zh'][index] - f_ent_vs_h[i]
                f_xm_bl_z_cldN[i,j] = f_df_cldN['xm_bl_z_m'][index] - f_xm_bl_zs[i]
                f_xm23_bl_z_cldN[i,j] = f_df_cldN['xm23_bl_z_m'][index] - f_xm23_bl_zs[i]
                f_xm_bl_zrat_cldN[i,j] = f_df_cldN['xm_bl_zrat_m'][index] - f_xm_bl_zrats[i]
                f_xm23_bl_zrat_cldN[i,j] = f_df_cldN['xm23_bl_zrat_m'][index] - f_xm23_bl_zrats[i]
                f_RH_maxz_cldN[i,j] = f_df_cldN['RH_maxz_m'][index] - f_RH_maxzs[i]
                f_xm23_maxz_cldN[i,j] = f_df_cldN['xm23_maxz_m'][index] - f_xm23_maxzs[i]
                f_h_xm_bl_z_cldN[i,j] = f_df_cldN['xm_bl_z_h'][index] - f_h_xm_bl_zs[i]
                f_h_xm23_bl_z_cldN[i,j] = f_df_cldN['xm23_bl_z_h'][index] - f_h_xm23_bl_zs[i]
                f_h_xm_bl_zrat_cldN[i,j] = f_df_cldN['xm_bl_zrat_h'][index] - f_h_xm_bl_zrats[i]
                f_h_xm23_bl_zrat_cldN[i,j] = f_df_cldN['xm23_bl_zrat_h'][index] - f_h_xm23_bl_zrats[i]
                f_h_RH_maxz_cldN[i,j] = f_df_cldN['RH_maxz_h'][index] - f_h_RH_maxzs[i]
                f_h_xm23_maxz_cldN[i,j] = f_df_cldN['xm23_maxz_h'][index] - f_h_xm23_maxzs[i]
                
            else:
                f_durs_cldN[i,j] = np.nan
                f_h_dzblt_cldN[i,j] = np.nan
                f_sedt_cldN[i,j] = np.nan
                f_ent_xm_cldN[i,j] = np.nan
                f_ent_sed_cldN[i,j] = np.nan
                f_ent_v_cldN[i,j] = np.nan
                f_sedt_cldN_h[i,j] = np.nan
                f_ent_xm_cldN_h[i,j] = np.nan
                f_ent_sed_cldN_h[i,j] = np.nan
                f_ent_v_cldN_h[i,j] = np.nan
                f_xm_bl_z_cldN[i,j] = np.nan
                f_xm23_bl_z_cldN[i,j] = np.nan
                f_xm_bl_zrat_cldN[i,j] = np.nan
                f_xm23_bl_zrat_cldN[i,j] = np.nan
                f_RH_maxz_cldN[i,j] = np.nan
                f_xm23_maxz_cldN[i,j] = np.nan
                f_h_xm_bl_z_cldN[i,j] = np.nan
                f_h_xm23_bl_z_cldN[i,j] = np.nan
                f_h_xm_bl_zrat_cldN[i,j] = np.nan
                f_h_xm23_bl_zrat_cldN[i,j] = np.nan
                f_h_RH_maxz_cldN[i,j] = np.nan
                f_h_xm23_maxz_cldN[i,j] = np.nan
                
    if np.shape(this_t_fog_cldN)[0] > 0:
        t_maxcldN = np.nanmax(this_t_fog_cldN)
        t_mincldN = np.nanmin(this_t_fog_cldN)
        t_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == t_maxcldN)[0] + 1])
        t_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == t_mincldN)[0] - 1])
        t_maxNcond = df_this_bg['cldN'] <= cldN_list[t_maxNi]
        t_minNcond = df_this_bg['cldN'] >= cldN_list[t_minNi]
        t_df_cldN = df_this_bg[t_maxNcond & t_minNcond].reset_index()
        
        t_durs[i] = np.nanmean(t_df_cldN['t_dur'])
        t_durweights[i] = np.shape(t_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in t_df_cldN['cldN'].to_numpy():
                index = np.where(t_df_cldN['cldN'] == cldN_list[j])[0]
                t_durs_cldN[i,j] = t_df_cldN['t_dur'][index] - t_durs[i]
            else:
                t_durs_cldN[i,j] = np.nan
                
    if np.shape(this_nl_fog_cldN)[0] > 0:
        nl_maxcldN = np.nanmax(this_nl_fog_cldN)
        nl_mincldN = np.nanmin(this_nl_fog_cldN)
        nl_maxNi = np.min([np.where(cldN_list == bg_maxcldN)[0], np.where(cldN_list == nl_maxcldN)[0] + 1])
        nl_minNi = np.max([np.where(cldN_list == bg_mincldN)[0], np.where(cldN_list == nl_mincldN)[0] - 1])
        nl_maxNcond = df_this_bg['cldN'] <= cldN_list[nl_maxNi]
        nl_minNcond = df_this_bg['cldN'] >= cldN_list[nl_minNi]
        nl_df_cldN = df_this_bg[nl_maxNcond & nl_minNcond].reset_index()
        
        nl_durs[i] = np.nanmean(nl_df_cldN['nl_dur'])
        nl_durweights[i] = np.shape(nl_df_cldN)[0]
        
        for j in range(10):
            if cldN_list[j] in nl_df_cldN['cldN'].to_numpy():
                index = np.where(nl_df_cldN['cldN'] == cldN_list[j])[0]
                nl_durs_cldN[i,j] = nl_df_cldN['nl_dur'][index] - nl_durs[i]
            else:
                nl_durs_cldN[i,j] = np.nan


durchange = np.nanmean(durs_cldN, axis = 0)                    
durstd = np.nanstd(durs_cldN, axis = 0)
dcount = np.count_nonzero(~np.isnan(durs_cldN), axis=0)
h_dzblt_change = np.nanmean(h_dzblt_cldN, axis = 0)                    
h_dzblt_std = np.nanstd(h_dzblt_cldN, axis = 0)
h_dzblt_dcount = np.count_nonzero(~np.isnan(h_dzblt_cldN), axis=0)
f_h_dzblt_change = np.nanmean(f_h_dzblt_cldN, axis = 0)                    
f_h_dzblt_std = np.nanstd(f_h_dzblt_cldN, axis = 0)
f_h_dzblt_dcount = np.count_nonzero(~np.isnan(f_h_dzblt_cldN), axis=0)

f_durchange = np.nanmean(f_durs_cldN, axis = 0)
t_durchange = np.nanmean(t_durs_cldN, axis = 0)
nl_durchange = np.nanmean(nl_durs_cldN, axis = 0)
f_durstd = np.nanstd(f_durs_cldN, axis = 0)
t_durstd = np.nanstd(t_durs_cldN, axis = 0)
nl_durstd = np.nanstd(nl_durs_cldN, axis = 0)
f_dcount = np.count_nonzero(~np.isnan(f_durs_cldN), axis=0)
t_dcount = np.count_nonzero(~np.isnan(t_durs_cldN), axis=0)
nl_dcount = np.count_nonzero(~np.isnan(nl_durs_cldN), axis=0)

sedt_change = np.nanmean(sedt_cldN, axis = 0)
ent_xm_change = np.nanmean(ent_xm_cldN, axis = 0)
ent_sed_change = np.nanmean(ent_sed_cldN, axis = 0)
ent_v_change = np.nanmean(ent_v_cldN, axis = 0)
sedt_change_h = np.nanmean(sedt_cldN_h, axis = 0)
ent_xm_change_h = np.nanmean(ent_xm_cldN_h, axis = 0)
ent_sed_change_h = np.nanmean(ent_sed_cldN_h, axis = 0)
ent_v_change_h = np.nanmean(ent_v_cldN_h, axis = 0)
xm_bl_z_change = np.nanmean(xm_bl_z_cldN, axis = 0)
xm23_bl_z_change = np.nanmean(xm23_bl_z_cldN, axis = 0)
xm_bl_zrat_change = np.nanmean(xm_bl_zrat_cldN, axis = 0)
xm23_bl_zrat_change = np.nanmean(xm23_bl_zrat_cldN, axis = 0)
RH_maxz_change = np.nanmean(RH_maxz_cldN, axis = 0)
xm23_maxz_change = np.nanmean(xm23_maxz_cldN, axis = 0)
h_xm_bl_z_change = np.nanmean(h_xm_bl_z_cldN, axis = 0)
h_xm23_bl_z_change = np.nanmean(h_xm23_bl_z_cldN, axis = 0)
h_xm_bl_zrat_change = np.nanmean(h_xm_bl_zrat_cldN, axis = 0)
h_xm23_bl_zrat_change = np.nanmean(h_xm23_bl_zrat_cldN, axis = 0)
h_RH_maxz_change = np.nanmean(h_RH_maxz_cldN, axis = 0)
h_xm23_maxz_change = np.nanmean(h_xm23_maxz_cldN, axis = 0)

f_sedt_change = np.nanmean(f_sedt_cldN, axis = 0)
f_ent_xm_change = np.nanmean(f_ent_xm_cldN, axis = 0)
f_ent_sed_change = np.nanmean(f_ent_sed_cldN, axis = 0)
f_ent_v_change = np.nanmean(f_ent_v_cldN, axis = 0)
f_sedt_change_h = np.nanmean(f_sedt_cldN_h, axis = 0)
f_ent_xm_change_h = np.nanmean(f_ent_xm_cldN_h, axis = 0)
f_ent_sed_change_h = np.nanmean(f_ent_sed_cldN_h, axis = 0)
f_ent_v_change_h = np.nanmean(f_ent_v_cldN_h, axis = 0)
f_xm_bl_z_change = np.nanmean(f_xm_bl_z_cldN, axis = 0)
f_xm23_bl_z_change = np.nanmean(f_xm23_bl_z_cldN, axis = 0)
f_xm_bl_zrat_change = np.nanmean(f_xm_bl_zrat_cldN, axis = 0)
f_xm23_bl_zrat_change = np.nanmean(f_xm23_bl_zrat_cldN, axis = 0)
f_RH_maxz_change = np.nanmean(f_RH_maxz_cldN, axis = 0)
f_xm23_maxz_change = np.nanmean(f_xm23_maxz_cldN, axis = 0)
f_h_xm_bl_z_change = np.nanmean(f_h_xm_bl_z_cldN, axis = 0)
f_h_xm23_bl_z_change = np.nanmean(f_h_xm23_bl_z_cldN, axis = 0)
f_h_xm_bl_zrat_change = np.nanmean(f_h_xm_bl_zrat_cldN, axis = 0)
f_h_xm23_bl_zrat_change = np.nanmean(f_h_xm23_bl_zrat_cldN, axis = 0)
f_h_RH_maxz_change = np.nanmean(f_h_RH_maxz_cldN, axis = 0)
f_h_xm23_maxz_change = np.nanmean(f_h_xm23_maxz_cldN, axis = 0)

sedt_std = np.nanstd(sedt_cldN, axis = 0)
ent_xm_std = np.nanstd(ent_xm_cldN, axis = 0)
ent_sed_std = np.nanstd(ent_sed_cldN, axis = 0)
ent_v_std = np.nanstd(ent_v_cldN, axis = 0)
sedt_std_h = np.nanstd(sedt_cldN_h, axis = 0)
ent_xm_std_h = np.nanstd(ent_xm_cldN_h, axis = 0)
ent_sed_std_h = np.nanstd(ent_sed_cldN_h, axis = 0)
ent_v_std_h = np.nanstd(ent_v_cldN_h, axis = 0)
xm_bl_z_std = np.nanstd(xm_bl_z_cldN, axis = 0)
xm23_bl_z_std = np.nanstd(xm23_bl_z_cldN, axis = 0)
xm_bl_zrat_std = np.nanstd(xm_bl_zrat_cldN, axis = 0)
xm23_bl_zrat_std = np.nanstd(xm23_bl_zrat_cldN, axis = 0)
RH_maxz_std = np.nanstd(RH_maxz_cldN, axis = 0)
xm23_maxz_std = np.nanstd(xm23_maxz_cldN, axis = 0)
h_xm_bl_z_std = np.nanstd(h_xm_bl_z_cldN, axis = 0)
h_xm23_bl_z_std = np.nanstd(h_xm23_bl_z_cldN, axis = 0)
h_xm_bl_zrat_std = np.nanstd(h_xm_bl_zrat_cldN, axis = 0)
h_xm23_bl_zrat_std = np.nanstd(h_xm23_bl_zrat_cldN, axis = 0)
h_RH_maxz_std = np.nanstd(h_RH_maxz_cldN, axis = 0)
h_xm23_maxz_std = np.nanstd(h_xm23_maxz_cldN, axis = 0)

f_sedt_std = np.nanstd(f_sedt_cldN, axis = 0)
f_ent_xm_std = np.nanstd(f_ent_xm_cldN, axis = 0)
f_ent_sed_std = np.nanstd(f_ent_sed_cldN, axis = 0)
f_ent_v_std = np.nanstd(f_ent_v_cldN, axis = 0)
f_sedt_std_h = np.nanstd(f_sedt_cldN_h, axis = 0)
f_ent_xm_std_h = np.nanstd(f_ent_xm_cldN_h, axis = 0)
f_ent_sed_std_h = np.nanstd(f_ent_sed_cldN_h, axis = 0)
f_ent_v_std_h = np.nanstd(f_ent_v_cldN_h, axis = 0)
f_xm_bl_z_std = np.nanstd(f_xm_bl_z_cldN, axis = 0)
f_xm23_bl_z_std = np.nanstd(f_xm23_bl_z_cldN, axis = 0)
f_xm_bl_zrat_std = np.nanstd(f_xm_bl_zrat_cldN, axis = 0)
f_xm23_bl_zrat_std = np.nanstd(f_xm23_bl_zrat_cldN, axis = 0)
f_RH_maxz_std = np.nanstd(f_RH_maxz_cldN, axis = 0)
f_xm23_maxz_std = np.nanstd(f_xm23_maxz_cldN, axis = 0)
f_h_xm_bl_z_std = np.nanstd(f_h_xm_bl_z_cldN, axis = 0)
f_h_xm23_bl_z_std = np.nanstd(f_h_xm23_bl_z_cldN, axis = 0)
f_h_xm_bl_zrat_std = np.nanstd(f_h_xm_bl_zrat_cldN, axis = 0)
f_h_xm23_bl_zrat_std = np.nanstd(f_h_xm23_bl_zrat_cldN, axis = 0)
f_h_RH_maxz_std = np.nanstd(f_h_RH_maxz_cldN, axis = 0)
f_h_xm23_maxz_std = np.nanstd(f_h_xm23_maxz_cldN, axis = 0)

sedt_dcount = np.count_nonzero(~np.isnan(sedt_cldN), axis=0)
ent_xm_dcount = np.count_nonzero(~np.isnan(ent_xm_cldN), axis=0)
ent_sed_dcount = np.count_nonzero(~np.isnan(ent_sed_cldN), axis=0)
ent_v_dcount = np.count_nonzero(~np.isnan(ent_v_cldN), axis=0)
sedt_dcount_h = np.count_nonzero(~np.isnan(sedt_cldN_h), axis=0)
ent_xm_dcount_h = np.count_nonzero(~np.isnan(ent_xm_cldN_h), axis=0)
ent_sed_dcount_h = np.count_nonzero(~np.isnan(ent_sed_cldN_h), axis=0)
ent_v_dcount_h = np.count_nonzero(~np.isnan(ent_v_cldN_h), axis=0)
xm_bl_z_dcount = np.count_nonzero(~np.isnan(xm_bl_z_cldN), axis=0)
xm23_bl_z_dcount = np.count_nonzero(~np.isnan(xm23_bl_z_cldN), axis=0)
xm_bl_zrat_dcount = np.count_nonzero(~np.isnan(xm_bl_zrat_cldN), axis=0)
xm23_bl_zrat_dcount = np.count_nonzero(~np.isnan(xm23_bl_zrat_cldN), axis=0)
RH_maxz_dcount = np.count_nonzero(~np.isnan(RH_maxz_cldN), axis=0)
xm23_maxz_dcount = np.count_nonzero(~np.isnan(xm23_maxz_cldN), axis=0)
h_xm_bl_z_dcount = np.count_nonzero(~np.isnan(h_xm_bl_z_cldN), axis=0)
h_xm23_bl_z_dcount = np.count_nonzero(~np.isnan(h_xm23_bl_z_cldN), axis=0)
h_xm_bl_zrat_dcount = np.count_nonzero(~np.isnan(h_xm_bl_zrat_cldN), axis=0)
h_xm23_bl_zrat_dcount = np.count_nonzero(~np.isnan(h_xm23_bl_zrat_cldN), axis=0)
h_RH_maxz_dcount = np.count_nonzero(~np.isnan(h_RH_maxz_cldN), axis=0)
h_xm23_maxz_dcount = np.count_nonzero(~np.isnan(h_xm23_maxz_cldN), axis=0)

f_sedt_dcount = np.count_nonzero(~np.isnan(f_sedt_cldN), axis=0)
f_ent_xm_dcount = np.count_nonzero(~np.isnan(f_ent_xm_cldN), axis=0)
f_ent_sed_dcount = np.count_nonzero(~np.isnan(f_ent_sed_cldN), axis=0)
f_ent_v_dcount = np.count_nonzero(~np.isnan(f_ent_v_cldN), axis=0)
f_sedt_dcount_h = np.count_nonzero(~np.isnan(f_sedt_cldN_h), axis=0)
f_ent_xm_dcount_h = np.count_nonzero(~np.isnan(f_ent_xm_cldN_h), axis=0)
f_ent_sed_dcount_h = np.count_nonzero(~np.isnan(f_ent_sed_cldN_h), axis=0)
f_ent_v_dcount_h = np.count_nonzero(~np.isnan(f_ent_v_cldN_h), axis=0)
f_xm_bl_z_dcount = np.count_nonzero(~np.isnan(f_xm_bl_z_cldN), axis=0)
f_xm23_bl_z_dcount = np.count_nonzero(~np.isnan(f_xm23_bl_z_cldN), axis=0)
f_xm_bl_zrat_dcount = np.count_nonzero(~np.isnan(f_xm_bl_zrat_cldN), axis=0)
f_xm23_bl_zrat_dcount = np.count_nonzero(~np.isnan(f_xm23_bl_zrat_cldN), axis=0)
f_RH_maxz_dcount = np.count_nonzero(~np.isnan(f_RH_maxz_cldN), axis=0)
f_xm23_maxz_dcount = np.count_nonzero(~np.isnan(f_xm23_maxz_cldN), axis=0)
f_h_xm_bl_z_dcount = np.count_nonzero(~np.isnan(f_h_xm_bl_z_cldN), axis=0)
f_h_xm23_bl_z_dcount = np.count_nonzero(~np.isnan(f_h_xm23_bl_z_cldN), axis=0)
f_h_xm_bl_zrat_dcount = np.count_nonzero(~np.isnan(f_h_xm_bl_zrat_cldN), axis=0)
f_h_xm23_bl_zrat_dcount = np.count_nonzero(~np.isnan(f_h_xm23_bl_zrat_cldN), axis=0)
f_h_RH_maxz_dcount = np.count_nonzero(~np.isnan(f_h_RH_maxz_cldN), axis=0)
f_h_xm23_maxz_dcount = np.count_nonzero(~np.isnan(f_h_xm23_maxz_cldN), axis=0)


wmean_dur = np.nansum(durs * bg_weights) / np.nansum(bg_weights)
wmean_h_dzblt = np.nansum(h_dzblt * bg_weights) / np.nansum(bg_weights)
wmean_f_h_dzblt = np.nansum(f_h_dzblt * bg_weights) / np.nansum(bg_weights)

wmean_f_dur = np.nansum(f_durs * f_durweights) / np.nansum(f_durweights)
wmean_t_dur = np.nansum(t_durs * t_durweights) / np.nansum(t_durweights)
wmean_nl_dur = np.nansum(nl_durs * nl_durweights) / np.nansum(nl_durweights)

wmean_sedt = np.nansum(sedts * bg_weights) / np.nansum(bg_weights)
wmean_ent_xm = np.nansum(ent_xms * bg_weights) / np.nansum(bg_weights)
wmean_ent_sed = np.nansum(ent_seds * bg_weights) / np.nansum(bg_weights)
wmean_ent_v = np.nansum(ent_vs * bg_weights) / np.nansum(bg_weights)
wmean_sedt_h = np.nansum(sedts_h * bg_weights) / np.nansum(bg_weights)
wmean_ent_xm_h = np.nansum(ent_xms_h * bg_weights) / np.nansum(bg_weights)
wmean_ent_sed_h = np.nansum(ent_seds_h * bg_weights) / np.nansum(bg_weights)
wmean_ent_v_h = np.nansum(ent_vs_h * bg_weights) / np.nansum(bg_weights)
wmean_xm_bl_z = np.nansum(xm_bl_zs * bg_weights) / np.nansum(bg_weights)
wmean_xm23_bl_z = np.nansum(xm23_bl_zs * bg_weights) / np.nansum(bg_weights)
wmean_xm_bl_zrat = np.nansum(xm_bl_zrats * bg_weights) / np.nansum(bg_weights)
wmean_xm23_bl_zrat = np.nansum(xm23_bl_zrats * bg_weights) / np.nansum(bg_weights)
wmean_RH_maxz = np.nansum(RH_maxzs * bg_weights) / np.nansum(bg_weights)
wmean_xm23_maxz = np.nansum(xm23_maxzs * bg_weights) / np.nansum(bg_weights)
h_wmean_xm_bl_z = np.nansum(h_xm_bl_zs * bg_weights) / np.nansum(bg_weights)
h_wmean_xm23_bl_z = np.nansum(h_xm23_bl_zs * bg_weights) / np.nansum(bg_weights)
h_wmean_xm_bl_zrat = np.nansum(h_xm_bl_zrats * bg_weights) / np.nansum(bg_weights)
h_wmean_xm23_bl_zrat = np.nansum(h_xm23_bl_zrats * bg_weights) / np.nansum(bg_weights)
h_wmean_RH_maxz = np.nansum(h_RH_maxzs * bg_weights) / np.nansum(bg_weights)
h_wmean_xm23_maxz = np.nansum(h_xm23_maxzs * bg_weights) / np.nansum(bg_weights)


f_wmean_sedt = np.nansum(f_sedts * f_durweights) / np.nansum(f_durweights)
f_wmean_ent_xm = np.nansum(f_ent_xms * f_durweights) / np.nansum(f_durweights)
f_wmean_ent_sed = np.nansum(f_ent_seds * f_durweights) / np.nansum(f_durweights)
f_wmean_ent_v = np.nansum(f_ent_vs * f_durweights) / np.nansum(f_durweights)
f_wmean_sedt_h = np.nansum(f_sedts_h * f_durweights) / np.nansum(f_durweights)
f_wmean_ent_xm_h = np.nansum(f_ent_xms_h * f_durweights) / np.nansum(f_durweights)
f_wmean_ent_sed_h = np.nansum(f_ent_seds_h * f_durweights) / np.nansum(f_durweights)
f_wmean_ent_v_h = np.nansum(f_ent_vs_h * f_durweights) / np.nansum(f_durweights)
f_wmean_xm_bl_z = np.nansum(f_xm_bl_zs * f_durweights) / np.nansum(f_durweights)
f_wmean_xm23_bl_z = np.nansum(f_xm23_bl_zs * f_durweights) / np.nansum(f_durweights)
f_wmean_xm_bl_zrat = np.nansum(f_xm_bl_zrats * f_durweights) / np.nansum(f_durweights)
f_wmean_xm23_bl_zrat = np.nansum(f_xm23_bl_zrats * f_durweights) / np.nansum(f_durweights)
f_wmean_RH_maxz = np.nansum(f_RH_maxzs * f_durweights) / np.nansum(f_durweights)
f_wmean_xm23_maxz = np.nansum(f_xm23_maxzs * f_durweights) / np.nansum(f_durweights)
f_h_wmean_xm_bl_z = np.nansum(f_h_xm_bl_zs * f_durweights) / np.nansum(f_durweights)
f_h_wmean_xm23_bl_z = np.nansum(f_h_xm23_bl_zs * f_durweights) / np.nansum(f_durweights)
f_h_wmean_xm_bl_zrat = np.nansum(f_h_xm_bl_zrats * f_durweights) / np.nansum(f_durweights)
f_h_wmean_xm23_bl_zrat = np.nansum(f_h_xm23_bl_zrats * f_durweights) / np.nansum(f_durweights)
f_h_wmean_RH_maxz = np.nansum(f_h_RH_maxzs * f_durweights) / np.nansum(f_durweights)
f_h_wmean_xm23_maxz = np.nansum(f_h_xm23_maxzs * f_durweights) / np.nansum(f_durweights)

con_95 = 2 * durstd / np.sqrt(dcount - 2)
h_dzblt_95con = 2 * h_dzblt_std / np.sqrt(h_dzblt_dcount - 2)
f_h_dzblt_95con = 2 * f_h_dzblt_std / np.sqrt(f_h_dzblt_dcount - 2)

f_95con = 2 * f_durstd / np.sqrt(f_dcount - 2)
t_95con = 2 * t_durstd / np.sqrt(t_dcount - 2)
nl_95con = 2 * nl_durstd / np.sqrt(nl_dcount - 2)

sedt_95con = 2 * sedt_std / np.sqrt(sedt_dcount - 2)
ent_xm_95con = 2 * ent_xm_std / np.sqrt(ent_xm_dcount - 2)
ent_sed_95con = 2 * ent_sed_std / np.sqrt(ent_sed_dcount - 2)
ent_v_95con = 2 * ent_v_std / np.sqrt(ent_v_dcount - 2)
sedt_95con_h = 2 * sedt_std_h / np.sqrt(sedt_dcount_h - 2)
ent_xm_95con_h = 2 * ent_xm_std_h / np.sqrt(ent_xm_dcount_h - 2)
ent_sed_95con_h = 2 * ent_sed_std_h / np.sqrt(ent_sed_dcount_h - 2)
ent_v_95con_h = 2 * ent_v_std_h / np.sqrt(ent_v_dcount_h - 2)
xm_bl_z_95con = 2 * xm_bl_z_std / np.sqrt(xm_bl_z_dcount - 2)
xm23_bl_z_95con = 2 * xm23_bl_z_std / np.sqrt(xm23_bl_z_dcount - 2)
xm_bl_zrat_95con = 2 * xm_bl_zrat_std / np.sqrt(xm_bl_zrat_dcount - 2)
xm23_bl_zrat_95con = 2 * xm23_bl_zrat_std / np.sqrt(xm23_bl_zrat_dcount - 2)
RH_maxz_95con = 2 * RH_maxz_std / np.sqrt(RH_maxz_dcount - 2)
xm23_maxz_95con = 2 * xm23_maxz_std / np.sqrt(xm23_maxz_dcount - 2)
h_xm_bl_z_95con = 2 * h_xm_bl_z_std / np.sqrt(h_xm_bl_z_dcount - 2)
h_xm23_bl_z_95con = 2 * h_xm23_bl_z_std / np.sqrt(h_xm23_bl_z_dcount - 2)
h_xm_bl_zrat_95con = 2 * h_xm_bl_zrat_std / np.sqrt(h_xm_bl_zrat_dcount - 2)
h_xm23_bl_zrat_95con = 2 * h_xm23_bl_zrat_std / np.sqrt(h_xm23_bl_zrat_dcount - 2)
h_RH_maxz_95con = 2 * h_RH_maxz_std / np.sqrt(h_RH_maxz_dcount - 2)
h_xm23_maxz_95con = 2 * h_xm23_maxz_std / np.sqrt(h_xm23_maxz_dcount - 2)

f_sedt_95con = 2 * f_sedt_std / np.sqrt(f_sedt_dcount - 2)
f_ent_xm_95con = 2 * f_ent_xm_std / np.sqrt(f_ent_xm_dcount - 2)
f_ent_sed_95con = 2 * f_ent_sed_std / np.sqrt(f_ent_sed_dcount - 2)
f_ent_v_95con = 2 * f_ent_v_std / np.sqrt(f_ent_v_dcount - 2)
f_sedt_95con_h = 2 * f_sedt_std_h / np.sqrt(f_sedt_dcount_h - 2)
f_ent_xm_95con_h = 2 * f_ent_xm_std_h / np.sqrt(f_ent_xm_dcount_h - 2)
f_ent_sed_95con_h = 2 * f_ent_sed_std_h / np.sqrt(f_ent_sed_dcount_h - 2)
f_ent_v_95con_h = 2 * f_ent_v_std_h / np.sqrt(f_ent_v_dcount_h - 2)
f_xm_bl_z_95con = 2 * f_xm_bl_z_std / np.sqrt(f_xm_bl_z_dcount - 2)
f_xm23_bl_z_95con = 2 * f_xm23_bl_z_std / np.sqrt(f_xm23_bl_z_dcount - 2)
f_xm_bl_zrat_95con = 2 * f_xm_bl_zrat_std / np.sqrt(f_xm_bl_zrat_dcount - 2)
f_xm23_bl_zrat_95con = 2 * f_xm23_bl_zrat_std / np.sqrt(f_xm23_bl_zrat_dcount - 2)
f_RH_maxz_95con = 2 * f_RH_maxz_std / np.sqrt(f_RH_maxz_dcount - 2)
f_xm23_maxz_95con = 2 * f_xm23_maxz_std / np.sqrt(f_xm23_maxz_dcount - 2)
f_h_xm_bl_z_95con = 2 * f_h_xm_bl_z_std / np.sqrt(f_h_xm_bl_z_dcount - 2)
f_h_xm23_bl_z_95con = 2 * f_h_xm23_bl_z_std / np.sqrt(f_h_xm23_bl_z_dcount - 2)
f_h_xm_bl_zrat_95con = 2 * f_h_xm_bl_zrat_std / np.sqrt(f_h_xm_bl_zrat_dcount - 2)
f_h_xm23_bl_zrat_95con = 2 * f_h_xm23_bl_zrat_std / np.sqrt(f_h_xm23_bl_zrat_dcount - 2)
f_h_RH_maxz_95con = 2 * f_h_RH_maxz_std / np.sqrt(f_h_RH_maxz_dcount - 2)
f_h_xm23_maxz_95con = 2 * f_h_xm23_maxz_std / np.sqrt(f_h_xm23_maxz_dcount - 2)


ent_factor = 86400
ent_factor_h = 86400 / 2


# In[14]:


fig, (ax1, ax2, ax3) = plt.subplots(1,3)
fig.set_size_inches(20,6)
plt.subplots_adjust(left=0.075, right=0.95, bottom=0.1, top=0.85, wspace=0.2, hspace=0.2)

ax1.plot(x, ent_xm_change_h * ent_factor_h + wmean_ent_xm_h * ent_factor_h, color = 'black', label = 'All Simulations')
ax1.plot(x, f_ent_xm_change_h * ent_factor_h + f_wmean_ent_xm_h * ent_factor_h, color = 'black', linestyle = '--', label = 'Foggy Simulations')
ax1.fill_between(cldN_list, 
                ent_xm_change_h * ent_factor_h + wmean_ent_xm_h * ent_factor_h - ent_xm_95con_h * ent_factor_h, 
                ent_xm_change_h * ent_factor_h + wmean_ent_xm_h * ent_factor_h + ent_xm_95con_h * ent_factor_h, 
                color='black', alpha=0.2)
ax1.fill_between(cldN_list, 
                f_ent_xm_change_h * ent_factor_h + f_wmean_ent_xm_h * ent_factor_h - f_ent_xm_95con_h * ent_factor_h, 
                f_ent_xm_change_h * ent_factor_h + f_wmean_ent_xm_h * ent_factor_h + f_ent_xm_95con_h * ent_factor_h, 
                color='black', alpha=0.1)

ax1.set_ylabel('Moisture Flux (kg/m^2)', fontsize = 20)
ax1.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 18)
ax1.set_title('Entrainment', fontsize = 22)
ax1.tick_params(axis='both', which='both', labelsize=16)
ax1.legend(fontsize = 18)


ax2.plot(x, -sedt_change_h - wmean_sedt_h, color = 'black', label = 'All Simulations')
ax2.plot(x, -f_sedt_change_h - f_wmean_sedt_h, color = 'black', linestyle = '--', label = 'Foggy Simulations')
ax2.fill_between(cldN_list, 
                -sedt_change_h - wmean_sedt_h - sedt_95con_h, 
                -sedt_change_h - wmean_sedt_h + sedt_95con_h, 
                color='black', alpha=0.2)
ax2.fill_between(cldN_list, 
                -f_sedt_change_h - f_wmean_sedt_h - f_sedt_95con_h, 
                -f_sedt_change_h - f_wmean_sedt_h + f_sedt_95con_h, 
                color='black', alpha=0.1)


ax2.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 18)
ax2.set_title('Surface Deposition', fontsize = 22)
ax2.tick_params(axis='both', which='both', labelsize=16)


ax3.plot(x, ent_sed_change_h + wmean_ent_sed_h, color = 'black', label = 'All Simulations')
ax3.plot(x, f_ent_sed_change_h + f_wmean_ent_sed_h, color = 'black', linestyle = '--', label = 'Foggy Simulations')
ax3.fill_between(cldN_list, 
                ent_sed_change_h + wmean_ent_sed_h - ent_sed_95con_h, 
                ent_sed_change_h + wmean_ent_sed_h + ent_sed_95con_h, 
                color='black', alpha=0.2)
ax3.fill_between(cldN_list, 
                f_ent_sed_change_h + f_wmean_ent_sed_h - f_ent_sed_95con_h, 
                f_ent_sed_change_h + f_wmean_ent_sed_h + f_ent_sed_95con_h, 
                color='black', alpha=0.1)


ax3.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 18)
ax3.set_title('Combined Moisture Flux', fontsize = 22)
ax3.tick_params(axis='both', which='both', labelsize=16)

fig.suptitle('12-Hour Surface Deposition and Entrainment Moisture Flux', fontsize = 24)

fname = desc + 'ent_sed_comb_h.png'
plt.savefig(fname)
plt.show()
plt.close()




fig, ax = plt.subplots()
fig.set_size_inches(6.5,4.5)
plt.subplots_adjust(left=0.175, right=0.925, bottom=0.135, top=0.95)

ax.plot(x, ent_sed_change_h, color = 'black', label = 'Combined, All Simulations')
ax.plot(x, f_ent_sed_change_h, color = 'black', linestyle = '--', label = 'Foggy Simulations')
ax.fill_between(cldN_list, 
                ent_sed_change_h - ent_sed_95con_h, 
                ent_sed_change_h + ent_sed_95con_h, 
                color='black', alpha=0.2)
ax.fill_between(cldN_list, 
                f_ent_sed_change_h - f_ent_sed_95con_h, 
                f_ent_sed_change_h + f_ent_sed_95con_h, 
                color='black', alpha=0.1)

ax.plot(x, - sedt_change_h, color = 'blue', label = 'Surface Deposition')
ax.plot(x, - f_sedt_change_h, color = 'blue', linestyle = '--')
ax.fill_between(cldN_list, 
                - sedt_change_h - sedt_95con_h, 
                - sedt_change_h + sedt_95con_h, 
                color='blue', alpha=0.2)
ax.fill_between(cldN_list, 
                - f_sedt_change_h - f_sedt_95con_h, 
                - f_sedt_change_h + f_sedt_95con_h, 
                color='blue', alpha=0.1)

ax.plot(x, ent_xm_change_h * ent_factor_h, color = 'red', label = 'Entrainment')
ax.plot(x, f_ent_xm_change_h * ent_factor_h, color = 'red', linestyle = '--')
ax.fill_between(cldN_list, 
                ent_xm_change_h * ent_factor_h - ent_xm_95con_h * ent_factor_h, 
                ent_xm_change_h * ent_factor_h + ent_xm_95con_h * ent_factor_h, 
                color='red', alpha=0.2)
ax.fill_between(cldN_list, 
                f_ent_xm_change_h * ent_factor_h - f_ent_xm_95con_h * ent_factor_h, 
                f_ent_xm_change_h * ent_factor_h + f_ent_xm_95con_h * ent_factor_h, 
                color='red', alpha=0.1)

ax.set_ylabel('Moisture Flux Delta (kg/m$^2$)', fontsize = 18)
ax.set_xlabel('Aerosol Number Concentration (#/cm$^3$)', fontsize = 18)
# ax.set_title('Moisture Flux Difference Vs. Na, First Half', fontsize = 20)
ax.tick_params(axis='both', which='both', labelsize=14)
ax.legend(fontsize = 16)

fname = desc + 'fig06_old.png'
plt.savefig(fname)
plt.show()
plt.close()


# In[15]:


th_max = 15
th_min = 5

df_gb_cldN = df_this_test.groupby(['cldN']).mean().reset_index()
df_gb_ug = df_this_test.groupby(['Ug']).mean().reset_index()
df_gb_wmax = df_this_test.groupby(['wmax']).mean().reset_index()
df_gb_dtsurf = df_this_test.groupby(['dtsurf']).mean().reset_index()

gb_cldN = (df_gb_cldN.index.to_numpy() - min(df_gb_cldN.index.to_numpy())) / (max(df_gb_cldN.index.to_numpy()) - min(df_gb_cldN.index.to_numpy()))
gb_ug = (df_gb_ug.index.to_numpy() - min(df_gb_ug.index.to_numpy())) / (max(df_gb_ug.index.to_numpy()) - min(df_gb_ug.index.to_numpy()))
gb_wmax = - (df_gb_wmax.index.to_numpy() - max(df_gb_wmax.index.to_numpy())) / (max(df_gb_wmax.index.to_numpy()) - min(df_gb_wmax.index.to_numpy()))
gb_dtsurf = (df_gb_dtsurf.index.to_numpy() - min(df_gb_dtsurf.index.to_numpy())) / (max(df_gb_dtsurf.index.to_numpy()) - min(df_gb_dtsurf.index.to_numpy()))
    
wmax_wmax = (gb_wmax * (th_max-th_min)) + th_min
ug_ug = (gb_ug * (th_max-th_min)) + th_min
dtsurf_dtsurf = (gb_dtsurf * (th_max-th_min)) + th_min
cldN_cldN = (gb_cldN * (th_max-th_min)) + th_min

factor = max([max(df_gb_ug['duration']), max(df_gb_wmax['duration']), max(df_gb_dtsurf['duration']), max(df_gb_cldN['duration'])])

fig, ax = plt.subplots()
fig.set_size_inches(6,4)
plt.subplots_adjust(left=0.15, right=0.925, bottom=0.135, top=0.925)

ax.scatter(df_gb_ug['duration'], df_gb_ug['max_sfc_xm2'] * 1000, color = 'green', s = ug_ug ** 2)
ax.plot(df_gb_ug['duration'], df_gb_ug['max_sfc_xm2'] * 1000, color = 'green', label = 'Geostrophic Wind')
ax.scatter(df_gb_wmax['duration'], df_gb_wmax['max_sfc_xm2'] * 1000, color = 'blue', s = wmax_wmax ** 2)
ax.plot(df_gb_wmax['duration'], df_gb_wmax['max_sfc_xm2'] * 1000, color = 'blue', label = 'Subsidence')
ax.scatter(df_gb_dtsurf['duration'], df_gb_dtsurf['max_sfc_xm2'] * 1000, color = 'red', s = dtsurf_dtsurf ** 2)
ax.plot(df_gb_dtsurf['duration'], df_gb_dtsurf['max_sfc_xm2'] * 1000, color = 'red', label = 'SFC Temp. Change')
ax.scatter(df_gb_cldN['duration'], df_gb_cldN['max_sfc_xm2'] * 1000, color = 'black', s = cldN_cldN ** 2)
ax.plot(df_gb_cldN['duration'], df_gb_cldN['max_sfc_xm2'] * 1000, color = 'black', label = 'Na')


ax.set_ylabel('Max Surface CWC (g/kg)', fontsize = 16)
ax.set_xlabel('Mean Fog Duration (Hours)', fontsize = 16)
# ax.set_title('Fog Duration and Max Sfc. CWC for Diff. Inputs', fontsize = 16)
ax.tick_params(axis='both', which='both', labelsize=12)
ax.legend(fontsize = 12)

fname = desc + 'fig04.png'
plt.savefig(fname)

plt.show()
plt.close()


# In[16]:


import warnings; warnings.simplefilter('ignore')

df_ExpDur = df_all[['Ug', 'wmax', 'dtsurf', 'cldN', 'duration', 't_dur', 'nl_dur']]

cldN_500 = df_all['cldN'] == 500.

df_cldN500 = df_all[cldN_500]
df_l_Ug_cldN500 = df_all[cldN_500 & low_Ug]
df_h_Ug_cldN500 = df_all[cldN_500 & high_Ug]
df_l_wmax_cldN500 = df_all[cldN_500 & low_wmax]
df_h_wmax_cldN500 = df_all[cldN_500 & high_wmax]
df_l_dtsurf_cldN500 = df_all[cldN_500 & low_dtsurf]
df_h_dtsurf_cldN500 = df_all[cldN_500 & high_dtsurf]

df_ED_vars = df_ExpDur[cond_outlier]
df_ED_h_Ug = df_ExpDur[cond_outlier & high_Ug]
df_ED_l_Ug = df_ExpDur[cond_outlier & low_Ug]
df_ED_h_wmax = df_ExpDur[cond_outlier & high_wmax]
df_ED_l_wmax = df_ExpDur[cond_outlier & low_wmax]
df_ED_h_dtsurf = df_ExpDur[cond_outlier & high_dtsurf]
df_ED_l_dtsurf = df_ExpDur[cond_outlier & low_dtsurf]

df_ED_bg = df_ED_vars.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_h_Ug = df_ED_h_Ug.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_l_Ug = df_ED_l_Ug.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_h_wmax = df_ED_h_wmax.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_l_wmax = df_ED_l_wmax.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_h_dtsurf = df_ED_h_dtsurf.groupby(['Ug', 'wmax', 'dtsurf']).mean()
df_ED_bg_l_dtsurf = df_ED_l_dtsurf.groupby(['Ug', 'wmax', 'dtsurf']).mean()

f_durs = np.empty(np.shape(df_ED_bg)[0])
f_durs[:] = np.nan
f_durweights = np.empty(np.shape(df_ED_bg)[0])
f_durs_h_Ug = np.empty(np.shape(df_ED_bg_h_Ug)[0])
f_durs_h_Ug[:] = np.nan
f_durweights_h_Ug = np.empty(np.shape(df_ED_bg_h_Ug)[0])
f_durs_l_Ug = np.empty(np.shape(df_ED_bg_l_Ug)[0])
f_durs_l_Ug[:] = np.nan
f_durweights_l_Ug = np.empty(np.shape(df_ED_bg_l_Ug)[0])
f_durs_h_wmax = np.empty(np.shape(df_ED_bg_h_wmax)[0])
f_durs_h_wmax[:] = np.nan
f_durweights_h_wmax = np.empty(np.shape(df_ED_bg_h_wmax)[0])
f_durs_l_wmax = np.empty(np.shape(df_ED_bg_l_wmax)[0])
f_durs_l_wmax[:] = np.nan
f_durweights_l_wmax = np.empty(np.shape(df_ED_bg_l_wmax)[0])
f_durs_h_dtsurf = np.empty(np.shape(df_ED_bg_h_dtsurf)[0])
f_durs_h_dtsurf[:] = np.nan
f_durweights_h_dtsurf = np.empty(np.shape(df_ED_bg_h_dtsurf)[0])
f_durs_l_dtsurf = np.empty(np.shape(df_ED_bg_l_dtsurf)[0])
f_durs_l_dtsurf[:] = np.nan
f_durweights_l_dtsurf = np.empty(np.shape(df_ED_bg_l_dtsurf)[0])

f_durs_cldN = np.empty((np.shape(df_ED_bg)[0], 10))
f_durs_cldN[:] = np.nan
f_durs_cldN_h_Ug = np.empty((np.shape(df_ED_bg_h_Ug)[0], 10))
f_durs_cldN_h_Ug[:] = np.nan
f_durs_cldN_l_Ug = np.empty((np.shape(df_ED_bg_l_Ug)[0], 10))
f_durs_cldN_l_Ug[:] = np.nan
f_durs_cldN_h_wmax = np.empty((np.shape(df_ED_bg_h_wmax)[0], 10))
f_durs_cldN_h_wmax[:] = np.nan
f_durs_cldN_l_wmax = np.empty((np.shape(df_ED_bg_l_wmax)[0], 10))
f_durs_cldN_l_wmax[:] = np.nan
f_durs_cldN_h_dtsurf = np.empty((np.shape(df_ED_bg_h_dtsurf)[0], 10))
f_durs_cldN_h_dtsurf[:] = np.nan
f_durs_cldN_l_dtsurf = np.empty((np.shape(df_ED_bg_l_dtsurf)[0], 10))
f_durs_cldN_l_dtsurf[:] = np.nan

t_durs = np.empty(np.shape(df_ED_bg)[0])
t_durweights = np.empty(np.shape(df_ED_bg)[0])

t_durs_cldN = np.empty((np.shape(df_ED_bg)[0], 10))
t_durs_cldN[:] = np.nan

nl_durs = np.empty(np.shape(df_ED_bg)[0])
nl_durweights = np.empty(np.shape(df_ED_bg)[0])

nl_durs_cldN = np.empty((np.shape(df_ED_bg)[0], 10))
nl_durs_cldN[:] = np.nan

for i in range(np.shape(df_ED_bg)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg.index[i]
    cond_Ug = df_ED_vars['Ug'] == Ug
    cond_wmax = df_ED_vars['wmax'] == wmax
    cond_dtsurf = df_ED_vars['dtsurf'] == dtsurf
    cond_fog = df_ED_vars['duration'] > 0.0
    cond_t = df_ED_vars['t_dur'] > 0.0
    cond_nl = df_ED_vars['t_dur'] > 0.0

    df_this_bg = df_ED_vars[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_vars[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    df_this_bg_t_fog = df_ED_vars[cond_Ug & cond_wmax & cond_dtsurf & cond_t].reset_index()
    df_this_bg_nl_fog = df_ED_vars[cond_Ug & cond_wmax & cond_dtsurf & cond_nl].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    this_t_fog_cldN = df_this_bg_t_fog['cldN'].to_numpy()
    this_nl_fog_cldN = df_this_bg_nl_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        
        cond500 = df_this_bg['cldN'] == 500
        f_df_cldN = df_this_bg[cond500]
        
        f_durs[i] = f_df_cldN['duration']
        f_durweights[i] = np.shape(df_this_bg)[0]
        
        for j in range(10):
            if cldN_list[j] in df_this_bg['cldN'].to_numpy():
                index = np.where(df_this_bg['cldN'] == cldN_list[j])[0][0]
                # print(index, df_this_bg['duration'][index])
                if df_this_bg['duration'][index].all() > 0:
                    f_durs_cldN[i,j] = df_this_bg['duration'][index] - f_durs[i]
                else:
                    f_durs_cldN[i,j] = np.nan
            else:
                f_durs_cldN[i,j] = np.nan
                
    if np.shape(this_t_fog_cldN)[0] > 0:
        
        cond500 = df_this_bg['cldN'] == 500
        t_df_cldN = df_this_bg[cond500]
        
        t_durs[i] = t_df_cldN['t_dur']
        t_durweights[i] = np.shape(df_this_bg)[0]
        
        for j in range(10):
            if cldN_list[j] in df_this_bg['cldN'].to_numpy():
                index = np.where(df_this_bg['cldN'] == cldN_list[j])[0][0]
                # print(index, df_this_bg['t_dur'][index])
                if df_this_bg['t_dur'][index].all() > 0:
                    t_durs_cldN[i,j] = df_this_bg['t_dur'][index] - t_durs[i]
                else:
                    t_durs_cldN[i,j] = np.nan
            else:
                t_durs_cldN[i,j] = np.nan
                
    if np.shape(this_nl_fog_cldN)[0] > 0:
        cond500 = df_this_bg['cldN'] == 500
        nl_df_cldN = df_this_bg[cond500]
        
        nl_durs[i] = nl_df_cldN['nl_dur']
        nl_durweights[i] = np.shape(df_this_bg)[0]
        
        for j in range(10):
            if cldN_list[j] in df_this_bg['cldN'].to_numpy():
                index = np.where(df_this_bg['cldN'] == cldN_list[j])[0]
                # print(index, df_this_bg['nl_dur'][index], df_this_bg['cldN'].to_numpy())
                nl_durs_cldN[i,j] = df_this_bg['nl_dur'][index] - nl_durs[i]
            else:
                nl_durs_cldN[i,j] = np.nan


for i in range(np.shape(df_ED_bg_h_Ug)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_h_Ug.index[i]
    cond_Ug = df_ED_h_Ug['Ug'] == Ug
    cond_wmax = df_ED_h_Ug['wmax'] == wmax
    cond_dtsurf = df_ED_h_Ug['dtsurf'] == dtsurf
    cond_fog = df_ED_h_Ug['duration'] > 0.0

    df_this_bg = df_ED_h_Ug[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_h_Ug[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        
        cond500 = df_this_bg['cldN'] == 500
        f_df_cldN = df_this_bg[cond500]
        
        f_durs_h_Ug[i] = f_df_cldN['duration']
        f_durweights_h_Ug[i] = np.shape(df_this_bg)[0]
        
        for j in range(10):
            if cldN_list[j] in df_this_bg['cldN'].to_numpy():
                index = np.where(df_this_bg['cldN'] == cldN_list[j])[0]
                if df_this_bg['duration'][index].all() > 0:
                    f_durs_cldN_h_Ug[i,j] = df_this_bg['duration'][index] - f_durs_h_Ug[i]
                else:
                    f_durs_cldN_h_Ug[i,j] = np.nan
            else:
                f_durs_cldN_h_Ug[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_l_Ug)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_l_Ug.index[i]
    cond_Ug = df_ED_l_Ug['Ug'] == Ug
    cond_wmax = df_ED_l_Ug['wmax'] == wmax
    cond_dtsurf = df_ED_l_Ug['dtsurf'] == dtsurf
    cond_fog = df_ED_l_Ug['duration'] > 0.0

    df_this_bg = df_ED_l_Ug[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_l_Ug[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        
        cond500 = df_this_bg['cldN'] == 500
        f_df_cldN = df_this_bg[cond500]
        
        f_durs_l_Ug[i] = f_df_cldN['duration']
        f_durweights_l_Ug[i] = np.shape(df_this_bg)[0]
        
        for j in range(10):
            if cldN_list[j] in df_this_bg['cldN'].to_numpy():
                index = np.where(df_this_bg['cldN'] == cldN_list[j])[0]
                if df_this_bg['duration'][index].all() > 0:
                    f_durs_cldN_l_Ug[i,j] = df_this_bg['duration'][index] - f_durs_l_Ug[i]
                else:
                    f_durs_cldN_l_Ug[i,j] = np.nan
            else:
                f_durs_cldN_l_Ug[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_h_wmax)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_h_wmax.index[i]
    cond_Ug = df_ED_h_wmax['Ug'] == Ug
    cond_wmax = df_ED_h_wmax['wmax'] == wmax
    cond_dtsurf = df_ED_h_wmax['dtsurf'] == dtsurf
    cond_fog = df_ED_h_wmax['duration'] > 0.0

    df_this_bg = df_ED_h_wmax[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_h_wmax[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        
        cond500 = df_this_bg['cldN'] == 500
        f_df_cldN = df_this_bg[cond500]
        
        f_durs_h_wmax[i] = f_df_cldN['duration']
        f_durweights_h_wmax[i] = np.shape(df_this_bg)[0]
        
        for j in range(10):
            if cldN_list[j] in df_this_bg['cldN'].to_numpy():
                index = np.where(df_this_bg['cldN'] == cldN_list[j])[0]
                if df_this_bg['duration'][index].all() > 0:
                    f_durs_cldN_h_wmax[i,j] = df_this_bg['duration'][index] - f_durs_h_wmax[i]
                else:
                    f_durs_cldN_h_wmax[i,j] = np.nan
            else:
                f_durs_cldN_h_wmax[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_l_wmax)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_l_wmax.index[i]
    cond_Ug = df_ED_l_wmax['Ug'] == Ug
    cond_wmax = df_ED_l_wmax['wmax'] == wmax
    cond_dtsurf = df_ED_l_wmax['dtsurf'] == dtsurf
    cond_fog = df_ED_l_wmax['duration'] > 0.0

    df_this_bg = df_ED_l_wmax[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_l_wmax[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        
        cond500 = df_this_bg['cldN'] == 500
        f_df_cldN = df_this_bg[cond500]
        
        f_durs_l_wmax[i] = f_df_cldN['duration']
        f_durweights_l_wmax[i] = np.shape(df_this_bg)[0]
        
        for j in range(10):
            if cldN_list[j] in df_this_bg['cldN'].to_numpy():
                index = np.where(df_this_bg['cldN'] == cldN_list[j])[0]
                if df_this_bg['duration'][index].all() > 0:
                    f_durs_cldN_l_wmax[i,j] = df_this_bg['duration'][index] - f_durs_l_wmax[i]
                else:
                    f_durs_cldN_l_wmax[i,j] = np.nan
            else:
                f_durs_cldN_l_wmax[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_h_dtsurf)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_h_dtsurf.index[i]
    cond_Ug = df_ED_h_dtsurf['Ug'] == Ug
    cond_wmax = df_ED_h_dtsurf['wmax'] == wmax
    cond_dtsurf = df_ED_h_dtsurf['dtsurf'] == dtsurf
    cond_fog = df_ED_h_dtsurf['duration'] > 0.0

    df_this_bg = df_ED_h_dtsurf[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_h_dtsurf[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        
        cond500 = df_this_bg['cldN'] == 500
        f_df_cldN = df_this_bg[cond500]
        
        f_durs_h_dtsurf[i] = f_df_cldN['duration']
        f_durweights_h_dtsurf[i] = np.shape(df_this_bg)[0]
        
        for j in range(10):
            if cldN_list[j] in df_this_bg['cldN'].to_numpy():
                index = np.where(df_this_bg['cldN'] == cldN_list[j])[0]
                if df_this_bg['duration'][index].all() > 0:
                    f_durs_cldN_h_dtsurf[i,j] = df_this_bg['duration'][index] - f_durs_h_dtsurf[i]
                else:
                    f_durs_cldN_h_dtsurf[i,j] = np.nan
            else:
                f_durs_cldN_h_dtsurf[i,j] = np.nan
                
for i in range(np.shape(df_ED_bg_l_dtsurf)[0]):
    (Ug, wmax, dtsurf) =  df_ED_bg_l_dtsurf.index[i]
    cond_Ug = df_ED_l_dtsurf['Ug'] == Ug
    cond_wmax = df_ED_l_dtsurf['wmax'] == wmax
    cond_dtsurf = df_ED_l_dtsurf['dtsurf'] == dtsurf
    cond_fog = df_ED_l_dtsurf['duration'] > 0.0

    df_this_bg = df_ED_l_dtsurf[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_this_bg_fog = df_ED_l_dtsurf[cond_Ug & cond_wmax & cond_dtsurf & cond_fog].reset_index()
    
    this_bg_cldN = df_this_bg['cldN'].to_numpy()
    this_fog_cldN = df_this_bg_fog['cldN'].to_numpy()
    
    bg_maxcldN = np.nanmax(this_bg_cldN)
    bg_mincldN = np.nanmin(this_bg_cldN)
    
    if np.shape(this_fog_cldN)[0] > 0:
        
        cond500 = df_this_bg['cldN'] == 500
        f_df_cldN = df_this_bg[cond500]
        
        f_durs_l_dtsurf[i] = f_df_cldN['duration']
        f_durweights_l_dtsurf[i] = np.shape(df_this_bg)[0]
        
        for j in range(10):
            if cldN_list[j] in df_this_bg['cldN'].to_numpy():
                index = np.where(df_this_bg['cldN'] == cldN_list[j])[0]
                if df_this_bg['duration'][index].all() > 0:
                    f_durs_cldN_l_dtsurf[i,j] = df_this_bg['duration'][index] - f_durs_l_dtsurf[i]
                else:
                    f_durs_cldN_l_dtsurf[i,j] = np.nan
            else:
                f_durs_cldN_l_dtsurf[i,j] = np.nan


                    
f_durchange = np.nanmean(f_durs_cldN, axis = 0)
t_durchange = np.nanmean(t_durs_cldN, axis = 0)
nl_durchange = np.nanmean(nl_durs_cldN, axis = 0)
f_durstd = np.nanstd(f_durs_cldN, axis = 0)
t_durstd = np.nanstd(t_durs_cldN, axis = 0)
nl_durstd = np.nanstd(nl_durs_cldN, axis = 0)
f_dcount = np.count_nonzero(~np.isnan(f_durs_cldN), axis=0)
t_dcount = np.count_nonzero(~np.isnan(t_durs_cldN), axis=0)
nl_dcount = np.count_nonzero(~np.isnan(nl_durs_cldN), axis=0)

f_durchange_h_Ug = np.nanmean(f_durs_cldN_h_Ug, axis = 0)
f_durchange_l_Ug = np.nanmean(f_durs_cldN_l_Ug, axis = 0)
f_durchange_h_wmax = np.nanmean(f_durs_cldN_h_wmax, axis = 0)
f_durchange_l_wmax = np.nanmean(f_durs_cldN_l_wmax, axis = 0)
f_durchange_h_dtsurf = np.nanmean(f_durs_cldN_h_dtsurf, axis = 0)
f_durchange_l_dtsurf = np.nanmean(f_durs_cldN_l_dtsurf, axis = 0)

f_durstd_h_Ug = np.nanstd(f_durs_cldN_h_Ug, axis = 0)
f_durstd_l_Ug = np.nanstd(f_durs_cldN_l_Ug, axis = 0)
f_durstd_h_wmax = np.nanstd(f_durs_cldN_h_wmax, axis = 0)
f_durstd_l_wmax = np.nanstd(f_durs_cldN_l_wmax, axis = 0)
f_durstd_h_dtsurf = np.nanstd(f_durs_cldN_h_dtsurf, axis = 0)
f_durstd_l_dtsurf = np.nanstd(f_durs_cldN_l_dtsurf, axis = 0)

f_dcount_h_Ug = np.count_nonzero(~np.isnan(f_durs_cldN_h_Ug), axis=0)
f_dcount_l_Ug = np.count_nonzero(~np.isnan(f_durs_cldN_l_Ug), axis=0)
f_dcount_h_wmax = np.count_nonzero(~np.isnan(f_durs_cldN_h_wmax), axis=0)
f_dcount_l_wmax = np.count_nonzero(~np.isnan(f_durs_cldN_l_wmax), axis=0)
f_dcount_h_dtsurf = np.count_nonzero(~np.isnan(f_durs_cldN_h_dtsurf), axis=0)
f_dcount_l_dtsurf = np.count_nonzero(~np.isnan(f_durs_cldN_l_dtsurf), axis=0)

# wmean_f_dur = np.nansum(f_durs * f_durweights) / np.nansum(f_durweights)
# wmean_t_dur = np.nansum(t_durs * t_durweights) / np.nansum(t_durweights)
# wmean_nl_dur = np.nansum(nl_durs * nl_durweights) / np.nansum(nl_durweights)

wmean_f_dur = np.nanmean(df_cldN500['duration'])
wmean_t_dur = np.nanmean(t_durs)
wmean_nl_dur = np.nanmean(nl_durs)

# wmean_f_dur_h_Ug = np.nansum(f_durs_h_Ug * f_durweights_h_Ug) / np.nansum(f_durweights_h_Ug)
# wmean_f_dur_l_Ug = np.nansum(f_durs_l_Ug * f_durweights_l_Ug) / np.nansum(f_durweights_l_Ug)
# wmean_f_dur_h_wmax = np.nansum(f_durs_h_wmax * f_durweights_h_wmax) / np.nansum(f_durweights_h_wmax)
# wmean_f_dur_l_wmax = np.nansum(f_durs_l_wmax * f_durweights_l_wmax) / np.nansum(f_durweights_l_wmax)
# wmean_f_dur_h_dtsurf = np.nansum(f_durs_h_dtsurf * f_durweights_h_dtsurf) / np.nansum(f_durweights_h_dtsurf)
# wmean_f_dur_l_dtsurf = np.nansum(f_durs_l_dtsurf * f_durweights_l_dtsurf) / np.nansum(f_durweights_l_dtsurf)

wmean_f_dur_h_Ug = np.nanmean(f_durs_h_Ug)
wmean_f_dur_l_Ug = np.nanmean(f_durs_l_Ug)
wmean_f_dur_h_wmax = np.nanmean(f_durs_h_wmax)
wmean_f_dur_l_wmax = np.nanmean(f_durs_l_wmax)
wmean_f_dur_h_dtsurf = np.nanmean(f_durs_h_dtsurf)
wmean_f_dur_l_dtsurf = np.nanmean(f_durs_l_dtsurf)

# wmean_f_dur_h_Ug = np.nanmean(df_h_Ug_cldN500['duration'])
# wmean_f_dur_l_Ug = np.nanmean(df_l_Ug_cldN500['duration'])
# wmean_f_dur_h_wmax = np.nanmean(df_h_wmax_cldN500['duration'])
# wmean_f_dur_l_wmax = np.nanmean(df_l_wmax_cldN500['duration'])
# wmean_f_dur_h_dtsurf = np.nanmean(df_h_dtsurf_cldN500['duration'])
# wmean_f_dur_l_dtsurf = np.nanmean(df_l_dtsurf_cldN500['duration'])


# In[17]:


df_vars_in = df_all[['cldN', 'Ug', 'wmax', 'dtsurf','duration', 't_dur', 'nl_dur']]
df_bg_in = df_vars_in.groupby(['Ug', 'wmax', 'dtsurf']).mean()

df_vars_in['bg_count_fog'] = np.empty(np.shape(df_vars_in)[0])
df_vars_in['bg_count_f'] = np.empty(np.shape(df_vars_in)[0])
df_vars_in['bg_count_t'] = np.empty(np.shape(df_vars_in)[0])
df_vars_in['bg_count_nl'] = np.empty(np.shape(df_vars_in)[0])
df_loop = df_vars_in.copy()
df_loop_bg = df_bg_in.copy()
for i in range(np.shape(df_loop_bg)[0]):
    
    # Create filter conditions
    (Ug, wmax, dtsurf) =  df_loop_bg.index[i]
    cond_Ug = df_loop['Ug'] == Ug
    cond_wmax = df_loop['wmax'] == wmax
    cond_dtsurf = df_loop['dtsurf'] == dtsurf
    cond_fog = df_loop['duration'] > 0.0
    cond_t = df_loop['t_dur'] > 0.0
    cond_nl = df_loop['t_dur'] > 0.0
    df_i = df_loop[cond_Ug & cond_wmax & cond_dtsurf]
    df_ic = df_loop[cond_Ug & cond_wmax & cond_dtsurf].reset_index()
    df_i_fog = df_loop[cond_Ug & cond_wmax & cond_dtsurf & cond_fog]
    df_i_t = df_loop[cond_Ug & cond_wmax & cond_dtsurf & cond_t]
    df_i_nl = df_loop[cond_Ug & cond_wmax & cond_dtsurf & cond_nl]
    
    count_fog = 0
    for j in range(10):
        print(i,j,df_ic['duration'][j])
        if df_ic['duration'][j] > 0:
            count_fog = j+1
        else:
            break
    print(count_fog)
    count_f = np.shape(df_i_f)[0]
    count_t = np.shape(df_i_t)[0]
    count_nl = np.shape(df_i_nl)[0]
    
    
    df_vars_in['bg_count_fog'][df_i.index] = count_fog
    df_vars_in['bg_count_f'][df_i.index] = count_f
    df_vars_in['bg_count_t'][df_i.index] = count_t
    df_vars_in['bg_count_nl'][df_i.index] = count_nl
    
cond_f_10 = df_vars_in['bg_count_fog'] == 10.
cond_f_9 = df_vars_in['bg_count_fog'] == 9.
cond_f_8 = df_vars_in['bg_count_fog'] == 8.
cond_f_7 = df_vars_in['bg_count_fog'] == 7.
cond_f_6 = df_vars_in['bg_count_fog'] == 6.
cond_f_5 = df_vars_in['bg_count_fog'] == 5.
cond_f_4 = df_vars_in['bg_count_fog'] == 4.
cond_f_3 = df_vars_in['bg_count_fog'] == 3.
cond_f_2 = df_vars_in['bg_count_fog'] == 2.
cond_f_1 = df_vars_in['bg_count_fog'] == 1.
cond_f_0 = df_vars_in['bg_count_fog'] == 0.
cond_500 = df_vars_in['cldN'] == 500.
c_fog = df_vars_in['duration'] > 0.

f10_500 = df_vars_in[cond_outlier & cond_500 & cond_f_10]
f9_500 = df_vars_in[cond_outlier & cond_500 & cond_f_9]
f8_500 = df_vars_in[cond_outlier & cond_500 & cond_f_8]
f7_500 = df_vars_in[cond_outlier & cond_500 & cond_f_7]
f6_500 = df_vars_in[cond_outlier & cond_500 & cond_f_6]
f5_500 = df_vars_in[cond_outlier & cond_500 & cond_f_5]
f4_500 = df_vars_in[cond_outlier & cond_500 & cond_f_4]
f3_500 = df_vars_in[cond_outlier & cond_500 & cond_f_3]
f2_500 = df_vars_in[cond_outlier & cond_500 & cond_f_2]
f1_500 = df_vars_in[cond_outlier & cond_500 & cond_f_1]

f10_500_f = df_vars_in[cond_outlier & cond_500 & cond_f_10 & c_fog]
f9_500_f = df_vars_in[cond_outlier & cond_500 & cond_f_9 & c_fog]
f8_500_f = df_vars_in[cond_outlier & cond_500 & cond_f_8 & c_fog]
f7_500_f = df_vars_in[cond_outlier & cond_500 & cond_f_7 & c_fog]
f6_500_f = df_vars_in[cond_outlier & cond_500 & cond_f_6 & c_fog]
f5_500_f = df_vars_in[cond_outlier & cond_500 & cond_f_5 & c_fog]
f4_500_f = df_vars_in[cond_outlier & cond_500 & cond_f_4 & c_fog]
f3_500_f = df_vars_in[cond_outlier & cond_500 & cond_f_3 & c_fog]
f2_500_f = df_vars_in[cond_outlier & cond_500 & cond_f_2 & c_fog]
f1_500_f = df_vars_in[cond_outlier & cond_500 & cond_f_1 & c_fog]

f10_dur_500 = np.mean(f10_500['duration'])
f9_dur_500 = np.mean(f9_500['duration'])
f8_dur_500 = np.mean(f8_500['duration'])
f7_dur_500 = np.mean(f7_500['duration'])
f6_dur_500 = np.mean(f6_500['duration'])
f5_dur_500 = np.mean(f5_500['duration'])
f4_dur_500 = np.mean(f4_500['duration'])
f3_dur_500 = np.mean(f3_500['duration'])
f2_dur_500 = np.mean(f2_500['duration'])
f1_dur_500 = np.mean(f1_500['duration'])

f10_dur_500_f = np.mean(f10_500_f['duration'])
f9_dur_500_f = np.mean(f9_500_f['duration'])
f8_dur_500_f = np.mean(f8_500_f['duration'])
f7_dur_500_f = np.mean(f7_500_f['duration'])
f6_dur_500_f = np.mean(f6_500_f['duration'])
f5_dur_500_f = np.mean(f5_500_f['duration'])
f4_dur_500_f = np.mean(f4_500_f['duration'])
f3_dur_500_f = np.mean(f3_500_f['duration'])
f2_dur_500_f = np.mean(f2_500_f['duration'])
f1_dur_500_f = np.mean(f1_500_f['duration'])

df_f10 = df_vars_in[cond_outlier & cond_f_10]
df_f9 = df_vars_in[cond_outlier & cond_f_9]
df_f8 = df_vars_in[cond_outlier & cond_f_8]
df_f7 = df_vars_in[cond_outlier & cond_f_7]
df_f6 = df_vars_in[cond_outlier & cond_f_6]
df_f5 = df_vars_in[cond_outlier & cond_f_5]
df_f4 = df_vars_in[cond_outlier & cond_f_4]
df_f3 = df_vars_in[cond_outlier & cond_f_3]
df_f2 = df_vars_in[cond_outlier & cond_f_2]
df_f1 = df_vars_in[cond_outlier & cond_f_1]
df_f0 = df_vars_in[cond_outlier & cond_f_0]

df_f10_f = df_vars_in[cond_outlier & cond_f_10 & c_fog]
df_f9_f = df_vars_in[cond_outlier & cond_f_9 & c_fog]
df_f8_f = df_vars_in[cond_outlier & cond_f_8 & c_fog]
df_f7_f = df_vars_in[cond_outlier & cond_f_7 & c_fog]
df_f6_f = df_vars_in[cond_outlier & cond_f_6 & c_fog]
df_f5_f = df_vars_in[cond_outlier & cond_f_5 & c_fog]
df_f4_f = df_vars_in[cond_outlier & cond_f_4 & c_fog]
df_f3_f = df_vars_in[cond_outlier & cond_f_3 & c_fog]
df_f2_f = df_vars_in[cond_outlier & cond_f_2 & c_fog]
df_f1_f = df_vars_in[cond_outlier & cond_f_1 & c_fog]

df_f10_cldN = df_f10.groupby(['cldN']).mean()
df_f9_cldN = df_f9.groupby(['cldN']).mean()
df_f8_cldN = df_f8.groupby(['cldN']).mean()
df_f7_cldN = df_f7.groupby(['cldN']).mean()
df_f6_cldN = df_f6.groupby(['cldN']).mean()
df_f5_cldN = df_f5.groupby(['cldN']).mean()
df_f4_cldN = df_f4.groupby(['cldN']).mean()
df_f3_cldN = df_f3.groupby(['cldN']).mean()
df_f2_cldN = df_f2.groupby(['cldN']).mean()
df_f1_cldN = df_f1.groupby(['cldN']).mean()

df_f10_cldN_f = df_f10_f.groupby(['cldN']).mean()
df_f9_cldN_f = df_f9_f.groupby(['cldN']).mean()
df_f8_cldN_f = df_f8_f.groupby(['cldN']).mean()
df_f7_cldN_f = df_f7_f.groupby(['cldN']).mean()
df_f6_cldN_f = df_f6_f.groupby(['cldN']).mean()
df_f5_cldN_f = df_f5_f.groupby(['cldN']).mean()
df_f4_cldN_f = df_f4_f.groupby(['cldN']).mean()
df_f3_cldN_f = df_f3_f.groupby(['cldN']).mean()
df_f2_cldN_f = df_f2_f.groupby(['cldN']).mean()
df_f1_cldN_f = df_f1_f.groupby(['cldN']).mean()


# In[18]:


df_cldN = df_ED_vars.groupby(['cldN']).mean()
df_cldN_h_Ug = df_ED_h_Ug.groupby(['cldN']).mean()
df_cldN_l_Ug = df_ED_l_Ug.groupby(['cldN']).mean()
df_cldN_h_wmax = df_ED_h_wmax.groupby(['cldN']).mean()
df_cldN_l_wmax = df_ED_l_wmax.groupby(['cldN']).mean()
df_cldN_h_dtsurf = df_ED_h_dtsurf.groupby(['cldN']).mean()
df_cldN_l_dtsurf = df_ED_l_dtsurf.groupby(['cldN']).mean()

fig, (ax2, ax3) = plt.subplots(1, 2)
fig.set_size_inches(12,5.5)
plt.subplots_adjust(left=0.075, right=0.95, bottom=0.125, top=0.925, wspace=0.25, hspace=0.2)

# ax1.plot(x, fogamount_cldN, color = 'black', label = 'Overall')
# ax1.plot(x, fogamount_cldN_high_wmax, color = 'blue')
# ax1.plot(x, fogamount_cldN_low_wmax, color = 'blue', linestyle = 'dashed')
# ax1.plot(x, fogamount_cldN_high_Ug, color = 'green')
# ax1.plot(x, fogamount_cldN_low_Ug, color = 'green', linestyle = 'dashed')
# ax1.plot(x, fogamount_cldN_high_dtsurf, color = 'red')
# ax1.plot(x, fogamount_cldN_low_dtsurf, color = 'red', linestyle = 'dashed')
# ax1.legend(fontsize = 14)

ax2.plot(x, fog_cldN, color = 'black', label = 'Overall')
ax2.plot(x, fog_cldN_high_wmax, color = 'blue', label = 'Weak Subsidence')
ax2.plot(x, fog_cldN_low_wmax, color = 'blue', linestyle = 'dashed', label = 'Strong Subsidence')
ax2.plot(x, fog_cldN_high_Ug, color = 'green', label = 'High Wind')
ax2.plot(x, fog_cldN_low_Ug, color = 'green', linestyle = 'dashed', label = 'Low Wind')
ax2.plot(x, fog_cldN_high_dtsurf, color = 'red', label = 'Warming Surface')
ax2.plot(x, fog_cldN_low_dtsurf, color = 'red', linestyle = 'dashed', label = 'Cooling Surface')
ax2.legend(fontsize = 14)

ax3.plot(x, wmean_f_dur + f_durchange, color = 'black', label = 'Overall')
ax3.plot(x, wmean_f_dur_h_wmax + f_durchange_h_wmax, color = 'blue', label = 'Weak Subsidence')
ax3.plot(x, wmean_f_dur_l_wmax + f_durchange_l_wmax, color = 'blue', linestyle = 'dashed', label = 'Strong Subsidence')
ax3.plot(x, wmean_f_dur_h_Ug + f_durchange_h_Ug, color = 'green', label = 'High Wind')
ax3.plot(x, wmean_f_dur_l_Ug + f_durchange_l_Ug, color = 'green', linestyle = 'dashed', label = 'Low Wind')
ax3.plot(x, wmean_f_dur_h_dtsurf + f_durchange_h_dtsurf, color = 'red', label = 'Warming Surface')
ax3.plot(x, wmean_f_dur_l_dtsurf + f_durchange_l_dtsurf, color = 'red', linestyle = 'dashed', label = 'Cooling Surface')
# ax3.legend(fontsize = 14)

# fig.suptitle('Fog Formation Fraction and Expected Duration vs Na Under Different Conditions, no Outliers', fontsize = 20)

# ax1.set_title('Simulation Fogginess', fontsize = 20)
# ax1.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 18)
# ax1.set_ylabel('Average Fog Proportion', fontsize = 18)
# ax1.tick_params(axis='both', which='both', labelsize=14)


# ax2.set_title('Fog Formation Fraction', fontsize = 20)
ax2.set_xlabel('Aerosol Number Concentration (#/cm$^3$)', fontsize = 17)
ax2.set_ylabel('Fraction That Form Fog', fontsize = 20)
ax2.tick_params(axis='both', which='both', labelsize=16)
ax2.text(0.9, 0.95, '(A)', transform = ax2.transAxes, fontsize = 18, fontweight = 'bold')

# ax3.set_title('Expected Fog Duration', fontsize = 20)
ax3.set_xlabel('Aerosol Number Concentration (#/cm$^3$)', fontsize = 17)
ax3.set_ylabel('Fog Duration Difference (hours)', fontsize = 20)
ax3.tick_params(axis='both', which='both', labelsize=16)
ax3.text(0.9, 0.95, '(B)', transform = ax3.transAxes, fontsize = 18, fontweight = 'bold')

fname = desc + 'fig03_old.png'
# plt.savefig(fname)
plt.show()
plt.close()

fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
fig.set_size_inches(20,6)
plt.subplots_adjust(left=0.075, right=0.95, bottom=0.1, top=0.925, wspace=0.25, hspace=0.2)


ax1.plot(x, fog_cldN, color = 'black', label = 'Overall')
ax1.plot(x, fog_cldN_high_wmax, color = 'blue', label = 'Weak Subsidence')
ax1.plot(x, fog_cldN_low_wmax, color = 'blue', linestyle = 'dashed', label = 'Strong Subsidence')
ax1.plot(x, fog_cldN_high_Ug, color = 'green', label = 'High Wind')
ax1.plot(x, fog_cldN_low_Ug, color = 'green', linestyle = 'dashed', label = 'Low Wind')
ax1.plot(x, fog_cldN_high_dtsurf, color = 'red', label = 'Warming Surface')
ax1.plot(x, fog_cldN_low_dtsurf, color = 'red', linestyle = 'dashed', label = 'Cooling Surface')
ax1.legend(fontsize = 14)
ax1.text(0.90, 0.95, '(A)', transform = ax1.transAxes, fontsize = 18, fontweight = 'bold')

ax2.plot(x, df_cldN['duration'], color = 'black', label = 'Overall')
ax2.plot(x, df_cldN_h_wmax['duration'], color = 'blue', label = 'Weak Subsidence')
ax2.plot(x, df_cldN_l_wmax['duration'], color = 'blue', linestyle = 'dashed', label = 'Strong Subsidence')
ax2.plot(x, df_cldN_h_Ug['duration'], color = 'green', label = 'High Wind')
ax2.plot(x, df_cldN_l_Ug['duration'], color = 'green', linestyle = 'dashed', label = 'Low Wind')
ax2.plot(x, df_cldN_h_dtsurf['duration'], color = 'red', label = 'Warming Surface')
ax2.plot(x, df_cldN_l_dtsurf['duration'], color = 'red', linestyle = 'dashed', label = 'Cooling Surface')
ax2.text(0.90, 0.95, '(B)', transform = ax2.transAxes, fontsize = 18, fontweight = 'bold')

ax3.plot(df_f10_cldN_f.index, df_f10_cldN_f['duration'], label = '10 sims')
ax3.plot(df_f9_cldN_f.index, df_f9_cldN_f['duration'], label = '9 sims')
ax3.plot(df_f8_cldN_f.index, df_f8_cldN_f['duration'], label = '8 sims')
ax3.plot(df_f7_cldN_f.index, df_f7_cldN_f['duration'], label = '7 sims')
ax3.plot(df_f6_cldN_f.index, df_f6_cldN_f['duration'], label = '6 sims')
ax3.text(0.90, 0.95, '(C)', transform = ax3.transAxes, fontsize = 18, fontweight = 'bold')
ax3.legend(fontsize = 14)

# fig.suptitle('Fog Formation Fraction and Expected Duration vs Na Under Different Conditions, no Outliers', fontsize = 20)


# ax2.set_title('Fog Formation Fraction', fontsize = 20)
ax1.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 16)
ax1.set_ylabel('Fraction That Form Fog', fontsize = 20)
ax1.tick_params(axis='both', which='both', labelsize=16)

# ax3.set_title('Expected Fog Duration', fontsize = 20)
ax2.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 16)
ax2.set_ylabel('Mean Fog Duration (hours)', fontsize = 20)
ax2.tick_params(axis='both', which='both', labelsize=16)

ax3.set_xlabel('Aerosol Number Concentration (#/cm^3)', fontsize = 16)
ax3.set_ylabel('Mean Fog Duration (hours)', fontsize = 20)
ax3.tick_params(axis='both', which='both', labelsize=16)

fname = desc + 'fig03.png'
plt.savefig(fname)
plt.show()
plt.close()


# In[19]:


fig, ax = plt.subplots()
fig.set_size_inches(7.5,5)
plt.subplots_adjust(left=0.1, right=0.925, bottom=0.135, top=0.925)

ax.plot(x, f_durchange, color = 'black', label = 'Default Fog Condition')
ax.plot(x, t_durchange, color = 'blue', label = '2 m RH Fog Condition')
ax.plot(x, nl_durchange, color = 'red', label = 'Mixed BL Fog Condition')
ax.fill_between(cldN_list, f_durchange - f_95con, f_durchange + f_95con, 
                 color='black', alpha=0.2)
ax.fill_between(cldN_list, t_durchange - t_95con, t_durchange + t_95con, 
                 color='blue', alpha=0.2)
ax.fill_between(cldN_list, nl_durchange - nl_95con, nl_durchange+nl_95con, 
                 color='red', alpha=0.2)

# ax.plot(x, f_durchange + wmean_f_dur, color = 'black', label = 'Default Fog Condition')
# ax.plot(x, t_durchange + wmean_t_dur, color = 'blue', label = '2 m RH Fog Condition')
# ax.plot(x, nl_durchange + wmean_nl_dur, color = 'red', label = 'Mixed BL Fog Condition')
# ax.fill_between(cldN_list, f_durchange + wmean_f_dur - f_95con, f_durchange + wmean_f_dur + f_95con, 
#                  color='black', alpha=0.2)
# ax.fill_between(cldN_list, t_durchange + wmean_t_dur - t_95con, t_durchange + wmean_t_dur + t_95con, 
#                  color='blue', alpha=0.2)
# ax.fill_between(cldN_list, nl_durchange + wmean_nl_dur - nl_95con, nl_durchange+wmean_nl_dur+nl_95con, 
#                  color='red', alpha=0.2)

ax.set_ylabel('Fog Duration Difference (hr)', fontsize = 18)
ax.set_xlabel('Aerosol Number Concentration (#/cm$^3$)', fontsize = 18)
ax.set_title('Fog Duration Using Different Conditions', fontsize = 20)
ax.tick_params(axis='both', which='both', labelsize=14)
ax.legend(fontsize = 16)

fname = desc + 'fig07_old.png'
# plt.savefig(fname)
plt.show()
plt.close()


# In[20]:


print(np.corrcoef(df_nout['cldN'], df_nout['duration'])[0,1], 
      np.corrcoef(df_nout['Ug'], df_nout['duration'])[0,1], 
      np.corrcoef(df_nout['wmax'], df_nout['duration'])[0,1], 
      np.corrcoef(df_nout['dtsurf'], df_nout['duration'])[0,1])


# In[21]:


cldN_100 = df_all['cldN'] == 100.
cldN_150 = df_all['cldN'] == 150.
cldN_200 = df_all['cldN'] == 200.
cldN_300 = df_all['cldN'] == 300.
cldN_400 = df_all['cldN'] == 400.
cldN_500 = df_all['cldN'] == 500.
cldN_750 = df_all['cldN'] == 750.
cldN_1000 = df_all['cldN'] == 1000.
cldN_1250 = df_all['cldN'] == 1250.
cldN_1500 = df_all['cldN'] == 1500.

print(np.shape(df_all[fog & cond_outlier & Ug_5])[0] / np.shape(df_all[cond_outlier & Ug_5])[0], 
      np.shape(df_all[fog & cond_outlier & Ug_7])[0] / np.shape(df_all[cond_outlier & Ug_7])[0], 
      np.shape(df_all[fog & cond_outlier & base_Ug])[0] / np.shape(df_all[cond_outlier & base_Ug])[0], 
      np.shape(df_all[fog & cond_outlier & Ug_12])[0] / np.shape(df_all[cond_outlier & Ug_12])[0], 
      np.shape(df_all[fog & cond_outlier & Ug_15])[0] / np.shape(df_all[cond_outlier & Ug_15])[0])

print(np.shape(df_all[fog & cond_outlier & wmax_35])[0] / np.shape(df_all[cond_outlier & wmax_35])[0], 
      np.shape(df_all[fog & cond_outlier & wmax_325])[0] / np.shape(df_all[cond_outlier & wmax_325])[0], 
      np.shape(df_all[fog & cond_outlier & base_wmax])[0] / np.shape(df_all[cond_outlier & base_wmax])[0], 
      np.shape(df_all[fog & cond_outlier & wmax_275])[0] / np.shape(df_all[cond_outlier & wmax_275])[0], 
      np.shape(df_all[fog & cond_outlier & wmax_25])[0] / np.shape(df_all[cond_outlier & wmax_25])[0])

print(np.shape(df_all[fog & cond_outlier & dtsurf__2])[0] / np.shape(df_all[cond_outlier & dtsurf__2])[0], 
      np.shape(df_all[fog & cond_outlier & dtsurf__1])[0] / np.shape(df_all[cond_outlier & dtsurf__1])[0], 
      np.shape(df_all[fog & cond_outlier & base_dtsurf])[0] / np.shape(df_all[cond_outlier & base_dtsurf])[0], 
      np.shape(df_all[fog & cond_outlier & dtsurf_1])[0] / np.shape(df_all[cond_outlier & dtsurf_1])[0], 
      np.shape(df_all[fog & cond_outlier & dtsurf_2])[0] / np.shape(df_all[cond_outlier & dtsurf_2])[0])

print(np.shape(df_all[fog & cond_outlier & cldN_100])[0] / np.shape(df_all[cond_outlier & cldN_100])[0], 
      np.shape(df_all[fog & cond_outlier & cldN_150])[0] / np.shape(df_all[cond_outlier & cldN_150])[0], 
      np.shape(df_all[fog & cond_outlier & cldN_200])[0] / np.shape(df_all[cond_outlier & cldN_200])[0], 
      np.shape(df_all[fog & cond_outlier & cldN_300])[0] / np.shape(df_all[cond_outlier & cldN_300])[0], 
      np.shape(df_all[fog & cond_outlier & cldN_400])[0] / np.shape(df_all[cond_outlier & cldN_400])[0], 
      np.shape(df_all[fog & cond_outlier & cldN_500])[0] / np.shape(df_all[cond_outlier & cldN_500])[0], 
      np.shape(df_all[fog & cond_outlier & cldN_750])[0] / np.shape(df_all[cond_outlier & cldN_750])[0], 
      np.shape(df_all[fog & cond_outlier & cldN_1000])[0] / np.shape(df_all[cond_outlier & cldN_1000])[0], 
      np.shape(df_all[fog & cond_outlier & cldN_1250])[0] / np.shape(df_all[cond_outlier & cldN_1250])[0], 
      np.shape(df_all[fog & cond_outlier & cldN_1500])[0] / np.shape(df_all[cond_outlier & cldN_1500])[0])


# In[22]:


dfv_no_out = df_vars_in[cond_outlier]

df_cldN_std = dfv_no_out.groupby(['cldN']).agg(['mean', 'std','count'])

fig, ax = plt.subplots()
fig.set_size_inches(7.5,5)
plt.subplots_adjust(left=0.1, right=0.925, bottom=0.135, top=0.925)

ax.plot(x, df_cldN_std['duration']['mean'] - df_cldN_std['duration']['mean'][500.], color = 'black', label = 'Default Fog Condition')
# ax.fill_between(x, df_cldN_std['duration']['mean'] - 2 * (df_cldN_std['duration']['std'] / np.sqrt(df_cldN_std['duration']['count'] - 2)) - df_cldN_std['duration']['mean'][500.], 
#                  df_cldN_std['duration']['mean'] + 2 * (df_cldN_std['duration']['std'] / np.sqrt(df_cldN_std['duration']['count'] - 2)) - df_cldN_std['duration']['mean'][500.], 
#                  color='black', alpha=0.2)
ax.plot(x, df_cldN_std['t_dur']['mean'] - df_cldN_std['t_dur']['mean'][500.], color = 'blue', label = '2 m RH Fog Condition')
# ax.fill_between(x, df_cldN_std['t_dur']['mean'] - 2 * (df_cldN_std['t_dur']['std'] / np.sqrt(df_cldN_std['t_dur']['count'] - 2)) - df_cldN_std['t_dur']['mean'][500.], 
#                  df_cldN_std['t_dur']['mean'] + 2 * (df_cldN_std['t_dur']['std'] / np.sqrt(df_cldN_std['t_dur']['count'] - 2)) - df_cldN_std['t_dur']['mean'][500.], 
#                  color='blue', alpha=0.2)
ax.plot(x, df_cldN_std['nl_dur']['mean'] - df_cldN_std['nl_dur']['mean'][500.], color = 'red', label = 'Mixed BL Fog Condition')
# ax.fill_between(x, df_cldN_std['nl_dur']['mean'] - 2 * (df_cldN_std['nl_dur']['std'] / np.sqrt(df_cldN_std['nl_dur']['count'] - 2)) - df_cldN_std['nl_dur']['mean'][500.], 
#                  df_cldN_std['nl_dur']['mean'] + 2 * (df_cldN_std['nl_dur']['std'] / np.sqrt(df_cldN_std['nl_dur']['count'] - 2)) - df_cldN_std['nl_dur']['mean'][500.], 
#                  color='red', alpha=0.2)

# ax.plot(x, f_durchange + wmean_f_dur, color = 'black', label = 'Default Fog Condition')
# ax.plot(x, t_durchange + wmean_t_dur, color = 'blue', label = '2 m RH Fog Condition')
# ax.plot(x, nl_durchange + wmean_nl_dur, color = 'red', label = 'Mixed BL Fog Condition')
# ax.fill_between(cldN_list, f_durchange + wmean_f_dur - f_95con, f_durchange + wmean_f_dur + f_95con, 
#                  color='black', alpha=0.2)
# ax.fill_between(cldN_list, t_durchange + wmean_t_dur - t_95con, t_durchange + wmean_t_dur + t_95con, 
#                  color='blue', alpha=0.2)
# ax.fill_between(cldN_list, nl_durchange + wmean_nl_dur - nl_95con, nl_durchange+wmean_nl_dur+nl_95con, 
#                  color='red', alpha=0.2)

ax.set_ylabel('Fog Duration Difference (hr)', fontsize = 18)
ax.set_xlabel('Aerosol Number Concentration (#/cm$^3$)', fontsize = 18)
ax.set_title('Fog Duration Using Different Conditions', fontsize = 20)
ax.tick_params(axis='both', which='both', labelsize=14)
ax.legend(fontsize = 16)

fname = desc + 'fig07.png'
plt.savefig(fname)
plt.show()
plt.close()


# In[ ]:





# In[ ]:





# In[ ]: