diff --git a/examples/dryland_maize/CheckOutput.ipynb b/examples/dryland_maize/CheckOutput.ipynb new file mode 100644 index 00000000..ac4dcc13 --- /dev/null +++ b/examples/dryland_maize/CheckOutput.ipynb @@ -0,0 +1,128 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEnCAYAAAB/kO72AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3zV1f3H8ddJbhYJCSOBbBJGEoYCMhwFREBBEVD81fWrtbUVtdpa2/6c9VdtrbaOWkerP1ct2ooDFK0ICrhl7xHCSAhkQiAhBLLv+f3xuTc3gYSVcW+Sz/PxOI/v937vOvcG8s453/M9x1hrUUoppVqTn7croJRSquPTsFFKKdXqNGyUUkq1Og0bpZRSrU7DRimlVKvTsFFKKdXqHN6uwKmKjIy0SUlJ3q6GUkqpJqxZs6bIWhvV2H3tJmySkpJYvXq1t6uhlFKqCcaY7Kbu81o3mjFmijEmwxiz0xhzr7fqoZRSqvV5JWyMMf7A34BLgUHAdcaYQd6oi1JKqdbnrW600cBOa20mgDFmDjAD2NoabzZ37ibWr88iODiA4GAHISEBhITItksX2Q8NbbgNDAwgICAAh8NRt3U4HBhjWqOKSinVoXkrbOKAvfVu5wDnHvsgY8wsYBZAYmLiGb/ZH+58jA25b53x8xtyYIwDP7+Auq2fnwN//wD8/ALw95d9f38HDkeAqzjqtgEBAQQEyDYwMICgIAm40FAHYWEBhITI/YGBgXTr1o2ePXvSrVs3/P3964qfn1+DrbzmqRcNTKVUW/NW2DT22+64GUGttS8BLwGMHDnyjGcMffnaEaT/Yy/lB0spx0F5fDIV/QdS0T+VytBwKiqqqayscW2rqaqqoaqq/raa6uoaqqtlW1NTTXV1NTU1sl9bW0NtrWyrqqqpra3G6azB6awGaoBq4Ei9/WO3jR2rauwraRHBwSFERvaid+/exMfHkZSUSGJiw9KrVy/8/HRkvFKqZRhvzPpsjDkfeMhaO9l1+z4Aa+1jTT1n5MiRtlmj0ayFTZvgo4+krFghx0ePhhtugGuugahGR+w1S00NVFZKqao6fv/IEThwoGEpKYGDBy0HDpRy4MABKioOUVXlpLKylqqqWqqra6mqclJVVYu1tXgC61TLEaAQKABygWzXMQ9jAgkKSiA0NJHw8DgiIqLo0aMXsbEJJCb2oV+/JAYMiCE62p/ISOjWDbTBpFTnZoxZY60d2eh9XgobB7AdmIj8tlsFXG+t3dLUc5odNsfKzYW33oI33oCNG8HhgClTJHimTYOQkJZ7r1ZUU+MJrqa29fcrKqC8HI4e9ZQjRyzFxSXs37+HAwf2UFy8h9LSPZSV7eHIkWyqqgqoqdnHsYEk40sigN4Yk0JIyACiolLo2zeFs85KYdiwGAYMMAwYAL16aRgp1dH5XNgAGGMuA/4K+AOvWWv/eKLHt3jY1LdxI7z5JvzrX5CXB2FhMHUqXHUVXHqp3FaUlh4hPX0v6enZbN++m927cygqKqGgIJeCgh0cPLiT2tqKes8IBVKB0QQHp5KUNIDU1FRSUpKIj3cQFwcxMRAbC9HREBzspQ+mlGoRPhk2p6tVw8atthY+/xzeeQc++AD275ffgFOmwMyZ8L3vQXKy/oneBKfTSU5ODtu3byc9fTtr1+5g/fpNZGSsprz8UL1HBgD9kCBylxTCwlLp1SuSmBhDXBzEx1O3jY+HlBSIjPTGJ1NKnQoNmzNRUwPffANz58K8edLiAejRA0aOlDJqlGzj4jSATsBay4EDB9i+fTsZGRls25bBpk2y3bt3JzU11XWPDQqKpEuXYcBwysqGU109HBiANIDl609NhX79oH9/2br3IyP1x6CUN2nYNJfTCRs2wKpVUlavlsEGtbVyvzHQpYt0t4WGHr8991y46SY5caEaqKmpITs7m4yMDDIyMtiyZQvr1q1j8+bNVFVVARAc3IWkpKH07Dkcf//hHDkynMLCweTmBlP/n2/Xro2HUL9+0jLSwXVKtS4Nm9ZQXi4BtGYNFBbKsLKysuO3hw5Bero8JyoKEhPlBEVUlHTRBQZ6SlCQbAMCGpbGjp3seDs/CVJVVUV6ejrr1q2rK+vXr+fw4cMA+Pn5kZSUTEJCGt27pxIUlEZ1dRolJWns3RvJ7t2Gak+DicBA6Nu3YQilpsrfARERXvqQSnUwGjbetmYNfPmlhE5uLhQUyPmg+sPGqqo8LaWWEBYGM2bAmDHyW/Xss6Fnz5Z7fS9wOp1kZWWxbt06Nm3a5OqS20ZGRgYVFZ6BCT169CAlJZWEhDS6dUsjICCNyso0DhxIJisrgJ075W8BkNbOsGEwbhxMmCBbDR+lzoyGTXtRWyuhU13dsJzusaoqWLZMzjUdPOh5/fh4Ocd0+eUwaZLc9vf33udtIU6nk71797Jt27bjSkFBQd3jHA4H/fv3Jy0tjYSENLp2PZuSkmFs2ZLCihX+VFRI+IwYARddBOPHy5iQ8HDvfTal2hMNm87K6ZSWVHq6DO9evx6+/hr27JH7HQ5ISICkJCl9+jTcxsVJt1w7dujQoboWkLtkZGSwY8cOql39bMHBwQwaNISYmGHU1g4lJ2cg27alUFMTh7+/H+eeK9k8aZJ0uwUGevlDKeWjNGyUh7USOitXQnY27N7t2ebn0+CMuzHS/Xb33XD11RJOHURVVRXbtm1j/fr1bNiwoW574MCBuscEBgYTEdGfmppUSkqGYe1wQkKGMX58LBdfbLj4Yhg8WEfAKeWmYaNOTWUl7N3rCZ89e+C992DrVhgwAP7rvyAtDS67rENe8GKtJS8vj4yMDLZv386OHTvYsWMHW7duZdeuXXWP8/ePorZ2BDCK7t1HMWnSKK67LpopU9rNxBNKtQoNG3XmnE65wPXxx2WgQ02NtHAGDoRbboHJk2V4Vwf/8760tJSNGzfWjYxbtmw1GRlbsNbpekQ8/v6jGDhwFDNmjOLWW0cSH9/Nq3VWqq1p2KiWUV0NmzfLha6LFsn1RiDDuMeNgx/9SKb56eDB43bkyBHWrVvH8uWr+Pjjlaxdu4rSUk8LKCwshaFDRzFt2ijGjh3F8OHDCdGmj+rANGxUy7MWtmyB776T8umncs7nxz+GV1/tNIFzrP37DzJ79mree28V69evoqJiFSCzT/j5+XPWWWdz4YVjGTduHGPHjqWXXuirOhANG9X6qqvhd7+Dxx6DP/4R7r/f2zXyOqdTcvi11/L44INVFBevws9vGcYso7a2HICBAwcyadIkJk6cyPjx44nQi3xUO6Zho9qGtXD99TKo4IMPpEtNARI8334rc7y++24VhYVrcTi+pGfPzykp+ZrKyqP4+fkxatQoJk6cyMSJE7ngggsIbsezQKjOR8NGtZ3iYhg7VrrYxo2Du+6S9YE6wMWjLaW21hM8770HhYWVOBzL6dNnCbCE3btXUFtbS3BwMGPGjGHixIlMnjyZYcOG6ZLeyqdp2Ki2VVkJf/sbPPOMDJ/u21cmIv3Rj+RCUVWntlYme/jgA3j/fcjMBD+/Us466ysiI5eQl7eY9PTNAMTFxXH55Zczbdo0JkyYoIMNlM/RsFHeUVMjv0H//nf44guZjeCnP4VHH5V1pFUD1spED+++C2+/DTt3SoPwvPMKSE7+hIMH/8NXX31KWVkZISEhTJo0iWnTpjF16lRiY2O9XX2lNGyUD9i1C/7yF3jxRejdG55/XhakU42yViYVf+89meIuPV0G+J17biVnn/0lFRUf8eWXH5GdnQ3AiBEj6lo9w4cPx0/XU1Be0GphY4x5ApgGVAG7gB9ba0tc990H/ASoBX5hrV3kOj4CeB0IARYAd9pTqISGTQexejXcfLNMmXPFFRI62rV2UunpEjpz58K6dXJs+HDLmDFbCAj4D8uXf8SyZcuw1hIZGcmECRPqSv/+/fVcj2oTrRk2lwBLrbU1xpg/A1hr7zHGDALeAkYDscBiIMVaW2uMWQncCSxHwuZZa+0nJ3svDZsOpLoann5ahkoHBMCf/gS33dZpr805XZmZEjzz5sn5HoBBg2Dy5P306LGQHTsWs2TJEnJzcwFISEhgwoQJTJw4kQkTJhCn4a5aSZt0oxljrgT+y1r7365WDdbax1z3LQIeAnYDn1tr01zHrwPGW2tvOdnra9h0QLt2wa23wuLF8MMfwiuvtPtZpttabq6cFps3T5ZMcjplPMaVV1pGj97B/v1L+OKLpXz++ed1k4ympqbWtXouuugierbzdY6U72irsPkIeNta+6Yx5nlgubX2Tdd9rwKfIGHzJ2vtJNfxscA91trLm3jNWcAsgMTExBHu/mnVgVgLf/iDtHImT5Z+otBQb9eqXdq/H+bPl+BZvFgakHFxcOWVcOWVTiIiNvLll0tZsmQJX331FWVlZRhjGDp0aF2rZ+zYsXTt2tXbH0W1U80KG2PMYiC6kbsesNbOdz3mAWAkMNNaa40xfwOWHRM2C4A9wGPHhM3d1tppJ/sQ2rLp4F55RSb2/N734D//0RXLmqmkRL7GefNg4UJZxTwqShZvveoqGDu2mg0bVrF0qYTPd999R1VVFQ6Hg/POO49JkyZx8cUXM2rUKAK0talOUau2bIwxNwK3AhOttUddx7QbTZ2+t9+GH/wAzjkHPvkEevTwdo06hCNHJHDmzpUAOnxYlr6eNk2CZ/JkgHK+/fZblixZwuLFi1mzZg3WWrp27cr48eO5+OKLmTRpEmlpaTrYQDWpNQcITAH+Alxord1f7/hg4N94BggsAQa4BgisAn4OrEBaO89Zaxec7L00bDqJDz+E739fljBYuBCiG2tUqzNVUQFLlkjwzJ8vq4Z36SJLFF11lWzDw+HgwYMsXbqUxYsXs3jx4rr1fBITE5k0aRKjRo1ixIgRnHXWWTqljqrTmmGzEwgC3MsbLrfW3uq67wHgJqAG+KV7xJkxZiSeoc+fAD/Xoc+qgU8/lWHRTqcs2HbzzTL1jf5F3aKqq+GrryR43n8fCgpkyetLLpHgmTYN3GMHsrKy+Oyzz1i4cCFffPEFxcXFADgcDgYPHsyIESM455xzOOeccxgwYAA9e/bUFlAnpBd1qvYnPV2mvHnzTTh0SFYKvflmGb2mJ7BbnNMpw6jnzpXzPNnZMnvBRRdJ8FxxhaeRaa0lOzubNWvWsGbNGtauXcuaNWsoKiqqe72uXbvSt29fkpOT6du3b11JTk4mKSlJW0MdlIaNar+OHpXL6F9+Gb75Rs5yP/gg/OxnOrlnK7EW1q6V4Jk7F7Zvl0bl974nkz5ceSUkJR37HMvevXtZv349mZmZZGVlkZmZWbdfXl7e4PFxcXH07duX1NRUUlNTSUtLIzU1leTkZBwOR9t9WNWiNGxUx7ByJdx3HyxdChdcAP/4B6SkeLtWHZq1sHWrp8WzYYMcHzhQRrb99KeyKviJX8NSWFjYIHwyMzPZuXMnGRkZ7N9fd7qXgIAA+vfvXxc+9YOoe/furfhJVUvQsFEdh7Xwr3/Bz38us0s/9pjs61xgbWLnTvjoIxksuHSpzFp98cUyan369DO7JvfgwYNkZGSQkZHBtm3b6ra7du2iurq67nHdu3cnJiaG6Ojo40pcXFxdCQsLa8FPrE6Hho3qePLyYNYs+PhjWT9n9uzj+3ZUq8rNlRXAX34ZcnLknM5PfiKn1vr0af7r19TUkJWVVRdAWVlZFBQU1JX8/PzjuucAIiIiSEhIqDs/lJSUVLefnJxMN51xvNVo2KiOyVr45z/hzjvlT+r33oPx471dq06ntlZaOi++CAtcFzFMniyLts6Y0XrX51prOXz4MPn5+eTm5taVnJwc9u7dS1ZWFllZWZSVlTV4XkRERIMAOjaMdAaFM6dhozq2HTukD2fnTnjuORmxprxizx6ZDOKf/5T9oCC5dufaa+Hyy+WanrZkraW4uJisrCx2797N7t27j9s/evRog+f06NGjLnji4+Pp2bNnXenRo0fdfq9evXQBu2No2KiO79Ah+VN6wQIZqfbXv+qknl7kdMLy5TIpxDvvyDU8oaFw6aUyjHrqVN9YP89aS1FRUZNBlJeXR2lpaaPPdTgcXHDBBSQlJREbG9ugxMTEEBMTQ1BQUBt/Iu/SsFGdQ20t3H8/PP64XCDy7rueqxKV19TWwtdfw5w5MmtBQQE4HHDhhRI+U6bIEgm+eg1odXU1Bw8e5MCBAw1KRkYGX3/9NXl5eeTn5zcYzODWs2dP4uPjSU1NZfjw4UyZMoWhQ4d22AteNWxU5/LGG3KWOjJSwuemm0AvIvQJTqeMYP/gAxnVtnWrHI+Pl/M8kyfDpEnQ3kY5O51ODhw4QF5eXl3Jz88nLy+P7Oxstm3bRmZmJgAxMTEMHjyYiIgIIiIiCA4OJigo6LSL+3nBwcEEBwcTEhJSt++tlVo1bFTns2oV3HUXfPut/CZ78EH48Y+1a83H7N0LixZJ+ewz6Q3184Nzz5UWzyWXwLBhHeNvhYKCAhYtWsTChQvJzs7m0KFDHDp0iIqKCiorK6msrGy0dXQm3CHkDiD3tn5IhYSEEBsbS3R0NF27diUsLIyUlBTGjRt3xu+rYaM6J2vh88/ht7+VuVj69oWHHpJzOzr7gM+pqZFWz8KFUlavlh+hwwGDB8tk4O4ydGjHXPbI6XRSVVVVFz6VlZXH3a6srKwLqIqKigalvLz8hNv6zzty5Ah5eXl189wBXHPNNcyZM+eM669hozo3a2XgwG9/C+vXw3nnyQUigwZ5u2bqBIqKZPXRtWulrFkjC8SBnN9JS2sYQMOG+cagg/amurqasrIyDh8+jMPhIDY29oxfS8NGKZATBm++Kd1rZWUSPvfeq11r7YS1ci2vO3zcJSfH85h+/TzhM3CgXOfbp4+GUFvRsFGqvn375ELQOXPg7LPh97+X63Q66Aihjm7fPli3rmELKCur4WMiIjzBk5TUcD86WgYtdrJRyq1Cw0apxnzwAfz615CZKScBHn5YQ6eDKC6Wa3yzs2H3bin194+ZVACQC0579PCUnj1lQKO7REU1vB0ZKc/Rfy4eGjZKNaWmBv79b3jkEZmJYOxYeOIJGQ6lOiRrJYzc4bNvn6xYeuCAbN2lqEiOFRVJD2xjgoObDqL6IdW9uwxo6NLFUwIDO15QadgodTI1NTLPyu9+J799pk6Fu++WFUJVp+Z0QkmJhE5jZf/+44+VlJz8df39G4aPuxwbSk0d69JFws5dQkIa3q5fgoLaZmL0Vg8bY8xvgCeAKGttkevYfcBPgFrgF9baRa7jI/AsC70AuFOXhVY+4/BheOYZKUVFcN11Mt+azkSgTkN1tad1VFQk+0ePHl+OHDm1Y+7jzfl1HRjoCZ76IRQQIMPLAwLkgtqHHz7z9zhR2DR7STxjTAJwMbCn3rFBwLXAYCAWWGyMSbHW1gIvALOA5UjYTAE+aW49lGoRXbvKKLVf/xqefBL+8Ae5Vue992SpSqVOQUAA9O4tpaVYK0s41Q+gigpPKS9vuF9eLo+vqPBsjy2VldKor66W0pqXn7XE+qtPA3cD8+sdmwHMsdZWAlnGmJ3AaGPMbiDcWrsMwBgzG7gCDRvla0JCZNaB6dPhv/5L5lr7+99laUqlvMAYT4ukRw9v1+b0NasXzxgzHci11m445q44YG+92zmuY3Gu/WOPN/X6s4wxq40xq+svHatUmxk6VC5rv+gimW/tjjvkT0Cl1Gk5adgYYxYbYzY3UmYADwD/29jTGjlmT3C8Udbal6y1I621I6Oiok5WVaVaR/fusiLob34Df/ubTNhVWOjtWinVrpw0bKy1k6y1Q44tQCaQDGxwdY/FA2uNMdFIiyWh3svEA3mu4/GNHFfKtzkcMiT6jTc886w99ph0eCulTuqMu9GstZustb2stUnW2iQkSM6x1hYAHwLXGmOCjDHJwABgpbU2HzhsjDnPyIIOP6ThuR6lfNsPfgAbNsiUxPffL9flbN/u7Vop5fNaZeS1tXYL8A6wFVgI3O4aiQZwG/AKsBPYhQ4OUO1NairMnQtvvQUZGTID5MsvN29cqlIdnF7UqVRz5OXBjTfC4sUyiOD553U2adVpneg6G+8s56ZURxEbKyt/vfCCLF8wdCjcc0/jk28p1Ylp2CjVXH5+cOut0qV2ww3w+OMym/Q333i7Zkr5DA0bpVpKVBS89hp8/bVcgTdunKyXU1np7Zop5XUaNkq1tDFjZMTaT38Kf/6zzCC9ebO3a6WUV2nYKNUawsLgpZfgo48gPx9GjICnnmp6rnqlOjgNG6Va0+WXS6vmsstkBoKJE2UhFaU6GQ0bpVpbVBTMmwf/+IesWXz22TB7tl6XozoVDRul2oIx8KMfybmcoUPl2pypU2HdOm/XTKk2oWGjVFtKTpb1cZ58EpYvh5Ej4Ze/lEXblOrANGyUamv+/rI4W1aWXJ/z7LMy68AHH3i7Zkq1Gg0bpbwlIkKWLFi2TFbDuvJKmDED9uw5+XOVamc0bJTytnPPhdWrZQmDxYullfP007p8gepQNGyU8gUBATI0eutWGD8efvUrGD0aVq3yds2UahEaNkr5kj595ELQ996DggJp9dx6Kxw44O2aKdUsGjZK+Rpj4KqrYNs2Gan2yiuQkiIzEtTWnvz5SvkgDRulfFV4OPzlL7J0wVlnwS23wHnnwcqV3q6ZUqet2WFjjPm5MSbDGLPFGPN4veP3GWN2uu6bXO/4CGPMJtd9z7qWh1ZKNWXIELk259//htxcCZybb4aiIm/XTKlT1qywMcZcBMwAzrbWDgaedB0fBFwLDAamAH83xvi7nvYCMAsY4CpTmlMHpToFY+C666Rr7Ve/gtdfl661F1/UrjXVLjS3ZXMb8CdrbSWAtXaf6/gMYI61ttJamwXsBEYbY2KAcGvtMivrUc8GrmhmHZTqPMLDZfaBDRtg2DC47TYZtbZ8ubdrptQJNTdsUoCxxpgVxpgvjTGjXMfjgL31HpfjOhbn2j/2uFLqdAwaBEuWwJw5Mmrt/PNlvjWdUVr5qJOGjTFmsTFmcyNlBuAAugPnAf8DvOM6B9PYeRh7guNNvfcsY8xqY8zq/fv3n9IHUqrTMAauuUaWo77nHnj7belau+su0P8vysecNGystZOstUMaKfORlsk8K1YCTiDSdTyh3svEA3mu4/GNHG/qvV+y1o601o6Mioo6/U+nVGcQFgZ/+hPs2AE33CBzrfXrB7//vU7wqXxGc7vRPgAmABhjUoBAoAj4ELjWGBNkjElGBgKstNbmA4eNMee5WkA/BOY3sw5KKYCEBLkmZ/NmuPhi+N3vJHSeew6qq71dO9XJNTdsXgP6GmM2A3OAG12tnC3AO8BWYCFwu7XWPWTmNuAVZNDALuCTZtZBKVXfwIEwdy6sWCHDpn/xC7lOZ8ECXbBNeY2x7eQf38iRI+3q1au9XQ2l2hdrJWR+9SvYvh0mT5YLRQcN8nbNVAdkjFljrR3Z2H06g4BSHZkxsiLopk0yk/SKFbIs9R136CAC1aY0bJTqDAIDZZ61HTtkYs8XX5RVQ++/Xyf5VG1Cw0apziQyEp5/XgYRTJsmo9iSk+HBB6Gw0Nu1Ux2Yho1SnVFaGrz1FmzcKOdxHnkEEhPhpptkdgKlWpiGjVKd2ZAh8O67MufaT38qF4YOGwZDh3oWc1OqBWjYKKUgNRX+9jfIyZHlqXv3hmeegcGDZSqcV17RC0RVs2jYKKU8uneXFs2nn8pyBk89BaWlsqRBTAz85Cfw3Xd6vY46bRo2SqnG9eol1+ds3gzLlskSB2+/Dd/7nrR4nnhCWkJKnQING6XUiRkjC7a9/LLMMP3qq9CtG9x9twwqGD9elqw+eNDbNVU+TMNGKXXqwsJkxNp338k1Ow8/LAF0yy3SzXbVVfD++1BR4e2aKh+jYaOUOjP9+8v1OenpsHo1/Oxn8M03MHMmREXB9dfDvHlw9Ki3a6p8gIaNUqp5jIERI2Q6nNxcWLRIzu989pm0dHr1knV33n0Xjhzxdm2Vl+hEnEqp1lFTA199Be+9Jy2cwkIICYFx42DCBJg4Ua7p8ff3dk1VCznRRJwaNkqp1ldbK11s8+bB4sWei0W7d5cBBhMnSgClpUlLSbVLJwobR1tXRinVCfn7w4UXSgEZVLB0KSxZIuX99+V4TIyEzpgxMgJuyBBw6K+pjkBbNkop78vMlNBZulTKvn1yvEsXGDVKgsddoqO9W1fVJO1GU0q1H9ZCVhYsXy5lxQpYt86ztHViIpx7rgTPuefCOefIuSDlda0WNsaYYcCLQDBQA/zMWrvSdd99wE+AWuAX1tpFruMjgNeBEGABcKc9hUpo2CjViVVUwNq1EjwrVkgIZWfLfQ6HDDS44AIYPVq63lJTITjYu3XuhFozbD4FnrbWfmKMuQy421o73hgzCHgLGA3EAouBFGttrTFmJXAnsBwJm2ettZ+c7L00bJRSDRQUwMqVMpXOsmWwapXnmh4/P+jXT6bVGTRItoMHawi1stYcIGCBcNd+BJDn2p8BzLHWVgJZxpidwGhjzG4g3Fq7zFWx2cAVwEnDRimlGoiOhunTpYB0s2VkyEi3LVukbN0K//mPDMMGTwilpcmice6SlCTb8PAm3041T3PD5pfAImPMk8gFohe4jschLRe3HNexatf+sceVUqp5AgKkC23IkIbHq6pkap36AZSRAZ9/DmVlDR/bowf07SslOVnODyUkeEqPHjo0+wydNGyMMYuBxoZ/PABMBO6y1s41xlwNvApMAhr7adgTHG/qvWcBswASExNPVlWllDpeYKCnG60+a2Xy0KwsT8nMlO3atTIc2z0owS0kREInMRH69JHi3k9MhPh4eT91nJOGjbV2UlP3ubrB7nTdfBd4xbWfAyTUe2g80sWW49o/9nhT7/0S8BLIOZuT1VUppU6ZMdCzp5SRjZxmcDpl1oO9e48ve/bAxx/LeaNjXzM2tvEwcu930q665naj5QEXAl8AE4AdruMfAv82xvwFGSAwAFjpGiBw2BhzHrAC+CHwXDProJRSLc/PTy4yjYmRUW6NqayU8MnOlgDKzvbsr1olMyZUVTV8TrduDUMoIUFmUujaVYLIXbp0gaAguS80tN133zU3bG4GntQwouQAACAASURBVDHGOIAKXF1e1totxph3gK3IkOjbrbW1rufchmfo8yfo4AClVHsVFCSzX/fv3/j97tZRY2GUnS1zxx06dPL3CQmRCU1PpURFyfkrH6MXdSqllDeVlcnS2/XLoUNQXi7XFx08KDMqNFaOPafk1r378SHUtasM+w4OlpCsvx8YKCUhQS6SPUM6N5pSSvmqsDApsbGn9zxrJZSaCiJ32boVvvhCQq2y8sSvefXVsvR3K9CwUUqp9sgYOf/TrRukpJzac6yVc0gVFVIqK2VbVSWlFQcvaNgopVRnYYx0mwUFQUREm761rtSplFKq1WnYKKWUanUaNkoppVpduxn6bIzZD2Q34yUigaIWqk5baY91hvZZb61z22mP9dY6n5o+1tqoxu5oN2HTXMaY1U2N//ZV7bHO0D7rrXVuO+2x3lrn5tNuNKWUUq1Ow0YppVSr60xh85K3K3AG2mOdoX3WW+vcdtpjvbXOzdRpztkopZTyns7UslFKKeUlnSJsjDFTjDEZxpidxph7vV2fphhjdhtjNhlj1htjVruO9TDGfGaM2eHadvdyHV8zxuwzxmyud6zJOhpj7nN97xnGmMneqXWT9X7IGJPr+r7XG2Muq3ef1+ttjEkwxnxujEk3xmwxxtzpOu6z3/cJ6uyz37UxJtgYs9IYs8FV54ddx335e26qzj77PWOt7dAF8Ad2AX2BQGADMMjb9WqirruByGOOPQ7c69q/F/izl+s4DjgH2HyyOgKDXN93EJDs+jn4+1C9HwJ+08hjfaLeQAxwjmu/K7DdVTef/b5PUGef/a6R5erDXPsByMKO5/n499xUnX32e+4MLZvRwE5rbaa1tgqYA8zwcp1Oxwzgn679fwJXeLEuWGu/Ag4ec7ipOs4A5lhrK621WcBO5OfR5pqod1N8ot7W2nxr7VrX/mEgHYjDh7/vE9S5Kb5QZ2utLXPdDHAVi29/z03VuSler3NnCJs4YG+92zmc+B+/N1ngU2PMGmPMLNex3tbafJD/yEAvr9WuaU3VsT1893cYYza6utnc3SQ+V29jTBIwHPkLtl1838fUGXz4uzbG+Btj1gP7gM+stT7/PTdRZ/DR77kzhE1jC3f76hC871lrzwEuBW43xozzdoWayde/+xeAfsAwIB94ynXcp+ptjAkD5gK/tNaWnuihjRzzSr0bqbNPf9fW2lpr7TAgHhhtjBlygof7cp199nvuDGGTAyTUux0P5HmpLidkrc1zbfcB7yPN3EJjTAyAa7vPezVsUlN19Onv3lpb6PoP6wRextOt4DP1NsYEIL+0/2Wtnec67NPfd2N1bg/fNYC1tgT4ApiCj3/PbvXr7Mvfc2cIm1XAAGNMsjEmELgW+NDLdTqOMSbUGNPVvQ9cAmxG6nqj62E3AvO9U8MTaqqOHwLXGmOCjDHJwABgpRfq1yj3LxKXK5HvG3yk3sYYA7wKpFtr/1LvLp/9vpuqsy9/18aYKGNMN9d+CDAJ2IZvf8+N1tmXv+c2G4ngzQJchoyK2QU84O36NFHHvshokQ3AFnc9gZ7AEmCHa9vDy/V8C2meVyN/Lf3kRHUEHnB97xnApT5W7zeATcBG5D9jjC/VGxiDdHVsBNa7ymW+/H2foM4++10DZwPrXHXbDPyv67gvf89N1dlnv2edQUAppVSr6wzdaEoppbxMw0YppVSr07BRSinV6jRslFJKtToNG6WUUq1Ow0YppVSr07BRSinV6jRslFJKtTqHtytwqiIjI21SUpK3q6GUUqoJa9asKbLWRjV2X7sJm6SkJFavXu3taiillGqCMSa7qfu0G00ppVSr07BRSinV6tpNN5pSSrkdOnSIffv2UVxcTElJCaWlpRw+fBiAuLg4Bg4cSEJCwkleRbUlDRullM/bvXs3CxYsYMGCBaxZs4aCgoKTPiclJYVrrrmGq6++miFDTrTwpmoL7WaJgZEjR1odIKBU57FhwwbeeOMNFixYQHp6OgD9+vVj7NixpKWlERsbS/fu3enevTvh4eF07doVp9NJTk4O69evZ/78+XzxxRc4nU4GDRrEzJkzGTp0KOeffz5xcXFe/nQdkzFmjbV2ZKP3adgopXzFrl27WLp0Kf/4xz9YtmwZgYGBXHjhhVx22WVMnTqVAQMGnNbrFRYWMm/ePN5++22++uor9yJiREdHEx8fT3x8PHFxccTFxdXtu7ehoaGt8RE7NA0bpZRPKyws5MEHH+TVV1/F6XSSmprKLbfcwo033kiPHj1a5D3KyspIT0/nm2++YdOmTeTm5taVkpKS4x4fERFBdHR0Xevp2BITE8OYMWOIiYnB4dAzEqBho5TyUSUlJTz33HM88cQTlJeXc8cddzBr1izS0tIwxrRZPY4cOVIXPDk5OXX7hYWFFBcXNyglJSXU/71pjCEyMpLevXs3KNHR0fTu3ZuYmBj69etHnz59OnwonShsOvYnV0r5pJqaGl588UUeeOABSktLueKKK/jzn/9MSkqKV+oTGhpKSkrKKb2/0+mktLSUzMxMVqxYQUFBAYWFhRQWFlJQUMCyZcsoLCzk6NGjDZ4XFhbGmDFjmDBhAikpKXXnmYKDgwkMDCQoKOi4rcPhaNPQbU3aslFKtZn9+/fz9NNP83//938cPHiQyZMn89hjjzF8+HBvV63FlZWVUVhYSG5uLrt27WL16tV8/vnndYMdTpU7fJoKJPc2ODiYoKCg40pTx0NDQ+nduzdRUVGEhoYSFhZGeHg4YWFhZ/yZtRtNKeU1JSUlzJ8/n7lz5/Lpp59SVVXFzJkzueGGG5g+fXqH+cv9VBUWFpKTk8Phw4cpLS2lsrKSqqqq47aNHTv2PvftioqKutvucuwxp9N50rpdddVVvPfee2f82bQbTSnV5tLT03nmmWeYPXs25eXlJCQkcNttt3HLLbeQlpbm7ep5jfucTlurqalpEEBHjhyhsLCQoqIijhw5QllZGX369Gm199ewUUq1qGXLlvHwww+zaNEigoKC+MEPfsDNN9/M6NGjO10rxpc4HA7CwsIadJOd7lDyZr1/m72TUqpD27NnD/feey9vvfUWvXv35pFHHmHWrFlERTU647zqZDRslFLNsmnTJp599lneeOMNjDE8+OCD3H333c060aw6Hg0bpdRpy8vLY/bs2bz11lts3LiR4OBgbrzxRu6///5W7fdX7ZeGjVLqlFlreeqpp3jwwQepqKjg/PPP59lnn+W6664jMjLS29VTPkzDRil1SioqKrj55pt58803ufLKK3n88cfp37+/t6ul2gkNG6XUSeXn53PllVeyYsUK/vCHP/DAAw/oyDJ1WjRslFIntGLFCmbOnElJSQlz585l5syZ3q6Saoc0bJRSDdTW1rJgwQLmz5/P0qVLycrKIikpie+++46hQ4e22vseOgRLlsAXX4Br0U2M8ZSwMOjfHwYMkJKcLMdV++C1sDHGTAGeAfyBV6y1f/JWXZRSEjJvvPEGjz76KDt27CAiIoIJEyZw66238pOf/ISePXue1utVVMDnn8Pu3RARAQkJkJgIvXtDfj5kZkrZsQO++w6WL4faWggNhZ49wVpPASgpgSNHPK+fmgoPPAD//d/g59dy34NqHV4JG2OMP/A34GIgB1hljPnQWrvVG/VRqrNbvHgxv/nNb9iwYQPnnHMOc+bMYebMmQQEBJzW6xQVwccfw4cfwqJFDcOhKYGBcPbZcO+9cMklcP750NjbWgsFBRJOW7bAyy/DD38Is2fD7bfDhAkQHn5a1VVtyFstm9HATmttJoAxZg4wA9CwUaqN1NbW8uGHH/Lkk0/y3XffkZSUxJw5c7j66qtP6+R/YaEEy2uvwddfg9MJsbFwww0wYwYMHQrFxbB3r5TCQoiJgb59pcTFgb//yd/HGHleTAyMGwe33AIvvQT33AOLF0toTZsGP/gBXHopBAU148tRLc5bYRMH7K13Owc4t7XebGTKzazZMRtwYIwDcODnJ/vGBODn5ziu+Ps78PcPwN/fgcPhqNu6S0BAgGvbsAQGOggMDHBtHQQFebZBQQ6CgwMIDnYQHOy+Lc+r/5ru/fDw8Lr11f39/fH398fPz69uW39fRwapU2WtZf78+dx9993s2LGD5ORknn32WWbNmkXQKf6GzsyEF16ABQtgq+tPxP79pVtrxgw455yG51NiYmDQoJb9HH5+cOutcNNNsGwZvP8+vPUWzJ0rXXHjx8PFF0trKS1Nz+94m7fCprEf+3FrHRhjZgGzABITE8/4zS6LCMTBBVRjqPYPpCY0lOrgYGqDg6hxOqmpqaG2tn6prjvmdNbgdB7F6azB2hqsrQZqGilNHW8bEjyOesW/QTg6HHI7JCSYqKhIevXqRVRUFL169aJPnz4MGDCAhIQEevfuTWBgYJvVW7WtDRs28Otf/5olS5YwaNAg3nnnHWbOnIn/qTQtgNxceOQReOUVuT1hgnRljR8Po0Z559xJYCBceKGUJ5+UVs5//gOffipdeiBh178/xMfLfmSknBeqX9zH9J9/6/DKejbGmPOBh6y1k1237wOw1j7W1HOavZ5NZiZ89pmUJUvkbGNAAEycCFddJX+OncKEgdbKSczKSqiq8mzrl8pKd7GUlzs5erSaiooajh6toaKihvJy2cq+3HfkSA1lZfIYKVUcOXKIo0eLKS8/TGVlLVVVTmprawEn4N46kVCrdZVjw+7YY+VAEcbsB/ZhbfFxnzEwsAehoTGEh8cTERFN167dCQ/vTvfuPYiOjiEuLp4+feIZMCCa3r396dGj8T525RucTieLFi3i6aef5rPPPqNbt278/ve/57bbbjvlZYqLiuDPf4bnn5d//7Nmwf33S3eZL9u9W/7Lf/UV7NkjYZmfD8csotlAWFjjIeQuPXpA9+6ebWSk7GvLyQcXTzPSl7UdmAjkAquA6621W5p6TosunlZbCytXSrt77lwJIj8/6Qi+6iq48krpSPZBNTVy0vXIESgrk/80jQXeiUp5uec1Skur2L8/m6KiHRw+nEdZWQHl5flUVuZRXZ2D01kIFANljdTGD+gBROFwxNKlSz8iIpLo3TuepKR4Bg3qz7Bh8aSkGPr1g+DgNv2qOr3q6mrmzp3Lo48+yqZNm4iNjeWOO+7glltuoUePHqf0GqtWwfz58Oyz8u/lhhvgd7+TYcftWUUFHDggpajIs3+iY8XH/11WJzhYWk29ekG3bjL6zr3t2lVKWJgU9379bWiovMYpNjB9ls+FDYAx5jLgr8jQ59estX880eNbbaVOa2HDBgmduXPBvWTr0KEwdSpcdhmcey6c4l+AHY218h+zuLiavLxisrPzyc7OYc+eHPbuzWXfvv0cPLifAwdyKCnZRWVl0TGvEAqkAKlERKQSF5dKfHwK/foNIDExnN69qSuRkVK6dNG/EpujrKyM559/nqeffpp9+/YxcOBA7rvvPq655ppT7iLduFFOvC9cKLdnzoQ//KHlz7u0JzU1EjjFxXDwoGe7b5+0mPbulVAqKZFrhtzb6upTfw9/fxnYEBgo2zMpwcHHl5CQ448FBUmPhMMhJSBAgq9XrzP/jnwybE5Xmy0LvXWrdPh+/DF8+620grp3lzONo0ZJCA0d2ryfSAd29OhRcnNzyc7ew7p121m7NoP09Ayys7dRUpJNw1Nz0UgQDXBtpQQF9aN376C6kUexsVKSkmT0Ur9+0uOpgeRRXV3N559/zjvvvMO8efMoLi5mypQp3H777Vx22WX4neLJlNxc+O1v4Z//lL/KH3gAbrzxlHqYVRMqK6UXoqxMLlZtbFtWVr/7vfml5gxPF3//+/DOO2f+WTVszlRJiXT4fvyxnHXMzfXcFx3tCZ6YGGkHu0uXLg1vx8drHxJQXl7Orl272L59O9u3b2fr1u2kp+9g167tFBfvq/dIQ3j4AEJCRmDt2Rw9OpCysoFAMiAnh0JC5ALBhATPxYLubWIi9OnTeb7yTz/9lJ/97Gfs2rWLsLAwpk+fzh133MH555/f5HOcTli6FLZvl3/WxcXyV/Wbb0rX7C9+AffdJ+ciVPtTUyOhU1Ehpbzcs1+/1NRIqa6WbUKCDPY4Uxo2LaWoSPoXNmzwlK1b5UTIifTqJeeC0tLkXFCvXnKmMThY2svuNrN7v7133J6BkpISduzYURdEGzduZPXq1eTk5NQ9xuFw0KtXMt27pxAUlIK1/Tl6tB/FxX3Zt68P0LCLKDZWzi0kJ0uLqP5+bGz7/5rz8/O56667ePvtt0lJSeHRRx9l6tSpBJ8kZRculCBZv15u+/tLK6amBkaPhr//XaaDUep0adi0ppoaaQu7z7gfOSJ/GnrOwMO770qX3KFDp/aa/v4SOvU7V0/WEXvsGce0NBkL2s6HiR06dIht27aRnp7eIIx27NhBeXl53eP8/PyIiUmgV6++RET0w+HoS3V1X0pL+7JvXz9yc7tTf8R9QIC0fvr3l6vXzz5bGqmpqb71lVl7fHdhdnY2Tz/9NK+//joVFRXcf//93HPPPQQFBbFzJ7z6qvzV2q2b/E0TFeW5wHHOHHj7bQnchx6CSZPkb5/2HrzKN2jY+AJrYf9+GXe5b5+cSTx2GFn92+7O14oKz7a8vGHb2F2OHvV0ANfWet4zIgKGD5eTHP36wVlnwYgR0u3XzjmdTvLz88nMzGTXrl1kZmbWlV27drFv374Gj4+IiCAuri89e/YjJKQvxvSlvLwf+/b1ZdeuBKqrJWECA2WwQvfu0vpxd8slJsrgBYdDvvLSUul6sFa6pEDua+yXttPpmdcrOlpCLi1NfgzHBklVlTSe582T8Sq7d8vjBgyw9OmzkaKiN1m48DmstUyfPp1HH32U6uoBvP++dIFt2yb1CA6WfxLHCgyUczL33KPXk6iWp2HTWbiHjpWWwooVMtBhyxYZ2l1Q4Hnc4MEwZoz8Bo2Lk3NK8fGy30HWjS8rKyMrK6vRINq9ezdV9bo+/f39iYnpQ7dufXE4+uLv3w9j+lJW1pf9+/ty4EC3VqljSIj0kQ8cKC2NDz+UH5W1Elrjx1fTs+fXbNgwn6ysD6mq2o20zv6bfv3+yIwZiRw8CK+/Lq937rlw/fUyciw+XhrdBw9K76+7fz46Wn7sSrUGDRslAbRpk8zrsXAhrFsnv4mOFRHhCZ74eGkZ/fCHHWqGw9raWvLy8hoNoszMTIqKGg7f7tatG716JdCzZzyRkfHExiaQmBhPfHw8cXHxxMTEExbWlZoaaVge21oxRrq0unTxzHackSHbPXvkx7J9O4wdaxkxIpdDhz4mP/8TvvlmKYcPHyY4OJhJkyZx6aUz6NPncjIyolm0CL78UgLkrrvgf/5HB0gq79OwUY0rL5ehSLm5kJPjKe7be/dKiygsTCagmjpVzhwnJXXoccelpaUNWkVZWVnk5uayd+9ecnJyjuuiAwgPDyc+Pp5evXoRHBzcoAQFBREcHFx3jYvT6cTpdGKtxel0UlpaSkbGdrZu3cIh13m9Pn36MGXKFCZPnswll1xCaGjoce/pvrC3d+/W/T6UOlUaNurMrVoll4+//bbn6rSUFBm2dMMN8L3vyfDuTqSyspK8vDxycnKOK/v376eyspKKiooGpbKyksrKSowxdROnuvfDwsLo168fgwcPJjU1lUmTJjFo0CCdXFW1Oxo2qvmKiuT8z6ZNMtXvypUyyMHPT046nHceXHONDG/SX5JKdUoaNqrllZfLhKarV0vrxz20+/rr5fLzTjq9j1Kd2YnCRn8jqDMTEgKXXy4FZEj2n/4kF29ERMiVgUop5aJho1pGUJBMB1xWJouKOJ0ysVZCgrdrppTyAV5Y6kh1aI89JgvCv/aaXEh61VVyvc/JpvRRSnVoGjaqZTkcssLWzp3w85/LovTTpsnVhLfdJpe4K6U6HQ0b1ToSE+Gpp+R6nY8+knWBXn8dhgyREDpwwNs1VEq1IQ0b1boCA2UQwZtvyuXyt9zimVb42WdPb2UppVS7pWGj2k5UFPztb7I0w4gRcOedMt3yggXerplSqpVp2Ki2N2QIfPqpLG5fWyvT4Pz859rKUaoD07BR3mEMTJ8OmzfDr34lgwouuaTxyUGVUu2eho3yrsBAGUgwezZ8953MtZaV5e1aKaVaWLPCxhjzhDFmmzFmozHmfWNMN9fxJGNMuTFmvau8WO85I4wxm4wxO40xzxqdbVCBTOr52WdQWCjzrK1c6e0aKaVaUHNbNp8BQ6y1ZwPbgfvq3bfLWjvMVW6td/wFYBYwwFWmNLMOqqMYN07W2wkNhfHj4ZlntFtNqQ6iWWFjrf3UWlvjurkciD/R440xMUC4tXaZlRlAZwNXNKcOqoNJTYXly2UJg1/+UtZE/v73ZRYCHUCgVLvVkudsbgI+qXc72RizzhjzpTFmrOtYHJBT7zE5rmNKefTqBV98IauJ/uxnsj9tmlwo+sQTsmqYUqpdOWnYGGMWG2M2N1Jm1HvMA0AN8C/XoXwg0Vo7HPgV8G9jTDiygPqxmlzjwBgzyxiz2hizev/+/afzuVRHMGwYPP20rBz6wQcyZPruu2XOtXfegXayPIZS6hTCxlo7yVo7pJEyH8AYcyNwOfDfrq4xrLWV1toDrv01wC4gBWnJ1O9qiwfyTvDeL1lrR1prR0ZFRZ3pZ1TtXWAgzJghAwi+/Rbi42WhtunTZelqpZTPa+5otCnAPcB0a+3ResejjDH+rv2+yECATGttPnDYGHOeaxTaD4H5zamD6mQuuEDO6Tz1FCxdCoMGwXPPycWhSimf1dxzNs8DXYHPjhniPA7YaIzZALwH3GqtdQ8rug14BdiJtHg+QanT4XDIhaCbN8t1Ob/4BVx8MRQUeLtmSqkm6LLQqn2zVmaTvv12CA+Hf/0LJk70dq2U6pROtCy0ziCg2jdj4Mc/hlWroEcPaeE8/LB2qynlYzRsVMcweLAEzg03wEMPweTJMhuBUsonaNiojiM0VLrUXn1VRq0NGwYLF+oQaaV8gIaN6liMgZtukrnVunWDSy+VQQQffghOp7drp1SnpWGjOqazzoK1a2Wxtvx8uU7nnHPg44+1paOUF2jYqI4rJESmu9mxQ5YwKCuTJarHjJFrdZRSbUbDRnV8DocMHEhPh//7P1kv5/zz4Uc/0mtzlGojGjaq8wgIgFmzICMD7rkH/v1vSEmBP/0JDh/2du2U6tA0bFTn07WrBMyWLbKGzn33yYzSDz6o6+co1Uo0bFTnNWCArJOzciVMmACPPAJJSfDAA3qNjlItTMNGqVGjYO5c2LgRpkyBxx6DhAS4/nrYudPbtVOqQ9CwUcrtrLNknZz0dJlr7T//kTV0HnpIVwlVqpk0bJQ6VmqqLNqWkQFXXSVzrY0bJ6PYlFJnRMNGqabExMgs0u7WzrBh8O673q6VUu2Sho1SJ/P978O6dTBwIFx9tQyfPnLE27VSql3RsFHqVCQnw9dfw733wssvy1DpJ5/UczlKnSING6VOVUCAjFT79ls47zz4n/+B4cPhq6+8XTOlfJ6GjVKn64ILZELP+fNlvrULL4RbbtGpb5Q6AQ0bpc7U9OmwdSv8+tfwyivQty/86lcaOko1QsNGqebo0kXO3WRkwDXXwLPPyvmd3/wGSkq8XTulfEazwsYY85AxJtcYs95VLqt3333GmJ3GmAxjzOR6x0cYYza57nvWGGOaUwelfEL//vCPf0joXHst/OUvMsnnSy9Bba23a6eU17VEy+Zpa+0wV1kAYIwZBFwLDAamAH83xvi7Hv8CMAsY4CpTWqAOSvmGfv0kdNasgbQ0OZczYgR8+aW3a6aUV7VWN9oMYI61ttJamwXsBEYbY2KAcGvtMmutBWYDV7RSHZTynuHDJWDefhuKi2H8eLlGZ/dub9dMKa9oibC5wxiz0RjzmjGmu+tYHLC33mNyXMfiXPvHHleq4zFGAmbbNvj972WutbQ0WcpALwpVncxJw8YYs9gYs7mRMgPpEusHDAPygafcT2vkpewJjjf13rOMMauNMav3799/0g+jlE8KCZGAcc+19sgjMnLtvvtg0yZv106pNnHSsLHWTrLWDmmkzLfWFlpra621TuBlYLTraTlAQr2XiQfyXMfjGzne1Hu/ZK0daa0dGRUVdbqfTSnfkpAgc619950sa/DEE3D22XDFFbBqlbdrp1Srau5otJh6N68ENrv2PwSuNcYEGWOSkYEAK621+cBhY8x5rlFoPwTmN6cOSrU7558vXWqFhTKj9FdfwejRMHmyzkagOqzmnrN53DWMeSNwEXAXgLV2C/AOsBVYCNxurXWP/7wNeAUZNLAL+KSZdVCqferZE/73fyE7G/78Z1i/XmYjGDcOPvkEnE5v11CpFmNkUJjvGzlypF29erW3q6FU6ykvl5kIHn8ccnKk2+3GG+FHP5Ih1Ur5OGPMGmvtyMbu0xkElPIVISHw85/LUtRz5sCgQfDHP8oFoxMmyPme8nJv11KpM6Jho5SvCQqSqW8WLoQ9e2T0WnY2/OAHEBsLd9wBa9dCO+mVUAo0bJTybfHx8MADsGMHLFkCl10mXW0jRsjKoX/9q1w0qpSP07BRqj3w8/N0peXnw9//Li2gu+6CuDhp9bhnK1DKB+kAAaXas/Xr4YUX4L334OBBmbVgxAiYOBEmTZLRbQEB3q6l6iRONEBAw0apjqC2FlasgM8+k+62ZcugpgZ69ICZM+Uc0Pjx4HB4u6aqA9OwUaqzKSuDpUvhnXc8K4pGRclsBVOnSssnLMzbtVQdjIaNUp1ZeTksWCDndD75RIInMFC62MaMgZEjZZh1YqKcG1LqDGnYKKVEVRV8842Ez8KFsqy1+3dAcDAMGCAzU8fESCBFRcnFpe4SG6vngFSTNGyUUo0rLYUNG2QZhIwM2W7bBkVFUFkJFRUNH+/n5wmdoCDo2lUCKTLSs42OlhIXJ7d1Md5O40Rho2cLlerMh7GjvwAACF5JREFUwsNh7FgpjSktlalz9u6Vkp0NWVkyiWhlpVx0unYt7N8vraZjdesGqamekpIi2/79ZcYE1Wlo2CilmhYeLudzBg068eOslXNB+/dDQYFcC7RnD2zfLi2mxYth9mzP442Rc0TuEBo8WM4dnXWWdN+pDkfDRinVfMZIl1rXrrIwXGPKyjzh495mZMj6PmVl8pigIBg6VNb7cZfUVPD3b7vPolqFnrNRSnmXtdI1t2qVp6xZ41k6OyxMLlQdOdITQMnJei7IB+kAAaVU+1JbKwMV6gfQhg2e80I9e0r4jB4tQ7gvuEDPAfkADRulVPtXVQWbNjUMoC1bZJG5oCAJnAsvlG64s86S1o9eN9SmNGyUUh3T4cOylPbSpTJNz8aNnuuGQkNhyBAJnvolMtK7de7ANGyUUp3DkSPS2tm4UVpBmzbJ/oEDnsfExHiCZ9AgGYbdv78c1/NAzaLX2SilOofQUDmPM3q055i1MhzbHTzuEHr+eblWyC0kRJbfdodP//6e2wkJOiKumTRslFIdmzHSaomJgUsu8RyvqZGLVHftkqW43WX7dplDrn4Q+fvLQnZ9+sj1QX36eEpiopQuXdr+s7UjzQobY8zbQKrrZjegxFo7zBiTBKQDGa77lltrb3U9ZwTwOhACLADutO2lL08p1XE4HNJy6devYQiBDDrIzfUE0O7dEkx79sg5otxcGTFXX1SUhE5cnCfcji29e3faZR6a9amttde4940xTwGH6t29y1o7rJGnvQDMApYjYTMF+KQ59VBKqRblngMuIQEuuuj4+2tqIC9PAsgdQu79rCy5ULWo6PjnGSOh1FQYRUfLAIagIJkY1b0NDm7355NaJGKNMQa4GphwksfFAOHW2mWu27OBK9CwUUq1Jw6Hp/usqXnlqqpkDrn8/IbFPZ1Pfr6cOyooOL6VdKygIAmhnj1Pfdu1q08FVEu158YChdbaHfWOJRtj1gGlwG+ttV8DcUBOvcfkuI4ppVTHEhjoaR2diNMprSB3ABUXe2bcrqyU9YiKi+UxBw7IdtMm2R48KM9vTEBA4yEUEiLhFRjoKe7bAwbAxRe3/HfBKYSNMWYxEN3IXQ9Ya+e79q8D3qp3Xz6QaK094DpH84ExZjDQWMw2eb7GGDML6XIjMTHxZFVVSqn2x88PevWSMnTo6T3X6YSSkoZBdOzWvb91q2wrKqTVVVnpuSbJ7eqrvRc21tpJJ7rfGOMAZgIj6j2nEqh07a8xxuwCUpCWTHy9p8cDeSd475eAl0CuszlZXZVSqlPx84MePaScLmul+66qylNacfBCS7zyJGCbtbaue8wYEwUctNbWGmP6AgOATGvtwf9v72xC46rCMPy8pGkVW9DaKqGIpuKmiNQsSkHJUttsgruurOBGUdCFi0hBqjsLdeNCUBSKiN1UsRvBIoo7a9QkTSmxqVb8KY1QxHbjT/1c3BM6xrmjcTpzzu28D1zumXMvk4d3mPnmnrk5R9JFSTuBT4GHgZevgoMxxpjVIFXFZc2avty2fTWKzR7+PoQGMA68IOkP4DLwWERcSMce58qtz+/jmwOMMeaap+tiExGPtOk7AhypOX8auLvbv2uMMaY5eEpUY4wxPacxE3FK+gn4toun2AS0+S+rommiMzTT2879o4nedv5v3B4Rm9sdaEyx6RZJ03WzkZZKE52hmd527h9N9LZz93gYzRhjTM9xsTHGGNNzBqnYvJpb4H/QRGdopred+0cTve3cJQPzm40xxph8DNKVjTHGmEwMRLGRtEvSgqRFSVO5feqQdFbSCUkzkqZT30ZJxySdTvubMju+IWlJ0nxLX62jpGdT7guSHsxjXeu9X9IPKe8ZSRMtx7J7S7pN0keSTkk6Kemp1F9s3h2ci81a0nWSjkuaTc7Pp/6Sc65zLjZnIuKa3oAh4AywFVgLzALbcnvVuJ4FNq3oOwBMpfYU8GJmx3FgDJj/N0dgW8p7HTCaXoehgrz3A8+0ObcIb2AEGEvtDcBXya3YvDs4F5s11Wz061N7mGrexp2F51znXGzOg3BlswNYjIivI+I34DAwmdlpNUwCh1L7ENVic9mIiE+ACyu66xwngcMR8WtEfAMsUr0efafGu44ivCPiXER8kdoXqZZa30LBeXdwrqME54iIS+nhcNqCsnOuc64ju/MgFJstwHctj0tesC2ADyR9ntbyAbg1Is5B9UYGbslmV0+dYxOyf1LSXBpmWx4mKc5b0h3AvVTfYBuR9wpnKDhrSUOSZoAl4FhEFJ9zjTMUmvMgFJtVLdiWmfsiYgzYDTwhaTy3UJeUnv0rwJ3AdqoF/w6m/qK8Ja2nmtj26Yj4pdOpbfqyeLdxLjrriLgcEdup1tjaIanTZMElOxeb8yAUm++B1nVZOy7YlpOI+DHtl4B3qS5zz0saAUj7pXyGtdQ5Fp19RJxPb9g/gde4MqxQjLekYaoP7bci4p3UXXTe7ZybkDVARPwMfAzsovCcl2l1LjnnQSg2nwF3SRqVtJZq/Z2jmZ3+gaQbJG1YbgMPAPNUrnvTaXuB99o/Q1bqHI8CeyStkzRKtYje8Qx+bVn+IEk8RJU3FOItScDrwKmIeKnlULF51zmXnLWkzZJuTO3rSQtCUnbObZ1LzrlvdyLk3IAJqrtizgD7cvvUOG6lultkFji57AncDHwInE77jZk936a6PP+d6tvSo50cgX0p9wVgd2HebwIngDmqN+NISd7A/VRDHXPATNomSs67g3OxWQP3AF8mt3ngudRfcs51zsXm7BkEjDHG9JxBGEYzxhiTGRcbY4wxPcfFxhhjTM9xsTHGGNNzXGyMMcb0HBcbY4wxPcfFxhhjTM9xsTHGGNNz/gJIJ3h92ig6LwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "#######\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "df_2001_c = pd.read_csv(\"./outputs/02001sccd1\", delim_whitespace=True, header=0)\n", + "df_2002_c = pd.read_csv(\"./outputs/02002sccd1\", delim_whitespace=True, header=0)\n", + "df_2003_c = pd.read_csv(\"./outputs/02003sccd1\", delim_whitespace=True, header=0)\n", + "\n", + "fig = plt.figure(figsize=(6.5, 7.5)) # Adjust the figure size as needed\n", + "ax1 = fig.add_subplot(311)\n", + "var='CO2_FLUX'\n", + "ax1.plot(np.arange(len(df_2001_c[var])), df_2001_c[var],color='red', label='2001')\n", + "ax1.plot(np.arange(len(df_2002_c[var])), df_2002_c[var],color='blue', label='2002')\n", + "ax1.plot(np.arange(len(df_2003_c[var])), df_2003_c[var],color='black', label='2003')\n", + "ax2 = fig.add_subplot(312)\n", + "var='NBP'\n", + "ax2.plot(np.arange(len(df_2001_c[var])), df_2001_c[var],color='red', label='2001')\n", + "ax2.plot(np.arange(len(df_2002_c[var])), df_2002_c[var],color='blue', label='2002')\n", + "ax2.plot(np.arange(len(df_2003_c[var])), df_2003_c[var],color='black', label='2003')\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEnCAYAAABsR64CAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXwV1f3/8dcnIQs7hH0PKi4IIoq4gIqKG1qxVisqFquWrxbqgm0Bta4/FFzQuot1wa24QAsirQp1K1g2QbaAbGEJIQlr2EOS8/vj3EDABCKXZOaG9/PxmMfcO3fuvZ8cIG/OnJk55pxDREQkGnFBFyAiIrFPYSIiIlFTmIiISNQUJiIiEjWFiYiIRE1hIiIiUasSdAEA9evXd6mpqUGXISIiBzBz5sx1zrkGJb0WijBJTU1lxowZQZchIiIHYGYrSntNh7lERCRqChMREYmawkREpBIrLCxk6dKl/OMf/2DdunXl9j2hGDMREZGfb/fu3axatYqcnBw2bNjA2rVrWb58Odu3b2fz5s3MnTuXefPmsW3bNgD++c9/0rNnz3KpRWEiIhICBQUF5Obmkpuby+bNm9m0aRMrVqxg6dKlrF69es9rW7Zs2fM4IyODgoKCfT4nLi6OpKQkqlevTvv27bnllls46aSTaN++Pe3bty+3+sscJmYWD8wAMpxzl5tZCvABkAqkA792zm2M7DsYuAUoAO5wzn12mOsWEQkl5xzLly9nzZo1ZGZmsnbtWtavX8/mzZv3BEVJ66Lew/7MjMaNG1O7dm1q1apFrVq1aNSoETVr1qRFixYcddRRNGrUiJSUFBo0aEDLli1JTEys4J/65/VM7gTSgFqR54OASc65oWY2KPJ8oJm1BXoBJwJNgYlmdqxzrqCkDxURqQwKCgoYMWIEQ4YMISMj4yev16xZk1q1au0JhTp16tCqVat9thV/XLt2bVq2bElqaipJSUkB/EQ/T5nCxMyaA5cBQ4ABkc09gW6RxyOBr4CBke2jnHO7gOVmtgToDHx32KoWEQmB7OxsxowZwz/+8Q8WLlzIypUrOeecc3jwwQdJTU2lcePGNG7cmJSUFOLj44Mut1yVtWfyLPBnoGaxbY2cc5kAzrlMM2sY2d4M+F+x/VZHtomIxLytW7fy8ssv89FHHzFz5kwKCws57rjjOPXUU3n66af51a9+hZlVSC2FhTB7NixaBEuWwKpVsGGDX7Ztg127IC/PL7t2wdtvw3nnlU8tBw0TM7scyHbOzTSzbmX4zJJa8SfTOZpZX6AvQMuWLcvwsSIiFaewsJBly5Yxe/ZsZsyYwaxZs0hPTycjI4Nt27Zx5plnMnjwYH7961/Tvn37Qw6Q7dshKwuys2HLFti92//yL77esQM2btwbFBs2wPr1PkSys/d+VsOGUK8epKRA3bqQlASJiXvXKSmHqXFKUJaeSRfgCjPrASQDtczsXSDLzJpEeiVNgKIfaTXQotj7mwNr9v9Q59wIYARAp06dNHewiFS4nTt3Mnr0aBYsWMCKFStYt24dOTk5ZGVlkZWVRX5+PgAJCQm0b9+ek08+me7du9O7d2/OPPPMMn1HYSGsWwerV0NGhl+ysuCbb2DaNNi6tez1Vq3qQyIlxS+XXALdu8Mpp0Dr1lCt2qG0wuFhP2cO+EjP5I+Rs7meBNYXG4BPcc792cxOBN7Hj5M0BSYBbQ40AN+pUyene3OJSHnLzMxk/PjxTJs2jfnz5/Pjjz+yfv164uPjadGiBQ0aNKB+/fo0btyYRo0acfTRR3PyySfTrl07kpOTS/zMHTsgPX3vsmIFrFzpw6MoQPLyfvq+k06Cs8+GFi18j6JhQ6hZ0/cgEhMhIWHvOjnZh0jVquXZOgdnZjOdc51Kei2a60yGAh+a2S3ASuAaAOfcfDP7EFgA5AP9dCaXiATBOUdaWhpjx45l7NixTJ06FYB69erRrl07Lr/8cm644Qa6detGQkJCqZ+TmQkLF8Ly5X5Ztmzv47Vr9903MRGaN/chcdZZ/nHz5tCsmV83bep7FUEHw+H2s3om5UU9ExE53D777DPuvPNOFi1aBECnTp3o2bMnPXv2pF27dgcc4ygshPnzfYB8/DF89BEU/aqMi/NB0bq1X446au/j1FRo3NjvUxmVV89ERCRUnHN8/fXXfPDBB7zyyiu0bduWl156iSuuuIJmzfY9qXT1an8GVFaW712sXQtpaTBzpt+2e7ffr3Zt+NOf4KKLfGC0aOEPPcm+FCYiUimsXbuWa6+9lm+++Ya4uDh+//vf89RTTxEXV5WpU2H06L3jGOnp8L//7e1tAFSpAq1a+VNnmzWD447z4xpt2/oxCzkwhYmIxLTc3Fz69OnDP//5T5KTk3n55Ze59tob+OKLmvziFzB5Muzc6fdNSto7fvHAA3DOOdCokT80Vbdu5T08VRE0ZiIioVZYWMjmzZtZt24dmZk5/Otf61i7dh3btm1k+/Z8vvtuJBs3LqZp03uoX783ycntWLHCH7Y66ii44gro1g3OOMOfMVVB1xNWShozEZFQy8zMZPr06cyaNYv58+eTnZ1NTk4O69atY/369T+5M25xZg1p2PBzjjvuPOLifFhccAFcdRX07AmV/C4moaEwEZHArF+/nltvvZWxY8finMPMOProo2natCnHH388NWrUZ/36BqxeXZ8ffqhPSkp97rmnPt26NaBOnbrUqFGFZs2SKv19r2KBwkREAjF37lx69uxJRkYG999/P5deeiknnXQS1atXJy8P3nsPBg/2Z1Y1agT33AMPPugv7JPwUZiISIXKzc3l9ddf5y9/+Qu1atXim2++4fTTTwegoMDfZuSxx+Czz6BDBxg3Dk47TWMdYacwEZEKUVhYyN/+9jfuv/9+cnJyOPfcc3n//fdp2rQpOTnw3HPwwQeweLEPjldegb59FSKxQmEiIuVu27ZtXHfddXzyySd07dqVN9/8hJUrT+f++2H6dFiwwO939tnw8MNw5pn+anKJHQoTESlXO3fu5LLLLuPbb7+lf//nyM3tT+/exqZNUL++P4T1y1/CddfBCScEXa0cKoWJiJQb5xy33XYbX3/9NR07vscLL1xPvXrQtSs8/jiceKIOY1UWChMRKRe7d+/m1ltv5e2336Zhw4dIS7ueYcPgjjt0e5LKSGEiIoddfn4+V199NePGjaNhw0fYvPl+xo/3EzlJ5aQwEZHDasOGAq699nYmThyH2fNAfz791F+VLpWXwkREopaRAWPHwujRW/nyy1449ykJCfdz8839eeyx8p17XMJBYSIih8Q5GD8eXngBPv+8EJhKUtLvgbkMGPAyjz9+G4mJQVcpFeWgN1w2s2Qzm2ZmP5jZfDN7OLI9xcy+MLPFkXXdYu8ZbGZLzGyRmV1cnj+AiFS8ggL4/e/hiiuW8t13D1KnztHAWSQmLmXChPE8/bSC5EhTlrv37wLOd851AE4GLjGzM4BBwCTnXBtgUuQ5ZtYW6AWcCFwCvGRmugubSCWxcWMuXbu+ziuvnAMcw9atj3LaaW14++23WblyJZdccknQJUoADnqYy/kJT7ZGniZEFgf0BLpFto8EvgIGRraPcs7tApab2RKgM/Dd4SxcRCrOjh07+OSTTxg9egyjR4+joGAH9esfx4ABj9G7d29atGgRdIkSsDKNmUR6FjOBY4AXnXNTzayRcy4TwDmXaWYNI7s3A/5X7O2rI9tEJMZs27aN559/nmeeeYbs7GySkxtSUNCH//u/m3j55c6YrjiUiDJNUumcK3DOnQw0BzqbWbsD7F7S366fTOdoZn3NbIaZzcjJySlbtSJSIZxzjBo1iuOOO47BgwfTsWNH+vefyM6da7j33pd55ZXTFSSyj58147FzbhP+cNYlQJaZNQGIrLMju60Givd5mwNrSvisEc65Ts65Tg0aNDiE0kWkPBQWFjJgwACuu+46GjVqxH//+1+uu+7fvPDCBVxxRTyPPhp0hRJGZTmbq4GZ1Yk8rgp0BxYC44A+kd36AGMjj8cBvcwsycxaA22AaYe7cBE5/PLy8rjxxht59tlnueOOO5g6dRpz5nTh5pvhwgv9LeLjftZ/QeVIUZYxkybAyMi4SRzwoXNuvJl9B3xoZrcAK4FrAJxz883sQ2ABkA/0c86VPoGziITC1q1bufrqq/nss8947LHHGDhwEIMHG088AZdeCh9/rHtqSenMn6wVrE6dOrkZM2YEXYbIEcs5x8UXX8ykSZN49dVXueGGW7nxRhg9Gm67DV58UT0SATOb6ZzrVNJr+ushInz88cd88cUXPPvss/z617dyww0wZgwMG6YgkbLR7VREjnAFBQXcd999nHTSSfTo8XtOPBFWr4bhw+Huu4OuTmKFwkTkCPfRRx+xePFi3nhjND16xLN1K3z3HZxxRtCVSSzRmInIEWzHjh20bduW6tVrUKfOD8ycGcfEidClS9CVSRhpzEREfsI5x1133UV6ejqnn/48kyfH8dprChI5NDrMJXIEWr58OX/4wx/49NNP+c1vBvL2293o1Qt69w66MjkkBQV+Upkff9y7ZGTAxo1+2bDBL6NHl9t0lwoTkSPIunXreOKJJ3j++eeJj4/njjue4aOP7qRJE3jppaCrkzJZswZmzICVK2HqVJg82T8uKHY5X/Xq0KIF1K0LjRvDiSf6GcqaNCm3shQmIkeALVu28OSTT/LMM8+wfft2brjhBurVG8Kzz7YgJQU++8z/3pGQycjwZ0MsWQLTpvklI2Pv6w0bwjnnwPXXQ8uW0KYNHHecD40KvneawkSkEsvPz+eNN97ggQceICsri2uuuYaHH36YjIwTuOgi/zvoxRehTp2gKxUAtm2Dr7+Gzz+HL76ABQv2vnbMMdCtG3TuDJ06QevWvtcRkhtuKkxEKiHnHGPGjOEvf/kLaWlpdO3alXHjxtG5c2dWrPBjIyecACNG+CMiUgE2bvThkJ4Omzf74Ni2DbZv9+sFC/whq927/X1rzj4bfvtbHyDHHBP6xFeYiFQiBQUFvPHGGzzzzDOkpaXRtm1bxowZw5VXXomZsWEDnHsu7NoFH354BAaJczB7NixfDvn5/hd3fn75P87Lgzlz/OP9JSf7P4iWLeGuu+Cii/wpdVWrVnz7REFhIlJJOOcYMGAAzz33HB07duTdd9+lV69exMfvnTX7kUdg1SqYMsWPyVZKGRmQne1/kRdf8vLgySf9YaSfKyEBqlTZuy7r48REqFbNP77gAt/LOOooPxherZpfKsm9ahQmIjFu586dvPfeewwfPpwFCxZw55138swzz/xk8qqpU/34yK23wumnB1RseVq82B8Wmjy59H2qVYPnnvOD1gkJZQuEYmEspVOYiMSw119/nXvvvZfs7Gw6dOjAyJEj6d27954gyc72vZBvv4WRI6FZM3/zxkqjsNBPsvL55/64XVKS730cc8zesCi+pKZCo0ZBV10pKUxEYtTUqVPp27cvZ511Fn//+98577zz9oTI2LHw6KMwc6bfNzHRH4ofOjT047hll5sLv/mN/2EbNPAX4730kk9MqXAKE5EYs27dOoYNG8YLL7xA06ZNGT9+PLVr197z+hNPwMCB0LYtPP64P6Jz6qn+P+0xqaDAj3msWwePPQb/+x/s3AmZmT5QnnkG7rij0ow9xCqFiUiMcM7x0ksvMXjwYLZu3Urv3r159NFH9wmSCRN8kFx7Lbz9tu+RhJZzPvnGjPGnzebl7V2KBszz8vyhrCKJiXD++VCzph/46d/fJ6UE7qBhYmYtgLeBxkAhMMI591czSwE+AFKBdODXzrmNkfcMBm4BCoA7nHOflUv1IkeIrKwsbr75ZiZMmMDFF1/M8OHDadu27T775OZC377Qrp0fHwltkHzwgT8TYO1aP2h+1ln+IryEBF900bL/88REP3/wCScE/RNICcrSM8kH7nHOfW9mNYGZZvYFcBMwyTk31MwGAYOAgWbWFugFnAg0BSaa2bGaB17k0EyYMIHf/va3bN68meeff55+/fr95Ewt8P/Jz8jw9/IL5SGtjAy45x4fJiec4FPvjjugX7/QXMUth+6gYeKcywQyI4+3mFka0AzoCXSL7DYS+AoYGNk+yjm3C1huZkuAzsB3h7t4kcps2bJl/OlPf2LMmDG0b9+eSZMm0a5dOwCysvxtmubNg7lz/bJgAVx3XQhP+83NhfffhyFD/OGsgQPh4YdDmnhyqH7WmImZpQIdgalAo0jQ4JzLNLOGkd2aAf8r9rbVkW0iUgZbt27l8ccf5+mnn6ZKlSoMGTKEAQMGkJyczDffwO2373vLplatoH176NkTBgwIrm4KCvz9pP77X3/LkKJlzRo/PnLyyfDJJ34tlU6Zw8TMagCjgbucc7kldbOLdi1h20+mczSzvkBfgJYtW5a1DJFKq7CwkPfee4+BAweSmZnJjTfeyOOPP06zyKmuixfDr34FtWv7a0W6dPEhUqtWwIWDvz/Ltdf603Tj4/3tz1NT4cIL/bpHDz8uosNZlVaZwsTMEvBB8p5zbkxkc5aZNYn0SpoA2ZHtq4EWxd7eHFiz/2c650YAI8BP23uI9YtUCunp6fzud79j4sSJdO7cmTFjxnBGsUnY09L8nOxVqsD48XD88QEWW2TtWvjhB1i40B/GmjbNXzDYv7+/35QcUcpyNpcBrwNpzrnhxV4aB/QBhkbWY4ttf9/MhuMH4NsA0w5n0SKVxdq1a3nqqad4/vnnqVKlCq+88gq/+93viCt2zYRz8Ic/+P/wz5zp/6MfiIUL4a23/I0SZ8/2AzdFGjaEjz6Cq68OqDgJWll6Jl2AG4G5ZjY7su1efIh8aGa3ACuBawCcc/PN7ENgAf5MsH46k0tkr4KCAr799lveeecd3n33XXbv3s1NN93EI488QvPmzffZd9Mmf7upSZP82bSBBElBge+BXHYZrF/v7xB56aXQoYNfTjzRX4GuQ1hHtLKczfVfSh4HAbiglPcMAYZEUZdIpfS3v/2NBx98kDVr1lCtWjV++9vfcs8999CmTZuf7Ltzp78EY/Fif5H37bdXcLFz5sCf/+xv7rVli79N+vff+1N6RfajK+BFKsisWbO47bbbOP300xk+fDiXX3451UuZUCQvD/7f//NjJePH+05BhXr3XX8FZO3acOONfsDm/PN13ysplcJEpJxlZ2czbNgwXnvtNRo0aMAnn3xCSkoKhYWwYgUsWgQ//uiXxYv9kp7ujy796lcVGCQbNsB//gPjxsE77/hZtD74QHfZlTJRmIiUI+ccvXv35quvvqJHjx4MHz6c9etTGDYMXn/dD0EUqVED2rTxt5rq1QtOO80PTZRTYX7ej//+1x+6mj/fd4Oc8xcT/vGP/qaKCQnlVIBUNgoTkcNo06ZN/Pjjj6SlpTFlyhS+/fZb0tLSeOGFF+jXrx9LlkDnzn4I4he/8GFx3HFw7LHQuHEFjmE/8QQMGuQfH3WUHwe55hp/n/rTTlOIyM+mMBE5DB544AFeffVVsrOz92yrVasWXbp04bbbbuO2225jxgx/2MrMn2V7zDEBFJqT40/hvf9+X8yrr0K9egEUIpWNwkTkEKSnp7Nw4UJWrFjB5MmTeeedd7jssss499xzOfbYY/csRfOvf/KJvwSjcWM/KWCFBMm6dfDddzBjBsya5ZfVq/1rF1zgj7MVu329SDQUJiI/08iRI7npppv2PK9Rowb9+vXjr3/9657wKG7LFvi///OTVX3xBdSvXw5FFRT4G3ZNmeIDZMoUP5IPftKo447zs2R17OgPY519tiaTksNKYSJSBnPmzGHixInMmTOH999/n27duvHoo4/SqlUrmjZtSnx8PIWF/vf3zJn+DK2iM7MWLICtW/0cUIclSHbs8LcLLroSffZsf03I1q3+9YYN4cwz4ZZb/PrUU/01IiLlSGEichAzZ86ka9eu7Ny5k9q1a3PTTTcxbNgw6tatC8D27fDGG/D00/6UXvDjIi1a+LOz+vTxwxPFbrVVNs75O+5++y2sWuUPW82bB5995nsi4Gcc7NDBf8npp/urHI86SlejS4VTmIgUk5GRwZAhQ8jKyiI3N5fc3FxmzZpF48aN+fbbb2nZsuWeialycuCll+CFF/zv+bPOgnvv9UeRjj/+EO51mJ/vuzUTJ/pl1izYvHnv64mJ0Lw53H2373F06ACtW+twlYSCwkQEmD17NrNmzWLo0KGsXLmSo48+mlq1alG3bl369+/PPffcs+dW8HPnwssvw5tv+lue/OIX/q4jXbuW8csKC+Hjj31YrF27d1m2zE8kBX7Oj+uv9zMSnnWWH/OoXl09DgkthYkc8T788EN69eqFc45mzZrx+eefc/bZZ++zz/LlPkDeesvfaT0x0d9l5J57SpmS3DmYMMEPmmRn77usXu2nsK1SxZ/e1bix73Gceaa/6vz88/2NE0ViiMJEjlg7d+6kd+/ejB49mjPPPJPXX399n9N5Fy2CZ5/1Z2AtXerf07atv+nijTce4PKMggI/gDJwoH8eH+8HxYuWc8+FSy6BG27QISqpNBQmckR67733eOihh1i6dCmPPPIIAwYM2Oemi/PmQffu/gSp88+HO++EC0/K4rj1U7B1OfD6Jj+esWnTT5dly/zxr1/+Ev72N6hTR6EhlZ7CRI4oGzZs4KGHHuL555/n1FNPZcyYMVx55ZV7Xl+wAH7/e3+ZRp06/nq/49sUwEMPQfehfpC8SJUq/qK/OnX8Uru2H9u45BI/OH7NNVC1asX/kCIBUJjIEWP+/PlcddVVLF26lL59+/Lcc8+RlJS05/Vx43yQ7N7tuOPabP583nQaTvgRRo/26fKb3/gpaZs29eFRrZoGxEUiFCZSaeTk5DB27Fi2bt3K9u3b9yw7duxg7dq1fPrpp9SuXZuvv/6aLl267PPe99/K44bfJtKmVR4fnzGUk959EN6NvNi6tb/1yM03V/wPJRIjyjIH/BvA5UC2c65dZFsK8AGQCqQDv3bObYy8Nhi4BSgA7nDOfVYulYtEFBYW8u9//5s77riDpUUj5UBcXBzVqlWjatWqVKtWk6uuuo1LL72fH35oyITxhaxK28rK1Ub6yjhW5FSnK9/ynxXnk7AiH+64w8+X27Il1K2rHojIQZhz7sA7mJ0DbAXeLhYmTwAbnHNDzWwQUNc5N9DM2gJ/BzoDTYGJwLEHmwO+U6dObsaMGdH/NHLEyM/PZ+jQoSxatIh58+Yxe/ZsGjRowIcffkiHDh2YM6caTz6ZyJw5RmbmvkMdAFUsnyZuDa1YQUtWcmqDVdw8uBF1GiVBkybQrZsCRGQ/ZjbTOdeppNfKMgf8N2aWut/mnkC3yOORwFfAwMj2Uc65XcByM1uCD5bvDqVwkZLs2rWL22+/nTfffJNWrVqRkpLCW2+9Ra9evUhKSmL2bOjRA2rVggsv9JdwNG8OzdbOpNlbQ2ieMZWGiZuIG/oYtGoF1RrAmVf4W5OIyCE51DGTRs65TADnXKaZNYxsbwb8r9h+qyPbRH62adOmMWrUKGbMmMGsWbMoKCjAOUd+fj75+fncf//9PProo/u8Z/lyuPJKSEmB6dOhcSMHo0bBgw/6U3aPPx763OxvgpiaGswPJlIJHe4B+JKOC5R4HM3M+gJ9AVq2bHmYy5BY9v333zNixAheffVVqlatSrt27bjppptITk7GzDAzLrjgAi666KK9byoo4N4/5/P48CSSkwr576gMGmdtgFvvg08/hU6d4E9/8tPRajIokcPuUMMky8yaRHolTYCi6eVWAy2K7dccWFPSBzjnRgAjwI+ZHGIdUols376dd999l9tvv53CwkL+8Ic/MGTIEGqWdPipoMDPUf7ZZ0xe0Zz+K/7IbDpyA+9y364hnPDLhX6/atX8Jet/+IO/El1EysWhhsk4oA8wNLIeW2z7+2Y2HD8A3waYFm2RUvllZ2fTpUsXlixZQteuXRkzZgwNIven2rYN8vL8/RG//RbmzS2kcMxY8mYXMr/uw3yeewaN62zjgdO+476eO0isMcgHR1wcdOnix0VEpFyV5dTgv+MH2+ub2WrgQXyIfGhmtwArgWsAnHPzzexDYAGQD/Q72JlccmTbvXs3Tz75JE8//TTbt2/n47fe4hdHH03i99/Dxo089dcEBk79JYWu+O1I4oCrgKtoUx+uuQKGDEmmWbMzgTOD+UFEjnAHPTW4IujU4CPTa6+9xt133822bdu4rGkLeqxrTvW8NmyiDlupwVd0YyIX8ouqE+lebQoWH0fr5EwurDWVKr2uJm7wQCxOp++KVJSoTg0WKQ+vvPIKt99+O12atKb3DsfzmZ/Sz7XdZ582LXfy/24p4E+DupOY2D2gSkWkLBQmUqF+/PFHnnvueV588QXgMiZnjmYySSQkON590888m5LiZymsVu3nTlUoIkFRmEiFcM5x772PMHToI0Ah0J9La/bk6kEbSGjRhHbtjI4dg65SRA6VwkTKTe4XU/n0ya+Zkj6X+RtX8uW6bzCu53rrwDlNdtBnekeSmuqaD5HKQGEih92uXbuYM/RNfvVQIqt4AVgFQJJdyVutW9LrkhXwwAPQSEEiUlkoTOSw2bFjByNefZU///GP5BX4M8Jr1kjhr899yamnduT442uRmKizr0QqI4WJlJlzjunTp5OWlrbnPllz585l/vxVLFq0mIyMBRQWFpBMN6rGX81jw1L57a1dqV27dtCli0g5U5hImWzcuJHf/OY3jB8/fp/tSUnJ5OUdhXOp+JtGn0arlM6891kjTu2kXojIkUJhIge1bNkyLrvsMpYuXcqvfvUUM2f2ZPfuBPLzISurAS2rb2fczotJqTKKGk8dT93+jYMuWUQqmMJEcA5ycvZOIOUcrFuXzbRpk/n++2mMHPkqu3ZBlSoTGT36HM46LY82efOpsnQRLUijT+FHpPbtBn/5i59YSkSOOAqTSmzzZvjyS/jf/2D7dn+jxPXr4fvvYcMGHxrOwc6dO9i+fSwwFz9jwErgY2AXflaBS+nc+RnOPvtYOtRO5/pXzyV+zSro3t1PbXvlTKhaNcCfVESCpjCpZNLTYeRI+Ne//ORQhYWQkADVq/ub6NasCR07+g7E+vU/sE/rFOcAABZtSURBVGTJ31m06B1gDWZxxMXFU716Ch06/IYuXW6hRYv2nHFGNTp2BPtkHFx3nb9EfcYMOOWUoH9cEQkJhUmMW77cz/00fTpMnQqLFvmpy884A+67Dy64AM46ywcK+ImnxowZw5dffsmUKVOoUqUK3bt3Z8CAtzjvvPOoUmW/vxLOweZN8PzbcPfdcOqpMG4cNNa4iIjspTAJubw8mD/fH5patcqPbaxbB5mZsHChfw7QtCl06AB9+8JVV+07I61zjkmT/sMbb7zBqFGjMDM6dOjAsGHDuOWWW6hXfObBjRv9cbHXX4cpU/yX7d7tX7v8cj8FbvXqFfbzi0hsUJiEwLx5flmzZu+SmenX6ek+UIqkpECDBtCwIfTsCSedBBdfDMceW/Jnb9myhZtvvpmPP/6YlJQU+vfvz8MPP0ydOnX27rRxI7z1FqSl+V5HVpYPjGuu8T2Q+vXhhBPgootg/56LiAgKk0Ds2AGrV8OHH/rf3dOKzUWZnAzNmvkxjZNPhiuv9EMTJ59cSHJyJqtXp5Oenk5GRgb5+fls3lzIBx84CgsLcc6vix7n5eUxfvx4fvzxRx5//HHuuusukqtU8aPxWVmwYgUMGwZffeVH5Bs29KHx9tt+YCUy06GIyMFocqwKkpkJL74I773nextFzjjDH5Y65ZTVrFw5hbS0mWzcuIHc3Nw9S3Z2NitWrGB30eGmAzAzDD8XYbwZR1etynPJyVywe7dPsf0/IyUFLr0U7rkH3bZXRA7kQJNjlVuYmNklwF+BeOBvzrmhpe1bGcPEOfjhBxg/HmbP9j2Q/Hzo0QNOO203+fmzKSycwpIlU5gyZQqrV68GIDExkXr16lGrVq09S/369UlNTSU1NZVWrVqRmppK8+bNSUpKwsyIi4vDnMNGjcIeeMCPyter54OienU/mFKnDlSr5k/hLVrXqOGDRD0QESmDCp9p0czigReBC4HVwHQzG+ecW1Ae31cutmzxh4MKCiA/H5dfQMGufAryCijIKyB/VwG7C+LI2pTEquwklqyK58cVm8jZtJlFywpYtDSfbTt2AztIqbWRE1tNpkrcdBbNymHSF+vYGRkIaVmvHl2POYazunblzFat6NCwIQmFhb4HUXwpGomfPds/z831I/Bm/nzftWth2TJ/bOzf//bjG6bbmYhIxSivMZPOwBLn3DIAMxuFv3FTuYTJ4MseZfv2HeQXFpBfkE9BYQG783dTWJiHuTwKdm5n5+at7Ny6nd35BezGsa0wgW2FVch3UEAceS6e7YXJbC9MpsDFUYjhAIdRiOEv3gN/UV/RegOwBsgEcvCTPv3UhlzYngud8MsVwBnAmUDz9ev9lYRTp5b8w8XF+fN6i5bERL+uWtWPb8TF+WA57jgYMgR+/Wu/TUSkApVXmDSjaBILbzVwevEdzKwv0BegZcuWUX3Z0AlPAbll2DMeSNizGFUw4v2hIhzx5oiPLyTewOLiwMAwzPxYhK977+NqSbVoULMuzeqewtENa9CmUTUa1ahO1cREkpOTqVqtGsnVqlG1enWObtGCxGrV9g2EkkJi/yU+Pqq2ERGpCOUVJiUdX9lncMY5NwIYAX7MJJovG/vWJ0AcCQkJJCYmkpCYQEJCAnEJyRRSlfhqNahZN5nq1eP2/H5OSdFZriIih0t5/TpdDbQo9rw5/nhQubiizznl9dEiIlIG5XVwfTrQxsxam1ki0AsYV07fJSIiASuXnolzLt/M+gOf4Qcq3nDOzS+P7xIRkeCV26iBc24CMKG8Pl9ERMIjFFfAm1kOsCLKj6kPrDsM5VQk1VxxYrFu1VxxYrHuIGpu5Zwr8SrnUITJ4WBmM0q7MjOsVHPFicW6VXPFicW6w1azrm4TEZGoKUxERCRqlSlMRgRdwCFQzRUnFutWzRUnFusOVc2VZsxERESCU5l6JiIiEpCYDxMzu8TMFpnZEjMbFHQ9B2Jm6WY218xmm9mMyLYUM/vCzBZH1nUDrvENM8s2s3nFtpVao5kNjrT9IjO7OEQ1P2RmGZG2nm1mPUJWcwsz+9LM0sxsvpndGdke2rY+QM1hb+tkM5tmZj9E6n44sj3MbV1azeFta+dczC74q+uXAkcBicAPQNug6zpAvelA/f22PQEMijweBAwLuMZzgFOAeQerEWgbafMkoHXkzyI+JDU/BPyxhH3DUnMT4JTI45rAj5HaQtvWB6g57G1tQI3I4wRgKn4WiDC3dWk1h7atY71nsmfeFOdcHlA0b0os6QmMjDweCVwZYC04577BT9RSXGk19gRGOed2OeeWA0vwfyYVqpSaSxOWmjOdc99HHm8B0vBTN4S2rQ9Qc2kCrxnAeVsjT4vmoHCEu61Lq7k0gdcc62FS0rwpB/rLHTQHfG5mMyPzuQA0cs5lgv/HCjQMrLrSlVZj2Nu/v5nNiRwGKzqEEbqazSwV6Ij/32dMtPV+NUPI29rM4s1sNpANfOGcC31bl1IzhLStYz1MDjpvSsh0cc6dAlwK9DOzWL93fpjb/2XgaOBk/FSYT0e2h6pmM6sBjAbucs4daIa30NRdQs2hb2vnXIFz7mT8dBidzazdAXYPRd2l1Bzato71MKnQeVOi5ZxbE1lnA//Ad0OzzKwJQGSdHVyFpSqtxtC2v3MuK/KPsRB4jb1d/tDUbGYJ+F/K7znnxkQ2h7qtS6o5Ftq6iHNuE/AVcAkhb+sixWsOc1vHepjEzLwpZlbdzGoWPQYuAubh6+0T2a0PMDaYCg+otBrHAb3MLMnMWgNtgGkB1PcTRb8kIn6Jb2sISc1mZsDrQJpzbnixl0Lb1qXVHANt3cDM6kQeVwW6AwsJd1uXWHOo27oiR/vLYwF64M8qWQrcF3Q9B6jzKPzZFj8A84tqBeoBk4DFkXVKwHX+Hd993o3/384tB6oRuC/S9ouAS0NU8zvAXGAO/h9ak5DV3BV/GGIOMDuy9AhzWx+g5rC39UnArEh984AHItvD3Nal1RzattYV8CIiErVYP8wlIiIhoDAREZGoKUxERCRqChMREYmawkRERKKmMBERkagpTEREJGoKExERiVqVoAsAqF+/vktNTQ26DBEROYCZM2euc841KOm1UIRJamoqM2bMCLoMERE5ADNbUdprOswlIiJRU5iIiEjUQnGYS+SQTJgAK0rtdYvI/nr0gFatyuWjFSYSm7Ky4LLLgq5CJLZ8+qnCRGQfkyf79YQJcMopwdYiEivq1Cm3j1aYSGyaPBmSkuD88/1aRAKlMJFgpafD229DYeHPe98//wmnnaYgEQkJhYkE64UX4OmnD+29/fsf3lpE5JApTCRYq1fDMcfA4sVBVyIiUdB1JhKsNWugWbOgqxCRKClMJFhr1kDTpkFXISJRUphIcJyDjAyFiUgloDCR4GzaBDt36jCXSCWgMJHgrFnj1+qZiMQ8hYkER2EiUmno1OBYtGED9OsH27YFXUl0MjL8Woe5RGKewiQW/etfMGoUtGsHCQlBVxOdX/wCWrYMugoRiZLCJBbNmAFVq8KsWVBFf4QiEjyNmcSi6dOhY0cFiYiEhsIk1uTn+x5Jp05BVyIisofCJNYsXgzbt8OppwZdiYjIHgqTWLN0qV8fe2ywdYiIFBN1mJhZvJnNMrPxkecpZvaFmS2OrOtGX6bssWyZXx91VLB1iIgUczh6JncCacWeDwImOefaAJMiz+VwWbYMqleHBg2CrkREZI+owsTMmgOXAX8rtrknMDLyeCRwZTTfIftZtsz3SsyCrkREZI9oeybPAn8Gis+52sg5lwkQWTcs6Y1m1tfMZpjZjJycnCjLOIIUhYmISIgccpiY2eVAtnNu5qG83zk3wjnXyTnXqYEO2ZSNcwoTEQmlaK566wJcYWY9gGSglpm9C2SZWRPnXKaZNQGyD0ehAqxfDzt2QKtWQVciIrKPQ+6ZOOcGO+eaO+dSgV7Af5xzvYFxQJ/Ibn2AsVFXKV52JJcbNQq2DhGR/ZTHdSZDgQvNbDFwYeS5HA5FYdKwxGEoEZHAHJabOznnvgK+ijxeD1xwOD5X9lN0ooLGmEQkZHQFfCxRz0REQkphEkuKeib16gVbh4jIfhQmsSQ72weJbj0vIiGjMIklOTk6xCUioaQwiSXZ2Rp8F5FQUpjEEvVMRCSkFCaxRD0TEQkphUmsyMjwt1M5+uigKxER+QmFSayYNMmvL9D1oCISPgqTWDFpEtSvDyedFHQlIiI/oTCJBc7BxIlw/vkQpz8yEQkf/WaKBQsXwpo10L170JWIiJRIYRILNF4iIiGnMIkFkyZB69aaYVFEQkthEgvS0uCUU4KuQkSkVAqTWLBmDTRrFnQVIiKlUpiE3ZYtfmnaNOhKRERKpTAJuzVr/Fo9ExEJMYVJ2GVk+LXCRERCTGESdkU9Ex3mEpEQU5iEXVHPRGEiIiGmMAm7NWugZk2/iIiElMIk7NasUa9EREJPYRJ2GzdCSkrQVYiIHJDCJOw2b4batYOuQkTkgBQmYZebqzARkdBTmISdeiYiEgMUJmG3eTPUqhV0FSIiB6QwCbO8PNi5Uz0TEQk9hUmY5eb6tcJEREJOYRJmmzf7tQ5ziUjIKUzCrChM1DMRkZBTmISZDnOJSIxQmISZeiYiEiMOOUzMrIWZfWlmaWY238zujGxPMbMvzGxxZF338JV7hNGYiYjEiGh6JvnAPc65E4AzgH5m1hYYBExyzrUBJkWey6FQz0REYsQhh4lzLtM5933k8RYgDWgG9ARGRnYbCVwZbZFHrKIxE/VMRCTkDsuYiZmlAh2BqUAj51wm+MABGh6O7zgibd4MSUl+EREJsajDxMxqAKOBu5xzuT/jfX3NbIaZzcjJyYm2jMppwwaoqyEnEQm/qMLEzBLwQfKec25MZHOWmTWJvN4EyC7pvc65Ec65Ts65Tg0aNIimjMpr1Spo0SLoKkREDiqas7kMeB1Ic84NL/bSOKBP5HEfYOyhl3eEW7ECWrYMugoRkYOKpmfSBbgRON/MZkeWHsBQ4EIzWwxcGHkuP5dzsHIltGoVdCUiIgdV5VDf6Jz7L2ClvHzBoX6uRKxfDzt2qGciIjFBV8CH1YoVfq0wEZEYoDAJq5Ur/VqHuUQkBihMwio93a/VMxGRGKAwCaOCAnjzTTjmGKhXL+hqREQO6pAH4OUwW7gQxo2Dd9/115ds2gSjRoGVdo6DiEh4KEzC4PPPoUcP3yPp0gVuuAE6dIBrrgm6MhGRMlGYBG3bNrj+emjb1vdMUlODrkhE5GdTmATtzTf9NSUKEhGJYRqAD9prr8Hpp8NZZwVdiYjIIVOYBGn7dpg3Dy6+OOhKRESiojAJ0pw5UFgIHTsGXYmISFQUJkGaNcuvFSYiEuM0AB+EggKYNAkeeMBPfqWr3EUkxqlnEoSHH/bjJOvWwQUX6MJEEYl56pkEYepUaNcO/vMf3S5FRCoF9UyCMG8enHIKNGgAcfojEJHYp99kFW3DBlizxvdMREQqCYVJRZs3z6/btw+2DhGRw0hhUtGmT/frE08Mtg4RkcNIA/DlzTk/0dX69bBlCzz9NJxxBjRvHnRlIiKHjcKkPGVnwznnwKJFe7fFx8PHH+t0YBGpVBQmh8P8+XD77b7nUVx2th9wf+EFf2FijRr+zsCtWwdSpohIeVGYRMs56N8f5s71vZDiWrWC3/0OLrssmNpERCqIwiRas2bBV1/Bs8/CnXcGXY2ISCB0Nle0li/3627dAi1DRCRICpNorV3r140bB1uHiEiAFCbRWrvW3xKlfv2gKxERCYzCJFpr10LDhv6UXxGRI5TCJFpr1+oQl4gc8RQm0VKYiIgoTKKmMBERUZhExTnIylKYiMgRT2ESjY0bYfduhYmIHPEUJtFIT/frZs0CLUNEJGgKk2jMnevXmuhKRI5w5RYmZnaJmS0ysyVmNqi8vidQc+ZAcjIcc0zQlYiIBKpcbvRoZvHAi8CFwGpgupmNc84tKI/vC8ycOX4ud12wKFKuFi70/9wkOmefDU2alM9nl9ddgzsDS5xzywDMbBTQEyiXMLm56yB25e0uj48+sB8cNGsLN7xS8d8tcoTYuRPGjYP8/KArKU9bgVWRZXu5fcvw4Y9w992dy+WzyytMmuFbpchq4PTiO5hZX6AvQMuWLaP6sjcnvwzkRvUZh2w5sPztYL5bRCqNatVq0KRJS6pXr1lu39G+ffklcnmFSUlz0rp9njg3AhgB0KlTJ1fC/mU25z+TKSyI6iMOTZV4qJtS8d8rcoSpVw+qVOLZl5KTk6lduzYWw9N5l9cfz2qgRbHnzYE15fRdtD+vXXl9tIiIlEF5nc01HWhjZq3NLBHoBYwrp+8SEZGAlUvPxDmXb2b9gc+AeOAN59z88vguEREJXrkdhXTOTQAmlNfni4hIeJhzAQxc71+EWQ6wIsqPqQ+sOwzlVCTVXHFisW7VXHFise4gam7lnGtQ0guhCJPDwcxmOOc6BV3Hz6GaK04s1q2aK04s1h22mnVvLhERiZrCREREolaZwmRE0AUcAtVccWKxbtVccWKx7lDVXGnGTEREJDiVqWciIiIBifkwiaV5U8ws3czmmtlsM5sR2ZZiZl+Y2eLIum7ANb5hZtlmNq/YtlJrNLPBkbZfZGYXh6jmh8wsI9LWs82sR8hqbmFmX5pZmpnNN7M7I9tD29YHqDnsbZ1sZtPM7IdI3Q9Htoe5rUurObxt7ZyL2QV/df1S4CggEfgBaBt0XQeoNx2ov9+2J4BBkceDgGEB13gOcAow72A1Am0jbZ4EtI78WcSHpOaHgD+WsG9Yam4CnBJ5XBP4MVJbaNv6ADWHva0NqBF5nABMBc4IeVuXVnNo2zrWeyZ75k1xzuUBRfOmxJKewMjI45HAlQHWgnPuG2DDfptLq7EnMMo5t8s5txxYgv8zqVCl1FyasNSc6Zz7PvJ4C5CGn7ohtG19gJpLE3jNAM7bGnmaEFkc4W7r0mouTeA1x3qYlDRvyoH+cgfNAZ+b2czIfC4AjZxzmeD/sQINA6uudKXVGPb2729mcyKHwYoOYYSuZjNLBTri//cZE229X80Q8rY2s3gzmw1kA18450Lf1qXUDCFt61gPk4POmxIyXZxzpwCXAv3M7JygC4pSmNv/ZeBo4GQgE3g6sj1UNZtZDWA0cJdz7kAzvIWm7hJqDn1bO+cKnHMn46fD6GxmB5q3IhR1l1JzaNs61sOkQudNiZZzbk1knQ38A98NzTKzJgCRdXZwFZaqtBpD2/7OuazIP8ZC4DX2dvlDU7OZJeB/Kb/nnBsT2Rzqti6p5lho6yLOuU3AV8AlhLytixSvOcxtHethEjPzpphZdTOrWfQYuAiYh6+3T2S3PsDYYCo8oNJqHAf0MrMkM2sNtAGmBVDfTxT9koj4Jb6tISQ1m5kBrwNpzrnhxV4KbVuXVnMMtHUDM6sTeVwV6A4sJNxtXWLNoW7rihztL48F6IE/q2QpcF/Q9RygzqPwZ1v8AMwvqhWoB0wCFkfWKQHX+Xd893k3/n87txyoRuC+SNsvAi4NUc3vAHOBOfh/aE1CVnNX/GGIOcDsyNIjzG19gJrD3tYnAbMi9c0DHohsD3Nbl1ZzaNtaV8CLiEjUYv0wl4iIhIDCREREoqYwERGRqClMREQkagoTERGJmsJERESipjAREZGoKUxERCRq/x8BcvJu0XlTBwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "df_2001_c = pd.read_csv(\"./outputs/02001scwd1\", delim_whitespace=True, header=0)\n", + "df_2002_c = pd.read_csv(\"./outputs/02002scwd1\", delim_whitespace=True, header=0)\n", + "df_2003_c = pd.read_csv(\"./outputs/02003scwd1\", delim_whitespace=True, header=0)\n", + "\n", + "fig = plt.figure(figsize=(6.5, 7.5)) # Adjust the figure size as needed\n", + "ax1 = fig.add_subplot(311)\n", + "var='ET'\n", + "ax1.plot(np.arange(len(df_2001_c[var])), df_2001_c[var],color='red', label='2001')\n", + "ax1.plot(np.arange(len(df_2002_c[var])), df_2002_c[var],color='blue', label='2002')\n", + "ax1.plot(np.arange(len(df_2003_c[var])), df_2003_c[var],color='black', label='2003')\n", + "ax2 = fig.add_subplot(312)\n", + "var='RUNOFF'\n", + "ax2.plot(np.arange(len(df_2001_c[var])), df_2001_c[var],color='red', label='2001')\n", + "ax2.plot(np.arange(len(df_2002_c[var])), df_2002_c[var],color='blue', label='2002')\n", + "ax2.plot(np.arange(len(df_2003_c[var])), df_2003_c[var],color='black', label='2003')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/dryland_maize/mesite b/examples/dryland_maize/mesite index f5e8b183..4e7419b7 100644 --- a/examples/dryland_maize/mesite +++ b/examples/dryland_maize/mesite @@ -1,4 +1,4 @@ -41.2 899.0 10.0 0 +41.2 899.0 10.0 0 0 2.1E+05 7.8E+05 370.0 1.8 0.3 0.005 31 0 0 3 50.0 50.0 0.0 0.0 0.0 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.0 diff --git a/examples/lake/vasite b/examples/lake/vasite index f6a2a36d..f8d9a910 100644 --- a/examples/lake/vasite +++ b/examples/lake/vasite @@ -1,4 +1,4 @@ -38.4 129.0 16.3 0 +38.4 129.0 16.3 0 0 2.1E+05 7.8E+05 275.0 1.8 0.3 0.002 35 0 0 1 50.0 50.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 diff --git a/examples/sample_run/st022852 b/examples/sample_run/st022852 index 54abd713..7a80329d 100755 --- a/examples/sample_run/st022852 +++ b/examples/sample_run/st022852 @@ -1,6 +1,6 @@ -69.1 130.0 -8.0 1.0 -210000.0 780000.0 282.9 1.8 0.3 0.005 -62 0 0 1 3 100.0 0.0 +69.1 130.0 -8.0 1.0 0 +210000.0 780000.0 326.0 1.8 0.3 0.005 +0 0 0 1 3 100.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0 diff --git a/f77src/ecosys_core/BLOCKDATA001.f b/f77src/ecosys_core/BLOCKDATA001.f deleted file mode 100644 index c7fe1cb3..00000000 --- a/f77src/ecosys_core/BLOCKDATA001.f +++ /dev/null @@ -1,4 +0,0 @@ - block data - include "blk17.h" - DATA ICOR/1,-1,0,0,1,1,2,3,3,4,4,5/ - END diff --git a/f77src/ecosys_core/CMakeLists.txt b/f77src/ecosys_core/CMakeLists.txt index da8a4688..23a98763 100644 --- a/f77src/ecosys_core/CMakeLists.txt +++ b/f77src/ecosys_core/CMakeLists.txt @@ -53,7 +53,7 @@ configure_file(filec.h filec.h COPYONLY) configure_file(files.h files.h COPYONLY) configure_file(parameters.h parameters.h COPYONLY) set(ECOSYS_SOURCES - BLOCKDATA001.f + bkinit.f day.f endrun.f erosion.f diff --git a/f77src/ecosys_core/bkinit.f b/f77src/ecosys_core/bkinit.f new file mode 100644 index 00000000..67e73d81 --- /dev/null +++ b/f77src/ecosys_core/bkinit.f @@ -0,0 +1,5 @@ + subroutine bkinit + include "blk17.h" + + ICOR=(/1,-1,0,0,1,1,2,3,3,4,4,5/) + END subroutine bkinit diff --git a/f77src/ecosys_core/blkc.h b/f77src/ecosys_core/blkc.h index b60a813e..40380275 100644 --- a/f77src/ecosys_core/blkc.h +++ b/f77src/ecosys_core/blkc.h @@ -12,4 +12,5 @@ 2,NPR,NPS,JOUT,IOUT,KOUT,IOLD,ILAST,IRUN,IBEGIN,ISTART,IEND 3,ISOIL(4,JZ,JY,JX),LYRX,LYRC,LSG(JZ,JY,JX),NP(JY,JX) 4,NP0(JY,JX),IFLGI(JP,JY,JX),IFLGC(JP,JY,JX),IETYP(JY,JX) + 5,iLAKE(JY,JX) diff --git a/f77src/ecosys_core/readi.f b/f77src/ecosys_core/readi.f index d70b25a3..b990f644 100644 --- a/f77src/ecosys_core/readi.f +++ b/f77src/ecosys_core/readi.f @@ -67,12 +67,13 @@ SUBROUTINE readi(NA,ND,NT,NE,NAX,NDX,NTX,NEX,NF,NFX,NTZ C :=2 means allowing freeze-thaw plus SOC accumulation to change elevation C :=3 means allowing freeze-thaw plus SOC accumulation, plus erosion to change elevation C :=-1 means no change in elevation. - READ(1,*)(datav(jj),jj=1,4) + READ(1,*)(datav(jj),jj=1,5) ALATG=datav(1) ALTIG=datav(2) ATCAG=datav(3) IDTBLG=int(datav(4)) -C READ(1,*)ALATG,ALTIG,ATCAG,IDTBLG + iLAKEG=int(datav(5)) +C READ(1,*)ALATG,ALTIG,ATCAG,IDTBLG,iLAKEG READ(1,*)OXYEG,Z2GEG,CO2EIG,CH4EG,Z2OEG,ZNH3EG READ(1,*)IETYPG,ISALTG,IERSNG,NCNG,DTBLIG,DTBLDIG,DTBLGG READ(1,*)RCHQNG,RCHQEG,RCHQSG,RCHQWG,RCHGNUG,RCHGEUG,RCHGSUG @@ -80,6 +81,7 @@ SUBROUTINE readi(NA,ND,NT,NE,NAX,NDX,NTX,NEX,NF,NFX,NTZ READ(1,*)(DHI(NX),NX=1,NHE) READ(1,*)(DVI(NY),NY=1,NVS) CLOSE(1) + DO 9895 NX=NHW,NHE DO 9890 NY=NVN,NVS ALAT(NY,NX)=ALATG @@ -93,7 +95,8 @@ SUBROUTINE readi(NA,ND,NT,NE,NAX,NDX,NTX,NEX,NF,NFX,NTZ Z2OE(NY,NX)=Z2OEG ZNH3E(NY,NX)=ZNH3EG IETYP(NY,NX)=IETYPG - NCN(NY,NX)=NCNG + NCN(NY,NX)=max(NCNG,1) + iLAKE(NY,NX)=iLAKEG DTBLI(NY,NX)=DTBLIG DTBLDI(NY,NX)=DTBLDIG DTBLG(NY,NX)=DTBLGG diff --git a/f77src/ecosys_core/reads.f b/f77src/ecosys_core/reads.f index 4602d0ce..51ce3eae 100644 --- a/f77src/ecosys_core/reads.f +++ b/f77src/ecosys_core/reads.f @@ -102,6 +102,7 @@ SUBROUTINE reads(NA,ND,NT,NE,NAX,NDX,NTX,NEX,NF,NFX,NTZ DCNOR(N)=DCNOR(N-1) 26 CONTINUE READ(4,*)NPX,NPY,JOUT,IOUT,KOUT,ICLM + C C INCREMENTS IN START AND END DATES FOR SUCCESSIVE SCENARIOS C FROM LOOPS FOR SCENES, SCENARIOS IN RUNSCRIPT SET IN MAIN.F diff --git a/f77src/ecosys_core/redist.f b/f77src/ecosys_core/redist.f index c2f5df43..d5016498 100644 --- a/f77src/ecosys_core/redist.f +++ b/f77src/ecosys_core/redist.f @@ -2663,6 +2663,9 @@ SUBROUTINE redist(I,J,NHW,NHE,NVN,NVS) TFLWX(N3,N2,N1)=TFLWX(N3,N2,N1)+FLWX(N,N3,N2,N1)-FLWXNU(N5,N4) TFLWH(N3,N2,N1)=TFLWH(N3,N2,N1)+FLWH(N,N3,N2,N1)-FLWHNU(N5,N4) THFLW(N3,N2,N1)=THFLW(N3,N2,N1)+HFLW(N,N3,N2,N1)-HFLWNU(N5,N4) +C if(n3==1)then +C write(*,*)'n3u',HFLW(N,N3,N2,N1),HFLWNU(N5,N4) +C endif ELSE TFLW(N3,N2,N1)=TFLW(N3,N2,N1)+FLW(N,N3,N2,N1)-FLW(N,N6,N5,N4) TFLWX(N3,N2,N1)=TFLWX(N3,N2,N1)+FLWX(N,N3,N2,N1) @@ -2671,7 +2674,13 @@ SUBROUTINE redist(I,J,NHW,NHE,NVN,NVS) 2-FLWH(N,N6,N5,N4) THFLW(N3,N2,N1)=THFLW(N3,N2,N1)+HFLW(N,N3,N2,N1) 2-HFLW(N,N6,N5,N4) +C if(N3==2)then +C write(*,*)'TKSFLW',HFLW(N,N3,N2,N1) +C 2,HFLW(N,N6,N5,N4),N3,N6 +C endif ENDIF +C write(*,*)TFLW(N3,N2,N1),FLW(N,N3,N2,N1),FLWNU(N5,N4) +C 2,FLW(N,N6,N5,N4) C IF(N1.EQ.1.AND.N3.EQ.1)THEN C WRITE(*,6632)'TFLW',I,J,N,N1,N2,N3,N4,N5,N6,NU(N2,N1) C 2,TFLW(N3,N2,N1),FLW(N,N3,N2,N1),FLW(N,N6,N5,N4),FLWNU(N5,N4) @@ -3456,8 +3465,11 @@ SUBROUTINE redist(I,J,NHW,NHE,NVN,NVS) IF(VHCP(NUM(NY,NX),NY,NX).GT.ZEROS(NY,NX))THEN TKS(NUM(NY,NX),NY,NX)=(ENGY+HFLWS) 2/VHCP(NUM(NY,NX),NY,NX) - - if(abs(TKS(NUM(NY,NX),NY,NX)/tksx-1.)>0.025)then + if(TKS(NUM(NY,NX),NY,NX)/=TKS(NUM(NY,NX),NY,NX))then + write(*,*)ENGY,HFLWS,VHCP(NUM(NY,NX),NY,NX) + pause + endif + if(abs(TKS(NUM(NY,NX),NY,NX)/tksx-1.)>0.25)then TKS(NUM(NY,NX),NY,NX)=TKSX endif ELSE @@ -4413,6 +4425,9 @@ SUBROUTINE redist(I,J,NHW,NHE,NVN,NVS) TVOLI=0.0 TVOLIH=0.0 TENGY=0.0 +C write(*,*)'NU=',NU(NY,NX),NL(NY,NX),NUM(NY,NX),NUI(NY,NX) +C write(*,*)TFLW(NU(NY,NX),NY,NX),iLAKE(NY,NX) +C if(NU(NY,NX)/=1)pause DO 125 L=NU(NY,NX),NL(NY,NX) C C WATER, ICE, HEAT, TEMPERATURE @@ -4421,8 +4436,14 @@ SUBROUTINE redist(I,J,NHW,NHE,NVN,NVS) VHCPX=VHCP(L,NY,NX) VOLWXX=VOLW(L,NY,NX) VOLIXX=VOLI(L,NY,NX) + VOLW0=VOLW(L,NY,NX) + if(iLAKE(NY,NX).EQ.1)then + VOLW(L,NY,NX)=VOLW(L,NY,NX)+TFLW(L,NY,NX)+FINH(L,NY,NX) + 2+TTHAW(L,NY,NX)+FLU(L,NY,NX) + else VOLW(L,NY,NX)=VOLW(L,NY,NX)+TFLW(L,NY,NX)+FINH(L,NY,NX) 2+TTHAW(L,NY,NX)+TUPWTR(L,NY,NX)+FLU(L,NY,NX) + endif C if(VOLW(L,NY,NX)<0. .and. L==NU(NY,NX))then C write(*,*)VOLWXX,VOLW(L,NY,NX),TFLW(L,NY,NX),FINH(L,NY,NX) C 2,TTHAW(L,NY,NX),TUPWTR(L,NY,NX),FLU(L,NY,NX) @@ -4478,14 +4499,26 @@ SUBROUTINE redist(I,J,NHW,NHE,NVN,NVS) IF(VHCP(L,NY,NX).GT.ZEROS(NY,NX))THEN TKS(L,NY,NX)=(ENGY+THFLW(L,NY,NX)+THTHAW(L,NY,NX) 2+TUPHT(L,NY,NX)+HWFLU(L,NY,NX))/VHCP(L,NY,NX) - if(L==1.and.abs(TKS(L,NY,NX)/TKS10-1.)>0.025)then + + IF(TKS(L,NY,NX)/=TKS(L,NY,NX) + 2.or.TKS(L,NY,NX)>400.)then + write(*,*)ENGY,THFLW(L,NY,NX),THTHAW(L,NY,NX) + 2,TUPHT(L,NY,NX),HWFLU(L,NY,NX) + endif + if(L==1.and.abs(TKS(L,NY,NX)/TKS10-1.)>0.25)then TKS(L,NY,NX)=TKS10 endif ELSE TKS(L,NY,NX)=TKS(NUM(NY,NX),NY,NX) ENDIF - +C if(L==2)then +C write(*,*)'watawTKS',TKS(L-1,NY,NX),TKS(L,NY,NX) +C 2,THFLW(L,NY,NX),THTHAW(L,NY,NX) +C 3,TUPHT(L,NY,NX),HWFLU(L,NY,NX) +C endif TCS(L,NY,NX)=TKS(L,NY,NX)-273.15 +C if(L.eq.1)write(*,*)TCS(L,NY,NX) + if(TCS(L,NY,NX)/=TCS(L,NY,NX))PAUSE TSMX(L,NY,NX)=AMAX1(TSMX(L,NY,NX),TCS(L,NY,NX)) TSMN(L,NY,NX)=AMIN1(TSMN(L,NY,NX),TCS(L,NY,NX)) UN2GS(NY,NX)=UN2GS(NY,NX)+XN2GS(L,NY,NX) @@ -5669,7 +5702,7 @@ SUBROUTINE redist(I,J,NHW,NHE,NVN,NVS) C C SOIL SUBSIDENCE C - IF(IERSNG.GE.0)THEN + IF(IERSNG.GE.0.and.iLAKE(NY,NX).eq.0)THEN IF(BKDS(NU(NY,NX),NY,NX).LE.ZERO)THEN ICHKLX=0 ELSE diff --git a/f77src/ecosys_core/watsub.f b/f77src/ecosys_core/watsub.f index f59ecdbf..89f07d54 100644 --- a/f77src/ecosys_core/watsub.f +++ b/f77src/ecosys_core/watsub.f @@ -262,6 +262,7 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) FLU1(L,NY,NX)=0.0 HWFLU1(L,NY,NX)=0.0 ENDIF + IF(CDPTH(L,NY,NX).GE.DTBLX(NY,NX))THEN AREAU(L,NY,NX)=AMIN1(1.0,AMAX1(0.0 2,(CDPTH(L,NY,NX)-DTBLX(NY,NX)) @@ -2041,6 +2042,7 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) ELSE VFLXG=EVAPG(NY,NX)*4.19*TKQ(NY,NX) ENDIF + if(iLAKE(NY,NX).EQ.1)VFLXG=0.0 VOLW2(NUM(NY,NX),NY,NX)=VOLW2(NUM(NY,NX),NY,NX)+EVAPG(NY,NX) C C SOLVE FOR SOIL SURFACE TEMPERATURE AT WHICH ENERGY @@ -2287,6 +2289,10 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) FLV2=0.0 HWFLV2=0.0 ENDIF + if(iLAKE(NY,NX).eq.1)then + FLV2=0.0 + HWFLV2=0.0 + endif TKXR=TKR1-HWFLV2/VHCPR2 TK1X=TKS1+HWFLV2/VHCP12 TKY=(TKXR*VHCPR2+TK1X*VHCP12)/(VHCPR2+VHCP12) @@ -2415,6 +2421,7 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) FLWLXG=FLQM+EVAPG(NY,NX)+FLV1 FLWHLG=FLHM HFLWLG=HWFLQM+HFLXG+HWFLV1+HFLCR1 +C write(*,*)'watx',HWFLQM,HFLXG,HWFLV1,HFLCR1 FLWRLG=FLYM+EVAPR(NY,NX)-FLV1 HFLWRLG=HWFLYM+HFLXR-HWFLV1-HFLCR1 FLWVLS=(VOLW1(NUM(NY,NX),NY,NX)-VOLWX1(NUM(NY,NX),NY,NX))*XNPH @@ -2487,6 +2494,9 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) FLWLX(3,NUM(NY,NX),NY,NX)=FLWLXW+FLWLXG FLWHL(3,NUM(NY,NX),NY,NX)=FLWHLW+FLWHLG HFLWL(3,NUM(NY,NX),NY,NX)=HFLWLW+HFLWLG + if(HFLWL(3,NUM(NY,NX),NY,NX)/=HFLWL(3,NUM(NY,NX),NY,NX))then + write(*,*)'watsub',HFLWLW,HFLWLG + endif FLWRL(NY,NX)=FLWRLW+FLWRLG HFLWRL(NY,NX)=HFLWRLW+HFLWRLG C IF(I.GT.350.AND.NX.EQ.1)THEN @@ -2603,7 +2613,9 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) FLQ2=FLQ2+AMIN1(0.0,AMAX1(-VOLW1(NUM(NY,NX),NY,NX)*XNPX 2,VOLP1Z(NUM(NY,NX),NY,NX))) ENDIF - + if(iLAKE(NY,NX).eq.1)then + FLQR=0.0 + endif IF(FLQR.GT.0.0)THEN C from layer 0 to layer 1 HFLQR=4.19*TK1(0,NY,NX)*FLQR @@ -2641,7 +2653,11 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) 4322 FORMAT(A8,8I4,40E12.4) C ENDIF ELSE + IF(iLAKE(NY,NX).EQ.1)then + FLQR=0.0 + else FLQR=XVOLW(NY,NX)*XNPX + endif HFLQR=4.19*TK1(0,NY,NX)*FLQR FLWL(3,NUM(NY,NX),NY,NX)=FLWL(3,NUM(NY,NX),NY,NX)+FLQR HFLWL(3,NUM(NY,NX),NY,NX)=HFLWL(3,NUM(NY,NX),NY,NX)+HFLQR @@ -2676,8 +2692,9 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) C FLWRL,HFLWRL=total litter water,heat flux C IF(VOLPH1(NUM(NY,NX),NY,NX).GT.0.0 - 2.AND.XVOLW(NY,NX).GT.0.0)THEN + 2.AND.XVOLW(NY,NX).GT.0.0)THEN FLQHR=AMIN1(XVOLW(NY,NX)*XNPX,VOLPH1(NUM(NY,NX),NY,NX)) + if(iLAKE(NY,NX).EQ.1)FLQHR=0.0 HFLQHR=FLQHR*4.19*TK1(0,NY,NX) FLWHL(3,NUM(NY,NX),NY,NX)=FLWHL(3,NUM(NY,NX),NY,NX)+FLQHR HFLWL(3,NUM(NY,NX),NY,NX)=HFLWL(3,NUM(NY,NX),NY,NX)+HFLQHR @@ -2691,6 +2708,7 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) 4357 FORMAT(A8,6I4,40E12.4) C ENDIF ENDIF + C C FREEZE-THAW IN RESIDUE SURFACE FROM NET CHANGE IN RESIDUE C SURFACE HEAT STORAGE @@ -2982,12 +3000,14 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) C THAWR(NY,NX)=THAWR(NY,NX)+WFLXR(NY,NX) HTHAWR(NY,NX)=HTHAWR(NY,NX)+TFLXR(NY,NX) + IF(iLAKE(NY,NX).EQ.0)THEN FLW(3,NUM(NY,NX),NY,NX)=FLW(3,NUM(NY,NX),NY,NX) 2+FLWL(3,NUM(NY,NX),NY,NX) FLWX(3,NUM(NY,NX),NY,NX)=FLWX(3,NUM(NY,NX),NY,NX) 2+FLWLX(3,NUM(NY,NX),NY,NX) FLWH(3,NUM(NY,NX),NY,NX)=FLWH(3,NUM(NY,NX),NY,NX) 2+FLWHL(3,NUM(NY,NX),NY,NX) + ENDIF HFLW(3,NUM(NY,NX),NY,NX)=HFLW(3,NUM(NY,NX),NY,NX) 2+HFLWL(3,NUM(NY,NX),NY,NX) FLWR(NY,NX)=FLWR(NY,NX)+FLWRL(NY,NX) @@ -3332,7 +3352,8 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) C AVCNDL=source-destination hydraulic conductance C DLYR=layer thickness C - IF(CND1.GT.ZERO.AND.CNDL.GT.ZERO)THEN + IF(CND1.GT.ZERO.AND.CNDL.GT.ZERO.and. + 2iLAKE(N2,N1).EQ.0)THEN AVCNDL=2.0*CND1*CNDL/(CND1*DLYR(N,N6,N5,N4) 2+CNDL*DLYR(N,N3,N2,N1)) ELSE @@ -3383,17 +3404,20 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) C FLQL3=FLQX+AMAX1(FLQL1,AMIN1(0.0,FLQL2))*XNPX C FLQL4=AMIN1(0.0,AMAX1(FLQL3,-VOLP1(N3,N2,N1)*XNPX)) ENDIF +C write(*,*)'FLQL0',FLQL,VOLP1Z(N6,N5,N4),N,VOLW2(N6,N5,N4) IF(N.EQ.3.AND.VOLP1Z(N6,N5,N4).LT.0.0)THEN FLQL=FLQL+AMIN1(0.0,AMAX1(-VOLW2(N6,N5,N4)*XNPX 2,VOLP1Z(N6,N5,N4))) FLQ2=FLQ2+AMIN1(0.0,AMAX1(-VOLW2(N6,N5,N4)*XNPX 2,VOLP1Z(N6,N5,N4))) ENDIF + if(iLAKE(N2,N1).EQ.1)FLQL=0.0 IF(FLQL.GT.0.0)THEN HWFLQL=4.19*TK1(N3,N2,N1)*FLQL ELSE HWFLQL=4.19*TK1(N6,N5,N4)*FLQL ENDIF +C WRITE(*,*)'FLQL',AVCNDL,FLQL,n1,n2,n3,n4,n5,n6 VOLW2(N3,N2,N1)=VOLW2(N3,N2,N1)-FLQL VOLW2(N6,N5,N4)=VOLW2(N6,N5,N4)+FLQL C @@ -3435,7 +3459,8 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) ENDIF IF(N.EQ.3)THEN FLWHL(N,N6,N5,N4)=FLWHL(N,N6,N5,N4)+AMIN1(0.0,VOLPH1Z(N6,N5,N4)) - ENDIF + if(iLAKE(N5,N4).eq.1)FLWHL(N,N6,N5,N4)=0.0 + ENDIF FLWHM(M,N,N6,N5,N4)=FLWHL(N,N6,N5,N4) C IF(N4.EQ.1)THEN C WRITE(*,5478)'FLWH',I,J,M,N1,N2,N3,IFLGH @@ -3511,6 +3536,12 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) FLVL=0.0 HWFLVL=0.0 ENDIF + if(iLAKE(N2,N1).EQ.1)then + FLVC=0.0 + FLVX=0.0 + HWFLVL=0.0 + ENDIF + C C FLWL=total water+vapor flux to destination C FLWLX=total unsaturated water+vapor flux to destination @@ -3628,6 +3659,13 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) ENDIF TKY=(VHCP1(N3,N2,N1)*TK1X+VHCP1(N6,N5,N4)*TKLX) 2/(VHCP1(N3,N2,N1)+VHCP1(N6,N5,N4)) + if(TKY>400.0 .OR. TK1X>400.0)THEN + WRITE(*,*)M,N1,N2,N3,N4,N5,N6,TK1X,TKLX,TKY + write(*,*)TKS(N3,N2,N1),TKS(N6,N5,N4) + WRITE(*,*)HWFLVL,HFLXG + write(*,*)'tk>400' + PAUSE + ENDIF HFLWX=(TK1X-TKY)*VHCP1(N3,N2,N1)*XNPX HFLWC=ATCNDL*(TK1X-TKLX)*AREA(N,N3,N2,N1)*XNPH IF(HFLWC.GE.0.0)THEN @@ -3636,6 +3674,11 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) HFLWSX=AMIN1(0.0,AMAX1(HFLWX,HFLWC)) ENDIF HFLWL(N,N6,N5,N4)=HWFLWL+HWFLHL+HFLWSX + if(N6==N6X(N2,N1))then +C write(*,*)'wataw',M,HWFLWL,HWFLHL,HFLWSX,tk1x,tky +C 2,TK1(N3,N2,N1),TK1(N6,N5,N4),TKS(N3,N2,N1),TKS(N6,N5,N4) + if(abs(HFLWSX)>1.e10)pause + endif C IF(N3.EQ.1)THEN C WRITE(*,8765)'HFLWL',I,J,M,N1,N2,N3,N4,N5,N6,N C 2,HFLWL(N,N3,N2,N1),HFLWL(N,N6,N5,N4) @@ -3659,10 +3702,14 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) C HFLW=total heat flux C FLWM=water flux used for solute flux calculations in trnsfr.f C + IF(iLAKE(NY,NX).EQ.0)THEN FLW(N,N6,N5,N4)=FLW(N,N6,N5,N4)+FLWL(N,N6,N5,N4) FLWX(N,N6,N5,N4)=FLWX(N,N6,N5,N4)+FLWLX(N,N6,N5,N4) FLWH(N,N6,N5,N4)=FLWH(N,N6,N5,N4)+FLWHL(N,N6,N5,N4) + ENDIF HFLW(N,N6,N5,N4)=HFLW(N,N6,N5,N4)+HFLWL(N,N6,N5,N4) +C if(N6==N6X(N2,N1).or.N6<3)write(*,*)'wwww',M,N6,HFLW(N,N6,N5,N4) +C 2,HFLWL(N,N6,N5,N4) FLWM(M,N,N6,N5,N4)=FLWL(N,N6,N5,N4) IF(N.EQ.3)THEN C IF(I.EQ.55)THEN @@ -3718,7 +3765,20 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) C 3,VOLX(N3,N2,N1),VOLX(N6,N5,N4) C ENDIF ENDIF +C if(N6<=2)then +C write(*,*)'in4320',M,N,N6,N5,N4,N1,N2,N3,HFLWL(N,N3,N2,N1) +C 2,HFLWL(N,N6,N5,N4) +C endif +C TFLWL(N3,N2,N1)=TFLWL(N3,N2,N1)+FLWL(N,N3,N2,N1) +C 2-FLWL(N,N6,N5,N4) +C TFLWLX(N3,N2,N1)=TFLWLX(N3,N2,N1)+FLWLX(N,N3,N2,N1) +C 2-FLWLX(N,N6,N5,N4) +C TFLWHL(N3,N2,N1)=TFLWHL(N3,N2,N1)+FLWHL(N,N3,N2,N1) +C 2-FLWHL(N,N6,N5,N4) +C THFLWL(N3,N2,N1)=THFLWL(N3,N2,N1)+HFLWL(N,N3,N2,N1) +C 2-HFLWL(N,N6,N5,N4) 4320 CONTINUE + 4400 CONTINUE 9890 CONTINUE 9895 CONTINUE @@ -3730,6 +3790,10 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) C DO 9595 NX=NHW,NHE DO 9590 NY=NVN,NVS + +C write(*,*)'in950',M,HFLWL(N,1,1,1) +C 2,HFLWL(N,2,1,1) + DO 9585 L=NUM(NY,NX),NL(NY,NX) VOLP2=VOLA1(L,NY,NX)-VOLW1(L,NY,NX)-VOLI1(L,NY,NX) VOLPX2=VOLP2 @@ -4129,6 +4193,10 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) FLWHL(N,M6,M5,M4)=AMIN1(VOLWH1(L,NY,NX)*XNPX 2,XN*0.0098*-ABS(SLOPE(N,N2,N1))*CNDH1(L,NY,NX)*AREA(3,N3,N2,N1)) 3*RCHGFU*RCHGFT*XNPH + if(iLAKE(n2,n1).eq.1)then + FLWL(N,M6,M5,M4)=0.0 + FLWHL(N,M6,M5,M4)=0.0 + endif HFLWL(N,M6,M5,M4)=4.19*TK1(N3,N2,N1) 2*(FLWL(N,M6,M5,M4)+FLWHL(N,M6,M5,M4)) C IF(I.EQ.336)THEN @@ -4169,6 +4237,7 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) IF(PSISWT.LT.0.0)PSISWT=PSISWT-PSISWD FLWT=PSISWT*HCND(N,1,N3,N2,N1)*AREA(N,N3,N2,N1) 2*(1.0-AREAU(N3,N2,N1))/(RCHGFU+1.0)*RCHGFT*XNPH + if(iLAKE(N2,N1).eq.1)FLWT=0.0 FLWL(N,M6,M5,M4)=XN*FLWT FLWLX(N,M6,M5,M4)=XN*FLWT HFLWL(N,M6,M5,M4)=4.19*TK1(N3,N2,N1)*XN*FLWT @@ -4256,6 +4325,7 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) IF(PSISWT.LT.0.0)PSISWT=PSISWT-PSISWD FLWT=PSISWT*HCND(N,1,N3,N2,N1)*AREA(N,N3,N2,N1) 2*(1.0-AREAUD(N3,N2,N1))/(RCHGFU+1.0)*RCHGFT*XNPH + if(iLAKE(N2,N1).eq.1)FLWT=0.0 FLWL(N,M6,M5,M4)=FLWL(N,M6,M5,M4)+XN*FLWT FLWLX(N,M6,M5,M4)=FLWLX(N,M6,M5,M4)+XN*FLWT HFLWL(N,M6,M5,M4)=HFLWL(N,M6,M5,M4)+4.19*TK1(N3,N2,N1)*XN*FLWT @@ -4353,6 +4423,10 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) FLWUL=FLWU FLWUX=FLWU ENDIF + if(iLAKE(N2,N1).EQ.1)THEN + FLWUL=0.0 + FLWUX=0.0 + ENDIF FLWL(N,M6,M5,M4)=FLWL(N,M6,M5,M4)+XN*FLWUL FLWLX(N,M6,M5,M4)=FLWLX(N,M6,M5,M4)+XN*FLWUX HFLWL(N,M6,M5,M4)=HFLWL(N,M6,M5,M4)+4.19*TK1(N3,N2,N1) @@ -4437,10 +4511,14 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) 2-TKSD(N2,N1))*TCNDG/(DPTHSK(N2,N1)-CDPTH(N3,N2,N1)) 3*AREA(N,N3,N2,N1)*XNPH ENDIF + IF(iLAKE(NY,NX).EQ.0)THEN FLW(N,M6,M5,M4)=FLW(N,M6,M5,M4)+FLWL(N,M6,M5,M4) FLWX(N,M6,M5,M4)=FLWX(N,M6,M5,M4)+FLWLX(N,M6,M5,M4) FLWH(N,M6,M5,M4)=FLWH(N,M6,M5,M4)+FLWHL(N,M6,M5,M4) + ENDIF HFLW(N,M6,M5,M4)=HFLW(N,M6,M5,M4)+HFLWL(N,M6,M5,M4) +C if(M6==N6X(N2,N1).or.M6<3)write(*,*)'mmmm',M6 +C 2,HFLW(N,M6,M5,M4),HFLWL(N,M6,M5,M4) FLWM(M,N,M6,M5,M4)=FLWL(N,M6,M5,M4) FLWHM(M,N,M6,M5,M4)=FLWHL(N,M6,M5,M4) ENDIF @@ -4452,7 +4530,15 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) FLWM(M,N,M6,M5,M4)=0.0 FLWHM(M,N,M6,M5,M4)=0.0 ENDIF +C if(N6<=2)then +C write(*,*)'n9575',M,M6,M6,M4,HFLWL(N,M6,M5,M4) +C endif 9575 CONTINUE +C if(N6<=2)then +C write(*,*)'in9575',M,N6,N5,N4,N1,N2,N3,HFLWL(N,N3,N2,N1) +C 2,HFLWL(N,N6,N5,N4) +C endif + C C NET WATER AND HEAT FLUXES IN RUNOFF AND SNOW DRIFT C @@ -4518,6 +4604,7 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) ENDIF 1200 CONTINUE 1201 CONTINUE +CCC IF(VOLX(N3,N2,N1).GT.ZEROS2(N2,N1))THEN TFLWL(N3,N2,N1)=TFLWL(N3,N2,N1)+FLWL(N,N3,N2,N1) 2-FLWL(N,N6,N5,N4) @@ -4527,6 +4614,11 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) 2-FLWHL(N,N6,N5,N4) THFLWL(N3,N2,N1)=THFLWL(N3,N2,N1)+HFLWL(N,N3,N2,N1) 2-HFLWL(N,N6,N5,N4) +C IF(N3<=2)THEN +C WRITE(*,*)'watawTHFLWL',M,N3,N,N1,N6,N5,N4 +C 2,THFLWL(N3,N2,N1),HFLWL(N,N3,N2,N1) +C 3,HFLWL(N,N6,N5,N4),HFLW(N,N3,N2,N1),HFLW(N,N6,N5,N4) +C ENDIF C IF(N3.EQ.1)THEN C WRITE(*,3378)'THFLW',I,J,M,N1,N2,N3,N4,N5,N6,N C 2,TFLWL(N3,N2,N1),FLWL(N,N3,N2,N1),FLWL(N,N6,N5,N4) @@ -4570,7 +4662,9 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) FINHL(N3,N2,N1)=AMIN1(0.0,AMAX1(FINHX,-VOLPH1X,-VOLW1X)) ENDIF FINHM(M,N3,N2,N1)=FINHL(N3,N2,N1) + IF(iLAKE(NY,NX).EQ.0)THEN FINH(N3,N2,N1)=FINH(N3,N2,N1)+FINHL(N3,N2,N1) + ENDIF C IF(NX.EQ.1.AND.NY.EQ.1)THEN C WRITE(*,3366)'FINHL',I,J,M,N4,N5,N6,IFLGH,FINHL(N3,N2,N1) C 3,FINHX,VOLWH1(N3,N2,N1),VOLPH1(N3,N2,N1),VOLP1(N3,N2,N1) @@ -4691,7 +4785,7 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) C C UPDATE STATE VARIABLES FROM FLUXES CALCULATED ABOVE C - IF(M.NE.NPH)THEN +C IF(M.NE.NPH)THEN DO 9795 NX=NHW,NHE DO 9790 NY=NVN,NVS C @@ -4896,11 +4990,25 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) C VHCP1A,VHCP1B=volumetric heat capacities of soil+micropore,macropore C TK1=soil temperature C +C FOR LAKE turnoff the update of water content + if(iLAKE(NY,NX).EQ.1)THEN + DO 9784 L=NUM(NY,NX),NL(NY,NX) + TFLWL(L,NY,NX)=0.0 + TFLWLX(L,NY,NX)=0.0 + FLU1(L,NY,NX)=0.0 + FINHL(L,NY,NX)=0.0 +9784 CONTINUE + ENDIF DO 9785 L=NUM(NY,NX),NL(NY,NX) IF(VOLT(L,NY,NX).GT.ZEROS2(NY,NX))THEN - VOLw10=VOLW1(L,NY,NX) + VOLw10=VOLW1(L,NY,NX) VOLW1(L,NY,NX)=VOLW1(L,NY,NX)+TFLWL(L,NY,NX) 2+FINHL(L,NY,NX)+TWFLXL(L,NY,NX)+FLU1(L,NY,NX) +C if(L<=1)then +C write(*,*)'volw',M,L,VOLW(L,NY,NX),VOLw10,VOLW1(L,NY,NX) +C 2,TFLWL(L,NY,NX),FINHL(L,NY,NX),TWFLXL(L,NY,NX) +C 3,FLU1(L,NY,NX) +C endif C if(VOLW1(L,NY,NX)/=VOLW1(L,NY,NX))then C print*,'volw',L,VOLw10 C print*,TFLWL(L,NY,NX), @@ -4982,12 +5090,15 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) C C END ARTIFICIAL SOIL WARMING C - IF(VHCP1(L,NY,NX).GT.ZEROS(NY,NX))THEN tk1l=TK1(L,NY,NX) + IF(VHCP1(L,NY,NX).GT.ZEROS(NY,NX))THEN TK1(L,NY,NX)=(ENGY1+THFLWL(L,NY,NX)+TTFLXL(L,NY,NX) 2+HWFLU1(L,NY,NX))/VHCP1(L,NY,NX) - if(abs(TK1(L,NY,NX)/tk1l-1.)>0.025)then - TK1(L,NY,NX)=tk1l + + if(TK1(L,NY,NX)>400.0)then + write(*,*)'watawTK1',M,L,TK1(L,NY,NX),tk1l + 2,ENGY1,THFLWL(L,NY,NX),TTFLXL(L,NY,NX) + 2,HWFLU1(L,NY,NX),VHCP1(L,NY,NX) endif ELSEIF(L.EQ.1)THEN TK1(L,NY,NX)=TKA(NY,NX) @@ -5072,9 +5183,11 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) 9971 CONTINUE 9790 CONTINUE 9795 CONTINUE - ELSE +C ENDIF + IF(M.EQ.NPH)THEN DO 9695 NX=NHW,NHE DO 9690 NY=NVN,NVS + IF(NUM(NY,NX).EQ.NU(NY,NX))THEN FLWNU(NY,NX)=FLW(3,N6X(NY,NX),NY,NX) FLWXNU(NY,NX)=FLWX(3,N6X(NY,NX),NY,NX) @@ -5086,6 +5199,11 @@ SUBROUTINE watsub(I,J,NHW,NHE,NVN,NVS) FLWHNU(NY,NX)=FLWHNX(NY,NX) HFLWNU(NY,NX)=HFLWNX(NY,NX) ENDIF + if(iLAKE(NY,NX).EQ.1)then + FLWNU(NY,NX)=0.0 + FLWXNU(NY,NX)=0.0 + FLWHNU(NY,NX)=0.0 + endif 9690 CONTINUE 9695 CONTINUE ENDIF diff --git a/f77src/ecosys_driver/main.f b/f77src/ecosys_driver/main.f index 3cf5c057..b177281e 100755 --- a/f77src/ecosys_driver/main.f +++ b/f77src/ecosys_driver/main.f @@ -33,6 +33,7 @@ PROGRAM main outdir='./outputs/' ENDIF call system('mkdir -p '//trim(outdir)) + call bkinit() C C READ INPUT FILES C diff --git a/f77src/ecosys_driver/soil.f b/f77src/ecosys_driver/soil.f index c12c33d3..a32a3c17 100644 --- a/f77src/ecosys_driver/soil.f +++ b/f77src/ecosys_driver/soil.f @@ -192,7 +192,7 @@ SUBROUTINE soil(NA,ND,NT,NE,NAX,NDX,NTX,NEX C IN 'WOUTS', 'WOUTP' AND 'WOUTQ' C IF((I/KOUT)*KOUT.EQ.I.OR.I.EQ.IFIN)THEN -C WRITE(*,333)'WOUTS' + WRITE(*,*)'WOUTS I KOUT IFIN',I,KOUT,IFIN CALL WOUTS(I,NHW,NHE,NVN,NVS) CALL WOUTP(I,NHW,NHE,NVN,NVS) CALL WOUTQ(I,NHW,NHE,NVN,NVS)