time_dependent_model_loggrid_multiple_zo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time-dependent boundary layer model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "import gsw\n",
    "from scipy.interpolate import interp1d\n",
    "from scipy.integrate import solve_ivp\n",
    "import datapath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ds = xr.open_dataset(datapath.adcpnc())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lat = 35+8.4585/60  # latitude (for f calculation)\n",
    "bini = 15 # index of top bin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Smooth ADCP time series\n",
    "* 30-min running median to remove spikes\n",
    "* 60-min running mean to smooth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'binheight' ()>\n",
       "array(17.14)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds['binheight'][bini]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ds['Eas_sm'] = ds['Eas'].rolling(time=6,center=True).median().rolling(time=12,center=True).mean()\n",
    "ds['Nor_sm'] = ds['Nor'].rolling(time=6,center=True).median().rolling(time=12,center=True).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate time derivative of ADCP velocity for model forcing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# add instrument height (18 inches)\n",
    "zf = ds['binheight'][bini]+ 0.4572 # upper boundary height\n",
    "zobs = np.array(ds['binheight'][0:(bini+1)]) + 0.4572"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "i = 1j\n",
    "kappa = 0.41 # Von Karman constant\n",
    "f = gsw.f(lat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gi = np.isfinite(ds['Eas_sm'][:,0])\n",
    "tobs = np.array((ds['time'][gi] - ds['time'][gi][0])/np.timedelta64(1,'s')) # array of times in seconds\n",
    "deltat = tobs[1]-tobs[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "wobs = np.array(ds['Eas_sm'][gi,0:(bini+1)]+i*ds['Nor_sm'][gi,0:(bini+1)])\n",
    "dwobsdt = np.gradient(wobs,deltat,axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define model equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dwdt_bbl(t, w_in):\n",
    "    \n",
    "    global zo\n",
    "    \n",
    "    N = int(len(w_in)/2)\n",
    "    w = w_in[:N]+i*w_in[N:]\n",
    "    \n",
    "    winf = np.interp(t,tobs,wobs[:,bini])\n",
    "    dwinfdt = np.interp(t,tobs,dwobsdt[:,bini])\n",
    "\n",
    "    F = dwinfdt + i*f*winf\n",
    "    ustar = kappa*zo*np.abs(w[0])/dz[0]\n",
    "    coef1 = 2*kappa*ustar/(dz[:-1]+dz[1:])\n",
    "\n",
    "    wall = np.concatenate([np.atleast_1d(0+i*0),w,np.atleast_1d(winf)])\n",
    "    \n",
    "    dwalldt = np.nan*np.ones(np.shape(wall))+i*np.nan*np.ones(np.shape(wall))\n",
    "    dwalldt[1:-1] = (F +\n",
    "                     coef1*(z[1:-1]/dz[1:])*wall[2:] +\n",
    "                     (-i*f - coef1*(z[1:-1]/dz[1:]+z[:-2]/dz[:-1]))*wall[1:-1] +\n",
    "                     coef1*(z[:-2]/dz[:-1])*wall[:-2])\n",
    "    \n",
    "    dwdt = dwalldt[1:-1]\n",
    "    \n",
    "    dwdt_out = np.concatenate([np.real(dwdt),np.imag(dwdt)])\n",
    "    \n",
    "    return dwdt_out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run model for varying $z_o$ values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zo: 0.0001 ( 1 / 7 )\n",
      "944 seconds\n",
      "zo: 0.00031622776601683794 ( 2 / 7 )\n",
      "734 seconds\n",
      "zo: 0.001 ( 3 / 7 )\n",
      "413 seconds\n",
      "zo: 0.0031622776601683794 ( 4 / 7 )\n",
      "449 seconds\n",
      "zo: 0.01 ( 5 / 7 )\n",
      "434 seconds\n",
      "zo: 0.03162277660168379 ( 6 / 7 )\n",
      "479 seconds\n",
      "zo: 0.1 ( 7 / 7 )\n",
      "516 seconds\n"
     ]
    }
   ],
   "source": [
    "nzo = 7\n",
    "zo_array = np.logspace(-4,-1,nzo) # roughness length array\n",
    "\n",
    "tf = 31*86400\n",
    "ntimes = len(tobs[tobs<=tf])\n",
    "\n",
    "nz = 201\n",
    "\n",
    "w_sol_all = np.nan*np.ones([ntimes,nz-2,nzo]) + i*np.nan*np.ones([ntimes,nz-2,nzo])\n",
    "w_sol_obs_all = np.nan*np.ones([ntimes,len(zobs),nzo]) + i*np.nan*np.ones([ntimes,len(zobs),nzo])\n",
    "\n",
    "for zi,zo in enumerate(zo_array):\n",
    "    print('zo:',zo,'(',zi+1,'/',nzo,')')\n",
    "\n",
    "    z = np.logspace(np.log10(zo),np.log10(zf),nz)\n",
    "    dz = np.diff(z)\n",
    "    \n",
    "    w0_obs = np.interp(z,zobs[1:bini],wobs[0,1:bini])[1:-1]\n",
    "    w0 = np.concatenate([np.real(w0_obs),\n",
    "                           np.imag(w0_obs)])\n",
    "\n",
    "    tic = np.datetime64('now')\n",
    "    a = solve_ivp(dwdt_bbl, t_span=[tobs[0],tf], y0 = w0, t_eval = tobs[tobs<=tf], method='BDF')\n",
    "    toc = np.datetime64('now')\n",
    "    print(toc-tic)\n",
    "    \n",
    "    N = int(np.shape(a.y)[0]/2)\n",
    "    t_sol = a.t\n",
    "    w_sol = a.y[:N,:]+i*a.y[N:,:]\n",
    "    \n",
    "    fint= interp1d(z[1:-1],w_sol,axis=0,\n",
    "                   bounds_error=False,fill_value = np.nan+i*np.nan)\n",
    "    w_sol_obs = fint(zobs)\n",
    "    \n",
    "    w_sol_all[:,:,zi] = w_sol.T\n",
    "    w_sol_obs_all[:,:,zi] = w_sol_obs.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save model runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bbl_sols = xr.Dataset(\n",
    "        data_vars={'uobs_subset':    (('tobs', 'zobs'), np.real(wobs[tobs<=tf,:])),\n",
    "                   'vobs_subset':    (('tobs', 'zobs'), np.imag(wobs[tobs<=tf,:])),\n",
    "                   'usol_all': (('tobs','zsol','zo'), np.real(w_sol_all)),\n",
    "                   'vsol_all': (('tobs','zsol','zo'), np.imag(w_sol_all)),\n",
    "                   'usol_obs_all': (('tobs','zobs','zo'), np.real(w_sol_obs_all)),\n",
    "                   'vsol_obs_all': (('tobs','zobs','zo'), np.imag(w_sol_obs_all)),                  },\n",
    "        coords={'zobs': zobs,\n",
    "                'tobs': tobs[tobs<=tf],\n",
    "                'zo':zo_array})\n",
    "bbl_sols.to_netcdf('data/bbl_model_solutions.nc',mode='w')\n",
    "bbl_sols.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py37]",
   "language": "python",
   "name": "conda-env-py37-py"
  },
  "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}