ocean_only/global_ALE/z/Surface flux analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MOM6 diagnostic boundary fluxes of scalars and their global budgets \n",
    "\n",
    "Results from this notebook: \n",
    "1. Maps of surface boundary fluxes of water, heat, and salt crossing into the liquid seawater in MOM6;\n",
    "2. Computation of self-consistency checks, including global heat, salt, and mass budgets to verify that the model is conserving scalars over the global domain.   \n",
    "\n",
    "Caveats regarding this notebook:\n",
    "1. This notebook is written for the MOM6-examples/ocean_only/global_ALE/z and \n",
    "   MOM6-examples/ocean_only/global_ALE/layer test cases.  \n",
    "   It is nearly the same as the notebook for MOM6-examples/ice_ocean_SIS/GOLD_SIS.\n",
    "2. It only considers tendencies over 12 hours.  \n",
    "\n",
    "Hopes for the use of this notebook: \n",
    "1. To provide a starting point to document boundary fluxes of scalar fields;\n",
    "2. To teach MOM6 users about the boundary fluxes, their patterns, units, and sign conventions;\n",
    "3. To perform self-consistency checks to ensure the model is conserving scalar fields;\n",
    "4. To illustrate a self-contained iPython notebook of use for MOM6 analysis.  \n",
    "\n",
    "This jupyter notebook was originally developed at NOAA/GFDL, and it is provided freely to the MOM6 community. GFDL scientists developing MOM6 make extensive use of Python for diagnostics. We solicit modifications/fixes that are useful to the MOM6 community."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import netCDF4\n",
    "import numpy\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path    = './'\n",
    "static  = netCDF4.Dataset(path+'ocean_static.nc')\n",
    "forcing = netCDF4.Dataset(path+'forcing.nc')\n",
    "surface = netCDF4.Dataset(path+'surface.nc'); tvar = 'Time'\n",
    "cp   = 3925.0\n",
    "rho0 = 1035.0\n",
    "n    = 1\n",
    "\n",
    "print('Forcing: ',[v for v in forcing.variables])\n",
    "print('Surface: ',[v for v in surface.variables])\n",
    "print('Static: ',[v for v in static.variables])    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This section fills the fields used in this notebook.\n",
    "\n",
    "#--------------------------------------------------------------\n",
    "# geometric factors \n",
    "lon  = static.variables['geolon'][:]\n",
    "lat  = static.variables['geolat'][:]\n",
    "wet  = static.variables['wet'][:]\n",
    "area = static.variables['area_t'][:]*wet\n",
    "\n",
    "#--------------------------------------------------------------\n",
    "# time in days \n",
    "time = surface.variables[tvar][:]*86400.\n",
    "\n",
    "#--------------------------------------------------------------\n",
    "# sea surface temperature \n",
    "sst = surface.variables['SST'][n]\n",
    "\n",
    "#--------------------------------------------------------------\n",
    "# \\int \\rho \\, dz \\, (1,Theta,S), with \\rho = \\rho_{0} for Bousssinesq.  \n",
    "mass_wt = surface.variables['mass_wt']\n",
    "tomint  = surface.variables['opottempmint']\n",
    "somint  = surface.variables['somint']\n",
    "\n",
    "\n",
    "#--------------------------------------------------------------\n",
    "# mass flux of water crossing ocean surface [kg/(m^2 s)]\n",
    "# positive values indicate mass entering ocean; \n",
    "# negative values indicate mass leaving ocean. \n",
    "\n",
    "# net mass flux entering ocean \n",
    "net_massin = forcing.variables['net_massin'][n]\n",
    "\n",
    "# net mass flux leaving ocean\n",
    "net_massout = forcing.variables['net_massout'][n]\n",
    "\n",
    "# evaporation (negative) and condensation (positive)\n",
    "evap = forcing.variables['evap'][n]\n",
    "\n",
    "# liquid runoff entering ocean (non-negative)\n",
    "lrunoff = forcing.variables['lrunoff'][n]\n",
    "\n",
    "# frozen runoff entering ocean (non-negative)\n",
    "frunoff = forcing.variables['frunoff'][n]\n",
    "\n",
    "# liquid precipitation entering ocean.\n",
    "# note: includes exchanges with sea-ice, with \n",
    "#       melt adding mass to ocean; formation removing mass.\n",
    "lprec = forcing.variables['lprec'][n]\n",
    "\n",
    "# frozen precipitation entering ocean.\n",
    "fprec = forcing.variables['fprec'][n]\n",
    "\n",
    "# virtual precipitation arising from conversion of salt restoring to water flux\n",
    "vprec = forcing.variables['vprec'][n]\n",
    "\n",
    "# net mass flux crossing surface (including exchange with sea-ice)\n",
    "PRCmE = forcing.variables['PRCmE'][n]\n",
    "\n",
    "#--------------------------------------------------------------\n",
    "# heat flux crossing ocean surface and bottom [Watt/m^2]\n",
    "# positive values indicate heat entering ocean;\n",
    "# negative values indicate heat leaving ocean.\n",
    "\n",
    "# geothermal heating at ocean bottom \n",
    "geothermal = forcing.variables['internal_heat'][n]\n",
    "\n",
    "# net heat crossing ocean surface due to all processes, except restoring \n",
    "net_heat_surface = forcing.variables['net_heat_surface'][n]\n",
    "\n",
    "# net heat passed through coupler from shortwave, longwave, latent, sensible.\n",
    "# note: latent includes heat to vaporize liquid and heat to melt ice/snow. \n",
    "# note: sensible includes air-sea and ice-sea sensible heat fluxes.  \n",
    "net_heat_coupler = forcing.variables['net_heat_coupler'][n]\n",
    "\n",
    "# sum of longwave + latent + sensible\n",
    "LwLatSens = forcing.variables['LwLatSens'][n]\n",
    "\n",
    "# net shortwave passing through ocean surface\n",
    "SW = forcing.variables['SW'][n]\n",
    "\n",
    "# heating of liquid seawater due to formation of frazil sea ice\n",
    "frazil = forcing.variables['frazil'][n]\n",
    "\n",
    "# heat flux due to restoring \n",
    "heat_restore = forcing.variables['heat_added'][n]\n",
    "\n",
    "# net heat content associated with transfer of mass across ocean surface, \n",
    "# computed relative to 0C. Both diagnostics should be the same, though \n",
    "# they are computed differently in MOM6.   \n",
    "heat_pme = forcing.variables['Heat_PmE'][n]\n",
    "heat_content_surfwater = forcing.variables['heat_content_surfwater'][n]\n",
    "\n",
    "# heat content associated with water mass leaving ocean\n",
    "heat_content_massout = forcing.variables['heat_content_massout'][n]\n",
    "\n",
    "# heat content associated with water mass entering ocean\n",
    "heat_content_massin = forcing.variables['heat_content_massin'][n]\n",
    "\n",
    "# heat content associated with liquid precipitation \n",
    "heat_content_lprec = forcing.variables['heat_content_lprec'][n]\n",
    "\n",
    "# heat content associated with frozen precipitation \n",
    "heat_content_fprec = forcing.variables['heat_content_fprec'][n]\n",
    "\n",
    "# heat content associated with virtual precipitation \n",
    "heat_content_vprec = forcing.variables['heat_content_vprec'][n]\n",
    "\n",
    "# heat content associated with liquid runoff \n",
    "heat_content_lrunoff = forcing.variables['heat_content_lrunoff'][n]\n",
    "\n",
    "# heat content associated with frozen runoff \n",
    "heat_content_frunoff = forcing.variables['heat_content_frunoff'][n]\n",
    "\n",
    "# heat content associated with liquid condensation \n",
    "heat_content_cond = forcing.variables['heat_content_cond'][n]\n",
    "\n",
    "#--------------------------------------------------------------\n",
    "# salt flux crossing ocean surface and bottom [kg/(m^2 s)]\n",
    "# positive values indicate salt entering ocean; \n",
    "# negative values indicate salt leaving ocean.\n",
    "\n",
    "# salt flux arising from ocean-ice interactions. \n",
    "# this term is zero for this test case, as there is no icea model. \n",
    "#salt_flux = forcing.variables['salt_flux'][n]\n",
    "salt_flux  = 0.0*forcing.variables['heat_content_cond'][n]\n",
    "\n",
    "# salt flux associated with surface restoring.  \n",
    "# this term is zero for this test case, as the SSS restoring \n",
    "# is converted to vprec.\n",
    "#salt_restore = forcing.variables['salt_restore'][n]\n",
    "salt_restore = 0.0*forcing.variables['heat_content_cond'][n]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# for easy setup of subplots\n",
    "def newSP(y,x):\n",
    "    global __spv, __spi ; __spv = (y,x) ; __spi = 1 ; plt.subplot(__spv[0], __spv[1], __spi)\n",
    "def nextSP():\n",
    "    global __spv, __spi ; __spi = __spi + 1 ; plt.subplot(__spv[0], __spv[1], __spi)\n",
    "    \n",
    "# to reduce the amount of code when plotting fields\n",
    "def make_plot(lon, lat, field, title, xmin=-280, xmax=80, ymin=-80, ymax=90, cmin=-200, cmax=200, xlabel=False):\n",
    "   '''\n",
    "   Uses pcolormesh to plot field as a function of lat and lon.\n",
    "   Writes the max, min and ave values of field into the plot.\n",
    "   '''\n",
    "   global area\n",
    "\n",
    "   field_min = numpy.amin(field)\n",
    "   field_max = numpy.amax(field)\n",
    "   field_ave = (field*area).sum() / area.sum()\n",
    "   ch = plt.pcolormesh(lon,lat,field)\n",
    "   cbax=plt.colorbar(ch, extend='both')\n",
    "   plt.title(r''+title)\n",
    "   if (cmin != 0.0 or cmax != 0.0):\n",
    "     plt.clim(cmin,cmax)\n",
    "\n",
    "   plt.xlim(xmin,xmax)\n",
    "   plt.ylim(ymin,ymax)\n",
    "   plt.ylabel(r'Latitude [$\\degree$N]')\n",
    "   if xlabel: plt.xlabel(r'Longitude')\n",
    "   axis = plt.gca()\n",
    "   axis.annotate('max=%5.2f\\nmin=%5.2f\\nave=%5.2f'%(field_max,field_min,field_ave),xy=(0.01,0.73),\n",
    "              xycoords='axes fraction', verticalalignment='bottom', fontsize=8, color='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<hr --- />\n",
    "\n",
    "<h1 align=\"center\">Mass fluxes and global seawater mass budget</h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Global seawater mass budget consistency check</h2>\n",
    "\n",
    "We compute the change in seawater mass over a given time period.  Two different methods are used, and the two methods should agree at the level of truncation error.  Note that \"truncation error\" precision is somewhat larger using offline diagnostics relative to online calculations, particularly if single precision output is saved rather than double precision.  \n",
    "\n",
    "The net mass per time of water (units of kg/s) entering through the ocean boundaries is given by the area integral\n",
    "$$\\begin{equation*}\n",
    "\\mbox{boundary water mass entering liquid seawater} = \\int Q_{W} \\, dA,\n",
    "\\end{equation*}$$\n",
    "where the net water flux (units of $\\mbox{kg}~\\mbox{m}^{-2}~\\mbox{s}^{-1}$) is given by \n",
    "$$\\begin{align*}\n",
    " Q_{W} &= {\\tt PRCmE}\n",
    "\\end{align*}$$\n",
    "A nonzero surface mass flux is associated with liquid and solid precipitation and runoff; evaporation and condensation; sea ice melt/formation; and surface restoring.  \n",
    "\n",
    "The time change of liquid seawater mass is computed according to \n",
    "$$\\begin{equation*}\n",
    "\\mbox{seawater mass change} = \n",
    "\\frac{1}{\\tau_{n+1} - \\tau_{n} } \\int dA \\left(\\int (\\rho_{n+1} - \\rho_{n}) \\, \\mathrm{d}z \\right) \n",
    "\\end{equation*}$$\n",
    "where $\\tau_{n+1} - \\tau_{n}$ is the time increment in seconds.  Note that we make use of the MOM6 diagnostic for depth integrated density \n",
    "$$\\begin{equation*}\n",
    " {\\tt mass\\_wt} =  \\int \\rho \\, \\mathrm{d}z.\n",
    "\\end{equation*}$$\n",
    "For a Boussinesq fluid, the in-situ $\\rho$ factor is set to $\\rho_{0}$, in which case the diagnostic field {\\tt mass\\_wt} measures the thickness of a fluid column, multiplied by $\\rho_{0}$.  For self-consistency, we should have the following equality holding to within truncation error  \n",
    "$$\\begin{equation*}\n",
    "\\boxed{\n",
    " \\mbox{boundary water mass entering liquid seawater} = \\mbox{seawater mass change}.\n",
    "}\n",
    "\\end{equation*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n0       = n-1\n",
    "dmass_wt = mass_wt[n] - mass_wt[n0]\n",
    "dt = time[n] - time[n0]\n",
    "lhs = area * dmass_wt / dt\n",
    "rhs = area * ( PRCmE )\n",
    "\n",
    "print ('Total seawater mass at time step n  [kg seawater]  =',(mass_wt[n]*area).sum())\n",
    "print ('Total seawater mass at time step n0 [kg seawater]  =',(mass_wt[n0]*area).sum())\n",
    "print ('Total seawater mass content change [kg seawater]   =',dt*lhs.sum())\n",
    "print ('Net water mass through boundaries [kg seawater]    =',dt*rhs.sum())\n",
    "print ('Residual [kg seawater]                             =',dt*lhs.sum() - dt*rhs.sum())\n",
    "print ('Non-dimensional residual (based on dmass_wt)       =',(  lhs.sum() - rhs.sum() )/lhs.sum())\n",
    "print ('Non-dimensional residual (based on mass_wt[n])     =',(  lhs.sum() - rhs.sum() )/(mass_wt[n]*area).sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Surface mass fluxes I: combined fields</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,12))\n",
    "newSP(2,2);\n",
    "\n",
    "field     = 86400.0*PRCmE\n",
    "make_plot(lon,lat,field, '$PRCmE$ [$kg/m^2/day$]',cmin=-20,cmax=20)\n",
    "\n",
    "nextSP()\n",
    "field     = 86400.0*net_massin\n",
    "make_plot(lon,lat,field, 'net_massin [$kg/m^2/day$]',cmin=-20,cmax=20)\n",
    "\n",
    "nextSP()\n",
    "field     = 86400.0*net_massout\n",
    "make_plot(lon,lat,field, 'net_massout [$kg/m^2/day$]',cmin=-20,cmax=20,xlabel=True)\n",
    "\n",
    "nextSP()\n",
    "# this should be within \"truncation error\" precision\n",
    "field     = 86400.0*(PRCmE - net_massout - net_massin) \n",
    "make_plot(lon,lat,field, '$PRCmE - M_{in} - M_{out}$ [$kg/m^2/day$]',cmin=0.0,cmax=0.0,xlabel=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Surface mass fluxes II: component fields</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,12))\n",
    "newSP(3,2);\n",
    "\n",
    "field     = 86400.0*lrunoff\n",
    "make_plot(lon,lat,field, 'lrunoff [$kg/m^2/day$]',cmin=-10,cmax=10)\n",
    "\n",
    "nextSP()\n",
    "field     = 86400.0*frunoff\n",
    "make_plot(lon,lat,field, 'frunoff [$kg/m^2/day$]',cmin=-10,cmax=10)\n",
    "\n",
    "nextSP()\n",
    "# lprec contains a contribution from sea ice melt/formation.\n",
    "field     = 86400.0*lprec\n",
    "make_plot(lon,lat,field, 'lprec [$kg/m^2/day$]',cmin=-10,cmax=10)\n",
    "\n",
    "nextSP()\n",
    "field     = 86400.0*fprec\n",
    "make_plot(lon,lat,field, 'fprec [$kg/m^2/day$]',cmin=-1,cmax=1)\n",
    "\n",
    "nextSP()\n",
    "# evaporation and condensation\n",
    "field     = 86400.0*evap\n",
    "make_plot(lon,lat,field, 'evap [$kg/m^2/day$]',cmin=-10,cmax=10,xlabel=True)\n",
    "\n",
    "nextSP()\n",
    "field     = 86400.0*vprec\n",
    "make_plot(lon,lat,field, 'vprec [$kg/m^2/day$]',cmin=-10,cmax=10,xlabel=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Surface mass flux self-consistency check</h2>\n",
    "\n",
    "We will now check if PRCmE = lprec + fprec + lrunoff + frunoff + vprec + evap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,12))\n",
    "newSP(2,1);\n",
    "# this should be within \"truncation error\" precision\n",
    "field     = 86400.0*(PRCmE -lprec -fprec -lrunoff -frunoff -vprec -evap)\n",
    "make_plot(lon,lat,field, 'PRCmE - lprec - fprec - lrunoff - frunoff - vprec - evap [$kg/m^2/day$]',cmin=0.0,cmax=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<hr --- />\n",
    "\n",
    "<h1 align=\"center\">Heat fluxes and global ocean heat budget</h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Global heat budget consistency check</h2>\n",
    "\n",
    "We compute the change in seawater heat content over a given time period.  Two different methods are used, and the two methods should agree at the level of truncation error.  If larger differences exist, then there is a bug.\n",
    "\n",
    "The net heat per time (units of Watts) entering through the ocean boundaries is given by the area integral\n",
    "$$\\begin{equation*}\n",
    "\\mbox{boundary heating of liquid seawater} = \\int Q \\, dA,\n",
    "\\end{equation*}$$\n",
    "where the net heat flux (units of $\\mbox{W}~\\mbox{m}^{-2}~\\mbox{s}^{-1}$) is given by \n",
    "$$\\begin{align*}\n",
    " Q &= {\\tt (net\\_heat\\_coupler + heat\\_pme + frazil) + internal\\_heat + heat\\_restore} \\\\\n",
    "   &= {\\tt net\\_heat\\_surface + internal\\_heat + heat\\_restore}\n",
    "\\end{align*}$$\n",
    "The time change of liquid seawater heat is computed according to \n",
    "$$\\begin{equation*}\n",
    "\\mbox{seawater heat content change} = \n",
    "\\frac{C_p }{\\tau_{n+1} - \\tau_{n} } \\int dA \\left(\\rho_0 \\int (\\Theta_{n+1} - \\Theta_{n}) \\, \\mathrm{d}z \\right) \n",
    "\\end{equation*}$$\n",
    " where $\\tau_{n+1} - \\tau_{n}$ is the time increment in seconds.  Note that we make use of the MOM6 diagnostic for depth integrated potential/conservative temperature  \n",
    "$$\\begin{equation*}\n",
    " {\\tt tomint} = \\rho_0 \\int \\Theta \\, \\mathrm{d}z,\n",
    "\\end{equation*}$$\n",
    " where the Boussinesq reference density, $\\rho_{0}$, is used since this test case makes the Boussinesq approximation. For self-consistency, we should have the following equality holding to within truncation error  \n",
    "$$\\begin{equation*}\n",
    "\\boxed{\n",
    " \\mbox{boundary heating of liquid seawater} = \\mbox{seawater heat content change}.\n",
    "}\n",
    "\\end{equation*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n0      = n-1\n",
    "dtomint = tomint[n] - tomint[n0]\n",
    "dt      = time[n] - time[n0]\n",
    "lhs     = cp * area * dtomint / dt\n",
    "rhs     = area * ( net_heat_coupler + heat_pme + geothermal + frazil + heat_restore )\n",
    "\n",
    "print ('Total seawater heat at time step n  [Joules]   =',cp * (area * tomint[n]).sum())\n",
    "print ('Total seawater heat at time step n0 [Joules]   =',cp * (area* tomint[n0]).sum())\n",
    "print ('Total seawater heat content change [Joules]    =',dt*lhs.sum())\n",
    "print ('Net heat through boundaries [Joules]           =',dt*rhs.sum())\n",
    "print ('Residual [Joules]                              =',dt*lhs.sum() - dt*rhs.sum())\n",
    "print ('Non-dimensional residual (based on difference) =',( lhs.sum() - rhs.sum() )/lhs.sum())\n",
    "print ('Non-dimensional residual (based on tomint[n])  =',( lhs.sum() - rhs.sum()) / (cp * (area * tomint[n]).sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Basic components to surface heat flux</h2> \n",
    "\n",
    "Self-consistency check: net_heat_surface = heat_pme + frazil + net_heat_coupler + heat_restore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(18,12))\n",
    "newSP(3,2);\n",
    "\n",
    "field     = net_heat_surface\n",
    "make_plot(lon,lat,field, 'net_heat_surface[$W/m^2$]')\n",
    "\n",
    "nextSP()\n",
    "field     = net_heat_coupler\n",
    "make_plot(lon,lat,field, 'net_heat_coupler[$W/m^2$]')\n",
    "\n",
    "nextSP()\n",
    "field     = frazil\n",
    "make_plot(lon,lat,field, 'frazil[$W/m^2$]',cmin=-20,cmax=20)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_pme\n",
    "make_plot(lon,lat,field, 'heat_pme[$W/m^2$]',cmin=-20,cmax=20)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_restore\n",
    "make_plot(lon,lat,field, 'heat_restore [$W/m^2$]',cmin=-20.0,cmax=20.0,xlabel=True)\n",
    "\n",
    "nextSP()\n",
    "# this should be within \"truncation error\" precision\n",
    "field = net_heat_surface-net_heat_coupler-frazil-heat_pme-heat_restore\n",
    "make_plot(lon,lat,field, 'Residual(error)[$W/m^2$]',cmin=0.0,cmax=0.0,xlabel=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Heat fluxes crossing ocean surface via the coupler</h2> \n",
    "\n",
    "Heat fluxes crossing ocean surface via the coupler: net_heat_coupler =  LwLatSens + SW, where LwLatSens = LW + Latent + Sensible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,10))\n",
    "newSP(2,2);\n",
    "\n",
    "field     = net_heat_coupler\n",
    "make_plot(lon,lat,field, 'net_heat_coupler[$W/m^2$]',cmin=-200,cmax=200)\n",
    "\n",
    "nextSP()\n",
    "field     = LwLatSens\n",
    "make_plot(lon,lat,field, 'LwLatSens [$W/m^2$]',cmin=-200,cmax=200)\n",
    "\n",
    "nextSP()\n",
    "field     = SW\n",
    "make_plot(lon,lat,field, 'SW [$W/m^2$]',cmin=-200,cmax=200, xlabel=True)\n",
    "\n",
    "nextSP()\n",
    "# this should be within \"truncation error\" precision\n",
    "field     = net_heat_coupler - SW - LwLatSens\n",
    "make_plot(lon,lat,field, 'Residual(error) [$W/m^2$]',cmin=0.0,cmax=0.0, xlabel=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Relation between heat_PmE, heat_content_massin, and heat_content_massout</h2>\n",
    "\n",
    "Alternative means to compute to heat_PmE via heat_content_massin and heat_content_massout, where heat_PmE = heat_content_massin + heat_content_massout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(18,12))\n",
    "newSP(3,2);\n",
    "\n",
    "field     = heat_content_massout\n",
    "make_plot(lon,lat,field, 'heat_content_massout [$W/m^2$]',cmin=-20,cmax=0)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_content_massin\n",
    "make_plot(lon,lat,field, 'heat_content_massin [$W/m^2$]',cmin=0,cmax=20)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_content_massin + heat_content_massout\n",
    "make_plot(lon,lat,field, 'heat_content_massin + heat_content_massout [$W/m^2$]',cmin=-20,cmax=20)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_pme\n",
    "make_plot(lon,lat,field, 'heat_pme [$W/m^2$]',cmin=-20,cmax=20,xlabel=True)\n",
    "\n",
    "nextSP()\n",
    "# this should be within \"truncation error\" precision\n",
    "field     = heat_content_massout + heat_content_massin - heat_pme\n",
    "make_plot(lon,lat,field, 'heat_massin + heat_massout - heat_pme [$W/m^2$]',cmin=0.0,cmax=0.0, xlabel=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Components of heat content from surface mass fluxes</h2>\n",
    "\n",
    "Components of heat content of surface mass fluxes: heat_PmE = heat_content_lprec + heat_content_fprec + heat_content_vprec + heat_content_lrunoff + heat_content_frunoff + heat_content_cond + heat_content_massout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,12))\n",
    "newSP(3,2);\n",
    "\n",
    "field     = heat_content_lprec\n",
    "make_plot(lon,lat,field, 'heat_content_lprec [$W/m^2$]',cmin=-20.0,cmax=20.0)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_content_lrunoff\n",
    "make_plot(lon,lat,field, 'heat_content_lrunoff [$W/m^2$]',cmin=-20.0,cmax=20.0)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_content_frunoff\n",
    "make_plot(lon,lat,field, 'heat_content_frunoff [$W/m^2$]',cmin=-20.0,cmax=20.0)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_content_cond\n",
    "make_plot(lon,lat,field, 'heat_content_cond [$W/m^2$]',cmin=-1.0,cmax=1.0)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_content_fprec\n",
    "make_plot(lon,lat,field, 'heat_content_fprec [$W/m^2$]',cmin=-1.0,cmax=1.0)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_content_vprec\n",
    "make_plot(lon,lat,field, 'heat_content_vprec [$W/m^2$]',cmin=-5.0,cmax=5.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Self-consistency of diagnosed heat content from mass entering ocean</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,12))\n",
    "newSP(3,2);\n",
    "\n",
    "heat_content_sum = ( heat_content_lprec + heat_content_fprec + heat_content_vprec +\n",
    "                     heat_content_lrunoff + heat_content_frunoff + heat_content_cond )\n",
    "\n",
    "field     = heat_content_massin\n",
    "make_plot(lon,lat,field, 'heat_content_massin [$W/m^2$]',cmin=-20.0,cmax=20.0)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_content_sum\n",
    "make_plot(lon,lat,field, '$\\Sigma$ of heat_contents entering ocean [$W/m^2$]',cmin=-20.0,cmax=20.0)\n",
    "\n",
    "nextSP()\n",
    "# this should be within \"truncation error\" precision\n",
    "field     = heat_content_massin - heat_content_sum\n",
    "make_plot(lon,lat,field, 'heat_content_massin - $\\Sigma$ of components [$W/m^2$]',cmin=0.0,cmax=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Self-consistency between heat_pme and heat_content_surfwater</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comp_sum = ( heat_content_lprec + heat_content_fprec + heat_content_vprec + heat_content_lrunoff\n",
    "           + heat_content_frunoff + heat_content_cond + heat_content_massout )\n",
    "\n",
    "plt.figure(figsize=(16,12))\n",
    "newSP(3,2);\n",
    "\n",
    "field     = heat_content_surfwater\n",
    "make_plot(lon,lat,field, 'heat_content_surfwater [$W/m^2$]',cmin=-20.0,cmax=20.0)\n",
    "\n",
    "nextSP()\n",
    "field     = comp_sum\n",
    "make_plot(lon,lat,field, '$\\Sigma$ of heat_content components [$W/m^2$]',cmin=-20.0,cmax=20.0)\n",
    "\n",
    "nextSP()\n",
    "# this should be within \"truncation error\" precision\n",
    "field     = comp_sum - heat_content_surfwater\n",
    "make_plot(lon,lat,field, '$\\Sigma$ - heat_content_surfwater [$W/m^2$]',cmin=0.0,cmax=0.0)\n",
    "\n",
    "nextSP()\n",
    "field     = heat_pme\n",
    "make_plot(lon,lat,field, 'heat_pme [$W/m^2$]',cmin=-20.0,cmax=20.0, xlabel=True)\n",
    "\n",
    "nextSP()\n",
    "# this should be within \"truncation error\" precision\n",
    "field     = heat_pme - heat_content_surfwater\n",
    "make_plot(lon,lat,field, 'heat_pme - heat_content_surfwater [$W/m^2$]',cmin=0.0,cmax=0.0, xlabel=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Map effective temperatures</h2>\n",
    "\n",
    "The following \"effective\" temperatures differ generally from the SST due to the means by which the water is exchanged across the ocean surface boundary. In particular, there are some occasions whereby layers deeper than k=1 need to be sampled in order to exchange water with the atmosphere and/or sea ice.  These maps provide us with a sense for where such occurrances take place.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TinEff  = heat_content_massin/(net_massin*cp)\n",
    "ToutEff = heat_content_massout/(net_massout*cp)\n",
    "TnetEff = heat_pme/(PRCmE*cp)\n",
    "\n",
    "plt.figure(figsize=(16,12))\n",
    "newSP(3,2);\n",
    "\n",
    "field     = TinEff\n",
    "make_plot(lon,lat,field, '$\\Theta_{in} [{\\degree}C$]',cmin=-2,cmax=30)\n",
    "\n",
    "nextSP()\n",
    "field     = TinEff - sst\n",
    "make_plot(lon,lat,field, '$\\Theta_{in} - SST [{\\degree}C$]',cmin=0.0,cmax=0.0)\n",
    "\n",
    "nextSP()\n",
    "field     = ToutEff\n",
    "make_plot(lon,lat,field, '$\\Theta_{out} [{\\degree}C$]',cmin=-2,cmax=30)\n",
    "\n",
    "nextSP()\n",
    "field     = ToutEff - sst\n",
    "make_plot(lon,lat,field, '$\\Theta_{out} - SST [{\\degree}C$]',cmin=0.0,cmax=0.0)\n",
    "\n",
    "nextSP()\n",
    "field     = TnetEff\n",
    "make_plot(lon,lat,field, '$\\Theta_{net} [{\\degree}C$]',cmin=-2,cmax=30, xlabel=True)\n",
    "\n",
    "nextSP()\n",
    "field     = TnetEff - sst\n",
    "make_plot(lon,lat,field, '$\\Theta_{net} - SST [{\\degree}C$]',cmin=0.0,cmax=0.0, xlabel=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Geothermal heat flux</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,12))\n",
    "newSP(2,2);\n",
    "\n",
    "field     = geothermal\n",
    "make_plot(lon,lat,field, 'Geothermal heat flux[$W/m^2$]',cmin=0.0,cmax=0.1,xlabel=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<hr --- />\n",
    "\n",
    "<h1 align=\"center\">Salt fluxes and global ocean salt budget</h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 align=\"center\">Global salt budget consistency check</h2>\n",
    "\n",
    "We compute the change in seawater salt content over a given time period.  Two different methods are used, and the two methods should agree at the level of truncation error.  If larger differences exist, then there is a bug.\n",
    "\n",
    "The net salt per time (units of kg/s) entering through the ocean boundaries is given by the area integral\n",
    "$$\\begin{equation*}\n",
    "\\mbox{boundary salt entering liquid seawater} = \\int Q_{S} \\, dA,\n",
    "\\end{equation*}$$\n",
    "where the net salt flux (units of $\\mbox{kg}~\\mbox{m}^{-2}~\\mbox{s}^{-1}$) is given by \n",
    "$$\\begin{align*}\n",
    " Q_{S} &= {\\tt salt\\_flux} + {\\tt salt\\_restore}.\n",
    "\\end{align*}$$\n",
    "A nonzero salt flux is associated with exchanges between liquid seawater and solid sea ice.  It also arises from simulations using a restoring boundary flux associated with damping to observed sea surface salinity.  Finally, there can be a salt flux when using sponges to damp the ocean interior back to an observed value.  \n",
    "\n",
    "The time change of liquid seawater salt content is computed according to \n",
    "$$\\begin{equation*}\n",
    "\\mbox{seawater salt content change} = \n",
    "\\frac{1}{\\tau_{n+1} - \\tau_{n} } \\int dA \\left(\\rho_0 \\int (S_{n+1} - S_{n}) \\, \\mathrm{d}z \\right) \n",
    "\\end{equation*}$$\n",
    "where $\\tau_{n+1} - \\tau_{n}$ is the time increment in seconds.  Note that we make use of the MOM6 diagnostic for depth integrated salinity \n",
    "$$\\begin{equation*}\n",
    " {\\tt somint} = \\rho_0 \\int S \\, \\mathrm{d}z,\n",
    "\\end{equation*}$$\n",
    "where the Boussinesq reference density, $\\rho_{0}$, is used since this test case makes the Boussinesq approximation. For self-consistency, we should have the following equality holding to within truncation error  \n",
    "$$\\begin{equation*}\n",
    "\\boxed{\n",
    " \\mbox{boundary salt entering liquid seawater} = \\mbox{seawater salt content change}.\n",
    "}\n",
    "\\end{equation*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n0      = n-1\n",
    "dsomint = somint[n] - somint[n0]\n",
    "time    = surface.variables[tvar][:]*86400.\n",
    "dt      = time[n] - time[n0]\n",
    "lhs     = area * dsomint / dt\n",
    "rhs     = area * ( salt_flux + salt_restore )\n",
    "\n",
    "print('\\n \\n')\n",
    "print ('Total salt at time step n  [kg salt]          =', (area * somint[n]).sum())\n",
    "print ('Total salt at time step n0 [kg salt]          =',(area* somint[n0]).sum())\n",
    "print ('Total seawater salt content change [kg salt]  =',dt*lhs.sum())\n",
    "print ('Net salt through boundaries [kg salt]         =',dt*rhs.sum())\n",
    "print ('Residual [kg salt]                            =', (area * somint[n]).sum() - (area* somint[n0]).sum())\n",
    "print ('Non-dimensional residual                      =', ((area * somint[n]).sum()- (area* somint[n0]).sum())/(area * somint[n]).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,12))\n",
    "newSP(2,2)\n",
    "\n",
    "field     = 86400.0*salt_flux\n",
    "make_plot(lon,lat,field, 'Surface salt flux from ice-ocean exchange [kg m$^{-2}$ day$^{-1}$]',cmin=0.0,cmax=0.0,xlabel=True)\n",
    "\n",
    "nextSP()\n",
    "field     = 86400.0*salt_restore\n",
    "make_plot(lon,lat,field, 'Surface salt flux from restoring [kg m$^{-2}$ day$^{-1}$]',cmin=0.0,cmax=0.0,xlabel=True)"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}