LayerThFlux.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Heat budget for the different\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import glob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xp1 = 51 # Number of cell corner points\n",
    "# Find the uppler and lower layer depth above the sill crest\n",
    "def Layer(u):\n",
    "    # Layer interface, zero velocity crossing point\n",
    "    z_th = np.empty(Xp1)\n",
    "    itf = np.empty(Xp1, dtype=int)\n",
    "    z_th[[0,1,-2,-1]] = 0\n",
    "    itf[[0,1,-2,-1]] = 0\n",
    "    for i in range(2,Xp1-2):\n",
    "        l = np.where(u[:, i] < 0)[-1][0]\n",
    "        itf[i] = l\n",
    "        z_th[i] = -w[l, i] * (pt[l-1:l+1,i].mean()-T0)\n",
    "    return (itf, z_th)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Cp = 4220  # Specific heat content J/kg/C\n",
    "rho = 1025 # Density\n",
    "T0 = 9 # Initial temperature\n",
    "ny = 9 # Number of grids in Y direction\n",
    "ga = np.repeat(np.array([[.5],[1],[1.5],[4],[8]]),10) * 200 # Grid area\n",
    "\n",
    "\n",
    "files = sorted(glob.glob('/Users/weiyangbao/Documents/Model_outputs/State700/Hs*.nc'))\n",
    "Hs = np.array([10, 12, 15, 110])\n",
    "\n",
    "H_f = np.empty(len(Hs)) # Heat flux through the ice front\n",
    "H_f1 = np.empty(len(Hs)) # 1: upper layer\n",
    "H_s = np.empty(len(Hs)) # Heat flux through the sill\n",
    "H_s1 = np.empty(len(Hs))\n",
    "H_z = np.empty(len(Hs)) # Vertical heat flux\n",
    "for i in range(len(Hs)):\n",
    "    state0 = xr.open_dataset(files[i])\n",
    "    t0 = 0 # Start from day 20\n",
    "    tn = len(state0.T) # End day\n",
    "    state = state0.isel(Y = range(1,10), T = range(t0,tn)) # Remove the boundary walls in y direction\n",
    "    pt = state.Temp.data.mean(axis=(0,2)) # Z, X\n",
    "    w = state.W.data.mean(axis=(0,2)) # Zl, X\n",
    "    u = state.U.data.mean(axis=(0,2)) \n",
    "    itf, z_th = Layer(u)\n",
    "    # sill crest: itf[27]\n",
    "    H_z[i] = z_th[2:29].sum() * 40000 * ny\n",
    "    H_f[i] = np.sum(u[:,2] * (pt[:,1]-T0) * ga) * ny\n",
    "    H_f1[i] = np.sum(u[:itf[2],2] * (pt[:itf[2],1]-T0) * ga[:itf[2]]) * ny\n",
    "    H_s[i] = -np.sum(u[:,28] * (pt[:,27]-T0) * ga) * ny\n",
    "    H_s1[i] = -np.sum(u[:itf[28],28] * (pt[:itf[28],27:29].mean(1)-T0) * ga[:itf[28]]) * ny\n",
    "itf[28]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    H_z[i] = z_th[2:29].sum() * 40000 * ny\n",
    "    Qf = u[:,2] * (pt[:,1]-T0) * Area * ny\n",
    "    H_f[i] = Qf.sum()\n",
    "    #H_f2[i] = Qf[np.int(itf[2]):].sum()\n",
    "    H_f2[i] = Qf[Qf>0].sum()\n",
    "    Qs = -u[:,28] * (pt[:,27]-T0) * Area * ny\n",
    "    H_s[i] = Qs.sum()\n",
    "    H_s2[i] = Qs[np.int(itf[28]):].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-4344.30991727, -2736.43736357,  -709.57352969,   408.75530743]),\n",
       " array([-5487.81237323, -5561.86994947, -5668.96521575, -5631.00153356]),\n",
       " array([3527.43614763, 3920.319888  , 4792.94689894, 6290.99725052]),\n",
       " array([-1960.3762256 , -1641.55006148,  -876.01831681,   659.99571696]))"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_z,H_f, H_s, H_f+H_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fig, ax =plt.subplots(figsize=(10,6))\n",
    "fig, (ax1, ax2) = plt.subplots(1,2,sharex=True, figsize=(16,16))\n",
    "labels = [10, 12, 15]\n",
    "#ax1.set_xticklabels(labels)\n",
    "ax1.grid(linestyle='--')\n",
    "plt.yticks(fontsize = 15)\n",
    "plt.xticks(np.arange(3), labels, fontsize=15)\n",
    "plt.xlabel('$H_s\\;(m)$', size=20)\n",
    "plt.ylabel('Heat flux ($\\circ C\\cdot m^3/s$)', size=20)\n",
    "plt.title('Hf', size=20)\n",
    "ax1.plot(np.arange(3), H_f, 'o-', lw=2)\n",
    "\n",
    "ax2.grid(linestyle='--')\n",
    "plt.yticks(fontsize = 15)\n",
    "plt.xticks(np.arange(3), labels, fontsize=15)\n",
    "plt.xlabel('$H_s\\;(m)$', size=20)\n",
    "#ax2.ylabel('Heat flux ($\\circ C\\cdot m^3/s$)', size=20)\n",
    "plt.title('Hs', size=20)\n",
    "ax2.plot(np.arange(3), H_s, 'o-', lw=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}