EUC_indexing_issue.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netCDF4 import Dataset\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following routines do essentially the same calculations, but one uses the full 3d `umo` field whereas the other just uses the subset saved for an individual strait"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_euc_from_umo(geolon_u, geolat_u, z_i, umo):\n",
    "  assert len(umo.shape)==4,'umo must have dimensions [time x depth x lat x lon]'\n",
    "  # Get the indices corresponding to the EUC region as defined in the OMIP paper (Griffies et al., 2018)\n",
    "  euc_idx = np.nonzero( (geolon_u >= -155) &\n",
    "                        (geolon_u <= -155) &\n",
    "                        (geolat_u >= -2)   &\n",
    "                        (geolat_u <= 2 )   ) \n",
    "  zlim = 350 # Define the depth to integrate to\n",
    "  \n",
    "  # Subset the data\n",
    "  umo_euc = umo[:,:,euc_idx[0],euc_idx[1]].copy()\n",
    "  zidx = np.sum(z_i[:]<=zlim)-1 # Find the last index to integrate down to\n",
    "  if z_i[zidx] < zlim: # This will not happen if the WOA05 levels are used\n",
    "    frac = (z_i[zidx+1] - zlim)/(z_i[zidx+1]-z_i[zidx])\n",
    "    umo_euc[:,zidx+1,:] = umo_euc[:,zidx+1,:]*frac\n",
    "    umo_euc[:,zidx+2:,:] = 0.\n",
    "  else: # Zero out anything below\n",
    "    umo_euc[:,zidx+1:,:] = 0.\n",
    "  transport = umo_euc.sum((1,2))\n",
    "  return transport\n",
    "\n",
    "def calc_euc_from_slice(z_i, umo):\n",
    "  zlim = 350 # Define the depth to integrate to\n",
    "  \n",
    "  # Subset the data\n",
    "  umo_euc = umo.copy()\n",
    "  zidx = np.sum(z_i[:]<=zlim)-1 # Find the last index to integrate down to\n",
    "  if z_i[zidx] < zlim: # This will not happen if the WOA05 levels are used\n",
    "    frac = (z_i[zidx+1] - zlim)/(z_i[zidx+1]-z_i[zidx])\n",
    "    umo_euc[:,zidx+1,...] = umo_euc[:,zidx+1,:]*frac\n",
    "    umo_euc[:,zidx+2:,...] = 0.\n",
    "  else: # Zero out anything below\n",
    "    umo_euc[:,zidx+1:,...] = 0.\n",
    "  transport = umo_euc.sum((1,2,3))\n",
    "  return transport"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the transports using both methods\n",
    "datapath = '/raid/ra40/data/ras/refine_diag/'\n",
    "geolon_u = Dataset(datapath + '06500101.ocean_static.nc')['geolon_u'][:,:]\n",
    "geolat_u = Dataset(datapath + '06500101.ocean_static.nc')['geolat_u'][:,:]\n",
    "\n",
    "umo = {}\n",
    "umo['month3d'] = Dataset(datapath+'17080101.ocean_month_z.nc')['umo'][:]\n",
    "z_i = Dataset(datapath+'17080101.ocean_month_z.nc')['z_i'][:]\n",
    "umo['month_strait'] = Dataset(datapath+'17080101.ocean_month_EUC.nc')['umo'][:]\n",
    "\n",
    "xq = Dataset(datapath+'17080101.ocean_month_z.nc')['xq'][:]\n",
    "yh = Dataset(datapath+'17080101.ocean_month_z.nc')['yh'][:]\n",
    "euc = {}\n",
    "euc['month3d'] = calc_euc_from_umo(geolon_u,geolat_u,z_i,umo['month3d'])\n",
    "euc['month_strait'] = calc_euc_from_slice(z_i,umo['month_strait'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "From 3D UMO:     7955284992.0\n",
      "From Strait UMO: 8426403840.0\n"
     ]
    }
   ],
   "source": [
    "print('From 3D UMO:     {}'.format(euc['month3d'][0]))\n",
    "print('From Strait UMO: {}'.format(euc['month_strait'][0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The transports are different by a small amount (~5%), but this is because the index that's used to find the 155W longitude is based on xq and NOT on the actual tripolar grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Index of 155W based on xq: 582\n"
     ]
    }
   ],
   "source": [
    "# Find the point in xq closes to 155W\n",
    "lonidx = np.argmin( np.abs(xq + 155.) )\n",
    "print('Index of 155W based on xq: {}'.format(lonidx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Longitude from geolon_u: -154.25\n",
      "Transport: 8426403840.0\n"
     ]
    }
   ],
   "source": [
    "latidx = np.nonzero( (yh<=2) & (yh>=-2)) # This corresponds to 2N/S in both geolat and yh\n",
    "print('Longitude from geolon_u: {}'.format(geolon_u[latidx[0][0],lonidx]))\n",
    "transport = calc_euc_from_slice(z_i,umo['month3d'][:,:,latidx,lonidx])\n",
    "print(f'Transport: {transport[0]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This transport is identical to the one put out from the Pacific Undercurrent history file. However, if the same correction is done 'correctly' based on geolon_u, the result is the same as when calculating from the global umo'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Longitude from geolon_u: -155.0\n",
      "Transport: 7955284992.0\n"
     ]
    }
   ],
   "source": [
    "lonidx = 579 # Found manually to in geolon_u\n",
    "print(f'Longitude from geolon_u: {geolon_u[latidx[0][0],lonidx]}')\n",
    "transport = calc_euc_from_slice(z_i,umo['month3d'][:,:,latidx,lonidx])\n",
    "print(f'Transport: {transport[0]}')"
   ]
  },
  {
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}