From 3a6508ff24c9f3db7fd8ec8d2260f128fa2a5253 Mon Sep 17 00:00:00 2001
From: Peter La Follette <peterlafollette@Peters-MacBook-Pro-2.local>
Date: Sun, 3 Nov 2024 13:09:12 -0800
Subject: [PATCH] adding notebooks for animating LGARTO synthetic tests

---
 tests/lgarto_profile_animation_synth_1.ipynb | 890 +++++++++++++++++++
 tests/lgarto_profile_animation_synth_2.ipynb | 885 ++++++++++++++++++
 tests/lgarto_profile_animation_synth_3.ipynb | 890 +++++++++++++++++++
 tests/lgarto_profile_animation_synth_4.ipynb | 890 +++++++++++++++++++
 4 files changed, 3555 insertions(+)
 create mode 100644 tests/lgarto_profile_animation_synth_1.ipynb
 create mode 100644 tests/lgarto_profile_animation_synth_2.ipynb
 create mode 100644 tests/lgarto_profile_animation_synth_3.ipynb
 create mode 100644 tests/lgarto_profile_animation_synth_4.ipynb

diff --git a/tests/lgarto_profile_animation_synth_1.ipynb b/tests/lgarto_profile_animation_synth_1.ipynb
new file mode 100644
index 0000000..1d90f8f
--- /dev/null
+++ b/tests/lgarto_profile_animation_synth_1.ipynb
@@ -0,0 +1,890 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8b26a1c1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from matplotlib import pyplot as plt\n",
+    "import matplotlib.animation as animation\n",
+    "from matplotlib.animation import FuncAnimation, PillowWriter \n",
+    "import sys, os\n",
+    "import imageio, re\n",
+    "import importlib as imp\n",
+    "import copy\n",
+    "import pandas as pd\n",
+    "\n",
+    "from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes\n",
+    "from mpl_toolkits.axes_grid1.inset_locator import mark_inset\n",
+    "import matplotlib.dates as mdates\n",
+    "from matplotlib.patches import Rectangle\n",
+    "\n",
+    "import pickle\n",
+    "\n",
+    "from datetime import datetime\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2ede7ad3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "layer_depth = [0.0, 10.0, 10.0+45.0, 10.0+45.0+45.0]\n",
+    "#51.0,99.0,53.0[cm]\n",
+    "time_step_for_plotting = 1\n",
+    "theta_e_vec = [0.41,0.43,0.45] #edit for your scenario \n",
+    "max_depth = layer_depth[-1]\n",
+    "# 36.0,111.0,107.0\n",
+    "\n",
+    "synth_case = str(1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3a9f5ba8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "c_path_layers = 'outputs/synthetic'+synth_case + '_LGARTO/data_layers.csv'\n",
+    "sim_case = 'outputs/synthetic'+synth_case+'_LGARTO/data_variables.csv' \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "72f8e45a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# hydrus_path_layers = '/Users/peterlafollette/Desktop/ngen_test_2/ngen/extern/LGARTO-C/tests/HYDRUS_outputs/Geneva/Nod_Inf.csv'\n",
+    "# data_sim_hydrus = pd.read_csv(hydrus_path_layers, sep='\\s+')\n",
+    "\n",
+    "# # data_sim_hydrus = data_sim_hydrus[0:100]\n",
+    "# # data_sim_hydrus = data_sim_hydrus[101:201]\n",
+    "# # data_sim_hydrus = data_sim_hydrus[202:302]\n",
+    "# data_sim_hydrus = data_sim_hydrus[202:]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "efb644ba",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5ef98a57",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data_sim_c = np.loadtxt(c_path_layers,dtype='str')\n",
+    "# data_sim_py = data_sim_c #setting these equal for now to get code working -- attempting to remove data_sim_py from code \n",
+    "\n",
+    "Depth_d = []\n",
+    "Theta_d = []\n",
+    "Head_d = []\n",
+    "Dat_u = []\n",
+    "Dat_v = []\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "29115339",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f820a6b0",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f9ff759c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "count = 0\n",
+    "for u in data_sim_c:\n",
+    "    u1 = u[1:-1].split('|')\n",
+    "    val_u = []\n",
+    "\n",
+    "    for u2 in u1:\n",
+    "        u3 = u2[1:-1].split(',')\n",
+    "        u4 = [float(u3[0]), float(u3[1]), float(u3[4])]\n",
+    "        val_u.append(u4)\n",
+    "        \n",
+    "    count = count + 1\n",
+    "    Dat_u.append(val_u)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "276f3f5c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def data_gen1():\n",
+    "    for u in data_sim_c[::12]:\n",
+    "        v1 = u[1:-1].split('|')\n",
+    "        val_c = []\n",
+    "        \n",
+    "        for v2 in v1:\n",
+    "            v3 = v2[1:-1].split(',')\n",
+    "            v4 = [float(v3[0]), float(v3[1]), float(v3[4])]\n",
+    "            val_c.append(v4)\n",
+    "        \n",
+    "        yield val_c\n",
+    "        \n",
+    "time = np.arange(0,len(Dat_u)/12.,1.) # 1/12 = 5/60 (min/hr)\n",
+    "time = [round (t,3) for t in time]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cd9c78e9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def data_gen2():\n",
+    "    c = -1\n",
+    "    #print (len(time))\n",
+    "    for u in data_sim_c[::12]:\n",
+    "        u1 = u[1:-1].split('|')\n",
+    "#         v1 = v[1:-1].split('],')\n",
+    "        val_u = []\n",
+    "#         val_v = []\n",
+    "        for u2 in u1:\n",
+    "            u3 = u2[1:-1].split(',')\n",
+    "            u4 = [float(u3[0]), float(u3[1]), float(u3[4])]\n",
+    "            val_u.append(u4)\n",
+    "\n",
+    "#         for v2 in v1:\n",
+    "#             v3 = v2[1:-1].split(',')\n",
+    "#             v3[0] = v3[0].replace('[', '')\n",
+    "#             depth = float(v3[0]) + layer_depth[int(v3[2])]\n",
+    "#             v4 = [depth, round(float(v3[1]),6), float(v3[2])]\n",
+    "#             #v4 = [float(v3[0]), round(float(v3[1]),6), float(v3[3])]\n",
+    "#             val_v.append(v4)\n",
+    "        c = c +1\n",
+    "        \n",
+    "        yield time[c], val_u\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b0e3bb2a",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b982e58a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %matplotlib notebook\n",
+    "fig, axs = plt.subplots(2,1,figsize=(6*1.3,4*1.3))\n",
+    "fig.tight_layout(pad=1.8)\n",
+    "len_simulation = len(Dat_u)\n",
+    "\n",
+    "# for i in range(time_step_for_plotting-1,time_step_for_plotting):\n",
+    "i = time_step_for_plotting-1\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "###plots theta vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "axs[0].set_ylim(0,layer_depth[-1])\n",
+    "axs[0].invert_yaxis()\n",
+    "axs[0].set_xlim(0,0.5)\n",
+    "\n",
+    "axs[0].set_ylabel('Depth (cm)',fontsize=12)\n",
+    "axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "        #print (k)\n",
+    "    if (k!=layer_depth[-1]):\n",
+    "        axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "    if ((count==1)|(count==2)):\n",
+    "        axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "max_val_psi = 0\n",
+    "for row in d:\n",
+    "    if (max_val_psi<max(row)):\n",
+    "        max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1\n",
+    "\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "###plots psi vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "axs[1].set_ylim(0,layer_depth[-1])\n",
+    "axs[1].invert_yaxis()\n",
+    "axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "axs[1].set_ylabel('Depth (cm)',fontsize=12)\n",
+    "axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", \"theta_e\", \"capillary head = -depth\"],loc=5)\n",
+    "\n",
+    "leg.legendHandles[1].set_color('blue')\n",
+    "leg.legendHandles[2].set_color('purple')\n",
+    "\n",
+    "plt.savefig(\"profile_xxx.pdf\")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1f5227f1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "i = time_step_for_plotting-1\n",
+    "i=200"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6e7f8504",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %matplotlib notebook\n",
+    "fig, axs = plt.subplots(4,1,figsize=(4.5*1.7,6*1.7))\n",
+    "fig.tight_layout(pad=2.5)\n",
+    "len_simulation = len(Dat_u)\n",
+    "\n",
+    "# for i in range(time_step_for_plotting-1,time_step_for_plotting):\n",
+    "# i = time_step_for_plotting-1\n",
+    "# i=3000\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "###plots theta vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "axs[0].set_ylim(0,layer_depth[-1])\n",
+    "axs[0].invert_yaxis()\n",
+    "axs[0].set_xlim(0,0.5)\n",
+    "\n",
+    "axs[0].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "        #print (k)\n",
+    "    if (k!=layer_depth[-1]):\n",
+    "        axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "    if ((count==1)|(count==2)):\n",
+    "        axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "max_val_psi = 0\n",
+    "for row in d:\n",
+    "    if (max_val_psi<max(row)):\n",
+    "        max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1\n",
+    "\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "###plots psi vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "axs[1].set_ylim(0,layer_depth[-1])\n",
+    "axs[1].invert_yaxis()\n",
+    "axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "axs[1].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", r\"$\\theta_{e}$\"],loc=3)\n",
+    "\n",
+    "leg.legendHandles[1].set_color('blue')\n",
+    "\n",
+    "\n",
+    "\n",
+    "leg2 = axs[1].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"capillary head = -depth\"],loc=3)\n",
+    "leg2.legendHandles[0].set_color('purple')\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "c_path_vars = sim_case\n",
+    "\n",
+    "\n",
+    "num_steps = 864 + 1\n",
+    "\n",
+    "data_c_vars = pd.read_csv(c_path_vars)\n",
+    "data_c_vars = data_c_vars[0:num_steps]\n",
+    "    \n",
+    "start = pd.Timestamp(\"2020-10-01 00:00:00\")\n",
+    "end = pd.Timestamp(\"2020-10-03 00:00:00\")\n",
+    "times = pd.date_range(start=start,end=end, periods=num_steps) #freq='min'\n",
+    "times_bmi = pd.date_range(start=start,end=end, periods=num_steps)\n",
+    "times_bmi = []\n",
+    "for item in data_c_vars['Time']:\n",
+    "    times_bmi.append(datetime.strptime(item, '%Y-%m-%d %H:%M:%S'))\n",
+    "#     datetime.datetime.strptime(data_c_vars_LGARTO['Time'][i], '%Y-%m-%d %H:%M:%S') - datetime.datetime.strptime(data_c_vars_LGARTO['Time'][0],'%Y-%m-%d %H:%M:%S')\n",
+    "# times_bmi = times_bmi[0:len(times_bmi)-1]\n",
+    "\n",
+    "xfmt = mdates.DateFormatter('%Y/%m/%d')\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "m_to_mm = 1000\n",
+    "\n",
+    "\n",
+    "vars_bmi = ['soil_storage', 'surface_runoff', 'infiltration', 'actual_evapotranspiration', 'mass_balance', 'percolation', 'precipitation', 'potential_evapotranspiration']\n",
+    "\n",
+    "# vars_hydrus = [')     Volume', 'sum(RunOff', 'sum(Infil', ')   sum(vRoot', ')    sum(vBot']\n",
+    "m_to_mm = 1000\n",
+    "m_to_cm = 100\n",
+    "timestep = 300/3600\n",
+    "y_labels = ['Storage (cm)', 'Runoff (cm)', 'Infiltration (cm)', 'Actual ET (cm)', 'Mass balance (cm)', 'Recharge (cm)']\n",
+    "\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm),'b',label='C (bmi)',color='orange')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[1]],'k',linestyle='dashed',label='Hydrus',color='orange')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm),'b',label='C (bmi)',color='blue')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[2]],'k',linestyle='dashed',label='Hydrus',color='blue')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm),'b',label='C (bmi)',color='green')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[3]],'k',linestyle='dashed',label='Hydrus',color='green')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm),'b',label='C (bmi)',color='purple')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[4]],'k',linestyle='dashed',label='Hydrus',color='purple')\n",
+    "\n",
+    "axs[2].vlines(x=times_bmi[i], ymin=-100, ymax=1000, colors='red')\n",
+    "\n",
+    "#axs[2].plot(times_bmi, data_c_vars[vars_bmi[5]] * m_to_cm*12, label='C (bmi)', color='green')\n",
+    "\n",
+    "\n",
+    "leg = axs[2].legend()#(handles=[runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS],labels=[\"LGARTO\", \"HYDRUS\", \"theta_e\", \"capillary head = -depth\"])\n",
+    "\n",
+    "mylabels = ['LGARTO runoff','LGARTO infiltration','LGARTO AET','LGARTO recharge', 'current time']\n",
+    "# axs[2].legend(labels=mylabels)\n",
+    "axs[2].legend(labels=mylabels, ncol=3, bbox_to_anchor=(0.5, -0.0), loc='lower center')\n",
+    "\n",
+    "\n",
+    "axs[2].text(x=times_bmi[550], y=50, s='synth '+synth_case, fontsize = 16)\n",
+    "\n",
+    "axs[2].set_ylabel(\"cumulative \\n output flux (cm)\",fontsize=16)\n",
+    "axs[2].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "axs[2].set_ylim(min(np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm))-6, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm))))\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm),'b',label='C (bmi)',color='black')\n",
+    "axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm),'b',label='C (bmi)',color='gold')\n",
+    "\n",
+    "axs[3].set_ylim(0, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm))))\n",
+    "\n",
+    "axs[3].set_ylabel(\"cumulative \\n forcing flux (cm)\",fontsize=16)\n",
+    "\n",
+    "axs[3].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "axs[3].vlines(x=times_bmi[i], ymin=-100, ymax=1000, colors='red')\n",
+    "\n",
+    "# axs[3].text(x=times_bmi[600], y=32+5, s='current time: '+str(round(i/12,2))+' h', fontsize = 12)\n",
+    "# axs[3].text(x=times_bmi[600], y=15-1, s=str(times_bmi[i]), fontsize = 12)\n",
+    "\n",
+    "\n",
+    "mylabels_forcing = ['precipitation', 'PET', 'current time']\n",
+    "axs[3].legend(labels=mylabels_forcing)\n",
+    "\n",
+    "fig.subplots_adjust(left=0.14)\n",
+    "\n",
+    "plt.savefig(\"profile_xxx.pdf\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "625de2c5",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "587eb41f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c4c1c574",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f4acea80",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4991dd06",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8d85ad43",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a3b54e38",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d = Dat_u[0]\n",
+    "max_val_psi = 0\n",
+    "\n",
+    "for d in Dat_u:\n",
+    "    for row in d:\n",
+    "        if (max_val_psi<max(row)):\n",
+    "            max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1 + 10\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ec71d705",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "max_val_psi"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4a9acecd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def init():\n",
+    "    for ax in axs:\n",
+    "        ax.clear()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dcd556ab",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def animate(frame):\n",
+    "    \n",
+    "    for ax in axs:\n",
+    "        ax.clear()\n",
+    "    \n",
+    "    d = Dat_u[frame]\n",
+    "    num_wf = len(d)\n",
+    "    \n",
+    "    fig.subplots_adjust(left=0.14)\n",
+    "\n",
+    "    lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "    ###plots theta vs depth\n",
+    "    for j in range(num_wf):\n",
+    "        if j ==0 :\n",
+    "            axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "        else:\n",
+    "            axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "            lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "    axs[0].set_ylim(0,layer_depth[-1])\n",
+    "    axs[0].invert_yaxis()\n",
+    "    axs[0].set_xlim(0,0.5)\n",
+    "\n",
+    "    axs[0].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "    axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "    count = 0\n",
+    "    for k in layer_depth:\n",
+    "            #print (k)\n",
+    "        if (k!=layer_depth[-1]):\n",
+    "            axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "        if ((count==1)|(count==2)):\n",
+    "            axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "        if (k==layer_depth[-1]):\n",
+    "            axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "        else:\n",
+    "            axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "        count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "#     max_val_psi = 0\n",
+    "#     for row in d:\n",
+    "#         if (max_val_psi<max(row)):\n",
+    "#             max_val_psi = max(row)\n",
+    "\n",
+    "#     max_val_psi = max_val_psi*0.1\n",
+    "#     max_val_psi = 400\n",
+    "\n",
+    "    d = Dat_u[frame]\n",
+    "    num_wf = len(d)\n",
+    "\n",
+    "    ###plots psi vs depth\n",
+    "    for j in range(num_wf):\n",
+    "        if j ==0 :\n",
+    "            axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "        else:\n",
+    "            axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "            axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "    axs[1].set_ylim(0,layer_depth[-1])\n",
+    "    axs[1].invert_yaxis()\n",
+    "    axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "    axs[1].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "    axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "    count = 0\n",
+    "    for k in layer_depth:\n",
+    "\n",
+    "        if (k==layer_depth[-1]):\n",
+    "            axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "        else:\n",
+    "            axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "        count = count + 1\n",
+    "\n",
+    "    hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "    hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "    hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "    leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", r\"$\\theta_{e}$\"],loc=3,)\n",
+    "\n",
+    "    leg.legendHandles[1].set_color('blue')\n",
+    "#     leg.legendHandles[2].set_color('purple')\n",
+    "\n",
+    "    leg2 = axs[1].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"capillary head = -depth\"],loc=3)\n",
+    "    leg2.legendHandles[0].set_color('purple')\n",
+    "    \n",
+    "    \n",
+    "\n",
+    "\n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm),'b',label='C (bmi)',color='orange')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[1]],'k',linestyle='dashed',label='Hydrus',color='orange')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm),'b',label='C (bmi)',color='blue')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[2]],'k',linestyle='dashed',label='Hydrus',color='blue')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm),'b',label='C (bmi)',color='green')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[3]],'k',linestyle='dashed',label='Hydrus',color='green')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm),'b',label='C (bmi)',color='purple')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[4]],'k',linestyle='dashed',label='Hydrus',color='purple')\n",
+    "\n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    axs[2].vlines(x=times_bmi[frame], ymin=-100, ymax=1000, colors='red')\n",
+    "    \n",
+    "    leg = axs[2].legend()#(handles=[runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS],labels=[\"LGARTO\", \"HYDRUS\", \"theta_e\", \"capillary head = -depth\"])\n",
+    "\n",
+    "    mylabels = ['LGARTO runoff','LGARTO infiltration','LGARTO AET','LGARTO recharge', 'current time']\n",
+    "    axs[2].legend(labels=mylabels, ncol=3, bbox_to_anchor=(0.5, -0.0), loc='lower center')\n",
+    "\n",
+    "    axs[2].text(x=times_bmi[550], y=50, s='synth '+synth_case, fontsize = 16)\n",
+    "\n",
+    "    axs[2].set_ylabel(\"cumulative \\n output flux (cm)\",fontsize=16)\n",
+    "    axs[2].set_xlabel(\"date\",fontsize=12)\n",
+    "    \n",
+    "#     axs[2].plot(times_bmi, data_c_vars[vars_bmi[5]] * m_to_cm*12, label='C (bmi)', color='green')\n",
+    "\n",
+    "\n",
+    "    axs[2].set_ylim(min(np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm))-6, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm))))\n",
+    "#     axs[2].set_ylim(-0.05,0.15)\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "    axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm),'b',label='C (bmi)',color='black')\n",
+    "    axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm),'b',label='C (bmi)',color='gold')\n",
+    "\n",
+    "    axs[3].set_ylim(0, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm))))\n",
+    "\n",
+    "    axs[3].vlines(x=times_bmi[frame], ymin=-100, ymax=1000, colors='red')\n",
+    "    \n",
+    "    axs[3].text(x=times_bmi[600], y=3, s='current time: '+str(round(frame/12,2))+' h', fontsize = 12)\n",
+    "#     axs[3].text(x=times_bmi[600], y=15-1, s=str(times_bmi[frame]), fontsize = 12)\n",
+    "    \n",
+    "    mylabels_forcing = ['precipitation', 'PET', 'current time']\n",
+    "    axs[3].legend(labels=mylabels_forcing, ncol=3)\n",
+    "    \n",
+    "    axs[3].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "    axs[3].set_ylabel(\"cumulative \\n forcing flux (cm)\",fontsize=16)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3af740cd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "init()\n",
+    "\n",
+    "start_move_writing_time = datetime.now()\n",
+    "\n",
+    "ani = FuncAnimation(fig, animate, frames = 10000,  init_func=init, blit=False,repeat=False)\n",
+    "# ani.save('LGARTO_animation.mp4', dpi=600, fps=50, extra_args=['-vcodec', 'libx264'])\n",
+    "ani.save('LGARTO_animation_synth_'+synth_case+'.mp4', dpi=400, fps=48, extra_args=['-vcodec', 'libx264', '-preset', 'veryslow', '-crf', '17'])\n",
+    "\n",
+    "# ani.save('LGARTO_animation.gif', fps=50)\n",
+    "\n",
+    "end_move_writing_time = datetime.now()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5b9c4b9b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#time it took to write movie. It can be a computationally intensive thing \n",
+    "end_move_writing_time - start_move_writing_time\n",
+    "#mp4 blit=False time: 174 seconds for 200 frames\n",
+    "#then edited with 'veryslow' and crf of 17, 188 seconds for 200 frames\n",
+    "#next trying with 30 fps and 300 dpi, using the slower extra_args ... 132 seconds for 200 frames\n",
+    "#dpi=300, fps=60, slow args: 135 seconds for 200 frames, and doesn't have the freezing problem, but does kind of feel jerky. Maybe too fast, 30 fps version looked smoother. But it's a bit slow, so 45 is a good medium \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c2b2537c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f960d57e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#extra_args=['-vcodec', 'libx264']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "553247b3",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2885dd2f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e4ba966c",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "36aefd84",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7f98f67c",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f1c884db",
+   "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.8.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/tests/lgarto_profile_animation_synth_2.ipynb b/tests/lgarto_profile_animation_synth_2.ipynb
new file mode 100644
index 0000000..f65cb33
--- /dev/null
+++ b/tests/lgarto_profile_animation_synth_2.ipynb
@@ -0,0 +1,885 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8b26a1c1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from matplotlib import pyplot as plt\n",
+    "import matplotlib.animation as animation\n",
+    "from matplotlib.animation import FuncAnimation, PillowWriter \n",
+    "import sys, os\n",
+    "import imageio, re\n",
+    "import importlib as imp\n",
+    "import copy\n",
+    "import pandas as pd\n",
+    "\n",
+    "from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes\n",
+    "from mpl_toolkits.axes_grid1.inset_locator import mark_inset\n",
+    "import matplotlib.dates as mdates\n",
+    "from matplotlib.patches import Rectangle\n",
+    "\n",
+    "import pickle\n",
+    "\n",
+    "from datetime import datetime\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2ede7ad3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "layer_depth = [0.0, 10.0, 10.0+45.0, 10.0+45.0+45.0]\n",
+    "#51.0,99.0,53.0[cm]\n",
+    "time_step_for_plotting = 1\n",
+    "theta_e_vec = [0.43,0.45,0.41] #edit for your scenario \n",
+    "max_depth = layer_depth[-1]\n",
+    "# 36.0,111.0,107.0"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3a9f5ba8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "c_path_layers = 'outputs/synthetic2_LGARTO/data_layers.csv'\n",
+    "sim_case = 'outputs/synthetic2_LGARTO/data_variables.csv' \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "72f8e45a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# hydrus_path_layers = '/Users/peterlafollette/Desktop/ngen_test_2/ngen/extern/LGARTO-C/tests/HYDRUS_outputs/Geneva/Nod_Inf.csv'\n",
+    "# data_sim_hydrus = pd.read_csv(hydrus_path_layers, sep='\\s+')\n",
+    "\n",
+    "# # data_sim_hydrus = data_sim_hydrus[0:100]\n",
+    "# # data_sim_hydrus = data_sim_hydrus[101:201]\n",
+    "# # data_sim_hydrus = data_sim_hydrus[202:302]\n",
+    "# data_sim_hydrus = data_sim_hydrus[202:]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "efb644ba",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5ef98a57",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data_sim_c = np.loadtxt(c_path_layers,dtype='str')\n",
+    "# data_sim_py = data_sim_c #setting these equal for now to get code working -- attempting to remove data_sim_py from code \n",
+    "\n",
+    "Depth_d = []\n",
+    "Theta_d = []\n",
+    "Head_d = []\n",
+    "Dat_u = []\n",
+    "Dat_v = []\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "29115339",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f820a6b0",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f9ff759c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "count = 0\n",
+    "for u in data_sim_c:\n",
+    "    u1 = u[1:-1].split('|')\n",
+    "    val_u = []\n",
+    "\n",
+    "    for u2 in u1:\n",
+    "        u3 = u2[1:-1].split(',')\n",
+    "        u4 = [float(u3[0]), float(u3[1]), float(u3[4])]\n",
+    "        val_u.append(u4)\n",
+    "        \n",
+    "    count = count + 1\n",
+    "    Dat_u.append(val_u)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "276f3f5c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def data_gen1():\n",
+    "    for u in data_sim_c[::12]:\n",
+    "        v1 = u[1:-1].split('|')\n",
+    "        val_c = []\n",
+    "        \n",
+    "        for v2 in v1:\n",
+    "            v3 = v2[1:-1].split(',')\n",
+    "            v4 = [float(v3[0]), float(v3[1]), float(v3[4])]\n",
+    "            val_c.append(v4)\n",
+    "        \n",
+    "        yield val_c\n",
+    "        \n",
+    "time = np.arange(0,len(Dat_u)/12.,1.) # 1/12 = 5/60 (min/hr)\n",
+    "time = [round (t,3) for t in time]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cd9c78e9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def data_gen2():\n",
+    "    c = -1\n",
+    "    #print (len(time))\n",
+    "    for u in data_sim_c[::12]:\n",
+    "        u1 = u[1:-1].split('|')\n",
+    "#         v1 = v[1:-1].split('],')\n",
+    "        val_u = []\n",
+    "#         val_v = []\n",
+    "        for u2 in u1:\n",
+    "            u3 = u2[1:-1].split(',')\n",
+    "            u4 = [float(u3[0]), float(u3[1]), float(u3[4])]\n",
+    "            val_u.append(u4)\n",
+    "\n",
+    "#         for v2 in v1:\n",
+    "#             v3 = v2[1:-1].split(',')\n",
+    "#             v3[0] = v3[0].replace('[', '')\n",
+    "#             depth = float(v3[0]) + layer_depth[int(v3[2])]\n",
+    "#             v4 = [depth, round(float(v3[1]),6), float(v3[2])]\n",
+    "#             #v4 = [float(v3[0]), round(float(v3[1]),6), float(v3[3])]\n",
+    "#             val_v.append(v4)\n",
+    "        c = c +1\n",
+    "        \n",
+    "        yield time[c], val_u\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b0e3bb2a",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b982e58a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %matplotlib notebook\n",
+    "fig, axs = plt.subplots(2,1,figsize=(6*1.3,4*1.3))\n",
+    "fig.tight_layout(pad=1.8)\n",
+    "len_simulation = len(Dat_u)\n",
+    "\n",
+    "# for i in range(time_step_for_plotting-1,time_step_for_plotting):\n",
+    "i = time_step_for_plotting-1\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "###plots theta vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "axs[0].set_ylim(0,layer_depth[-1])\n",
+    "axs[0].invert_yaxis()\n",
+    "axs[0].set_xlim(0,0.5)\n",
+    "\n",
+    "axs[0].set_ylabel('Depth (cm)',fontsize=12)\n",
+    "axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "        #print (k)\n",
+    "    if (k!=layer_depth[-1]):\n",
+    "        axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "    if ((count==1)|(count==2)):\n",
+    "        axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "max_val_psi = 0\n",
+    "for row in d:\n",
+    "    if (max_val_psi<max(row)):\n",
+    "        max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1\n",
+    "\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "###plots psi vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "axs[1].set_ylim(0,layer_depth[-1])\n",
+    "axs[1].invert_yaxis()\n",
+    "axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "axs[1].set_ylabel('Depth (cm)',fontsize=12)\n",
+    "axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", \"theta_e\", \"capillary head = -depth\"],loc=5)\n",
+    "\n",
+    "leg.legendHandles[1].set_color('blue')\n",
+    "leg.legendHandles[2].set_color('purple')\n",
+    "\n",
+    "plt.savefig(\"profile_xxx.pdf\")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1f5227f1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "i = time_step_for_plotting-1\n",
+    "i=200"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6e7f8504",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %matplotlib notebook\n",
+    "fig, axs = plt.subplots(4,1,figsize=(4.5*1.7,6*1.7))\n",
+    "fig.tight_layout(pad=2.5)\n",
+    "len_simulation = len(Dat_u)\n",
+    "\n",
+    "# for i in range(time_step_for_plotting-1,time_step_for_plotting):\n",
+    "# i = time_step_for_plotting-1\n",
+    "# i=3000\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "###plots theta vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "axs[0].set_ylim(0,layer_depth[-1])\n",
+    "axs[0].invert_yaxis()\n",
+    "axs[0].set_xlim(0,0.5)\n",
+    "\n",
+    "axs[0].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "        #print (k)\n",
+    "    if (k!=layer_depth[-1]):\n",
+    "        axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "    if ((count==1)|(count==2)):\n",
+    "        axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "max_val_psi = 0\n",
+    "for row in d:\n",
+    "    if (max_val_psi<max(row)):\n",
+    "        max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1\n",
+    "\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "###plots psi vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "axs[1].set_ylim(0,layer_depth[-1])\n",
+    "axs[1].invert_yaxis()\n",
+    "axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "axs[1].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", r\"$\\theta_{e}$\"],loc=3)\n",
+    "\n",
+    "leg.legendHandles[1].set_color('blue')\n",
+    "\n",
+    "\n",
+    "\n",
+    "leg2 = axs[1].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"capillary head = -depth\"],loc=3)\n",
+    "leg2.legendHandles[0].set_color('purple')\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "c_path_vars = sim_case\n",
+    "num_steps = 864 + 1\n",
+    "data_c_vars = pd.read_csv(c_path_vars)\n",
+    "data_c_vars = data_c_vars[0:num_steps]\n",
+    "\n",
+    "\n",
+    "    \n",
+    "start = pd.Timestamp(\"2020-10-01 00:00:00\")\n",
+    "end = pd.Timestamp(\"2020-10-03 00:00:00\")\n",
+    "times = pd.date_range(start=start,end=end, periods=num_steps) #freq='min'\n",
+    "times_bmi = pd.date_range(start=start,end=end, periods=num_steps)\n",
+    "times_bmi = []\n",
+    "for item in data_c_vars['Time']:\n",
+    "    times_bmi.append(datetime.strptime(item, '%Y-%m-%d %H:%M:%S'))\n",
+    "#     datetime.datetime.strptime(data_c_vars_LGARTO['Time'][i], '%Y-%m-%d %H:%M:%S') - datetime.datetime.strptime(data_c_vars_LGARTO['Time'][0],'%Y-%m-%d %H:%M:%S')\n",
+    "# times_bmi = times_bmi[0:len(times_bmi)-1]\n",
+    "\n",
+    "xfmt = mdates.DateFormatter('%Y/%m/%d')\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "m_to_mm = 1000\n",
+    "\n",
+    "\n",
+    "vars_bmi = ['soil_storage', 'surface_runoff', 'infiltration', 'actual_evapotranspiration', 'mass_balance', 'percolation', 'precipitation', 'potential_evapotranspiration']\n",
+    "\n",
+    "# vars_hydrus = [')     Volume', 'sum(RunOff', 'sum(Infil', ')   sum(vRoot', ')    sum(vBot']\n",
+    "m_to_mm = 1000\n",
+    "m_to_cm = 100\n",
+    "timestep = 300/3600\n",
+    "y_labels = ['Storage (cm)', 'Runoff (cm)', 'Infiltration (cm)', 'Actual ET (cm)', 'Mass balance (cm)', 'Recharge (cm)']\n",
+    "\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm),'b',label='C (bmi)',color='orange')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[1]],'k',linestyle='dashed',label='Hydrus',color='orange')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm),'b',label='C (bmi)',color='blue')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[2]],'k',linestyle='dashed',label='Hydrus',color='blue')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm),'b',label='C (bmi)',color='green')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[3]],'k',linestyle='dashed',label='Hydrus',color='green')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm),'b',label='C (bmi)',color='purple')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[4]],'k',linestyle='dashed',label='Hydrus',color='purple')\n",
+    "\n",
+    "axs[2].vlines(x=times_bmi[i], ymin=-100, ymax=1000, colors='red')\n",
+    "\n",
+    "#axs[2].plot(times_bmi, data_c_vars[vars_bmi[5]] * m_to_cm*12, label='C (bmi)', color='green')\n",
+    "\n",
+    "\n",
+    "leg = axs[2].legend()#(handles=[runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS],labels=[\"LGARTO\", \"HYDRUS\", \"theta_e\", \"capillary head = -depth\"])\n",
+    "\n",
+    "mylabels = ['LGARTO runoff','LGARTO infiltration','LGARTO AET','LGARTO recharge', 'current time']\n",
+    "axs[2].legend(labels=mylabels)\n",
+    "\n",
+    "axs[2].text(x=times_bmi[550], y=50, s='synth 2', fontsize = 16)\n",
+    "\n",
+    "axs[2].set_ylabel(\"cumulative \\n output flux (cm)\",fontsize=16)\n",
+    "axs[2].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "axs[2].set_ylim(min(np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm))-5, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm))))\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm),'b',label='C (bmi)',color='black')\n",
+    "axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm),'b',label='C (bmi)',color='gold')\n",
+    "\n",
+    "axs[3].set_ylim(0, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm))))\n",
+    "\n",
+    "axs[3].set_ylabel(\"cumulative \\n forcing flux (cm)\",fontsize=16)\n",
+    "\n",
+    "axs[3].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "axs[3].vlines(x=times_bmi[i], ymin=-100, ymax=1000, colors='red')\n",
+    "\n",
+    "# axs[3].text(x=times_bmi[600], y=32+5, s='current time: '+str(round(i/12,2))+' h', fontsize = 12)\n",
+    "# axs[3].text(x=times_bmi[600], y=15-1, s=str(times_bmi[i]), fontsize = 12)\n",
+    "\n",
+    "\n",
+    "mylabels_forcing = ['precipitation', 'PET', 'current time']\n",
+    "axs[3].legend(labels=mylabels_forcing)\n",
+    "\n",
+    "fig.subplots_adjust(left=0.14)\n",
+    "\n",
+    "plt.savefig(\"profile_xxx.pdf\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "625de2c5",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "587eb41f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c4c1c574",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f4acea80",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4991dd06",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8d85ad43",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a3b54e38",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d = Dat_u[0]\n",
+    "max_val_psi = 0\n",
+    "\n",
+    "for d in Dat_u:\n",
+    "    for row in d:\n",
+    "        if (max_val_psi<max(row)):\n",
+    "            max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1 + 10\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ec71d705",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "max_val_psi"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4a9acecd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def init():\n",
+    "    for ax in axs:\n",
+    "        ax.clear()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dcd556ab",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def animate(frame):\n",
+    "    \n",
+    "    for ax in axs:\n",
+    "        ax.clear()\n",
+    "    \n",
+    "    d = Dat_u[frame]\n",
+    "    num_wf = len(d)\n",
+    "    \n",
+    "    fig.subplots_adjust(left=0.14)\n",
+    "\n",
+    "    lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "    ###plots theta vs depth\n",
+    "    for j in range(num_wf):\n",
+    "        if j ==0 :\n",
+    "            axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "        else:\n",
+    "            axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "            lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "    axs[0].set_ylim(0,layer_depth[-1])\n",
+    "    axs[0].invert_yaxis()\n",
+    "    axs[0].set_xlim(0,0.5)\n",
+    "\n",
+    "    axs[0].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "    axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "    count = 0\n",
+    "    for k in layer_depth:\n",
+    "            #print (k)\n",
+    "        if (k!=layer_depth[-1]):\n",
+    "            axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "        if ((count==1)|(count==2)):\n",
+    "            axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "        if (k==layer_depth[-1]):\n",
+    "            axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "        else:\n",
+    "            axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "        count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "#     max_val_psi = 0\n",
+    "#     for row in d:\n",
+    "#         if (max_val_psi<max(row)):\n",
+    "#             max_val_psi = max(row)\n",
+    "\n",
+    "#     max_val_psi = max_val_psi*0.1\n",
+    "#     max_val_psi = 400\n",
+    "\n",
+    "    d = Dat_u[frame]\n",
+    "    num_wf = len(d)\n",
+    "\n",
+    "    ###plots psi vs depth\n",
+    "    for j in range(num_wf):\n",
+    "        if j ==0 :\n",
+    "            axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "        else:\n",
+    "            axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "            axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "    axs[1].set_ylim(0,layer_depth[-1])\n",
+    "    axs[1].invert_yaxis()\n",
+    "    axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "    axs[1].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "    axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "    count = 0\n",
+    "    for k in layer_depth:\n",
+    "\n",
+    "        if (k==layer_depth[-1]):\n",
+    "            axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "        else:\n",
+    "            axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "        count = count + 1\n",
+    "\n",
+    "    hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "    hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "    hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "    leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", r\"$\\theta_{e}$\"],loc=3,)\n",
+    "\n",
+    "    leg.legendHandles[1].set_color('blue')\n",
+    "#     leg.legendHandles[2].set_color('purple')\n",
+    "\n",
+    "    leg2 = axs[1].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"capillary head = -depth\"],loc=3)\n",
+    "    leg2.legendHandles[0].set_color('purple')\n",
+    "    \n",
+    "    \n",
+    "\n",
+    "\n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm),'b',label='C (bmi)',color='orange')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[1]],'k',linestyle='dashed',label='Hydrus',color='orange')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm),'b',label='C (bmi)',color='blue')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[2]],'k',linestyle='dashed',label='Hydrus',color='blue')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm),'b',label='C (bmi)',color='green')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[3]],'k',linestyle='dashed',label='Hydrus',color='green')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm),'b',label='C (bmi)',color='purple')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[4]],'k',linestyle='dashed',label='Hydrus',color='purple')\n",
+    "\n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    axs[2].vlines(x=times_bmi[frame], ymin=-100, ymax=1000, colors='red')\n",
+    "    \n",
+    "    leg = axs[2].legend()#(handles=[runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS],labels=[\"LGARTO\", \"HYDRUS\", \"theta_e\", \"capillary head = -depth\"])\n",
+    "\n",
+    "    mylabels = ['LGARTO runoff','LGARTO infiltration','LGARTO AET','LGARTO recharge', 'current time']\n",
+    "    axs[2].legend(labels=mylabels, ncol=3)\n",
+    "\n",
+    "    axs[2].text(x=times_bmi[550], y=50, s='synth 2', fontsize = 16)\n",
+    "\n",
+    "    axs[2].set_ylabel(\"cumulative \\n output flux (cm)\",fontsize=16)\n",
+    "    axs[2].set_xlabel(\"date\",fontsize=12)\n",
+    "    \n",
+    "#     axs[2].plot(times_bmi, data_c_vars[vars_bmi[5]] * m_to_cm*12, label='C (bmi)', color='green')\n",
+    "\n",
+    "\n",
+    "    axs[2].set_ylim(min(np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm))-5, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm))))\n",
+    "#     axs[2].set_ylim(-0.05,0.15)\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "    axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm),'b',label='C (bmi)',color='black')\n",
+    "    axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm),'b',label='C (bmi)',color='gold')\n",
+    "\n",
+    "    axs[3].set_ylim(0, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm))))\n",
+    "\n",
+    "    axs[3].vlines(x=times_bmi[frame], ymin=-100, ymax=1000, colors='red')\n",
+    "    \n",
+    "    axs[3].text(x=times_bmi[600], y=3, s='current time: '+str(round(frame/12,2))+' h', fontsize = 12)\n",
+    "#     axs[3].text(x=times_bmi[600], y=15-1, s=str(times_bmi[frame]), fontsize = 12)\n",
+    "    \n",
+    "    mylabels_forcing = ['precipitation', 'PET', 'current time']\n",
+    "    axs[3].legend(labels=mylabels_forcing, ncol=3)\n",
+    "    \n",
+    "    axs[3].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "    axs[3].set_ylabel(\"cumulative \\n forcing flux (cm)\",fontsize=16)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3af740cd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "init()\n",
+    "\n",
+    "start_move_writing_time = datetime.now()\n",
+    "\n",
+    "ani = FuncAnimation(fig, animate, frames = 10000,  init_func=init, blit=False,repeat=False)\n",
+    "# ani.save('LGARTO_animation.mp4', dpi=600, fps=50, extra_args=['-vcodec', 'libx264'])\n",
+    "ani.save('LGARTO_animation_synth_2.mp4', dpi=400, fps=48, extra_args=['-vcodec', 'libx264', '-preset', 'veryslow', '-crf', '17'])\n",
+    "\n",
+    "# ani.save('LGARTO_animation.gif', fps=50)\n",
+    "\n",
+    "end_move_writing_time = datetime.now()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5b9c4b9b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#time it took to write movie. It can be a computationally intensive thing \n",
+    "end_move_writing_time - start_move_writing_time\n",
+    "#mp4 blit=False time: 174 seconds for 200 frames\n",
+    "#then edited with 'veryslow' and crf of 17, 188 seconds for 200 frames\n",
+    "#next trying with 30 fps and 300 dpi, using the slower extra_args ... 132 seconds for 200 frames\n",
+    "#dpi=300, fps=60, slow args: 135 seconds for 200 frames, and doesn't have the freezing problem, but does kind of feel jerky. Maybe too fast, 30 fps version looked smoother. But it's a bit slow, so 45 is a good medium \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c2b2537c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f960d57e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#extra_args=['-vcodec', 'libx264']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "553247b3",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2885dd2f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e4ba966c",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "36aefd84",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7f98f67c",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f1c884db",
+   "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.8.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/tests/lgarto_profile_animation_synth_3.ipynb b/tests/lgarto_profile_animation_synth_3.ipynb
new file mode 100644
index 0000000..4e3b3d5
--- /dev/null
+++ b/tests/lgarto_profile_animation_synth_3.ipynb
@@ -0,0 +1,890 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8b26a1c1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from matplotlib import pyplot as plt\n",
+    "import matplotlib.animation as animation\n",
+    "from matplotlib.animation import FuncAnimation, PillowWriter \n",
+    "import sys, os\n",
+    "import imageio, re\n",
+    "import importlib as imp\n",
+    "import copy\n",
+    "import pandas as pd\n",
+    "\n",
+    "from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes\n",
+    "from mpl_toolkits.axes_grid1.inset_locator import mark_inset\n",
+    "import matplotlib.dates as mdates\n",
+    "from matplotlib.patches import Rectangle\n",
+    "\n",
+    "import pickle\n",
+    "\n",
+    "from datetime import datetime\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2ede7ad3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "layer_depth = [0.0, 10.0, 10.0+45.0, 10.0+45.0+45.0]\n",
+    "#51.0,99.0,53.0[cm]\n",
+    "time_step_for_plotting = 1\n",
+    "theta_e_vec = [0.45,0.41,0.38] #edit for your scenario \n",
+    "max_depth = layer_depth[-1]\n",
+    "# 36.0,111.0,107.0\n",
+    "\n",
+    "synth_case = str(3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3a9f5ba8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "c_path_layers = 'outputs/synthetic'+synth_case + '_LGARTO/data_layers.csv'\n",
+    "sim_case = 'outputs/synthetic'+synth_case+'_LGARTO/data_variables.csv' \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "72f8e45a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# hydrus_path_layers = '/Users/peterlafollette/Desktop/ngen_test_2/ngen/extern/LGARTO-C/tests/HYDRUS_outputs/Geneva/Nod_Inf.csv'\n",
+    "# data_sim_hydrus = pd.read_csv(hydrus_path_layers, sep='\\s+')\n",
+    "\n",
+    "# # data_sim_hydrus = data_sim_hydrus[0:100]\n",
+    "# # data_sim_hydrus = data_sim_hydrus[101:201]\n",
+    "# # data_sim_hydrus = data_sim_hydrus[202:302]\n",
+    "# data_sim_hydrus = data_sim_hydrus[202:]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "efb644ba",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5ef98a57",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data_sim_c = np.loadtxt(c_path_layers,dtype='str')\n",
+    "# data_sim_py = data_sim_c #setting these equal for now to get code working -- attempting to remove data_sim_py from code \n",
+    "\n",
+    "Depth_d = []\n",
+    "Theta_d = []\n",
+    "Head_d = []\n",
+    "Dat_u = []\n",
+    "Dat_v = []\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "29115339",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f820a6b0",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f9ff759c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "count = 0\n",
+    "for u in data_sim_c:\n",
+    "    u1 = u[1:-1].split('|')\n",
+    "    val_u = []\n",
+    "\n",
+    "    for u2 in u1:\n",
+    "        u3 = u2[1:-1].split(',')\n",
+    "        u4 = [float(u3[0]), float(u3[1]), float(u3[4])]\n",
+    "        val_u.append(u4)\n",
+    "        \n",
+    "    count = count + 1\n",
+    "    Dat_u.append(val_u)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "276f3f5c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def data_gen1():\n",
+    "    for u in data_sim_c[::12]:\n",
+    "        v1 = u[1:-1].split('|')\n",
+    "        val_c = []\n",
+    "        \n",
+    "        for v2 in v1:\n",
+    "            v3 = v2[1:-1].split(',')\n",
+    "            v4 = [float(v3[0]), float(v3[1]), float(v3[4])]\n",
+    "            val_c.append(v4)\n",
+    "        \n",
+    "        yield val_c\n",
+    "        \n",
+    "time = np.arange(0,len(Dat_u)/12.,1.) # 1/12 = 5/60 (min/hr)\n",
+    "time = [round (t,3) for t in time]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cd9c78e9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def data_gen2():\n",
+    "    c = -1\n",
+    "    #print (len(time))\n",
+    "    for u in data_sim_c[::12]:\n",
+    "        u1 = u[1:-1].split('|')\n",
+    "#         v1 = v[1:-1].split('],')\n",
+    "        val_u = []\n",
+    "#         val_v = []\n",
+    "        for u2 in u1:\n",
+    "            u3 = u2[1:-1].split(',')\n",
+    "            u4 = [float(u3[0]), float(u3[1]), float(u3[4])]\n",
+    "            val_u.append(u4)\n",
+    "\n",
+    "#         for v2 in v1:\n",
+    "#             v3 = v2[1:-1].split(',')\n",
+    "#             v3[0] = v3[0].replace('[', '')\n",
+    "#             depth = float(v3[0]) + layer_depth[int(v3[2])]\n",
+    "#             v4 = [depth, round(float(v3[1]),6), float(v3[2])]\n",
+    "#             #v4 = [float(v3[0]), round(float(v3[1]),6), float(v3[3])]\n",
+    "#             val_v.append(v4)\n",
+    "        c = c +1\n",
+    "        \n",
+    "        yield time[c], val_u\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b0e3bb2a",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b982e58a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %matplotlib notebook\n",
+    "fig, axs = plt.subplots(2,1,figsize=(6*1.3,4*1.3))\n",
+    "fig.tight_layout(pad=1.8)\n",
+    "len_simulation = len(Dat_u)\n",
+    "\n",
+    "# for i in range(time_step_for_plotting-1,time_step_for_plotting):\n",
+    "i = time_step_for_plotting-1\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "###plots theta vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "axs[0].set_ylim(0,layer_depth[-1])\n",
+    "axs[0].invert_yaxis()\n",
+    "axs[0].set_xlim(0,0.5)\n",
+    "\n",
+    "axs[0].set_ylabel('Depth (cm)',fontsize=12)\n",
+    "axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "        #print (k)\n",
+    "    if (k!=layer_depth[-1]):\n",
+    "        axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "    if ((count==1)|(count==2)):\n",
+    "        axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "max_val_psi = 0\n",
+    "for row in d:\n",
+    "    if (max_val_psi<max(row)):\n",
+    "        max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1\n",
+    "\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "###plots psi vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "axs[1].set_ylim(0,layer_depth[-1])\n",
+    "axs[1].invert_yaxis()\n",
+    "axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "axs[1].set_ylabel('Depth (cm)',fontsize=12)\n",
+    "axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", \"theta_e\", \"capillary head = -depth\"],loc=5)\n",
+    "\n",
+    "leg.legendHandles[1].set_color('blue')\n",
+    "leg.legendHandles[2].set_color('purple')\n",
+    "\n",
+    "plt.savefig(\"profile_xxx.pdf\")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1f5227f1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "i = time_step_for_plotting-1\n",
+    "i=685+4"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6e7f8504",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %matplotlib notebook\n",
+    "fig, axs = plt.subplots(4,1,figsize=(4.5*1.7,6*1.7))\n",
+    "fig.tight_layout(pad=2.5)\n",
+    "len_simulation = len(Dat_u)\n",
+    "\n",
+    "# for i in range(time_step_for_plotting-1,time_step_for_plotting):\n",
+    "# i = time_step_for_plotting-1\n",
+    "# i=3000\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "###plots theta vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "axs[0].set_ylim(0,layer_depth[-1])\n",
+    "axs[0].invert_yaxis()\n",
+    "axs[0].set_xlim(0.3,0.5)\n",
+    "\n",
+    "axs[0].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "        #print (k)\n",
+    "    if (k!=layer_depth[-1]):\n",
+    "        axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "    if ((count==1)|(count==2)):\n",
+    "        axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "max_val_psi = 0\n",
+    "for row in d:\n",
+    "    if (max_val_psi<max(row)):\n",
+    "        max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1\n",
+    "\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "###plots psi vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "axs[1].set_ylim(0,layer_depth[-1])\n",
+    "axs[1].invert_yaxis()\n",
+    "axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "axs[1].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", r\"$\\theta_{e}$\"],loc=3)\n",
+    "\n",
+    "leg.legendHandles[1].set_color('blue')\n",
+    "\n",
+    "\n",
+    "\n",
+    "leg2 = axs[1].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"capillary head = -depth\"],loc=3)\n",
+    "leg2.legendHandles[0].set_color('purple')\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "c_path_vars = sim_case\n",
+    "num_steps = 864 + 1\n",
+    "data_c_vars = pd.read_csv(c_path_vars)\n",
+    "data_c_vars = data_c_vars[0:num_steps]\n",
+    "\n",
+    "\n",
+    "\n",
+    "    \n",
+    "start = pd.Timestamp(\"2020-10-01 00:00:00\")\n",
+    "end = pd.Timestamp(\"2020-10-03 00:00:00\")\n",
+    "times = pd.date_range(start=start,end=end, periods=num_steps) #freq='min'\n",
+    "times_bmi = pd.date_range(start=start,end=end, periods=num_steps)\n",
+    "times_bmi = []\n",
+    "for item in data_c_vars['Time']:\n",
+    "    times_bmi.append(datetime.strptime(item, '%Y-%m-%d %H:%M:%S'))\n",
+    "#     datetime.datetime.strptime(data_c_vars_LGARTO['Time'][i], '%Y-%m-%d %H:%M:%S') - datetime.datetime.strptime(data_c_vars_LGARTO['Time'][0],'%Y-%m-%d %H:%M:%S')\n",
+    "# times_bmi = times_bmi[0:len(times_bmi)-1]\n",
+    "\n",
+    "xfmt = mdates.DateFormatter('%Y/%m/%d')\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "m_to_mm = 1000\n",
+    "\n",
+    "\n",
+    "vars_bmi = ['soil_storage', 'surface_runoff', 'infiltration', 'actual_evapotranspiration', 'mass_balance', 'percolation', 'precipitation', 'potential_evapotranspiration']\n",
+    "\n",
+    "# vars_hydrus = [')     Volume', 'sum(RunOff', 'sum(Infil', ')   sum(vRoot', ')    sum(vBot']\n",
+    "m_to_mm = 1000\n",
+    "m_to_cm = 100\n",
+    "timestep = 300/3600\n",
+    "y_labels = ['Storage (cm)', 'Runoff (cm)', 'Infiltration (cm)', 'Actual ET (cm)', 'Mass balance (cm)', 'Recharge (cm)']\n",
+    "\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm),'b',label='C (bmi)',color='orange')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[1]],'k',linestyle='dashed',label='Hydrus',color='orange')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm),'b',label='C (bmi)',color='blue')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[2]],'k',linestyle='dashed',label='Hydrus',color='blue')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm),'b',label='C (bmi)',color='green')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[3]],'k',linestyle='dashed',label='Hydrus',color='green')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm),'b',label='C (bmi)',color='purple')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[4]],'k',linestyle='dashed',label='Hydrus',color='purple')\n",
+    "\n",
+    "axs[2].vlines(x=times_bmi[i], ymin=-100, ymax=1000, colors='red')\n",
+    "\n",
+    "#axs[2].plot(times_bmi, data_c_vars[vars_bmi[5]] * m_to_cm*12, label='C (bmi)', color='green')\n",
+    "\n",
+    "\n",
+    "leg = axs[2].legend()#(handles=[runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS],labels=[\"LGARTO\", \"HYDRUS\", \"theta_e\", \"capillary head = -depth\"])\n",
+    "\n",
+    "mylabels = ['LGARTO runoff','LGARTO infiltration','LGARTO AET','LGARTO recharge', 'current time']\n",
+    "# axs[2].legend(labels=mylabels)\n",
+    "axs[2].legend(labels=mylabels, ncol=3, bbox_to_anchor=(0.5, -0.0), loc='lower center')\n",
+    "\n",
+    "\n",
+    "axs[2].text(x=times_bmi[550], y=50, s='synth '+synth_case, fontsize = 16)\n",
+    "\n",
+    "axs[2].set_ylabel(\"cumulative \\n output flux (cm)\",fontsize=16)\n",
+    "axs[2].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "axs[2].set_ylim(min(np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm))-3, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm))))\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm),'b',label='C (bmi)',color='black')\n",
+    "axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm),'b',label='C (bmi)',color='gold')\n",
+    "\n",
+    "axs[3].set_ylim(0, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm))))\n",
+    "\n",
+    "axs[3].set_ylabel(\"cumulative \\n forcing flux (cm)\",fontsize=16)\n",
+    "\n",
+    "axs[3].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "axs[3].vlines(x=times_bmi[i], ymin=-100, ymax=1000, colors='red')\n",
+    "\n",
+    "# axs[3].text(x=times_bmi[600], y=32+5, s='current time: '+str(round(i/12,2))+' h', fontsize = 12)\n",
+    "# axs[3].text(x=times_bmi[600], y=15-1, s=str(times_bmi[i]), fontsize = 12)\n",
+    "\n",
+    "\n",
+    "mylabels_forcing = ['precipitation', 'PET', 'current time']\n",
+    "axs[3].legend(labels=mylabels_forcing)\n",
+    "\n",
+    "fig.subplots_adjust(left=0.14)\n",
+    "\n",
+    "plt.savefig(\"profile_xxx.pdf\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "625de2c5",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "587eb41f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c4c1c574",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f4acea80",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4991dd06",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8d85ad43",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a3b54e38",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d = Dat_u[0]\n",
+    "max_val_psi = 0\n",
+    "\n",
+    "for d in Dat_u:\n",
+    "    for row in d:\n",
+    "        if (max_val_psi<max(row)):\n",
+    "            max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1 + 10\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ec71d705",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "max_val_psi"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4a9acecd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def init():\n",
+    "    for ax in axs:\n",
+    "        ax.clear()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dcd556ab",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def animate(frame):\n",
+    "    \n",
+    "    for ax in axs:\n",
+    "        ax.clear()\n",
+    "    \n",
+    "    d = Dat_u[frame]\n",
+    "    num_wf = len(d)\n",
+    "    \n",
+    "    fig.subplots_adjust(left=0.14)\n",
+    "\n",
+    "    lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "    ###plots theta vs depth\n",
+    "    for j in range(num_wf):\n",
+    "        if j ==0 :\n",
+    "            axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "        else:\n",
+    "            axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "            lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "    axs[0].set_ylim(0,layer_depth[-1])\n",
+    "    axs[0].invert_yaxis()\n",
+    "    axs[0].set_xlim(0,0.5)\n",
+    "\n",
+    "    axs[0].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "    axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "    count = 0\n",
+    "    for k in layer_depth:\n",
+    "            #print (k)\n",
+    "        if (k!=layer_depth[-1]):\n",
+    "            axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "        if ((count==1)|(count==2)):\n",
+    "            axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "        if (k==layer_depth[-1]):\n",
+    "            axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "        else:\n",
+    "            axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "        count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "#     max_val_psi = 0\n",
+    "#     for row in d:\n",
+    "#         if (max_val_psi<max(row)):\n",
+    "#             max_val_psi = max(row)\n",
+    "\n",
+    "#     max_val_psi = max_val_psi*0.1\n",
+    "#     max_val_psi = 400\n",
+    "\n",
+    "    d = Dat_u[frame]\n",
+    "    num_wf = len(d)\n",
+    "\n",
+    "    ###plots psi vs depth\n",
+    "    for j in range(num_wf):\n",
+    "        if j ==0 :\n",
+    "            axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "        else:\n",
+    "            axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "            axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "    axs[1].set_ylim(0,layer_depth[-1])\n",
+    "    axs[1].invert_yaxis()\n",
+    "    axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "    axs[1].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "    axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "    count = 0\n",
+    "    for k in layer_depth:\n",
+    "\n",
+    "        if (k==layer_depth[-1]):\n",
+    "            axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "        else:\n",
+    "            axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "        count = count + 1\n",
+    "\n",
+    "    hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "    hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "    hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "    leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", r\"$\\theta_{e}$\"],loc=3,)\n",
+    "\n",
+    "    leg.legendHandles[1].set_color('blue')\n",
+    "#     leg.legendHandles[2].set_color('purple')\n",
+    "\n",
+    "    leg2 = axs[1].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"capillary head = -depth\"],loc=3)\n",
+    "    leg2.legendHandles[0].set_color('purple')\n",
+    "    \n",
+    "    \n",
+    "\n",
+    "\n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm),'b',label='C (bmi)',color='orange')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[1]],'k',linestyle='dashed',label='Hydrus',color='orange')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm),'b',label='C (bmi)',color='blue')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[2]],'k',linestyle='dashed',label='Hydrus',color='blue')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm),'b',label='C (bmi)',color='green')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[3]],'k',linestyle='dashed',label='Hydrus',color='green')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm),'b',label='C (bmi)',color='purple')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[4]],'k',linestyle='dashed',label='Hydrus',color='purple')\n",
+    "\n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    axs[2].vlines(x=times_bmi[frame], ymin=-100, ymax=1000, colors='red')\n",
+    "    \n",
+    "    leg = axs[2].legend()#(handles=[runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS],labels=[\"LGARTO\", \"HYDRUS\", \"theta_e\", \"capillary head = -depth\"])\n",
+    "\n",
+    "    mylabels = ['LGARTO runoff','LGARTO infiltration','LGARTO AET','LGARTO recharge', 'current time']\n",
+    "    axs[2].legend(labels=mylabels, ncol=3, bbox_to_anchor=(0.5, -0.0), loc='lower center')\n",
+    "\n",
+    "    axs[2].text(x=times_bmi[550], y=50, s='synth '+synth_case, fontsize = 16)\n",
+    "\n",
+    "    axs[2].set_ylabel(\"cumulative \\n output flux (cm)\",fontsize=16)\n",
+    "    axs[2].set_xlabel(\"date\",fontsize=12)\n",
+    "    \n",
+    "#     axs[2].plot(times_bmi, data_c_vars[vars_bmi[5]] * m_to_cm*12, label='C (bmi)', color='green')\n",
+    "\n",
+    "\n",
+    "    axs[2].set_ylim(min(np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm))-3, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm))))\n",
+    "#     axs[2].set_ylim(-0.05,0.15)\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "    axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm),'b',label='C (bmi)',color='black')\n",
+    "    axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm),'b',label='C (bmi)',color='gold')\n",
+    "\n",
+    "    axs[3].set_ylim(0, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm))))\n",
+    "\n",
+    "    axs[3].vlines(x=times_bmi[frame], ymin=-100, ymax=1000, colors='red')\n",
+    "    \n",
+    "    axs[3].text(x=times_bmi[600], y=3, s='current time: '+str(round(frame/12,2))+' h', fontsize = 12)\n",
+    "#     axs[3].text(x=times_bmi[600], y=15-1, s=str(times_bmi[frame]), fontsize = 12)\n",
+    "    \n",
+    "    mylabels_forcing = ['precipitation', 'PET', 'current time']\n",
+    "    axs[3].legend(labels=mylabels_forcing, ncol=3)\n",
+    "    \n",
+    "    axs[3].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "    axs[3].set_ylabel(\"cumulative \\n forcing flux (cm)\",fontsize=16)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3af740cd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "init()\n",
+    "\n",
+    "start_move_writing_time = datetime.now()\n",
+    "\n",
+    "ani = FuncAnimation(fig, animate, frames = 10000,  init_func=init, blit=False,repeat=False)\n",
+    "# ani.save('LGARTO_animation.mp4', dpi=600, fps=50, extra_args=['-vcodec', 'libx264'])\n",
+    "ani.save('LGARTO_animation_synth_'+synth_case+'.mp4', dpi=400, fps=48, extra_args=['-vcodec', 'libx264', '-preset', 'veryslow', '-crf', '17'])\n",
+    "\n",
+    "# ani.save('LGARTO_animation.gif', fps=50)\n",
+    "\n",
+    "end_move_writing_time = datetime.now()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5b9c4b9b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#time it took to write movie. It can be a computationally intensive thing \n",
+    "end_move_writing_time - start_move_writing_time\n",
+    "#mp4 blit=False time: 174 seconds for 200 frames\n",
+    "#then edited with 'veryslow' and crf of 17, 188 seconds for 200 frames\n",
+    "#next trying with 30 fps and 300 dpi, using the slower extra_args ... 132 seconds for 200 frames\n",
+    "#dpi=300, fps=60, slow args: 135 seconds for 200 frames, and doesn't have the freezing problem, but does kind of feel jerky. Maybe too fast, 30 fps version looked smoother. But it's a bit slow, so 45 is a good medium \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c2b2537c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f960d57e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#extra_args=['-vcodec', 'libx264']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "553247b3",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2885dd2f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e4ba966c",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "36aefd84",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7f98f67c",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f1c884db",
+   "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.8.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/tests/lgarto_profile_animation_synth_4.ipynb b/tests/lgarto_profile_animation_synth_4.ipynb
new file mode 100644
index 0000000..f74598e
--- /dev/null
+++ b/tests/lgarto_profile_animation_synth_4.ipynb
@@ -0,0 +1,890 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8b26a1c1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from matplotlib import pyplot as plt\n",
+    "import matplotlib.animation as animation\n",
+    "from matplotlib.animation import FuncAnimation, PillowWriter \n",
+    "import sys, os\n",
+    "import imageio, re\n",
+    "import importlib as imp\n",
+    "import copy\n",
+    "import pandas as pd\n",
+    "\n",
+    "from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes\n",
+    "from mpl_toolkits.axes_grid1.inset_locator import mark_inset\n",
+    "import matplotlib.dates as mdates\n",
+    "from matplotlib.patches import Rectangle\n",
+    "\n",
+    "import pickle\n",
+    "\n",
+    "from datetime import datetime\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2ede7ad3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "layer_depth = [0.0, 10.0, 10.0+45.0, 10.0+45.0+45.0]\n",
+    "#51.0,99.0,53.0[cm]\n",
+    "time_step_for_plotting = 1\n",
+    "theta_e_vec = [0.43,0.46,0.45] #edit for your scenario \n",
+    "max_depth = layer_depth[-1]\n",
+    "# 36.0,111.0,107.0\n",
+    "\n",
+    "synth_case = str(4)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3a9f5ba8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "c_path_layers = 'outputs/synthetic'+synth_case + '_LGARTO/data_layers.csv'\n",
+    "sim_case = 'outputs/synthetic'+synth_case+'_LGARTO/data_variables.csv' \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "72f8e45a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# hydrus_path_layers = '/Users/peterlafollette/Desktop/ngen_test_2/ngen/extern/LGARTO-C/tests/HYDRUS_outputs/Geneva/Nod_Inf.csv'\n",
+    "# data_sim_hydrus = pd.read_csv(hydrus_path_layers, sep='\\s+')\n",
+    "\n",
+    "# # data_sim_hydrus = data_sim_hydrus[0:100]\n",
+    "# # data_sim_hydrus = data_sim_hydrus[101:201]\n",
+    "# # data_sim_hydrus = data_sim_hydrus[202:302]\n",
+    "# data_sim_hydrus = data_sim_hydrus[202:]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "efb644ba",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5ef98a57",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data_sim_c = np.loadtxt(c_path_layers,dtype='str')\n",
+    "# data_sim_py = data_sim_c #setting these equal for now to get code working -- attempting to remove data_sim_py from code \n",
+    "\n",
+    "Depth_d = []\n",
+    "Theta_d = []\n",
+    "Head_d = []\n",
+    "Dat_u = []\n",
+    "Dat_v = []\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "29115339",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f820a6b0",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f9ff759c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "count = 0\n",
+    "for u in data_sim_c:\n",
+    "    u1 = u[1:-1].split('|')\n",
+    "    val_u = []\n",
+    "\n",
+    "    for u2 in u1:\n",
+    "        u3 = u2[1:-1].split(',')\n",
+    "        u4 = [float(u3[0]), float(u3[1]), float(u3[4])]\n",
+    "        val_u.append(u4)\n",
+    "        \n",
+    "    count = count + 1\n",
+    "    Dat_u.append(val_u)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "276f3f5c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def data_gen1():\n",
+    "    for u in data_sim_c[::12]:\n",
+    "        v1 = u[1:-1].split('|')\n",
+    "        val_c = []\n",
+    "        \n",
+    "        for v2 in v1:\n",
+    "            v3 = v2[1:-1].split(',')\n",
+    "            v4 = [float(v3[0]), float(v3[1]), float(v3[4])]\n",
+    "            val_c.append(v4)\n",
+    "        \n",
+    "        yield val_c\n",
+    "        \n",
+    "time = np.arange(0,len(Dat_u)/12.,1.) # 1/12 = 5/60 (min/hr)\n",
+    "time = [round (t,3) for t in time]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cd9c78e9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def data_gen2():\n",
+    "    c = -1\n",
+    "    #print (len(time))\n",
+    "    for u in data_sim_c[::12]:\n",
+    "        u1 = u[1:-1].split('|')\n",
+    "#         v1 = v[1:-1].split('],')\n",
+    "        val_u = []\n",
+    "#         val_v = []\n",
+    "        for u2 in u1:\n",
+    "            u3 = u2[1:-1].split(',')\n",
+    "            u4 = [float(u3[0]), float(u3[1]), float(u3[4])]\n",
+    "            val_u.append(u4)\n",
+    "\n",
+    "#         for v2 in v1:\n",
+    "#             v3 = v2[1:-1].split(',')\n",
+    "#             v3[0] = v3[0].replace('[', '')\n",
+    "#             depth = float(v3[0]) + layer_depth[int(v3[2])]\n",
+    "#             v4 = [depth, round(float(v3[1]),6), float(v3[2])]\n",
+    "#             #v4 = [float(v3[0]), round(float(v3[1]),6), float(v3[3])]\n",
+    "#             val_v.append(v4)\n",
+    "        c = c +1\n",
+    "        \n",
+    "        yield time[c], val_u\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b0e3bb2a",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b982e58a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %matplotlib notebook\n",
+    "fig, axs = plt.subplots(2,1,figsize=(6*1.3,4*1.3))\n",
+    "fig.tight_layout(pad=1.8)\n",
+    "len_simulation = len(Dat_u)\n",
+    "\n",
+    "# for i in range(time_step_for_plotting-1,time_step_for_plotting):\n",
+    "i = time_step_for_plotting-1\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "###plots theta vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "axs[0].set_ylim(0,layer_depth[-1])\n",
+    "axs[0].invert_yaxis()\n",
+    "axs[0].set_xlim(0,0.5)\n",
+    "\n",
+    "axs[0].set_ylabel('Depth (cm)',fontsize=12)\n",
+    "axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "        #print (k)\n",
+    "    if (k!=layer_depth[-1]):\n",
+    "        axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "    if ((count==1)|(count==2)):\n",
+    "        axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "max_val_psi = 0\n",
+    "for row in d:\n",
+    "    if (max_val_psi<max(row)):\n",
+    "        max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1\n",
+    "\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "###plots psi vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "axs[1].set_ylim(0,layer_depth[-1])\n",
+    "axs[1].invert_yaxis()\n",
+    "axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "axs[1].set_ylabel('Depth (cm)',fontsize=12)\n",
+    "axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", \"theta_e\", \"capillary head = -depth\"],loc=5)\n",
+    "\n",
+    "leg.legendHandles[1].set_color('blue')\n",
+    "leg.legendHandles[2].set_color('purple')\n",
+    "\n",
+    "plt.savefig(\"profile_xxx.pdf\")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1f5227f1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "i = time_step_for_plotting-1\n",
+    "i=200"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6e7f8504",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %matplotlib notebook\n",
+    "fig, axs = plt.subplots(4,1,figsize=(4.5*1.7,6*1.7))\n",
+    "fig.tight_layout(pad=2.5)\n",
+    "len_simulation = len(Dat_u)\n",
+    "\n",
+    "# for i in range(time_step_for_plotting-1,time_step_for_plotting):\n",
+    "# i = time_step_for_plotting-1\n",
+    "# i=3000\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "###plots theta vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "axs[0].set_ylim(0,layer_depth[-1])\n",
+    "axs[0].invert_yaxis()\n",
+    "axs[0].set_xlim(0,0.5)\n",
+    "\n",
+    "axs[0].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "        #print (k)\n",
+    "    if (k!=layer_depth[-1]):\n",
+    "        axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "    if ((count==1)|(count==2)):\n",
+    "        axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "max_val_psi = 0\n",
+    "for row in d:\n",
+    "    if (max_val_psi<max(row)):\n",
+    "        max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1\n",
+    "\n",
+    "d = Dat_u[i]\n",
+    "num_wf = len(d)\n",
+    "\n",
+    "###plots psi vs depth\n",
+    "for j in range(num_wf):\n",
+    "    if j ==0 :\n",
+    "        axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "    else:\n",
+    "        axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "        axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "axs[1].set_ylim(0,layer_depth[-1])\n",
+    "axs[1].invert_yaxis()\n",
+    "axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "axs[1].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "count = 0\n",
+    "for k in layer_depth:\n",
+    "\n",
+    "    if (k==layer_depth[-1]):\n",
+    "        axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "    else:\n",
+    "        axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "    count = count + 1\n",
+    "\n",
+    "hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", r\"$\\theta_{e}$\"],loc=3)\n",
+    "\n",
+    "leg.legendHandles[1].set_color('blue')\n",
+    "\n",
+    "\n",
+    "\n",
+    "leg2 = axs[1].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"capillary head = -depth\"],loc=3)\n",
+    "leg2.legendHandles[0].set_color('purple')\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "c_path_vars = sim_case\n",
+    "num_steps = 864 + 1\n",
+    "data_c_vars = pd.read_csv(c_path_vars)\n",
+    "data_c_vars = data_c_vars[0:num_steps]\n",
+    "\n",
+    "\n",
+    "\n",
+    "    \n",
+    "start = pd.Timestamp(\"2020-10-01 00:00:00\")\n",
+    "end = pd.Timestamp(\"2020-10-03 00:00:00\")\n",
+    "times = pd.date_range(start=start,end=end, periods=num_steps) #freq='min'\n",
+    "times_bmi = pd.date_range(start=start,end=end, periods=num_steps)\n",
+    "times_bmi = []\n",
+    "for item in data_c_vars['Time']:\n",
+    "    times_bmi.append(datetime.strptime(item, '%Y-%m-%d %H:%M:%S'))\n",
+    "#     datetime.datetime.strptime(data_c_vars_LGARTO['Time'][i], '%Y-%m-%d %H:%M:%S') - datetime.datetime.strptime(data_c_vars_LGARTO['Time'][0],'%Y-%m-%d %H:%M:%S')\n",
+    "# times_bmi = times_bmi[0:len(times_bmi)-1]\n",
+    "\n",
+    "xfmt = mdates.DateFormatter('%Y/%m/%d')\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "m_to_mm = 1000\n",
+    "\n",
+    "\n",
+    "vars_bmi = ['soil_storage', 'surface_runoff', 'infiltration', 'actual_evapotranspiration', 'mass_balance', 'percolation', 'precipitation', 'potential_evapotranspiration']\n",
+    "\n",
+    "# vars_hydrus = [')     Volume', 'sum(RunOff', 'sum(Infil', ')   sum(vRoot', ')    sum(vBot']\n",
+    "m_to_mm = 1000\n",
+    "m_to_cm = 100\n",
+    "timestep = 300/3600\n",
+    "y_labels = ['Storage (cm)', 'Runoff (cm)', 'Infiltration (cm)', 'Actual ET (cm)', 'Mass balance (cm)', 'Recharge (cm)']\n",
+    "\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm),'b',label='C (bmi)',color='orange')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[1]],'k',linestyle='dashed',label='Hydrus',color='orange')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm),'b',label='C (bmi)',color='blue')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[2]],'k',linestyle='dashed',label='Hydrus',color='blue')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm),'b',label='C (bmi)',color='green')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[3]],'k',linestyle='dashed',label='Hydrus',color='green')\n",
+    "\n",
+    "axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm),'b',label='C (bmi)',color='purple')\n",
+    "# plt.plot(HYDRUS_output[vars_hydrus[4]],'k',linestyle='dashed',label='Hydrus',color='purple')\n",
+    "\n",
+    "axs[2].vlines(x=times_bmi[i], ymin=-100, ymax=1000, colors='red')\n",
+    "\n",
+    "#axs[2].plot(times_bmi, data_c_vars[vars_bmi[5]] * m_to_cm*12, label='C (bmi)', color='green')\n",
+    "\n",
+    "\n",
+    "leg = axs[2].legend()#(handles=[runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS],labels=[\"LGARTO\", \"HYDRUS\", \"theta_e\", \"capillary head = -depth\"])\n",
+    "\n",
+    "mylabels = ['LGARTO runoff','LGARTO infiltration','LGARTO AET','LGARTO recharge', 'current time']\n",
+    "# axs[2].legend(labels=mylabels)\n",
+    "axs[2].legend(labels=mylabels, ncol=3, bbox_to_anchor=(0.5, -0.0), loc='lower center')\n",
+    "\n",
+    "\n",
+    "axs[2].text(x=times_bmi[550], y=50, s='synth '+synth_case, fontsize = 16)\n",
+    "\n",
+    "axs[2].set_ylabel(\"cumulative \\n output flux (cm)\",fontsize=16)\n",
+    "axs[2].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "axs[2].set_ylim(min(np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm))-3, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm))))\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm),'b',label='C (bmi)',color='black')\n",
+    "axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm),'b',label='C (bmi)',color='gold')\n",
+    "\n",
+    "axs[3].set_ylim(0, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm))))\n",
+    "\n",
+    "axs[3].set_ylabel(\"cumulative \\n forcing flux (cm)\",fontsize=16)\n",
+    "\n",
+    "axs[3].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "axs[3].vlines(x=times_bmi[i], ymin=-100, ymax=1000, colors='red')\n",
+    "\n",
+    "# axs[3].text(x=times_bmi[600], y=32+5, s='current time: '+str(round(i/12,2))+' h', fontsize = 12)\n",
+    "# axs[3].text(x=times_bmi[600], y=15-1, s=str(times_bmi[i]), fontsize = 12)\n",
+    "\n",
+    "\n",
+    "mylabels_forcing = ['precipitation', 'PET', 'current time']\n",
+    "axs[3].legend(labels=mylabels_forcing)\n",
+    "\n",
+    "fig.subplots_adjust(left=0.14)\n",
+    "\n",
+    "plt.savefig(\"profile_xxx.pdf\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "625de2c5",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "587eb41f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c4c1c574",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f4acea80",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4991dd06",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8d85ad43",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a3b54e38",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d = Dat_u[0]\n",
+    "max_val_psi = 0\n",
+    "\n",
+    "for d in Dat_u:\n",
+    "    for row in d:\n",
+    "        if (max_val_psi<max(row)):\n",
+    "            max_val_psi = max(row)\n",
+    "\n",
+    "max_val_psi = max_val_psi*0.1 + 10\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ec71d705",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "max_val_psi"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4a9acecd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def init():\n",
+    "    for ax in axs:\n",
+    "        ax.clear()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dcd556ab",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def animate(frame):\n",
+    "    \n",
+    "    for ax in axs:\n",
+    "        ax.clear()\n",
+    "    \n",
+    "    d = Dat_u[frame]\n",
+    "    num_wf = len(d)\n",
+    "    \n",
+    "    fig.subplots_adjust(left=0.14)\n",
+    "\n",
+    "    lines_LGARTO = axs[0].hlines(0,0,0,color='k',lw=1)\n",
+    "    ###plots theta vs depth\n",
+    "    for j in range(num_wf):\n",
+    "        if j ==0 :\n",
+    "            axs[0].vlines(d[j][1],0,d[j][0]*0.1,color='k',lw=1)\n",
+    "        else:\n",
+    "            axs[0].vlines(d[j][1],d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "            lines_LGARTO = axs[0].hlines(d[j-1][0]*0.1,d[j-1][1],d[j][1],color='k',lw=1)\n",
+    "\n",
+    "    axs[0].set_ylim(0,layer_depth[-1])\n",
+    "    axs[0].invert_yaxis()\n",
+    "    axs[0].set_xlim(0,0.5)\n",
+    "\n",
+    "    axs[0].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "    axs[0].set_xlabel('Volumetric water content (-)',fontsize=12)\n",
+    "\n",
+    "    count = 0\n",
+    "    for k in layer_depth:\n",
+    "            #print (k)\n",
+    "        if (k!=layer_depth[-1]):\n",
+    "            axs[0].vlines(theta_e_vec[count], ymin=k, ymax=layer_depth[count+1],color='blue')\n",
+    "\n",
+    "        if ((count==1)|(count==2)):\n",
+    "            axs[0].hlines(k, xmin=theta_e_vec[count], xmax=theta_e_vec[count-1],color='blue')\n",
+    "\n",
+    "        if (k==layer_depth[-1]):\n",
+    "            axs[0].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "        else:\n",
+    "            axs[0].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "        count = count + 1\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "#     max_val_psi = 0\n",
+    "#     for row in d:\n",
+    "#         if (max_val_psi<max(row)):\n",
+    "#             max_val_psi = max(row)\n",
+    "\n",
+    "#     max_val_psi = max_val_psi*0.1\n",
+    "#     max_val_psi = 400\n",
+    "\n",
+    "    d = Dat_u[frame]\n",
+    "    num_wf = len(d)\n",
+    "\n",
+    "    ###plots psi vs depth\n",
+    "    for j in range(num_wf):\n",
+    "        if j ==0 :\n",
+    "            axs[1].vlines(-d[j][2]*0.1,0,d[j][0]*0.1,color='k',lw=1)\n",
+    "        else:\n",
+    "            axs[1].vlines(-d[j][2]*0.1,d[j-1][0]*0.1,d[j][0]*0.1,color='k',lw=1)\n",
+    "            axs[1].hlines(d[j-1][0]*0.1,-d[j-1][2]*0.1,-d[j][2]*0.1,color='k',lw=1)\n",
+    "\n",
+    "    axs[1].set_ylim(0,layer_depth[-1])\n",
+    "    axs[1].invert_yaxis()\n",
+    "    axs[1].set_xlim(-max_val_psi*1.05,0)\n",
+    "\n",
+    "    axs[1].set_ylabel('Depth (cm)',fontsize=16)\n",
+    "    axs[1].set_xlabel('Capillary head (cm)',fontsize=12)\n",
+    "\n",
+    "    count = 0\n",
+    "    for k in layer_depth:\n",
+    "\n",
+    "        if (k==layer_depth[-1]):\n",
+    "            axs[1].axhline(k, ls='solid', color='black', lw=1, zorder=10)\n",
+    "        else:\n",
+    "            axs[1].axhline(k, ls='dashed', color='lightgrey', lw=1, zorder=10)\n",
+    "\n",
+    "        count = count + 1\n",
+    "\n",
+    "    hydrostatic_line_data_x = [-1*layer_depth[-1], 0]\n",
+    "    hydrostatic_line_data_y = [0, layer_depth[-1]]\n",
+    "    hydrostatic_line = axs[1].plot(np.array(hydrostatic_line_data_x),np.array(hydrostatic_line_data_y),color='purple')\n",
+    "\n",
+    "    leg = axs[0].legend(handles=[lines_LGARTO,lines_LGARTO],labels=[\"LGARTO\", r\"$\\theta_{e}$\"],loc=3,)\n",
+    "\n",
+    "    leg.legendHandles[1].set_color('blue')\n",
+    "#     leg.legendHandles[2].set_color('purple')\n",
+    "\n",
+    "    leg2 = axs[1].legend(handles=[lines_LGARTO,lines_LGARTO,lines_LGARTO],labels=[\"capillary head = -depth\"],loc=3)\n",
+    "    leg2.legendHandles[0].set_color('purple')\n",
+    "    \n",
+    "    \n",
+    "\n",
+    "\n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm),'b',label='C (bmi)',color='orange')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[1]],'k',linestyle='dashed',label='Hydrus',color='orange')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm),'b',label='C (bmi)',color='blue')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[2]],'k',linestyle='dashed',label='Hydrus',color='blue')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm),'b',label='C (bmi)',color='green')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[3]],'k',linestyle='dashed',label='Hydrus',color='green')\n",
+    "\n",
+    "    axs[2].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm),'b',label='C (bmi)',color='purple')\n",
+    "    # plt.plot(HYDRUS_output[vars_hydrus[4]],'k',linestyle='dashed',label='Hydrus',color='purple')\n",
+    "\n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    \n",
+    "    axs[2].vlines(x=times_bmi[frame], ymin=-100, ymax=1000, colors='red')\n",
+    "    \n",
+    "    leg = axs[2].legend()#(handles=[runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS,runoff_line_HYDRUS],labels=[\"LGARTO\", \"HYDRUS\", \"theta_e\", \"capillary head = -depth\"])\n",
+    "\n",
+    "    mylabels = ['LGARTO runoff','LGARTO infiltration','LGARTO AET','LGARTO recharge', 'current time']\n",
+    "    axs[2].legend(labels=mylabels, ncol=3, bbox_to_anchor=(0.5, -0.0), loc='lower center')\n",
+    "\n",
+    "    axs[2].text(x=times_bmi[550], y=50, s='synth '+synth_case, fontsize = 16)\n",
+    "\n",
+    "    axs[2].set_ylabel(\"cumulative \\n output flux (cm)\",fontsize=16)\n",
+    "    axs[2].set_xlabel(\"date\",fontsize=12)\n",
+    "    \n",
+    "#     axs[2].plot(times_bmi, data_c_vars[vars_bmi[5]] * m_to_cm*12, label='C (bmi)', color='green')\n",
+    "\n",
+    "\n",
+    "    axs[2].set_ylim(min(np.cumsum(data_c_vars[vars_bmi[5]]*m_to_cm))-3, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[1]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[2]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[3]]*m_to_cm))))\n",
+    "#     axs[2].set_ylim(-0.05,0.15)\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "    axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm),'b',label='C (bmi)',color='black')\n",
+    "    axs[3].plot(times_bmi, np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm),'b',label='C (bmi)',color='gold')\n",
+    "\n",
+    "    axs[3].set_ylim(0, 1.05*max(max(np.cumsum(data_c_vars[vars_bmi[7]]*m_to_cm)), max(np.cumsum(data_c_vars[vars_bmi[6]]*m_to_cm))))\n",
+    "\n",
+    "    axs[3].vlines(x=times_bmi[frame], ymin=-100, ymax=1000, colors='red')\n",
+    "    \n",
+    "    axs[3].text(x=times_bmi[600], y=3, s='current time: '+str(round(frame/12,2))+' h', fontsize = 12)\n",
+    "#     axs[3].text(x=times_bmi[600], y=15-1, s=str(times_bmi[frame]), fontsize = 12)\n",
+    "    \n",
+    "    mylabels_forcing = ['precipitation', 'PET', 'current time']\n",
+    "    axs[3].legend(labels=mylabels_forcing, ncol=3)\n",
+    "    \n",
+    "    axs[3].set_xlabel(\"date\",fontsize=12)\n",
+    "\n",
+    "    axs[3].set_ylabel(\"cumulative \\n forcing flux (cm)\",fontsize=16)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3af740cd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "init()\n",
+    "\n",
+    "start_move_writing_time = datetime.now()\n",
+    "\n",
+    "ani = FuncAnimation(fig, animate, frames = 10000,  init_func=init, blit=False,repeat=False)\n",
+    "# ani.save('LGARTO_animation.mp4', dpi=600, fps=50, extra_args=['-vcodec', 'libx264'])\n",
+    "ani.save('LGARTO_animation_synth_'+synth_case+'.mp4', dpi=400, fps=48, extra_args=['-vcodec', 'libx264', '-preset', 'veryslow', '-crf', '17'])\n",
+    "\n",
+    "# ani.save('LGARTO_animation.gif', fps=50)\n",
+    "\n",
+    "end_move_writing_time = datetime.now()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5b9c4b9b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#time it took to write movie. It can be a computationally intensive thing \n",
+    "end_move_writing_time - start_move_writing_time\n",
+    "#mp4 blit=False time: 174 seconds for 200 frames\n",
+    "#then edited with 'veryslow' and crf of 17, 188 seconds for 200 frames\n",
+    "#next trying with 30 fps and 300 dpi, using the slower extra_args ... 132 seconds for 200 frames\n",
+    "#dpi=300, fps=60, slow args: 135 seconds for 200 frames, and doesn't have the freezing problem, but does kind of feel jerky. Maybe too fast, 30 fps version looked smoother. But it's a bit slow, so 45 is a good medium \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c2b2537c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f960d57e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#extra_args=['-vcodec', 'libx264']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "553247b3",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2885dd2f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e4ba966c",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "36aefd84",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7f98f67c",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f1c884db",
+   "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.8.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}