{
 "metadata": {
  "name": "",
  "signature": "sha256:6e21f1f3eff7c286e5b4bde8bc58193ea6bbdd789971c5a7d9c97b1f0f3bab1c"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook generates a tidal forcing for 2D domain."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plan: Use M2 and K1 consituents at JdF from 3D domain.Only U. Average across boundary.\n",
      "\n",
      "Increase currents by factor of 2. Horizontal constriction in sill area is missing so need to artificially increase currents.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import netCDF4 as nc\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import os\n",
      "\n",
      "from salishsea_tools import tidetools,nc_tools\n",
      "\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Tide Data"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First, quickly look at tide data for 3D model"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T={}; U={}; V={}\n",
      "\n",
      "path = '/data/nsoontie/MEOPAR/NEMO-forcing/open_boundaries/west/tides/'\n",
      "\n",
      "for t in ['M2','K1']:\n",
      "    fname = 'SalishSea2_Flux_RC13_west_tide_{}_grid_{}.nc'.format(t,'T')\n",
      "    T[t]=nc.Dataset(os.path.join(path,fname))\n",
      "    fname = 'SalishSea2_Flux_RC13_west_tide_{}_grid_{}.nc'.format(t,'U')\n",
      "    U[t]=nc.Dataset(os.path.join(path,fname))\n",
      "    fname = 'SalishSea2_Flux_RC13_west_tide_{}_grid_{}.nc'.format(t,'V')\n",
      "    V[t]=nc.Dataset(os.path.join(path,fname))\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      " nc_tools.show_variables(U['K1'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[u'xb', u'yb', u'nbidta', u'nbjdta', u'nbrdta', u'u1', u'u2']\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Look at amp of U/V"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "u1={}; u2={}; v1={}; v2={}; t1={}; t2={}\n",
      "\n",
      "for t in ['K1','M2']:\n",
      "    u1[t] = U[t].variables['u1']\n",
      "    u2[t] = U[t].variables['u2']\n",
      "    \n",
      "    v1[t] = V[t].variables['v1']\n",
      "    v2[t] = V[t].variables['v2']\n",
      "    \n",
      "    t1[t] = T[t].variables['z1']\n",
      "    t2[t] = T[t].variables['z2']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def average_amp_phase(x1,x2):\n",
      "    temp = np.sqrt(x1[:]**2 + x2[:]**2); temp = np.ma.array(temp,mask=temp==0)\n",
      "    amp = np.mean(temp)\n",
      "    pha=[]\n",
      "    for i in range(0,x1.shape[1]):\n",
      "       pha.append(math.atan2(x2[0,i],x1[0,i]))\n",
      "    pha=np.ma.array(pha,mask=pha==0)\n",
      "    pha=np.mean(pha)\n",
      "    \n",
      "    return amp, pha"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Uamp={}; Vamp={}; Tamp={}; Upha={}; Vpha={}; Tpha={};\n",
      "for t in ['K1','M2']:\n",
      "    Uamp[t], Upha[t]= average_amp_phase(u1[t],u2[t])\n",
      "    Vamp[t], Vpha[t]= average_amp_phase(v1[t],v2[t])\n",
      "    Tamp[t], Tpha[t]= average_amp_phase(t1[t],t2[t])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "u1['K1'][:].shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "(1, 87)"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for t in ['K1','M2']:\n",
      "    print t\n",
      "    print Uamp[t], Vamp[t], Tamp[t]\n",
      "    print Upha[t], Vpha[t], Tpha[t]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "K1\n",
        "0.186448183927 0.124745784662 0.402432726575\n",
        "0.565605422043 0.183389151205 -1.87524389342\n",
        "M2\n",
        "0.293222130119 0.181446491144 0.822501739898\n",
        "-1.38122753203 1.69071899544 -2.09972927195\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print t1['K1'].shape\n",
      "print u1['K1'].shape\n",
      "print v1['K1'].shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(1, 87)\n",
        "(1, 87)\n",
        "(1, 88)\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I think I will just use the U components for the velcoity in this very simplified model. Average along the whole boundary."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Construction of 2D forcing"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Average the phase and amplitude and write to the tide components."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "u_factor = 3\n",
      "\n",
      "u2dZ1={}; u2dZ2={}; v2dZ1={}; v2dZ2={}; t2dZ1={}; t2dZ2={}\n",
      "\n",
      "for t in ['K1','M2']:\n",
      "    u2dZ1[t] = u_factor*Uamp[t]*math.cos(Upha[t])*np.ones((1,10))\n",
      "    u2dZ2[t] = u_factor*Uamp[t]*math.sin(Upha[t])*np.ones((1,10))\n",
      "    \n",
      "    v2dZ1[t] = np.zeros((1,10))\n",
      "    v2dZ2[t] = np.zeros((1,10))\n",
      "    \n",
      "    t2dZ1[t] = Tamp[t]*np.cos(Tpha[t])*np.ones((1,10))\n",
      "    t2dZ2[t] = Tamp[t]*np.sin(Tpha[t])*np.ones((1,10))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Create netCDF"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def creat_netcdf(constituent, tidevar, Z1, Z2):   \n",
      "    nemo = nc.Dataset('../Salish2D_'+constituent+'_grid_'+tidevar+'.nc', 'w', zlib=True)\n",
      "\n",
      "    #start and end points\n",
      "    lengthi=10\n",
      "\n",
      "    # dataset attributes\n",
      "    nc_tools.init_dataset_attrs(\n",
      "        nemo,\n",
      "        title='Tidal Boundary Conditions 2D domain',\n",
      "        notebook_name='Generate_2D_tides',\n",
      "        nc_filepath='/data/nsoontie/MEOPAR/2Ddomain/Salish2D_'+constituent+'_grid_'+tidevar+'.nc',\n",
      "        comment='based on average values across mouth of JdF')\n",
      "\n",
      "    # dimensions (only need x and y, don't need depth or time_counter)\n",
      "    nemo.createDimension('xb', lengthi)\n",
      "    nemo.createDimension('yb', 1)\n",
      "\n",
      "    # variables\n",
      "    # nbidta, ndjdta, ndrdta\n",
      "    nbidta = nemo.createVariable('nbidta', 'int32' , ('yb','xb'))\n",
      "    nbidta.long_name = 'i grid position'\n",
      "    nbidta.units = 1\n",
      "    nbjdta = nemo.createVariable('nbjdta', 'int32' , ('yb','xb'))\n",
      "    nbjdta.long_name = 'j grid position'\n",
      "    nbjdta.units = 1\n",
      "    nbrdta = nemo.createVariable('nbrdta', 'int32' , ('yb','xb'))\n",
      "    nbrdta.long_name = 'position from boundary'\n",
      "    nbrdta.units = 1\n",
      "    print nbidta.shape\n",
      "    # add in the counter around the boundary (taken from Susan's code in Prepare Tide Files)\n",
      "    xb = nemo.createVariable('xb', 'int32', ('xb',),zlib=True)\n",
      "    xb.units = 'non dim'\n",
      "    xb.long_name = 'counter around boundary'\n",
      "    yb = nemo.createVariable('yb', 'int32', ('yb',),zlib=True)\n",
      "    yb.units = 'non dim'\n",
      "    yb.long_name = 'counter along boundary'\n",
      "    yb[0] = 1\n",
      "    xb[:] = np.arange(0,lengthi) +1\n",
      "\n",
      "    # values\n",
      "    # nbidta, nbjdta\n",
      "    nbidta[:] = 1\n",
      "    nbjdta[:] = np.arange(0,lengthi) +1\n",
      "    nbrdta[:] = 1\n",
      "\n",
      "    if tidevar=='T':\n",
      "        z1 = nemo.createVariable('z1','float32',('yb','xb'),zlib=True)\n",
      "        z1.units = 'm'\n",
      "        z1.long_name = 'tidal elevation: cosine'\n",
      "        z2 = nemo.createVariable('z2','float32',('yb','xb'),zlib=True)\n",
      "        z2.units = 'm'\n",
      "        z2.long_name = 'tidal elevation: sine'\n",
      "        z1[0,:] = Z1[0,:]\n",
      "        z2[0,:] = Z2[0,:]\n",
      "\n",
      "    if tidevar=='U':\n",
      "        u1 = nemo.createVariable('u1','float32',('yb','xb'),zlib=True)\n",
      "        u1.units = 'm'\n",
      "        u1.long_name = 'tidal x-velocity: cosine'\n",
      "        u2 = nemo.createVariable('u2','float32',('yb','xb'),zlib=True)\n",
      "        u2.units = 'm'\n",
      "        u2.long_name = 'tidal x-velocity: sine'\n",
      "        u1[0,:] = Z1[0,:]\n",
      "        u2[0,:] = Z2[0,:]\n",
      "\n",
      "    if tidevar=='V':\n",
      "        v1 = nemo.createVariable('v1','float32',('yb','xb'),zlib=True)\n",
      "        v1.units = 'm'\n",
      "        v1.long_name = 'tidal y-velocity: cosine'\n",
      "        v2 = nemo.createVariable('v2','float32',('yb','xb'),zlib=True)\n",
      "        v2.units = 'm'\n",
      "        v2.long_name = 'tidal y-velocity: sine'\n",
      "        v1[0,:] = Z1[0,:]\n",
      "        v2[0,:] = Z2[0,:]\n",
      "\n",
      "    nc_tools.check_dataset_attrs(nemo)\n",
      "    nemo.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "varis=['T','U','V']\n",
      "for t in['M2','K1']:\n",
      "    Z1s=[t2dZ1[t],u2dZ1[t],v2dZ1[t]];\n",
      "    Z2s=[t2dZ2[t],u2dZ2[t],v2dZ2[t]];\n",
      "    for var,Z1,Z2 in zip(varis,Z1s,Z2s):\n",
      "        creat_netcdf(t, var, Z1, Z2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "file format: NETCDF4\n",
        "Conventions: CF-1.6\n",
        "title: Tidal Boundary Conditions 2D domain\n",
        "institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia\n",
        "source: REQUIRED\n",
        "references: REQUIRED\n",
        "history: [2015-03-23 09:36:39] Created netCDF4 zlib=True dataset.\n",
        "comment: based on average values across mouth of JdF\n",
        "(1, 10)\n",
        "Missing value for dataset attribute: source\n",
        "Missing value for dataset attribute: references\n",
        "file format: NETCDF4\n",
        "Conventions: CF-1.6\n",
        "title: Tidal Boundary Conditions 2D domain\n",
        "institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia\n",
        "source: REQUIRED\n",
        "references: REQUIRED\n",
        "history: [2015-03-23 09:36:39] Created netCDF4 zlib=True dataset.\n",
        "comment: based on average values across mouth of JdF\n",
        "(1, 10)\n",
        "Missing value for dataset attribute: source\n",
        "Missing value for dataset attribute: references\n",
        "file format: NETCDF4"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Conventions: CF-1.6\n",
        "title: Tidal Boundary Conditions 2D domain\n",
        "institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia\n",
        "source: REQUIRED\n",
        "references: REQUIRED\n",
        "history: [2015-03-23 09:36:39] Created netCDF4 zlib=True dataset.\n",
        "comment: based on average values across mouth of JdF\n",
        "(1, 10)\n",
        "Missing value for dataset attribute: source\n",
        "Missing value for dataset attribute: references\n",
        "file format: NETCDF4\n",
        "Conventions: CF-1.6\n",
        "title: Tidal Boundary Conditions 2D domain\n",
        "institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia\n",
        "source: REQUIRED\n",
        "references: REQUIRED\n",
        "history: [2015-03-23 09:36:39] Created netCDF4 zlib=True dataset.\n",
        "comment: based on average values across mouth of JdF\n",
        "(1, 10)\n",
        "Missing value for dataset attribute: source\n",
        "Missing value for dataset attribute: references\n",
        "file format: NETCDF4"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Conventions: CF-1.6\n",
        "title: Tidal Boundary Conditions 2D domain\n",
        "institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia\n",
        "source: REQUIRED\n",
        "references: REQUIRED\n",
        "history: [2015-03-23 09:36:39] Created netCDF4 zlib=True dataset.\n",
        "comment: based on average values across mouth of JdF\n",
        "(1, 10)\n",
        "Missing value for dataset attribute: source\n",
        "Missing value for dataset attribute: references\n",
        "file format: NETCDF4\n",
        "Conventions: CF-1.6\n",
        "title: Tidal Boundary Conditions 2D domain\n",
        "institution: Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia\n",
        "source: REQUIRED\n",
        "references: REQUIRED\n",
        "history: [2015-03-23 09:36:39] Created netCDF4 zlib=True dataset.\n",
        "comment: based on average values across mouth of JdF\n",
        "(1, 10)\n",
        "Missing value for dataset attribute: source\n",
        "Missing value for dataset attribute: references\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    }
   ],
   "metadata": {}
  }
 ]
}