{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a491730b",
   "metadata": {},
   "source": [
    "# Energy Balance Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26fa988d",
   "metadata": {},
   "source": [
    "This code is based on that created by Dr. David Archer and adapted by Tatsam Garg for TROP ICSU, https://tropicsu.org/lesson-plan-create-climate-model-with-py/."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f858c79c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#This imports the desired packages.\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11d3654f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Setting up the main parameters\n",
    "#Define your time interval size\n",
    "timeStep = 10                      # years\n",
    "waterDepth = 4000                  # meters\n",
    "L = 1350                           # Watts/m^2 (Solar,Constant)\n",
    "albedo = 0.3\n",
    "epsilon = 1\n",
    "sigma = 5.67*10**(-8)              # The Stefan-Boltzmann constant in units of W/m2 K^4\n",
    "SpecificHeatCapacity = 4200        # (J/kg K) (for water)\n",
    "Density = 1000                     # (kg/m^3) (for water)\n",
    "Volume = waterDepth                # m^3 (of a water column with unit cross sectional area)\n",
    "WaterColumnMass = Density * Volume #(kg/m^2)\n",
    "HeatCapacity = SpecificHeatCapacity * WaterColumnMass #(J/m^2 K)\n",
    "print(\"Heat capacity is : \", HeatCapacity, \"J/m^2 K\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01d67852",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define the number of timesteps you want to take.\n",
    "nSteps = 100\n",
    "#Create empty arrays with as many spaces as the number of time steps you want to take.\n",
    "#The np.zeros(size) command creates empty arrays\n",
    "Heatrate = np.zeros(nSteps) #Heat content per unit time (as defined in equation 3)\n",
    "HeatContent = np.zeros(nSteps) #Total heat content at each time step\n",
    "T = np.zeros(nSteps) #Temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe4c5789",
   "metadata": {},
   "outputs": [],
   "source": [
    "T[0] = 0 #Initial Temperature\n",
    "HeatContent[0] = HeatCapacity * T[0] #Initial heat content (Heat content and Temperature are related by the heat capacity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4583c55d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define a linear array whose elements are consecutive integers and range is the same as the number of time steps.\n",
    "t = np.linspace(0, nSteps, nSteps)\n",
    "#Now the time elapsed can simply be:\n",
    "Time = t*timeStep # Time elapsed is nothing but integral multiples of the timeStep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3abb4fe",
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# A 'for' loop that iteratively determines different quantities of interest corresponding\n",
    "# to each time step using the previous step's data.\n",
    "for t in range (1, nSteps): # 't' has become our index\n",
    "    Heatrate[t] = L*(1-albedo)/4 - epsilon * sigma * T[t-1]**4\n",
    "    HeatContent[t] = HeatContent[t-1] + Heatrate[t] * 24*365*3600*timeStep\n",
    "    T[t] = HeatContent[t] / HeatCapacity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c522d2e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Time, T)\n",
    "plt.rcParams['figure.figsize']=[14,10] # Function to control Figure size\n",
    "plt.title(\"Temperature change over time\", fontsize = '20')\n",
    "plt.ylabel(\"Temperature (Kelvins)\", fontsize='15')\n",
    "#plt.ylabel(\"Temperature (Degrees Celcius)\", fontsize='15') \n",
    "plt.xlabel(\"Time (Years)\", fontsize='15')\n",
    "plt.show()\n"
   ]
  }
 ],
 "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.9.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}