{
 "metadata": {
  "name": "Lesson 1 - Simple ocean mixed layer"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Lesson 1 - Ocean mixed layer and radiative forcing without sea ice and atmosphere"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2>Scenario 1</h2>\n",
      "\n",
      " * Ocean mixed layer forced by shortwave radiation only\n",
      " * No atmosphere\n",
      " * No exchange with deeper ocean layers, immediate mixing\n",
      " * Heat balance at the sea surface: Short wave incoming radiation + long wave outgoing radiation\n",
      "\n",
      "<h2>Model equations</h2>\n",
      "    \n",
      "<h3>Surface energy balance</h3>\n",
      "$$ Q^\\downarrow + Q_{LW}^{\\uparrow} = Q_{srf} $$\n",
      "\n",
      " * $Q_{srf}$ energy balance at the surface\n",
      " * $Q^\\downarrow$ incoming energy flux\n",
      " * $Q_{LW}^{\\uparrow}$ outgoing longwave radiation\n",
      " \n",
      "<h3>Differential equation for the prognostic variable mixed layer temperature</h3>\n",
      "\n",
      "$\\frac{\\partial T_{ml}} {\\partial t} \\rho_w c_w h_{ml}=Q_{srf}$\n",
      "\n",
      " * $T_{ml}$ mixed layer temperature\n",
      " * $t$ time \n",
      " * $\\rho_w$ density of sea water\n",
      " * $c_w$ specific heat capacity of sea water \n",
      "\n",
      "<h3>Numerics</h3>\n",
      "forward-in-time integration, finite differences derived from Taylor series expansion\n",
      "\n",
      "$$f(t+\\Delta t)=f(t)+\\frac{\\Delta t}{1!} f'(t)+O(\\Delta t^2)$$\n",
      "\n",
      "\n",
      "\n",
      "<h2>Research questions</h2>\n",
      "\n",
      "Compute the time evolution of the ocean mixed layer temperature $T_{ml}(t)$ for different \n",
      "\n",
      "* $h_{ml}$\n",
      "* initial temperatures $T_{ml}(t=0)$\n",
      "* incoming energy fluxed $Q^{\\downarrow}$\n",
      "\n",
      "Estimate the typical time scale for stationarity and select appropriate time step $\\Delta t$ for model integration."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Main equations\n",
      "\n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We rearrange the Taylor series\n",
      "$$f(t+\\Delta t)=f(t)+\\frac{\\Delta t}{1!} f'(t)+O(\\Delta t^2)$$\n",
      "to obtain the time derivative of the function f\n",
      "$$f'(t)=\\frac{f(t+\\Delta t)-f(t)}{\\Delta t}$$\n",
      "whereas the truncation error was neglected. \n",
      "\n",
      "We identify the mixed layer temperature $T_{ml}$ as the function $f$. Thus it follows for the differential equation\n",
      "\n",
      "$$\\frac{\\partial T_{ml}} {\\partial t} \\rho_w c_w h_{ml}=Q_{srf}$$\n",
      "\n",
      "$$\\frac{ T_{ml}(t+\\Delta t)-T_{ml}(t)}{\\Delta t}\\rho_w c_w h_{ml}=Q_{srf}$$\n",
      "\n",
      "$$\\frac{ T_{ml}(t+\\Delta t)-T_{ml}(t)}{\\Delta t}\\rho_w c_w h_{ml}=Q^\\downarrow + Q_{LW}^{\\uparrow}$$\n",
      "\n",
      "The following equation allows to calculate the temperature $T_{ml}(t+\\Delta t)$ from the previous time step $T_{ml}(t)$:\n",
      "\n",
      "$$T_{ml}(t+\\Delta t)=\\frac{\\Delta t}{\\rho_w c_w h_{ml}} (Q^\\downarrow + Q_{LW}^{\\uparrow})+T_{ml}(t)$$\n",
      "\n",
      "Thus, it is straightforward to do the integration within one time loop.\n",
      "\n",
      "\n",
      "\n",
      "\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Model implementation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "<h2>Constants and functions</h2>\n",
      "<h3>Sea water</h3>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import seawater\n",
      "S=34 # salinityy [PSU}\n",
      "T=0.0 # Temperature [deg C]\n",
      "P=0 # Pressure in [dbar]\n",
      "rho0=seawater.csiro.dens(S,T,P) # Density of sea water\n",
      "cp=seawater.csiro.cp(S,T,P) # Heat capacity of sea water\n",
      "\n",
      "print 'Density [kg/m^3]', rho0\n",
      "print 'Heat capacity [J/kg/K]',cp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " Density [kg/m^3] 1027.29893846\n",
        "Heat capacity [J/kg/K] 3992.61688157\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h3>Radiation</h3>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.constants\n",
      "\n",
      "def outgoing(T,epsilon=0.97):\n",
      "    \"\"\"\"Calculate outgoing radiation from Stefan-Boltzmann law\n",
      "        \n",
      "        Input parameters:\n",
      "        T temperature in K\n",
      "        epsilon infrared emissvity, default 0.97\n",
      "\n",
      "        return radiation in W/m^2\n",
      "        \"\"\"\n",
      "    return epsilon*scipy.constants.sigma*T**4"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 24
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<h2>Variables</h2>\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Important variables (to play around with)\n",
      "h_ml=10 # [m] mixed layer depth\n",
      "T_ml_0=20 # [deg C] initial temperature\n",
      "dt=60*60*24*10 # [s] time step, 10 days\n",
      "\n",
      "# Not so important variables (treated as constants)\n",
      "import seawater\n",
      "rho0=seawater.csiro.dens(34.0,T_ml_0,0) # Density of sea water, salinity 34 PSU, pressur 0 dbar\n",
      "cp=seawater.csiro.cp(34.0,T_ml_0,0) # Heat capacity of sea water\n",
      "epsilon=0.97 # Infrared emissivity\n",
      "\n",
      "# Constants\n",
      "import scipy.constants\n",
      "sig=scipy.constants.sigma #Stefan-Boltzmann constant 5.670373e-08 W m^-2 K^-4"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    }
   ],
   "metadata": {}
  }
 ]
}