ParameterComp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5e92fc74",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Calculate and save selected properties in different parameter cases\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os, sys\n",
    "sys.path.append(os.path.abspath('./TEF'))\n",
    "import TEF_Variables as tef"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "495e491b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obtain reflux coefficients and entrainment flux Qe\n",
    "def reflux(case_id, xi):\n",
    "\n",
    "#     State = xr.open_dataset('/Volumes/Extreme SSD/MITgcm outputs/Qsg' + str(format(Qsg,'03d')) + '_tide003/state_' + str(format(hs,'03d')) + '.nc')\n",
    "#     Grid = xr.open_dataset('/Volumes/Extreme SSD/MITgcm outputs/Qsg' + str(format(Qsg,'03d')) + '_tide003/grid_' + str(format(hs,'03d')) + '.nc')\n",
    "    State0 = xr.open_dataset('/work/oceans/wbao/MITgcm_results/iceplume/weak_Tide_maxhs/state_' + str(format(case_id,'03d')) + '.nc')\n",
    "    Grid = xr.open_dataset('/work/oceans/wbao/MITgcm_results/iceplume/weak_Tide_maxhs/grid_' + str(format(case_id,'03d')) + '.nc')\n",
    "    State = State0.isel(T=~State0.get_index(\"T\").duplicated())\n",
    "    # Confine to the range of fjord\n",
    "    state = State.isel(X=range(260), Xp1=range(261), Y=range(35,45), T=range(1,len(State.T)))\n",
    "    grid = Grid.isel(X=range(260), Xp1=range(261), Y=range(35,45))\n",
    "\n",
    "\n",
    "    s = state.S.data.mean(2) # Salinity in T,Z,X\n",
    "    u = (state.U.data[:,:,:,1:].mean(2) + state.U.data[:,:,:,:-1].mean(2)) / 2 # Along channel velocity\n",
    "    ot = state.T.data # Time in seconds\n",
    "\n",
    "\n",
    "    HFacC1 = grid.HFacC.data.mean(1)\n",
    "    dyF1 = grid.dyF.data.mean(0)\n",
    "    drF1 = grid.drF.data\n",
    "    gridA = np.broadcast_to(drF1[:, np.newaxis], HFacC1.shape) * np.broadcast_to(dyF1[np.newaxis, :], HFacC1.shape) * HFacC1 # Grid Area\n",
    "    da = np.broadcast_to(gridA[np.newaxis,:,:], u.shape)\n",
    "\n",
    "    S = state.S.data\n",
    "    U = (state.U.data[:,:,:,1:] + state.U.data[:,:,:,:-1]) / 2\n",
    "    drF = np.broadcast_to(grid.drF.data[np.newaxis, :, np.newaxis, np.newaxis], U.shape)\n",
    "    dyF = np.broadcast_to(grid.dyF.data[np.newaxis, np.newaxis, :, :], U.shape)\n",
    "    HFacC = np.broadcast_to(grid.HFacC.data[np.newaxis, :, :, :], U.shape)\n",
    "    DA = drF * dyF * HFacC\n",
    "    \n",
    "    \n",
    "    Qin = np.empty(len(xi))\n",
    "    Qout = np.empty(len(xi))\n",
    "    Sin = np.empty(len(xi))\n",
    "    Sout = np.empty(len(xi))\n",
    "    t0 = 89 # 274\n",
    "#    Exclude the upper (surface) layer for TEF calculation\n",
    "#     for t in range(U.shape[0]):\n",
    "        \n",
    "#         Uzx = u[t,:,:] # Steady-state along-channel velocity\n",
    "#         for j in range(Uzx.shape[1]):\n",
    "#             if any(Uzx[:,j]>0)==False:\n",
    "#                 pass\n",
    "#             else:\n",
    "#                 mid = np.where(Uzx[:,j]>0)[0][0]\n",
    "#                 if mid>0:\n",
    "#                     S[t,:mid,:,j] = 0\n",
    "#                     U[t,:mid,:,j] = 0\n",
    "#                 else:\n",
    "#                     pass\n",
    "    \n",
    "    Uzx = u[t0:,:,:].mean(0) # Steady-state along-channel velocity\n",
    "    for j in range(Uzx.shape[1]):\n",
    "        if any(Uzx[:,j]>0)==False:\n",
    "            pass\n",
    "        else:\n",
    "            mid = np.where(Uzx[:,j]>0)[0][0]\n",
    "            if mid>0:\n",
    "                S[:,:mid,:,j] = 0\n",
    "                U[:,:mid,:,j] = 0\n",
    "            else:\n",
    "                pass\n",
    "    \n",
    "    for i in range(len(xi)):\n",
    "    \n",
    "        tef_q1, tef_vel1, tef_da1, tef_qs1, tef_qs21, sbins1 = tef.process_section(U,S,DA,ot,xi[i],23,testing=False)\n",
    "\n",
    "        qin1, qout1, sin1, sout1 = tef.bulk_calc(tef_q1, tef_vel1, tef_da1, tef_qs1, tef_qs21, sbins1, ot)\n",
    "    \n",
    "        Qin[i] = qin1[t0:].mean() # Time averaging from 1+t0 hour\n",
    "        Qout[i] = qout1[t0:].mean()\n",
    "        Sin[i] = sin1[t0:].mean() # Time averaging from 1+t0 hour\n",
    "        Sout[i] = sout1[t0:].mean()\n",
    "        \n",
    "    a002 = (Sin[0]/Sout[0])*(Sout[1]-Sout[0])/(Sout[1]-Sin[0])\n",
    "    a001 = (Qin[0]-Qin[1]) / Qin[0]    \n",
    "    a11 = (Sout[1]/Sin[1])*(Sin[1]-Sin[0])/(Sout[1]-Sin[0])    \n",
    "    a01 = (Sout[1]/Sout[0])*(Sout[0]-Sin[0])/(Sout[1]-Sin[0])\n",
    "    \n",
    "    a10 = (Sin[0]/Sin[1])*(Sout[1]-Sin[1])/(Sout[1]-Sin[0])\n",
    "    \n",
    "    #a00 = -(Qout[0]/Qin[0])*(Sout[1]-Sout[0])/(Sout[1]-Sin[0])\n",
    "    #a11 = -(Qin[1]/Qout[1])*(Sin[1]-Sin[0])/(Sout[1]-Sin[0])\n",
    "    #a01 = (Qin[1]/Qin[0])*(Sout[1]-Sin[1])/(Sout[1]-Sin[0])\n",
    "    #a10 = (Qout[0]/Qout[1])*(Sout[0]-Sin[0])/(Sout[1]-Sin[0])\n",
    "    \n",
    "    q0 = Qin[0]\n",
    "    q1 = -Qout[-1]\n",
    "    Q0 = -Qout[0]\n",
    "    Q1 = Qin[-1]\n",
    "    f0 = Sin[0]*q0\n",
    "    f1 = Sout[-1]*q1\n",
    "    F0 = Sout[0]*Q0\n",
    "    F1 = Sin[-1]*Q1\n",
    "    \n",
    "    A = np.array([[q1, q0, 0, 0], [f1, f0, 0, 0], [0, 0, q1, q0], [0, 0, f1, f0]])\n",
    "    B = np.array([Q1, F1, Q0, F0])\n",
    "    X = np.linalg.solve(A,B) # a11, a01, a10, a00\n",
    "    \n",
    "    Qe = Q0 - q1\n",
    "    \n",
    "#     if (a11 < 0.1) or (a11 > 1):\n",
    "#         a22 = a001\n",
    "#     else:\n",
    "#         a22 = a002\n",
    "    if (X[-1] < 0) or (X[-1] > 1):\n",
    "        a22 = a001\n",
    "        a11 = 0\n",
    "    else:\n",
    "        a22 = X[-1]\n",
    "        a11 = X[0]\n",
    "    \n",
    "    return a11,a22,q1,q0,Q1,Q0\n",
    "    #return a11,a22,q1,q0,Q1,Q0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "eef3318c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def vol_temp(case_id, xi):\n",
    "    \n",
    "#     State = xr.open_dataset('/Volumes/Extreme SSD/MITgcm outputs/Qsg' + str(format(Qsg,'03d')) + '_tide003/state_' + str(format(hs,'03d')) + '.nc')\n",
    "#     Grid = xr.open_dataset('/Volumes/Extreme SSD/MITgcm outputs/Qsg' + str(format(Qsg,'03d')) + '_tide003/grid_' + str(format(hs,'03d')) + '.nc')\n",
    "    State0 = xr.open_dataset('/work/oceans/wbao/MITgcm_results/iceplume/1_BaseCase/state_' + str(format(case_id,'03d')) + '.nc')\n",
    "    Grid = xr.open_dataset('/work/oceans/wbao/MITgcm_results/iceplume/1_BaseCase/grid_' + str(format(case_id,'03d')) + '.nc')\n",
    "    State = State0.isel(T=~State0.get_index(\"T\").duplicated())    \n",
    "    # Confine to the range of fjord\n",
    "    state = State.isel(X=range(260), Xp1=range(261), Y=range(35,45), T=range(1,len(State.T)))\n",
    "    grid = Grid.isel(X=range(260), Xp1=range(261), Y=range(35,45))\n",
    "\n",
    "    U = (state.U.data[:,:,:,1:] + state.U.data[:,:,:,:-1]) / 2 # Along-channel velocity\n",
    "    \n",
    "    drF = np.broadcast_to(grid.drF.data[np.newaxis, :, np.newaxis, np.newaxis], U.shape)\n",
    "    dyF = np.broadcast_to(grid.dyF.data[np.newaxis, np.newaxis, :, :], U.shape)\n",
    "    HFacC = np.broadcast_to(grid.HFacC.data[np.newaxis, :, :, :], U.shape)\n",
    "    DA = drF * dyF * HFacC\n",
    "    \n",
    "    rA = np.broadcast_to(grid.rA.data[np.newaxis, np.newaxis, :, :], U.shape)\n",
    "    CV = rA * drF * HFacC\n",
    "\n",
    "    \n",
    "    t0 = 89 # 274\n",
    "    da = DA[t0:,:,:,xi[0]+1:xi[1]+1].mean(axis=(0,2))\n",
    "    cv = CV[t0:,:,:,xi[0]+1:xi[1]+1].mean(axis=(0,2))\n",
    "    s = state.S.data[t0:,:,:,xi[0]+1:xi[1]+1].mean(axis=(0,2))\n",
    "    temp = state.Temp.data[t0:,:,:,xi[0]+1:xi[1]+1].mean(axis=(0,2))\n",
    "    u = U[t0:,:,:,xi[0]+1:xi[1]+1].mean(axis=(0,2))\n",
    "    sma = np.ma.masked_where(s==0, s)\n",
    "    topo = np.ma.getmask(sma) # Masked Topography\n",
    "    uma = np.ma.MaskedArray(u, mask=topo)\n",
    "    tma = np.ma.MaskedArray(temp, mask=topo)\n",
    "    \n",
    "    \n",
    "    Tv1 = np.empty(uma.shape[1])\n",
    "    Tv2 = np.empty(uma.shape[1])\n",
    "    Vol1 = np.empty(uma.shape[1])\n",
    "    Vol2 = np.empty(uma.shape[1])\n",
    "    Tvol1 = np.empty(uma.shape[1])\n",
    "    Tvol2 = np.empty(uma.shape[1])\n",
    "    L = np.empty(uma.shape[1])\n",
    "    \n",
    "    for i in range(uma.shape[1]):\n",
    "\n",
    "        # Exclude the surface layer \n",
    "        if any(uma[:,i]>0)==False:\n",
    "            pass\n",
    "        else:\n",
    "            mid = np.where(uma[:,i]>0)[0][0]\n",
    "            if mid>0:\n",
    "                uma[:mid,i] = 0\n",
    "            else:\n",
    "                pass\n",
    "            \n",
    "        if any(uma[:,i]<0)==False:\n",
    "            pass\n",
    "            \n",
    "        else:\n",
    "            #meanS = np.abs(sma[:,i] - (sma[:,i].max() + sma[:,i].min()) / 2)\n",
    "            #l = np.argmin(meanS)\n",
    "            l = np.where(uma[:,i] < 0)[-1][0]    \n",
    "            Tv1[i] = np.sum(tma[:l,i]*da[:l,i]*uma[:l,i]) / np.sum(da[:l,i]*uma[:l,i])\n",
    "            Tv2[i] = np.sum(tma[l:,i]*da[l:,i]*uma[l:,i]) / np.sum(da[l:,i]*uma[l:,i])\n",
    "            Vol1[i] = cv[:l,i].sum()\n",
    "            Vol2[i] = cv[l:,i].sum()\n",
    "            Tvol1[i] = np.sum(tma[:l,i]*cv[:l,i])\n",
    "            Tvol2[i] = np.sum(tma[l:,i]*cv[l:,i])\n",
    "            L[i] = l\n",
    "            \n",
    "        Tf = Tv1[:-1].mean()\n",
    "        Ts = Tv2[:-1].mean()\n",
    "        #Tf = Tvol1[:-1].sum()/Vol1[:-1].sum()\n",
    "        #Ts = Tvol2[:-1].sum()/Vol2[:-1].sum()\n",
    "        Ts_in = Tv2[-1]\n",
    "        \n",
    "    return Tf, Ts, Ts_in\n",
    "    #return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d1c99c8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obtain submarine melting Qsm\n",
    "def IFA(case_id):\n",
    "    # Grid areas\n",
    "#     Area = np.empty([50, 9])\n",
    "#     Area[:10,:] = 100\n",
    "#     Area[10:20,:] = 200\n",
    "#     Area[20:30,:] = 300\n",
    "#     Area[30:40,:] = 800\n",
    "#     Area[40:,:] = 1600\n",
    "    Area = np.empty([90, 10])\n",
    "    Area[:20,:] = 400\n",
    "    Area[20:50,:] = 800\n",
    "    Area[50:,:] = 1200\n",
    "    file0 = xr.open_dataset('/work/oceans/wbao/MITgcm_results/iceplume/1_BaseCase/icefrntA_' + str(format(case_id,'03d')) + '.nc')\n",
    "    file = file0.isel(T=~file0.get_index(\"T\").duplicated()) \n",
    "    t0 = 89\n",
    "    tn = len(file.T)\n",
    "    state = file.isel(Y=range(35,45), T=range(t0,tn))\n",
    "    MR = state.icefrntA.isel(X=1).data.mean(0) # Melt rate at the icefront\n",
    "    Qsm = (MR*Area).sum()/(24*3600)\n",
    "    return Qsm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "64667b01",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  time 0 out of 120\n",
      "  time 0 out of 120\n",
      "  time 0 out of 120\n",
      "  time 0 out of 120\n"
     ]
    }
   ],
   "source": [
    "\n",
    "case = np.array([1,5])\n",
    "#hsr = np.array([0.04,0.06,0.08,0.10,0.12])\n",
    "Qsg = np.array([25,50,100,250,500,1000])\n",
    "T0 = np.array([2,4,6,8,12])\n",
    "\n",
    "\n",
    "# xrange1 = np.array([5,47]) # X index range for the fjord interior\n",
    "# xrange2 = np.array([47,97]) # X index range for the sill segment\n",
    "xrange = np.array([120,235]) # X index range for the sill segment\n",
    "\n",
    "alp11 = np.empty(len(case))\n",
    "alp22 = np.empty(len(case))\n",
    "Qin0 = np.empty(len(case))\n",
    "Qin1 = np.empty(len(case))\n",
    "Qout0 = np.empty(len(case))\n",
    "Qout1 = np.empty(len(case))\n",
    "Ts = np.empty(len(case))\n",
    "Tf = np.empty(len(case))\n",
    "Ts_in = np.empty(len(case))\n",
    "\n",
    "\n",
    "\n",
    "# alp11_1 = np.empty(len(hs))\n",
    "# alp00_1 = np.empty(len(hs))\n",
    "# Qin0_1 = np.empty(len(hs))\n",
    "# Qin1_1 = np.empty(len(hs))\n",
    "# Qout0_1 = np.empty(len(hs))\n",
    "# Ts_1 = np.empty(len(hs))\n",
    "# Tf_1 = np.empty(len(hs))\n",
    "# Ts_in_1 = np.empty(len(hs))\n",
    "\n",
    "# alp11_2 = np.empty(len(hs))\n",
    "# alp00_2 = np.empty(len(hs))\n",
    "# Qin0_2 = np.empty(len(hs))\n",
    "# Qin1_2 = np.empty(len(hs))\n",
    "# Qout0_2 = np.empty(len(hs))\n",
    "# Ts_2 = np.empty(len(hs))\n",
    "# Tf_2 = np.empty(len(hs))\n",
    "# Ts_in_2 = np.empty(len(hs))\n",
    "\n",
    "Qsm = np.empty(len(case))\n",
    "\n",
    "for j in range(len(case)):\n",
    "#     alp11_1[j], alp00_1[j], Qin1_1[j], Qin0_1[j], Qout1_1, Qout0_1[j] = reflux(qsg,hs[j],xrange1)\n",
    "#     Tf_1[j], Ts_1[j], Ts_in_1[j] = vol_temp(qsg,hs[j],xrange1)\n",
    "    \n",
    "#     alp11_2[j], alp00_2[j], Qin1_2[j], Qin0_2[j], Qout1_2, Qout0_2[j] = reflux(qsg,hs[j],xrange2)\n",
    "#     Tf_2[j], Ts_2[j], Ts_in_2[j] = vol_temp(qsg,hs[j],xrange2)\n",
    "    alp11[j], alp22[j], Qin1[j], Qin0[j], Qout1[j], Qout0[j] = reflux(case[j],xrange)\n",
    "    Tf[j], Ts[j], Ts_in[j] = vol_temp(case[j],xrange)\n",
    "\n",
    "    Qsm[j] = IFA(case[j])\n",
    "\n",
    "# Ts inside fjord\n",
    "#Tse_1 = (Ts_in_1*Qin1_1*(1-alp11_1) + Tf_1*Qin0_1*alp00_1) / Qout0_1\n",
    "\n",
    "# Ts at sill\n",
    "#Tse_2 = (Ts_in_2*Qin1_2*(1-alp11_2) + Tf_2*Qin0_2*alp00_2) / Qout0_2\n",
    "Tse = (Ts_in*Qin1*(1-alp11) + Tf*Qin0*alp22) / Qout0\n",
    "Tr = Tse / Ts_in\n",
    "\n",
    "Qe = Qout0 - Qin1\n",
    "#Qe2 = Qout1 - Qin0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "070049dc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.95936279, 0.99274798])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Qout0, Qin0, Qin1, Qout1\n",
    "Ts / Ts_in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4a0da3a2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0.06345722, -0.02452123]),\n",
       " array([0.70398433, 0.17783125]),\n",
       " array([0.90468664, 0.98712014]),\n",
       " array([6353.05710586, 1629.93682006]),\n",
       " array([10.50586356, 11.23190402]))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#alp11_1.shape, Qe.shape, Tse_1.shape\n",
    "#Tf_1, Tf_2, Qin0_1, Qin0_2\n",
    "#ll = vol_temp(qsg, 5, xrange1)\n",
    "alp11,alp22,Tr, Qe, Qsm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d075fdd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save outputs\n",
    "Cd = np.array([1, 2.5e-3, 10e-3, 25e-3])\n",
    "U0 = np.array([1e-3, 2e-3, 3e-3, 4e-3, 5e-3])\n",
    "TEF_ds = xr.Dataset(\n",
    "    data_vars={'Cd' : Cd,\n",
    "    'alpha11' : alp11,\n",
    "    'alpha22' : alp22,\n",
    "    'Qe' : Qe,\n",
    "    'Ts_in' : Ts_in,\n",
    "    'Ts_e' : Tse,\n",
    "    'Ts' : Ts,\n",
    "    'Qsm' : Qsm})\n",
    "\n",
    "#outdir = \"/Volumes/Extreme SSD/MITgcm outputs/Qsg_comp/\"\n",
    "#outdir = \"/Volumes/Extreme SSD/MITgcm outputs/Strat_comp/\"\n",
    "#outdir = \"/Volumes/Extreme SSD/MITgcm outputs/Cd_comp/\"\n",
    "outdir = \"/home/1959/Parameters/\"\n",
    "if not os.path.exists(outdir):\n",
    "    os.makedirs(outdir)\n",
    "\n",
    "#TEF_ds.to_netcdf(outdir + 'Qsg_' + str(format(qsg,'04d')) + '.nc')\n",
    "#TEF_ds.to_netcdf(outdir + 'Cd_004.nc')\n",
    "TEF_ds.to_netcdf(outdir + 'Cd_025.nc')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7487a700",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quicklook of the output/input file\n",
    "h = 500\n",
    "\n",
    "file = xr.open_dataset('/Volumes/Extreme SSD/MITgcm outputs/Qsg_comp/Qsg_' + str(format(h,'04d')) + '.nc')\n",
    "\n",
    "file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15904381",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1_Plot for Qsg cases\n",
    "Qsg = np.array([100, 250, 500, 1000])\n",
    "\n",
    "fig1 = plt.figure(figsize=(10,10))\n",
    "#fig1.tight_layout()\n",
    "#plt.grid(linestyle='--', linewidth=0.5)    \n",
    "ax1 = fig1.add_subplot(221)\n",
    "ax2 = fig1.add_subplot(222)\n",
    "ax3 = fig1.add_subplot(223)\n",
    "ax4 = fig1.add_subplot(224)\n",
    "\n",
    "\n",
    "from matplotlib.legend_handler import HandlerBase\n",
    "\n",
    "\n",
    "# class AnyObjectHandler(HandlerBase):\n",
    "#     def create_artists(self, legend, orig_handle,\n",
    "#                        x0, y0, width, height, fontsize, trans):\n",
    "#         l1 = plt.Line2D([x0,y0+width], [0.7*height,0.7*height],\n",
    "#                            linestyle=orig_handle[1],color=orig_handle[0])\n",
    "#         l2 = plt.Line2D([x0,y0+width], [0.3*height,0.3*height], \n",
    "#                            color=orig_handle[0])\n",
    "#         return [l1, l2]\n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "for i in range(len(Qsg)):\n",
    "    \n",
    "    file1 = xr.open_dataset('/Volumes/Extreme SSD/MITgcm outputs/Qsg_comp/Qsg_' + str(format(Qsg[i],'04d')) + '.nc')\n",
    "    hsr = file1.hsr.data\n",
    "    alp11 = file1.alpha11.data\n",
    "    alp11[alp11<0] = 0\n",
    "    alp00 = file1.alpha00.data\n",
    "    Ts_f = file1.Ts_fjord.data\n",
    "    Ts_s = file1.Ts_sill.data\n",
    "    Qsm = file1.Qsm.data   \n",
    "    Qe = file1.Qe.data\n",
    "\n",
    "    # Subplot of reflux coefficients\n",
    "    if Qsg[i] == 100:\n",
    "        l1,l2 = ax1.plot(hsr,alp00,'o-', hsr,alp11,'^-', ms=7, alpha=0.7, color='#1f77b4',lw=2)\n",
    "    #ax1.set_ylim(0,.7)\n",
    "        \n",
    "    elif Qsg[i] == 250:\n",
    "        \n",
    "        l3,l4 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#ff7f0e',lw=2)\n",
    "        \n",
    "    elif Qsg[i] == 500:\n",
    "        \n",
    "        l5,l6 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#2ca02c',lw=2)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        l7,l8 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#d62728',lw=2)\n",
    "    \n",
    "#     ax1.legend((l1,l2,l3,l4,l5,l6,l7,l8),('100 $T_S$','100 $T_S2$','250 $T_S1$','250 $T_S2$','500 $T_S1$','500 $T_S2$','1000 $T_S1$','1000 $T_S2$'),\n",
    "#                title=r'$Q_{sg}\\;[m^3/s]$', loc='best', fontsize=12)\n",
    "\n",
    "\n",
    "\n",
    "#     #Plot analytic solution\n",
    "#     ax1.plot(hsr,alp11,'o-',lw=2,label=str(Qsg[i]))\n",
    "\n",
    "#     #Plot simulation\n",
    "#     ax1.plot(hsr,alp00,'^-',lw=2,label=str(Qsg[i]))\n",
    "\n",
    "    \n",
    "#     #Get artists and labels for legend and chose which ones to display\n",
    "#     handles, labels = ax1.get_legend_handles_labels()\n",
    "#     display = (0,1,2,3,4)\n",
    "\n",
    "#     #Create custom artists\n",
    "    \n",
    "#     simArtist = plt.Line2D((0,1),(0,0), color='k', marker='o', linestyle='')\n",
    "#     anyArtist = plt.Line2D((0,1),(0,0), color='k', marker='^', linestyle='')\n",
    "\n",
    "#     #Create legend from custom artist/label lists\n",
    "#     ax1.legend([handle for i,handle in enumerate(handles) if i in display]+[simArtist,anyArtist],\n",
    "#           [label for i,label in enumerate(labels) if i in display]+['Simulation', 'Analytic'], title=r'$Q_{sg}\\;[m^3/s]$',loc='best', fontsize=12)\n",
    "    \n",
    "    #line = Line2D([0,1],[0,1],linestyle='-', color='r')\n",
    "    line1 = plt.Line2D((0,1),(0,0), color='#1f77b4', lw=2, linestyle='-')\n",
    "    line2 = plt.Line2D((0,1),(0,0), color='#ff7f0e', lw=2, linestyle='-')\n",
    "    line3 = plt.Line2D((0,1),(0,0), color='#2ca02c', lw=2, linestyle='-')\n",
    "    line4 = plt.Line2D((0,1),(0,0), color='#d62728', lw=2, linestyle='-')\n",
    "    marker1 = plt.Line2D((0,1),(0,0), color='k', marker='o', linestyle='')\n",
    "    marker2 = plt.Line2D((0,1),(0,0), color='k', marker='^', linestyle='')    \n",
    "    \n",
    "    ax1.legend([line1,line2,line3,line4,marker1,marker2],['100','250','500','1000',r'$\\alpha_{11}$', r'$\\alpha_{00}$'],\n",
    "              title=r'$Q_{sg}\\;[m^3/s]$',loc='best', fontsize=10)\n",
    "\n",
    "    ax1.grid(linestyle='--', linewidth=0.5)\n",
    "    ax1.set_ylabel(r'Reflux Coefficient',size = 12)\n",
    "    #ax1.xaxis.set_visible(False)\n",
    "    ax1.set_title(r'(a)', loc='left', size=12, weight='bold')\n",
    "\n",
    "    \n",
    "    # Subplot of entrainment flux Qe\n",
    "    ax2.plot(hsr,Qe,'o-',lw=2,label=str(Qsg[i]))\n",
    "    #ax2.set_ylim(0,4500)\n",
    "    ax2.legend(title=r'$Q_{sg}\\;[m^3/s]$',loc='best', fontsize=10)\n",
    "    ax2.grid(linestyle='--', linewidth=0.5)\n",
    "    ax2.set_ylabel(r'$Q_e\\;[m^3/s]$',size = 12)\n",
    "    #ax2.xaxis.set_visible(False)\n",
    "    ax2.set_title(r'(b)', loc='left', size=12, weight='bold')\n",
    "       \n",
    "    \n",
    "    # Subplot of deep-layer temperature Ts\n",
    "    if Qsg[i] == 100:\n",
    "           \n",
    "        ax3.plot(hsr,Ts_s,'o-', hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#1f77b4',lw=2)\n",
    "    #ax1.set_ylim(0,.7)\n",
    "        \n",
    "    elif Qsg[i] == 250:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#ff7f0e',lw=2)\n",
    "        \n",
    "    elif Qsg[i] == 500:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#2ca02c',lw=2)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#d62728',lw=2)   \n",
    "    \n",
    "#     ax3.plot(hsr,Ts_s,'o-',lw=2,label=r'$Q_{sg}$')\n",
    "#     ax3.plot(hsr,Ts_f,'^-',lw=2,label=r'$Q_{sg}$')\n",
    "\n",
    "    #ax3.set_ylim(0,.7)\n",
    "    #ax3.legend(loc='best', fontsize=12)\n",
    "    ax3.legend([line1,line2,line3,line4,marker1,marker2],['100','250','500','1000','Sill', 'Fjord Interior'],\n",
    "              title=r'$Q_{sg}\\;[m^3/s]$',loc='best', fontsize=10)\n",
    "    \n",
    "    ax3.grid(linestyle='--', linewidth=0.5)\n",
    "    ax3.set_ylabel(r'$T_s\\;[^\\circ{C}]$',size = 12)\n",
    "    ax3.set_xlabel(r'$h_s/h$', size = 12)\n",
    "    ax3.set_title(r'(c)', loc='left', size=12, weight='bold')    \n",
    "    \n",
    "    \n",
    "    # Subplot of submarine melting Qsm\n",
    "    ax4.plot(hsr,Qsm,'o-',lw=2,label=str(Qsg[i]))\n",
    "    #ax4.set_ylim(0,4500)\n",
    "    ax4.legend(title=r'$Q_{sg}\\;[m^3/s]$',loc='best', fontsize=10)\n",
    "    ax4.grid(linestyle='--', linewidth=0.5)\n",
    "    ax4.set_ylabel(r'$Q_{sm}\\;[m^3/s]$',size = 12)\n",
    "    ax4.set_xlabel(r'$h_s/h$', size = 12)\n",
    "    ax4.set_title(r'(d)', loc='left', size=12, weight='bold')\n",
    "\n",
    "plt.subplots_adjust(left=0.1,\n",
    "                    bottom=0.1, \n",
    "                    right=0.9, \n",
    "                    top=0.9, \n",
    "                    wspace=0.25, \n",
    "                    hspace=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03d5ec19",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2_Plot for Stratification cases\n",
    "Strat = np.array([1, 2, 3])\n",
    "Sz = [r'$\\frac{\\partial s_0}{\\partial z}$', r'$2\\frac{\\partial s_0}{\\partial z}$', r'$4\\frac{\\partial s_0}{\\partial z}$']\n",
    "\n",
    "fig2 = plt.figure(figsize=(10,10))\n",
    "#fig1.tight_layout()\n",
    "#plt.grid(linestyle='--', linewidth=0.5)    \n",
    "ax1 = fig2.add_subplot(221)\n",
    "ax2 = fig2.add_subplot(222)\n",
    "ax3 = fig2.add_subplot(223)\n",
    "ax4 = fig2.add_subplot(224)\n",
    "\n",
    "\n",
    "from matplotlib.legend_handler import HandlerBase\n",
    "\n",
    "\n",
    "# class AnyObjectHandler(HandlerBase):\n",
    "#     def create_artists(self, legend, orig_handle,\n",
    "#                        x0, y0, width, height, fontsize, trans):\n",
    "#         l1 = plt.Line2D([x0,y0+width], [0.7*height,0.7*height],\n",
    "#                            linestyle=orig_handle[1],color=orig_handle[0])\n",
    "#         l2 = plt.Line2D([x0,y0+width], [0.3*height,0.3*height], \n",
    "#                            color=orig_handle[0])\n",
    "#         return [l1, l2]\n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "for i in range(len(Strat)):\n",
    "    \n",
    "    file1 = xr.open_dataset('/Volumes/Extreme SSD/MITgcm outputs/Strat_comp/Strat_' + str(format(Strat[i],'03d')) + '.nc')\n",
    "    hsr = file1.hsr.data\n",
    "    alp11 = file1.alpha11.data\n",
    "    alp11[alp11<0] = 0\n",
    "    alp00 = file1.alpha00.data\n",
    "    Ts_f = file1.Ts_fjord.data\n",
    "    Ts_s = file1.Ts_sill.data\n",
    "    Qsm = file1.Qsm.data   \n",
    "    Qe = file1.Qe.data\n",
    "\n",
    "    # Subplot of reflux coefficients\n",
    "    if Strat[i] == 1:\n",
    "        l1,l2 = ax1.plot(hsr,alp00,'o-', hsr,alp11,'^-', ms=7, alpha=0.7, color='#1f77b4',lw=2)\n",
    "    #ax1.set_ylim(0,.7)\n",
    "        \n",
    "    elif Strat[i] == 2:\n",
    "        \n",
    "        l3,l4 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#ff7f0e',lw=2)\n",
    "        \n",
    "    elif Strat[i] == 3:\n",
    "        \n",
    "        l5,l6 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#2ca02c',lw=2)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        l7,l8 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#d62728',lw=2)\n",
    "    \n",
    "#     ax1.legend((l1,l2,l3,l4,l5,l6,l7,l8),('100 $T_S$','100 $T_S2$','250 $T_S1$','250 $T_S2$','500 $T_S1$','500 $T_S2$','1000 $T_S1$','1000 $T_S2$'),\n",
    "#                title=r'$Q_{sg}\\;[m^3/s]$', loc='best', fontsize=12)\n",
    "\n",
    "\n",
    "\n",
    "#     #Plot analytic solution\n",
    "#     ax1.plot(hsr,alp11,'o-',lw=2,label=str(Qsg[i]))\n",
    "\n",
    "#     #Plot simulation\n",
    "#     ax1.plot(hsr,alp00,'^-',lw=2,label=str(Qsg[i]))\n",
    "\n",
    "    \n",
    "#     #Get artists and labels for legend and chose which ones to display\n",
    "#     handles, labels = ax1.get_legend_handles_labels()\n",
    "#     display = (0,1,2,3,4)\n",
    "\n",
    "#     #Create custom artists\n",
    "    \n",
    "#     simArtist = plt.Line2D((0,1),(0,0), color='k', marker='o', linestyle='')\n",
    "#     anyArtist = plt.Line2D((0,1),(0,0), color='k', marker='^', linestyle='')\n",
    "\n",
    "#     #Create legend from custom artist/label lists\n",
    "#     ax1.legend([handle for i,handle in enumerate(handles) if i in display]+[simArtist,anyArtist],\n",
    "#           [label for i,label in enumerate(labels) if i in display]+['Simulation', 'Analytic'], title=r'$Q_{sg}\\;[m^3/s]$',loc='best', fontsize=12)\n",
    "    \n",
    "    #line = Line2D([0,1],[0,1],linestyle='-', color='r')\n",
    "    line1 = plt.Line2D((0,1),(0,0), color='#1f77b4', lw=2, linestyle='-')\n",
    "    line2 = plt.Line2D((0,1),(0,0), color='#ff7f0e', lw=2, linestyle='-')\n",
    "    line3 = plt.Line2D((0,1),(0,0), color='#2ca02c', lw=2, linestyle='-')\n",
    "    line4 = plt.Line2D((0,1),(0,0), color='#d62728', lw=2, linestyle='-')\n",
    "    marker1 = plt.Line2D((0,1),(0,0), color='k', marker='o', linestyle='')\n",
    "    marker2 = plt.Line2D((0,1),(0,0), color='k', marker='^', linestyle='')    \n",
    "    \n",
    "    ax1.legend([line1,line2,line3,marker1,marker2],[r'$\\frac{\\partial s_0}{\\partial z}$',r'$2\\frac{\\partial s_0}{\\partial z}$',r'$4\\frac{\\partial s_0}{\\partial z}$',r'$\\alpha_{11}$', r'$\\alpha_{00}$'],\n",
    "              loc='best', fontsize=10)\n",
    "\n",
    "    ax1.grid(linestyle='--', linewidth=0.5)\n",
    "    ax1.set_ylabel(r'Reflux Coefficient',size = 12)\n",
    "    #ax1.xaxis.set_visible(False)\n",
    "    ax1.set_title(r'(a)', loc='left', size=12, weight='bold')\n",
    "\n",
    "    \n",
    "    # Subplot of entrainment flux Qe\n",
    "    ax2.plot(hsr,Qe,'o-',lw=2,label=Sz[i])\n",
    "    #ax2.set_ylim(0,4500)\n",
    "    ax2.legend(loc='best', fontsize=10)\n",
    "    ax2.grid(linestyle='--', linewidth=0.5)\n",
    "    ax2.set_ylabel(r'$Q_e\\;[m^3/s]$',size = 12)\n",
    "    #ax2.xaxis.set_visible(False)\n",
    "    ax2.set_title(r'(b)', loc='left', size=12, weight='bold')\n",
    "       \n",
    "    \n",
    "    # Subplot of deep-layer temperature Ts\n",
    "    if Strat[i] == 1:\n",
    "           \n",
    "        ax3.plot(hsr,Ts_s,'o-', hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#1f77b4',lw=2)\n",
    "    #ax1.set_ylim(0,.7)\n",
    "        \n",
    "    elif Strat[i] == 2:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#ff7f0e',lw=2)\n",
    "        \n",
    "    elif Strat[i] == 3:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#2ca02c',lw=2)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#d62728',lw=2)   \n",
    "    \n",
    "#     ax3.plot(hsr,Ts_s,'o-',lw=2,label=r'$Q_{sg}$')\n",
    "#     ax3.plot(hsr,Ts_f,'^-',lw=2,label=r'$Q_{sg}$')\n",
    "\n",
    "    #ax3.set_ylim(0,.7)\n",
    "    #ax3.legend(loc='best', fontsize=12)\n",
    "    ax3.legend([line1,line2,line3,marker1,marker2],[r'$\\frac{\\partial s_0}{\\partial z}$',r'$2\\frac{\\partial s_0}{\\partial z}$',r'$4\\frac{\\partial s_0}{\\partial z}$','Sill', 'Fjord Interior'],\n",
    "              loc='best', fontsize=10)\n",
    "    \n",
    "    ax3.grid(linestyle='--', linewidth=0.5)\n",
    "    ax3.set_ylabel(r'$T_s\\;[^\\circ{C}]$',size = 12)\n",
    "    ax3.set_xlabel(r'$h_s/h$', size = 12)\n",
    "    ax3.set_title(r'(c)', loc='left', size=12, weight='bold')    \n",
    "    \n",
    "    \n",
    "    # Subplot of submarine melting Qsm\n",
    "    ax4.plot(hsr,Qsm,'o-',lw=2,label=Sz[i])\n",
    "    #ax4.set_ylim(0,4500)\n",
    "    ax4.legend(loc='best', fontsize=10)\n",
    "    ax4.grid(linestyle='--', linewidth=0.5)\n",
    "    ax4.set_ylabel(r'$Q_{sm}\\;[m^3/s]$',size = 12)\n",
    "    ax4.set_xlabel(r'$h_s/h$', size = 12)\n",
    "    ax4.set_title(r'(d)', loc='left', size=12, weight='bold')\n",
    "\n",
    "plt.subplots_adjust(left=0.1,\n",
    "                    bottom=0.1, \n",
    "                    right=0.9, \n",
    "                    top=0.9, \n",
    "                    wspace=0.25, \n",
    "                    hspace=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e57ae83",
   "metadata": {},
   "outputs": [],
   "source": [
    "Sz[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a916c780",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3_Plot for Cd cases\n",
    "Id = np.array([1, 2, 3, 4])\n",
    "Cd = np.array([0e-3, 1.25e-3, 2.5e-3, 2.5e-2])\n",
    "\n",
    "fig3 = plt.figure(figsize=(10,10))\n",
    "#fig1.tight_layout()\n",
    "#plt.grid(linestyle='--', linewidth=0.5)    \n",
    "ax1 = fig3.add_subplot(221)\n",
    "ax2 = fig3.add_subplot(222)\n",
    "ax3 = fig3.add_subplot(223)\n",
    "ax4 = fig3.add_subplot(224)\n",
    "\n",
    "\n",
    "from matplotlib.legend_handler import HandlerBase\n",
    "\n",
    "\n",
    "# class AnyObjectHandler(HandlerBase):\n",
    "#     def create_artists(self, legend, orig_handle,\n",
    "#                        x0, y0, width, height, fontsize, trans):\n",
    "#         l1 = plt.Line2D([x0,y0+width], [0.7*height,0.7*height],\n",
    "#                            linestyle=orig_handle[1],color=orig_handle[0])\n",
    "#         l2 = plt.Line2D([x0,y0+width], [0.3*height,0.3*height], \n",
    "#                            color=orig_handle[0])\n",
    "#         return [l1, l2]\n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "for i in range(len(Cd)):\n",
    "    \n",
    "    file1 = xr.open_dataset('/Volumes/Extreme SSD/MITgcm outputs/Cd_comp/Cd_' + str(format(Id[i],'03d')) + '.nc')\n",
    "    hsr = file1.hsr.data\n",
    "    alp11 = file1.alpha11.data\n",
    "    alp11[alp11<0] = 0\n",
    "    alp00 = file1.alpha00.data\n",
    "    Ts_f = file1.Ts_fjord.data\n",
    "    Ts_s = file1.Ts_sill.data\n",
    "    Qsm = file1.Qsm.data   \n",
    "    Qe = file1.Qe.data\n",
    "\n",
    "    # Subplot of reflux coefficients\n",
    "    if Id[i] == 1:\n",
    "        l1,l2 = ax1.plot(hsr,alp00,'o-', hsr,alp11,'^-', ms=7, alpha=0.7, color='#1f77b4',lw=2)\n",
    "    #ax1.set_ylim(0,.7)\n",
    "        \n",
    "    elif Id[i] == 2:\n",
    "        \n",
    "        l3,l4 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#ff7f0e',lw=2)\n",
    "        \n",
    "    elif Id[i] == 3:\n",
    "        \n",
    "        l5,l6 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#2ca02c',lw=2)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        l7,l8 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#d62728',lw=2)\n",
    "    \n",
    "#     ax1.legend((l1,l2,l3,l4,l5,l6,l7,l8),('100 $T_S$','100 $T_S2$','250 $T_S1$','250 $T_S2$','500 $T_S1$','500 $T_S2$','1000 $T_S1$','1000 $T_S2$'),\n",
    "#                title=r'$Q_{sg}\\;[m^3/s]$', loc='best', fontsize=12)\n",
    "\n",
    "\n",
    "\n",
    "#     #Plot analytic solution\n",
    "#     ax1.plot(hsr,alp11,'o-',lw=2,label=str(Qsg[i]))\n",
    "\n",
    "#     #Plot simulation\n",
    "#     ax1.plot(hsr,alp00,'^-',lw=2,label=str(Qsg[i]))\n",
    "\n",
    "    \n",
    "#     #Get artists and labels for legend and chose which ones to display\n",
    "#     handles, labels = ax1.get_legend_handles_labels()\n",
    "#     display = (0,1,2,3,4)\n",
    "\n",
    "#     #Create custom artists\n",
    "    \n",
    "#     simArtist = plt.Line2D((0,1),(0,0), color='k', marker='o', linestyle='')\n",
    "#     anyArtist = plt.Line2D((0,1),(0,0), color='k', marker='^', linestyle='')\n",
    "\n",
    "#     #Create legend from custom artist/label lists\n",
    "#     ax1.legend([handle for i,handle in enumerate(handles) if i in display]+[simArtist,anyArtist],\n",
    "#           [label for i,label in enumerate(labels) if i in display]+['Simulation', 'Analytic'], title=r'$Q_{sg}\\;[m^3/s]$',loc='best', fontsize=12)\n",
    "    \n",
    "    #line = Line2D([0,1],[0,1],linestyle='-', color='r')\n",
    "    line1 = plt.Line2D((0,1),(0,0), color='#1f77b4', lw=2, linestyle='-')\n",
    "    line2 = plt.Line2D((0,1),(0,0), color='#ff7f0e', lw=2, linestyle='-')\n",
    "    line3 = plt.Line2D((0,1),(0,0), color='#2ca02c', lw=2, linestyle='-')\n",
    "    line4 = plt.Line2D((0,1),(0,0), color='#d62728', lw=2, linestyle='-')\n",
    "    marker1 = plt.Line2D((0,1),(0,0), color='k', marker='o', linestyle='')\n",
    "    marker2 = plt.Line2D((0,1),(0,0), color='k', marker='^', linestyle='')    \n",
    "    \n",
    "    ax1.legend([line1,line2,line3,line4,marker1,marker2],[str(Cd[0]),str(Cd[1]),str(Cd[2]),str(Cd[3]),r'$\\alpha_{11}$', r'$\\alpha_{00}$'],\n",
    "              title=r'$C_d$',loc='best', fontsize=10)\n",
    "\n",
    "    ax1.grid(linestyle='--', linewidth=0.5)\n",
    "    ax1.set_ylabel(r'Reflux Coefficient',size = 12)\n",
    "    #ax1.xaxis.set_visible(False)\n",
    "    ax1.set_title(r'(a)', loc='left', size=12, weight='bold')\n",
    "\n",
    "    \n",
    "    # Subplot of entrainment flux Qe\n",
    "    ax2.plot(hsr,Qe,'o-',lw=2,label=str(Cd[i]))\n",
    "    #ax2.set_ylim(0,4500)\n",
    "    ax2.legend(title=r'$C_d$',loc='best', fontsize=10)\n",
    "    ax2.grid(linestyle='--', linewidth=0.5)\n",
    "    ax2.set_ylabel(r'$Q_e\\;[m^3/s]$',size = 12)\n",
    "    #ax2.xaxis.set_visible(False)\n",
    "    ax2.set_title(r'(b)', loc='left', size=12, weight='bold')\n",
    "       \n",
    "    \n",
    "    # Subplot of deep-layer temperature Ts\n",
    "    if Id[i] == 1:\n",
    "           \n",
    "        ax3.plot(hsr,Ts_s,'o-', hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#1f77b4',lw=2)\n",
    "    #ax1.set_ylim(0,.7)\n",
    "        \n",
    "    elif Id[i] == 2:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#ff7f0e',lw=2)\n",
    "        \n",
    "    elif Id[i] == 3:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#2ca02c',lw=2)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#d62728',lw=2)   \n",
    "    \n",
    "#     ax3.plot(hsr,Ts_s,'o-',lw=2,label=r'$Q_{sg}$')\n",
    "#     ax3.plot(hsr,Ts_f,'^-',lw=2,label=r'$Q_{sg}$')\n",
    "\n",
    "    #ax3.set_ylim(0,.7)\n",
    "    #ax3.legend(loc='best', fontsize=12)\n",
    "    ax3.legend([line1,line2,line3,line4,marker1,marker2],[str(Cd[0]),str(Cd[1]),str(Cd[2]),str(Cd[3]),'Sill', 'Fjord Interior'],\n",
    "              title=r'$C_d$',loc='best', fontsize=10)\n",
    "    \n",
    "    ax3.grid(linestyle='--', linewidth=0.5)\n",
    "    ax3.set_ylabel(r'$T_s\\;[^\\circ{C}]$',size = 12)\n",
    "    ax3.set_xlabel(r'$h_s/h$', size = 12)\n",
    "    ax3.set_title(r'(c)', loc='left', size=12, weight='bold')    \n",
    "    \n",
    "    \n",
    "    # Subplot of submarine melting Qsm\n",
    "    ax4.plot(hsr,Qsm,'o-',lw=2,label=str(Cd[i]))\n",
    "    #ax4.set_ylim(0,4500)\n",
    "    ax4.legend(title=r'$C_d$',loc='best', fontsize=10)\n",
    "    ax4.grid(linestyle='--', linewidth=0.5)\n",
    "    ax4.set_ylabel(r'$Q_{sm}\\;[m^3/s]$',size = 12)\n",
    "    ax4.set_xlabel(r'$h_s/h$', size = 12)\n",
    "    ax4.set_title(r'(d)', loc='left', size=12, weight='bold')\n",
    "\n",
    "plt.subplots_adjust(left=0.1,\n",
    "                    bottom=0.1, \n",
    "                    right=0.9, \n",
    "                    top=0.9, \n",
    "                    wspace=0.25, \n",
    "                    hspace=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03969cc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4_Plot for Tide cases\n",
    "Id = np.array([0, 1, 3])\n",
    "Ue = 0.20 # Exchange velocity\n",
    "U0 = np.array([0, 0.01, 0.04]) # Tidal amplitude\n",
    "Ur = ['0', '1/20', '1/5']\n",
    "\n",
    "fig4 = plt.figure(figsize=(10,10))\n",
    "#fig1.tight_layout()\n",
    "#plt.grid(linestyle='--', linewidth=0.5)    \n",
    "ax1 = fig4.add_subplot(221)\n",
    "ax2 = fig4.add_subplot(222)\n",
    "ax3 = fig4.add_subplot(223)\n",
    "ax4 = fig4.add_subplot(224)\n",
    "\n",
    "\n",
    "from matplotlib.legend_handler import HandlerBase\n",
    "\n",
    "\n",
    "# class AnyObjectHandler(HandlerBase):\n",
    "#     def create_artists(self, legend, orig_handle,\n",
    "#                        x0, y0, width, height, fontsize, trans):\n",
    "#         l1 = plt.Line2D([x0,y0+width], [0.7*height,0.7*height],\n",
    "#                            linestyle=orig_handle[1],color=orig_handle[0])\n",
    "#         l2 = plt.Line2D([x0,y0+width], [0.3*height,0.3*height], \n",
    "#                            color=orig_handle[0])\n",
    "#         return [l1, l2]\n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "for i in range(len(Ur)):\n",
    "    \n",
    "    file1 = xr.open_dataset('/Volumes/Extreme SSD/MITgcm outputs/Tide_comp/Tide_' + str(format(Id[i],'03d')) + '.nc')\n",
    "    hsr = file1.hsr.data\n",
    "    alp11 = file1.alpha11.data\n",
    "    alp11[alp11<0] = 0\n",
    "    alp00 = file1.alpha00.data\n",
    "    Ts_f = file1.Ts_fjord.data\n",
    "    Ts_s = file1.Ts_sill.data\n",
    "    Qsm = file1.Qsm.data   \n",
    "    Qe = file1.Qe.data\n",
    "\n",
    "    # Subplot of reflux coefficients\n",
    "    if Id[i] == 0:\n",
    "        l1,l2 = ax1.plot(hsr,alp00,'o-', hsr,alp11,'^-', ms=7, alpha=0.7, color='#1f77b4',lw=2)\n",
    "    #ax1.set_ylim(0,.7)\n",
    "        \n",
    "    elif Id[i] == 1:\n",
    "        \n",
    "        l3,l4 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#ff7f0e',lw=2)\n",
    "        \n",
    "    elif Id[i] == 3:\n",
    "        \n",
    "        l5,l6 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#2ca02c',lw=2)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        l7,l8 = ax1.plot(hsr,alp00,'o-',hsr,alp11,'^-', ms=7, alpha=0.7, color='#d62728',lw=2)\n",
    "    \n",
    "#     ax1.legend((l1,l2,l3,l4,l5,l6,l7,l8),('100 $T_S$','100 $T_S2$','250 $T_S1$','250 $T_S2$','500 $T_S1$','500 $T_S2$','1000 $T_S1$','1000 $T_S2$'),\n",
    "#                title=r'$Q_{sg}\\;[m^3/s]$', loc='best', fontsize=12)\n",
    "\n",
    "\n",
    "\n",
    "#     #Plot analytic solution\n",
    "#     ax1.plot(hsr,alp11,'o-',lw=2,label=str(Qsg[i]))\n",
    "\n",
    "#     #Plot simulation\n",
    "#     ax1.plot(hsr,alp00,'^-',lw=2,label=str(Qsg[i]))\n",
    "\n",
    "    \n",
    "#     #Get artists and labels for legend and chose which ones to display\n",
    "#     handles, labels = ax1.get_legend_handles_labels()\n",
    "#     display = (0,1,2,3,4)\n",
    "\n",
    "#     #Create custom artists\n",
    "    \n",
    "#     simArtist = plt.Line2D((0,1),(0,0), color='k', marker='o', linestyle='')\n",
    "#     anyArtist = plt.Line2D((0,1),(0,0), color='k', marker='^', linestyle='')\n",
    "\n",
    "#     #Create legend from custom artist/label lists\n",
    "#     ax1.legend([handle for i,handle in enumerate(handles) if i in display]+[simArtist,anyArtist],\n",
    "#           [label for i,label in enumerate(labels) if i in display]+['Simulation', 'Analytic'], title=r'$Q_{sg}\\;[m^3/s]$',loc='best', fontsize=12)\n",
    "    \n",
    "    #line = Line2D([0,1],[0,1],linestyle='-', color='r')\n",
    "    line1 = plt.Line2D((0,1),(0,0), color='#1f77b4', lw=2, linestyle='-')\n",
    "    line2 = plt.Line2D((0,1),(0,0), color='#ff7f0e', lw=2, linestyle='-')\n",
    "    line3 = plt.Line2D((0,1),(0,0), color='#2ca02c', lw=2, linestyle='-')\n",
    "    line4 = plt.Line2D((0,1),(0,0), color='#d62728', lw=2, linestyle='-')\n",
    "    marker1 = plt.Line2D((0,1),(0,0), color='k', marker='o', linestyle='')\n",
    "    marker2 = plt.Line2D((0,1),(0,0), color='k', marker='^', linestyle='')    \n",
    "    \n",
    "    ax1.legend([line1,line2,line3,marker1,marker2],[str(Ur[0]),str(Ur[1]),str(Ur[2]),r'$\\alpha_{11}$', r'$\\alpha_{00}$'],\n",
    "              title=r'$U_0/U_e$',loc='best', fontsize=10)\n",
    "\n",
    "    ax1.grid(linestyle='--', linewidth=0.5)\n",
    "    ax1.set_ylabel(r'Reflux Coefficient',size = 12)\n",
    "    #ax1.xaxis.set_visible(False)\n",
    "    ax1.set_title(r'(a)', loc='left', size=12, weight='bold')\n",
    "\n",
    "    \n",
    "    # Subplot of entrainment flux Qe\n",
    "    ax2.plot(hsr,Qe,'o-',lw=2,label=str(Ur[i]))\n",
    "    #ax2.set_ylim(0,4500)\n",
    "    ax2.legend(title=r'$U_0/U_e$',loc='best', fontsize=10)\n",
    "    ax2.grid(linestyle='--', linewidth=0.5)\n",
    "    ax2.set_ylabel(r'$Q_e\\;[m^3/s]$',size = 12)\n",
    "    #ax2.xaxis.set_visible(False)\n",
    "    ax2.set_title(r'(b)', loc='left', size=12, weight='bold')\n",
    "       \n",
    "    \n",
    "    # Subplot of deep-layer temperature Ts\n",
    "    if Id[i] == 0:\n",
    "           \n",
    "        ax3.plot(hsr,Ts_s,'o-', hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#1f77b4',lw=2)\n",
    "    #ax1.set_ylim(0,.7)\n",
    "        \n",
    "    elif Id[i] == 1:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#ff7f0e',lw=2)\n",
    "        \n",
    "    elif Id[i] == 3:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#2ca02c',lw=2)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        ax3.plot(hsr,Ts_s,'o-',hsr,Ts_f,'^-', ms=7, alpha=0.7, color='#d62728',lw=2)   \n",
    "    \n",
    "#     ax3.plot(hsr,Ts_s,'o-',lw=2,label=r'$Q_{sg}$')\n",
    "#     ax3.plot(hsr,Ts_f,'^-',lw=2,label=r'$Q_{sg}$')\n",
    "\n",
    "    #ax3.set_ylim(0,.7)\n",
    "    #ax3.legend(loc='best', fontsize=12)\n",
    "    ax3.legend([line1,line2,line3,marker1,marker2],[str(Ur[0]),str(Ur[1]),str(Ur[2]),'Sill', 'Fjord Interior'],\n",
    "              title=r'$U_0/U_e$',loc='best', fontsize=10)\n",
    "    \n",
    "    ax3.grid(linestyle='--', linewidth=0.5)\n",
    "    ax3.set_ylabel(r'$T_s\\;[^\\circ{C}]$',size = 12)\n",
    "    ax3.set_xlabel(r'$h_s/h$', size = 12)\n",
    "    ax3.set_title(r'(c)', loc='left', size=12, weight='bold')    \n",
    "    \n",
    "    \n",
    "    # Subplot of submarine melting Qsm\n",
    "    ax4.plot(hsr,Qsm,'o-',lw=2,label=str(Ur[i]))\n",
    "    #ax4.set_ylim(0,4500)\n",
    "    ax4.legend(title=r'$U_0/U_e$',loc='best', fontsize=10)\n",
    "    ax4.grid(linestyle='--', linewidth=0.5)\n",
    "    ax4.set_ylabel(r'$Q_{sm}\\;[m^3/s]$',size = 12)\n",
    "    ax4.set_xlabel(r'$h_s/h$', size = 12)\n",
    "    ax4.set_title(r'(d)', loc='left', size=12, weight='bold')\n",
    "\n",
    "plt.subplots_adjust(left=0.1,\n",
    "                    bottom=0.1, \n",
    "                    right=0.9, \n",
    "                    top=0.9, \n",
    "                    wspace=0.25, \n",
    "                    hspace=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e030826",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Ts_f\n",
    "from fractions import Fraction\n",
    "from sympy import Rational\n",
    "A = Fraction(1,3)\n",
    "B = Fraction(1.5,3)\n",
    "print(A,B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c0e8857",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_path = \"/Users/weiyangbao/Documents/Results_Figs/Comp_cases/\"\n",
    "if not os.path.exists(fig_path):\n",
    "    os.makedirs(fig_path)\n",
    "\n",
    "#fig1.savefig(fig_path + '/Qsg_comp.png', dpi=200, bbox_inches='tight')\n",
    "#fig2.savefig(fig_path + '/Strat_comp.png', dpi=200, bbox_inches='tight')\n",
    "#fig3.savefig(fig_path + '/Cd_comp.png', dpi=200, bbox_inches='tight')\n",
    "fig4.savefig(fig_path + '/Tide_comp.png', dpi=200, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0d88d90",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}