trajectories_CM1_2nd_order.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# About\n",
    "In recent years, Cloud Model 1 (CM1; http://www2.mmm.ucar.edu/people/bryan/cm1/) has become a very popular tool for performing idealized studies of atmospheric phenomena. There exists very little support for computing trajectories using CM1 output, which are usually necessary to understand the processes of the atmospheric phenomena of interest. Natively, CM1 only supports 'online' forward trajectories in 2D simulations and in 3D simulation without terrain. I wrote this script because there are no adequate tools available to compute highly customizable 'offline' trajectories in simulations with terrain. This script is intended to be easily customizable.\n",
    "\n",
    "Notes:\n",
    "\n",
    "* Can compute backward or forward trajectories (Default is backward, but can be forward with simple changes to \"Calculate Trajectories\" block)\n",
    "* Written to work with 3D model output (can be modified to work with 2D output)\n",
    "* Will work with or without terrain\n",
    "* Initial location, number, and density of parcels can be easily specified in \"Initialize Parcels\" block\n",
    "* Uses xarray and Dask to distribute memory and calculation across multiple processors\n",
    "* With modifications, can be used with WRF output (several others have already done so)\n",
    "* Comments that say \"set by user\" are specific to model output and desired trajectories"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "from scipy import interpolate\n",
    "import time\n",
    "\n",
    "%config IPCompleter.greedy=True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in CM1 Output\n",
    "\n",
    "* User must insert path to data\n",
    "    * If model output is one file use ***xr.open_dataset***\n",
    "    * If model output is in multiple files use ***xr.open_mfdataset***\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Use xarray to open model output and specify chunking if data set is large (set by user)\n",
    "ds = xr.open_dataset('/uufs/chpc.utah.edu/common/home/steenburgh-group8/tom/cm1/output/12ms_2000m_tug.nc', chunks={'nk': 4})\n",
    "\n",
    "#Get model output dimensions\n",
    "num_x = ds.nx #Number of gridpoints in x\n",
    "num_y = ds.ny #Number of gridpoints in y\n",
    "num_z = ds.nz #Number of gridpoints in z\n",
    "\n",
    "x = np.arange(0,num_x,1)\n",
    "y = np.arange(0,num_y,1)\n",
    "z = np.arange(0,num_z,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create Dask Cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Option to use multiple processors and threads (set by user)\n",
    "from dask.distributed import Client, LocalCluster\n",
    "c = LocalCluster(n_workers=10, threads_per_worker=4)\n",
    "client = Client(c)\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initialize Parcels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "User must enter desired trajectory characteristics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Number of parcels in vertical (can be more than number of vertical levels; set by user) \n",
    "num_seeds_z = 10\n",
    "\n",
    "#Number of parcels in y (set by user) \n",
    "num_seeds_y = 10\n",
    "\n",
    "#Number of time steps to run trajectories back (set by user) \n",
    "time_steps = 20\n",
    "\n",
    "#Time step to start backward trajectories at (set by user) \n",
    "start_time_step = 210\n",
    "\n",
    "#Variable to record at each parcel's location throughout trajectory (code can be easily modified to add more; set by user) \n",
    "var_name1 = 'th'\n",
    "\n",
    "#Set as 'Y' or 'N' for 'yes' or 'no' if the u, v, and w model output is on the staggered grid \n",
    "#(unless you have interpolated u, v, and w to the scalar grid, they are most likely on the staggered grid (set by user)\n",
    "staggered = 'N'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Model output info**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Horizontal resolution of model output (meters)\n",
    "hor_resolution = (ds.xf[1].values-ds.xf[0].values)*1000\n",
    "\n",
    "#Vertical resolution of model output (meters). Changes in x and y, if there is terrain, and z, if grid is stretched.\n",
    "vert_resolution = ds.zh[0,1:,:,:].values-ds.zh[0,:-1,:,:].values \n",
    "                  \n",
    "#Model output time step length (seconds)\n",
    "time_step_length = (ds.time[1].values - ds.time[0].values)/np.timedelta64(1, 's')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create empty arrays to store x, y, and z positions of parcels**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xpos = np.zeros((time_steps, num_seeds_z, num_seeds_y)) #x-location (grid points on staggered grid)\n",
    "ypos = np.zeros((time_steps, num_seeds_z, num_seeds_y)) #y-location (grid points on staggered grid)\n",
    "zpos = np.zeros((time_steps, num_seeds_z, num_seeds_y)) #z-location (grid points on staggered grid)\n",
    "zpos_heightASL = np.zeros((time_steps, num_seeds_z, num_seeds_y)) #Height above sea level (meters)\n",
    "zpos_vert_res = np.zeros((time_steps, num_seeds_z, num_seeds_y)) #Vertical grid spacing at parcel location (meters)\n",
    "variable1 = np.zeros((time_steps, num_seeds_z, num_seeds_y)) #User specified variable to track"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Initial location of parcels in gridpoints, specifically on the scalar grid (set by user). Initializes an array of parcels in the the y-z domain (modification necessary for x-dimension or 3D array of parcels)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#x-position\n",
    "xpos[0,:,:] = 1000 #This example initializes all seeds at same x-position (1000th x-grpt, set by user)\n",
    "\n",
    "#y-position   \n",
    "for i in range(num_seeds_y):\n",
    "    ypos[0,:,i] = 4*i #This example initializes seeds evenly in y-dimension (0th, 4th, 8th, etc. y-grpt; set by user)\n",
    "\n",
    "#z-position\n",
    "for i in range(num_seeds_z):\n",
    "    zpos[0,i,:] = i #This example initializes seeds evenly starting in z-dimension (0th, 1st, 2nd, etc., z-grpt; set by user)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determine Initial Height of Parcels Above Sea Level\n",
    "Use the height of the models levels (meters above sea level) to convert from terrain following grid points to height above seal level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get height of surface\n",
    "try:\n",
    "    zs = ds.zs[0,:,:].values\n",
    "    print('Output has terrain')\n",
    "except: \n",
    "    zs = np.zeros((ds.ny,ds.nx))\n",
    "    print('Output does not have terrain')\n",
    "\n",
    "#Get height of vertical coordinates (scalar grid)\n",
    "try:\n",
    "    zh = ds.zh[0,:,:,:].values\n",
    "    print('Output has terrain')\n",
    "except:\n",
    "    zh1d = (ds.z[:].values)*1000\n",
    "    zh2d = np.repeat(zh1d,ds.ny, axis = 0).reshape(ds.nz, ds.ny)\n",
    "    zh = np.repeat(zh2d,ds.nx, axis = 0).reshape(ds.nz, ds.ny, ds.nx)\n",
    "    print('Output does not have terrain')\n",
    "\n",
    "#Create list of initial coordinates to get height\n",
    "xloc = (xpos[0,:,:]).flatten()\n",
    "yloc = (ypos[0,:,:]).flatten()\n",
    "zloc = (zpos[0,:,:]).flatten()\n",
    "coord_height = []\n",
    "for i in range(len(xloc)):\n",
    "    coord_height.append((zloc[i], yloc[i], xloc[i]))\n",
    "\n",
    "#Get the actual inital height of the parcels in meters above sea level\n",
    "zpos_heightASL[0,:,:] = np.reshape(interpolate.interpn((z,y,x), zh, coord_height, method='linear', bounds_error=False, fill_value= 0), (num_seeds_z, num_seeds_y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate Trajectories\n",
    "User must change model output varibale names\n",
    "\n",
    "The method used to calculate trajectories is a second-order semi-implicit discretization in space and time. The method is described in depth in section 2.1 of Miltenberger et al. 2013 (https://www.geosci-model-dev.net/6/1989/2013/gmd-6-1989-2013.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Loop over all time steps and compute trajectory\n",
    "for t in range(time_steps-1):\n",
    "    \n",
    "    start = time.time() #Timer\n",
    "    \n",
    "    ##########################################################################################################\n",
    "    ##########################################################################################################   \n",
    "    ##################### Get data for 'first guess' step of integration scheme ##############################\n",
    "    ##########################################################################################################\n",
    "    ##########################################################################################################\n",
    "    \n",
    "    #Get model data (set by user)\n",
    "    u = ds.uinterp[start_time_step-t,:,:,:].values\n",
    "    v = ds.vinterp[start_time_step-t,:,:,:].values\n",
    "    w = ds.winterp[start_time_step-t,:,:,:].values\n",
    "    var1 = getattr(ds,var_name1)[start_time_step-t,:,:,:].values\n",
    "        \n",
    "        \n",
    "        \n",
    "    ############## Generate coordinates for interpolations ###############\n",
    "\n",
    "    #x, y, and z on staggered and scalar grids\n",
    "    xloc = np.copy(xpos[t,:,:]).flatten()\n",
    "    xloc_stag = np.copy(xpos[t,:,:]+0.5).flatten()\n",
    "    yloc = np.copy(ypos[t,:,:]).flatten()\n",
    "    yloc_stag = np.copy(ypos[t,:,:]+0.5).flatten()\n",
    "    zloc = np.copy(zpos[t,:,:]).flatten()\n",
    "    zloc_stag = np.copy(zpos[t,:,:]+0.5).flatten()\n",
    "\n",
    "    #If u, v, and w are staggered, generate three staggered sets of coordinates:\n",
    "    #    1) u-grid (staggered in x)\n",
    "    #    2) v-grid (staggered in y)\n",
    "    #    3) w-grid (staggered in z)\n",
    "    \n",
    "    if staggered == 'Y':\n",
    "        coord_u = []\n",
    "        coord_v = []\n",
    "        coord_w = []\n",
    "        for i in range(len(xloc)):\n",
    "            coord_u.append((zloc[i], yloc[i], xloc_stag[i])) \n",
    "            coord_v.append((zloc[i], yloc_stag[i], xloc[i])) \n",
    "            coord_w.append((zloc_stag[i], yloc[i], xloc[i])) \n",
    "    \n",
    "    #If not, generate scalar coordinates\n",
    "    else: \n",
    "        coord_u = []\n",
    "        coord_v = []\n",
    "        coord_w = []\n",
    "        for i in range(len(xloc)):\n",
    "            coord_u.append((zloc[i], yloc[i], xloc[i])) \n",
    "            coord_v.append((zloc[i], yloc[i], xloc[i])) \n",
    "            coord_w.append((zloc[i], yloc[i], xloc[i])) \n",
    "    \n",
    "    #Scalar coordinates for all other variables\n",
    "    coord = []\n",
    "    for i in range(len(xloc)):\n",
    "        coord.append((zloc[i], yloc[i], xloc[i])) \n",
    "    \n",
    "    ##########################################################################################################   \n",
    "    ########################## Integrate 'first guess' of parcel's new location ##############################\n",
    "    ##########################################################################################################   \n",
    "\n",
    "    \n",
    "    #########################   Calc 'first guess' new xpos in grdpts   #######################################\n",
    "    dx_0 = np.reshape(interpolate.interpn((z,y,x), u, coord_u, method='linear', bounds_error=False, fill_value=np.nan)*time_step_length/hor_resolution, (num_seeds_z, num_seeds_y))\n",
    "    xpos_1 = xpos[t,:,:] - dx_0\n",
    "\n",
    "    #########################   Calc 'first guess' new ypos in grdpts   #######################################\n",
    "    dy_0 = np.reshape(interpolate.interpn((z,y,x), v, coord_v, method='linear', bounds_error=False, fill_value=np.nan)*time_step_length/hor_resolution, (num_seeds_z, num_seeds_y))\n",
    "    ypos_1 = ypos[t,:,:] - dy_0\n",
    "\n",
    "    #########################   Calc 'first guess' new zpos in meters above sea level ######################################\n",
    "    dz_0 = np.reshape(interpolate.interpn((z,y,x), w, coord_w, method='linear', bounds_error=False, fill_value= 0)*time_step_length, (num_seeds_z, num_seeds_y))\n",
    "    zpos_heightASL_1 = zpos_heightASL[t,:,:] - dz_0\n",
    "    \n",
    "    ############# Convert zpos from meters above sea level to gridpts abve surface for interpolation #########\n",
    "    #Get vertical grid spacing at each parcel's location\n",
    "    zpos_vert_res[t,:,:] = np.reshape(interpolate.interpn((z[:-1],y,x), vert_resolution, coord, method='linear', bounds_error=False, fill_value= np.nan), (num_seeds_z, num_seeds_y))\n",
    "\n",
    "\n",
    "    #Calculate change in surface height and change in parcel height\n",
    "    xloc = np.copy(xpos[t,:,:]).flatten()\n",
    "    yloc = np.copy(ypos[t,:,:]).flatten()\n",
    "    coord_zs1 = []\n",
    "    for i in range(len(xloc)):\n",
    "        coord_zs1.append((yloc[i], xloc[i]))\n",
    "        \n",
    "    xloc = np.copy(xpos_1).flatten()\n",
    "    yloc = np.copy(ypos_1).flatten()\n",
    "    coord_zs2 = []\n",
    "    for i in range(len(xloc)):\n",
    "        coord_zs2.append((yloc[i], xloc[i]))\n",
    "    \n",
    "    #Change in surface height over last timestep\n",
    "    zs1 = interpolate.interpn((y,x), zs, coord_zs1, method='linear', bounds_error=False, fill_value= np.nan)\n",
    "    zs2 = interpolate.interpn((y,x), zs, coord_zs2, method='linear', bounds_error=False, fill_value= np.nan)\n",
    "    zs_change = zs2-zs1\n",
    "    \n",
    "    #Change in parcel height over last times step\n",
    "    zpos_heightASL_change = zpos_heightASL_1.flatten()-zpos_heightASL[t,:,:].flatten()\n",
    "    \n",
    "    #Calculate zpos in grdpts above surface\n",
    "    zpos_1 = zpos[t,:,:] + np.reshape((zpos_heightASL_change - zs_change)/zpos_vert_res[t,:,:].flatten(), (num_seeds_z, num_seeds_y))\n",
    "    \n",
    "    #Prevent parcels from going into the ground\n",
    "    zpos_heightASL_1 = zpos_heightASL_1.clip(min=0)\n",
    "    zpos_1 = zpos_1.clip(min=0)\n",
    "    \n",
    "    \n",
    "    #Calculate value of variable at each parcel's location\n",
    "    variable1[t,:,:] = np.reshape(interpolate.interpn((z,y,x), var1, coord, method = 'linear', bounds_error=False, fill_value= np.nan), (num_seeds_z, num_seeds_y)) \n",
    "\n",
    "\n",
    "    ##########################################################################################################\n",
    "    ##########################################################################################################   \n",
    "    ##################### Get data for 'correction' step of integration scheme ###############################\n",
    "    ##########################################################################################################\n",
    "    ##########################################################################################################\n",
    "    \n",
    "    \n",
    "    \n",
    "    #Get model data for next time step (set by user)\n",
    "    t = t + 1\n",
    "    u = ds.uinterp[start_time_step-t,:,:,:].values\n",
    "    v = ds.vinterp[start_time_step-t,:,:,:].values\n",
    "    w = ds.winterp[start_time_step-t,:,:,:].values\n",
    "    t = t - 1\n",
    "        \n",
    "        \n",
    "    ############## Generate coordinates for interpolations ###############\n",
    "\n",
    "    #x, y, and z on staggered and scalar grids\n",
    "    xloc = np.copy(xpos_1).flatten()\n",
    "    xloc_stag = np.copy(xpos_1+0.5).flatten()\n",
    "    yloc = np.copy(ypos_1).flatten()\n",
    "    yloc_stag = np.copy(ypos_1+0.5).flatten()\n",
    "    zloc = np.copy(zpos_1).flatten()\n",
    "    zloc_stag = np.copy(zpos_1+0.5).flatten()\n",
    "\n",
    "    #If u, v, and w are staggered, generate three staggered sets of coordinates:\n",
    "    #    1) u-grid (staggered in x)\n",
    "    #    2) v-grid (staggered in y)\n",
    "    #    3) w-grid (staggered in z)\n",
    "    \n",
    "    if staggered == 'Y':\n",
    "        coord_u = []\n",
    "        coord_v = []\n",
    "        coord_w = []\n",
    "        for i in range(len(xloc)):\n",
    "            coord_u.append((zloc[i], yloc[i], xloc_stag[i])) \n",
    "            coord_v.append((zloc[i], yloc_stag[i], xloc[i])) \n",
    "            coord_w.append((zloc_stag[i], yloc[i], xloc[i])) \n",
    "    \n",
    "    #If not, generate scalar coordinates\n",
    "    else: \n",
    "        coord_u = []\n",
    "        coord_v = []\n",
    "        coord_w = []\n",
    "        for i in range(len(xloc)):\n",
    "            coord_u.append((zloc[i], yloc[i], xloc[i])) \n",
    "            coord_v.append((zloc[i], yloc[i], xloc[i])) \n",
    "            coord_w.append((zloc[i], yloc[i], xloc[i])) \n",
    "    \n",
    "    #Scalar coordinates for all other variables\n",
    "    coord = []\n",
    "    for i in range(len(xloc)):\n",
    "        coord.append((zloc[i], yloc[i], xloc[i])) \n",
    "        \n",
    "    \n",
    "    ##########################################################################################################   \n",
    "    ########################## Integrate 'correction' of parcel's new location ###############################\n",
    "    ##########################################################################################################   \n",
    "\n",
    "    \n",
    "    #########################   Calc 'correction' new xpos in grdpts   #######################################\n",
    "    dx_1 = np.reshape(interpolate.interpn((z,y,x), u, coord_u, method='linear', bounds_error=False, fill_value=np.nan)*time_step_length/hor_resolution, (num_seeds_z, num_seeds_y))\n",
    "    xpos[t+1,:,:] = xpos[t,:,:] - (dx_0 + dx_1)/2\n",
    "\n",
    "    #########################   Calc 'correction' new ypos in grdpts   #######################################\n",
    "    dy_1 = np.reshape(interpolate.interpn((z,y,x), v, coord_v, method='linear', bounds_error=False, fill_value=np.nan)*time_step_length/hor_resolution, (num_seeds_z, num_seeds_y))\n",
    "    ypos[t+1,:,:] = ypos[t,:,:] - (dy_0 + dy_1)/2\n",
    "\n",
    "    #########################   Calc 'correction' new zpos in meters above sea level ######################################\n",
    "    dz_1 = np.reshape(interpolate.interpn((z,y,x), w, coord_w, method='linear', bounds_error=False, fill_value= 0)*time_step_length, (num_seeds_z, num_seeds_y))\n",
    "    zpos_heightASL[t+1,:,:] = zpos_heightASL[t,:,:] - (dz_0 + dz_1)/2\n",
    "    \n",
    "    \n",
    "    \n",
    "    ############# Convert zpos from meters above sea level to gridpts abve surface for interpolation #########\n",
    "    #Get vertical grid spacing at each parcel's location\n",
    "    zpos_vert_res[t,:,:] = np.reshape(interpolate.interpn((z[:-1],y,x), vert_resolution, coord, method='linear', bounds_error=False, fill_value= np.nan), (num_seeds_z, num_seeds_y))\n",
    "\n",
    "    \n",
    "    #Calculate change in surface height and change in parcel height\n",
    "    xloc = np.copy(xpos[t,:,:]).flatten()\n",
    "    yloc = np.copy(ypos[t,:,:]).flatten()\n",
    "    coord_zs1 = []\n",
    "    for i in range(len(xloc)):\n",
    "        coord_zs1.append((yloc[i], xloc[i]))\n",
    "        \n",
    "    xloc = np.copy(xpos[t+1,:,:]).flatten()\n",
    "    yloc = np.copy(ypos[t+1,:,:]).flatten()\n",
    "    coord_zs2 = []\n",
    "    for i in range(len(xloc)):\n",
    "        coord_zs2.append((yloc[i], xloc[i]))\n",
    "    \n",
    "    #Change in surface height over last timestep\n",
    "    zs1 = interpolate.interpn((y,x), zs, coord_zs1, method='linear', bounds_error=False, fill_value= np.nan)\n",
    "    zs2 = interpolate.interpn((y,x), zs, coord_zs2, method='linear', bounds_error=False, fill_value= np.nan)\n",
    "    zs_change = zs2-zs1\n",
    "    \n",
    "    #Change in parcel height over last timestep\n",
    "    zpos_heightASL_change = zpos_heightASL[t+1,:,:].flatten()-zpos_heightASL[t,:,:].flatten()\n",
    "    \n",
    "    #Calculate zpos in grdpts above surface\n",
    "    zpos[t+1,:,:] = zpos[t,:,:] + np.reshape((zpos_heightASL_change - zs_change)/zpos_vert_res[t,:,:].flatten(), (num_seeds_z, num_seeds_y))\n",
    "    ##########################################################################################################\n",
    "    \n",
    "    \n",
    "    #Prevent parcels from going into the ground\n",
    "    zpos = zpos.clip(min=0)\n",
    "    zpos_heightASL = zpos_heightASL.clip(min=0)    \n",
    "        \n",
    "\n",
    "    #Timer\n",
    "    stop = time.time()\n",
    "    print(\"Integration {:01d} took {:.2f} seconds\".format(t, stop-start))\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get variable data for final time step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load variable data\n",
    "t = time_steps-1\n",
    "var1 = getattr(ds,var_name1)[start_time_step-t,:,:,:].values\n",
    "\n",
    "#Get get x, y, and z positions from scalar grid\n",
    "xloc = np.copy(xpos[t,:,:]).flatten()\n",
    "yloc = np.copy(ypos[t,:,:]).flatten()\n",
    "zloc = np.copy(zpos[t,:,:]).flatten()\n",
    "coord = []\n",
    "for i in range(len(xloc)):\n",
    "    coord.append((zloc[i], yloc[i], xloc[i])) \n",
    "\n",
    "#Variables\n",
    "variable1[t,:,:] = np.reshape(interpolate.interpn((z,y,x), var1, coord, method = 'linear', bounds_error=False, fill_value= np.nan), (num_seeds_z, num_seeds_y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save Trajectory Data\n",
    "The x, y, and z positions and user-specified variable values are saved in 3D numpy arrays. The first dimension is time and the other two are the positions and values of variables of all the parcels at that specifc time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('trajectory_arrays/xpos', xpos)\n",
    "np.save('trajectory_arrays/ypos', ypos)\n",
    "np.save('trajectory_arrays/zpos', zpos_heightASL)\n",
    "np.save('trajectory_arrays/%s' %var_name1, variable1)"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}