From 539faf45777d956116a15c105b2d65f3fe454a68 Mon Sep 17 00:00:00 2001
From: Jeremy Auclair <jeremy.auclair@cesbio.cnes.fr>
Date: Wed, 19 Jul 2023 17:51:48 +0200
Subject: [PATCH] Development of samir xarray

---
 DEV_inputs_test/land_cover.nc              |  Bin 0 -> 13158 bytes
 DEV_inputs_test/ndvi.nc                    |  Bin 0 -> 56198 bytes
 DEV_inputs_test/soil.nc                    |  Bin 0 -> 14251 bytes
 DEV_inputs_test/weather.nc                 |  Bin 0 -> 167106 bytes
 code/modspa_samir.py                       |  228 +++-
 dev_samir_xarray.ipynb                     | 1257 ++++++++++++++++++++
 input/calculate_ndvi.py                    |    3 +-
 input/lib_era5_land_pixel.py               |    3 +-
 parameters/csv_files/params_samir_test.csv |   24 +
 9 files changed, 1484 insertions(+), 31 deletions(-)
 create mode 100644 DEV_inputs_test/land_cover.nc
 create mode 100644 DEV_inputs_test/ndvi.nc
 create mode 100644 DEV_inputs_test/soil.nc
 create mode 100644 DEV_inputs_test/weather.nc
 create mode 100644 dev_samir_xarray.ipynb
 create mode 100644 parameters/csv_files/params_samir_test.csv

diff --git a/DEV_inputs_test/land_cover.nc b/DEV_inputs_test/land_cover.nc
new file mode 100644
index 0000000000000000000000000000000000000000..381acc9d9e01422e5725274e38e090a7357e6b3c
GIT binary patch
literal 13158
zcmeHOeQaCR6~8YpUrpUMl+cYXtzO+mM9dA2?Zol4s&yPYcAD71c1l>9J>Kl+<bmzy
zJ-?K+($r1@RHXeeDgpZ_2pAenjDKLFfF|)drm35jY2`zVbuiebiH*XfNlfe1opbMf
zj$h<X4a>wp??(2y@7{CIJ@<F+$GPV@x5r0f8#=ajG_|)AD_!8J;AbUXUhz=u=kKQC
zBk3mBkDA=HOl;v*R+E*qqlpN}(%FgMhd1FjAyG!i72Tv6vZ6p%wzHmDyH`ZkK<4dQ
zyD%d(H}PlT@Clw}{BLYv6&_yCNIJ6rI8}ec-GV?RGzo&x!FenycA>_XNqY$cWFHBv
zCxHu~LCZ$05wBhb>cgX~zcW52<OG*xjjQ#d<8p<ORU}=+1s@_ORWp|_1R}a_<h8P5
zs#Ven{UADoNCXvyzG<>3KDa&)Fo2j=P!@P2<Kaw(Tn-foH$pO!n207biBwX~M8e}y
z;_Ak7k+iS|$?!y!@*^*P7O6-PI5vW2AkxxoL@Q}|yH3E4$u)v{iex<u1cwsSip2v;
zaY1be;>AWn$cIM>32&?XjhU!>jZHM*V%mO}k=oSzpPC3wwdN)|3__<Cb}Tir2X?sg
zjW6y>MdJHi*WysokaTXy$;hOkS5(8)U~=Ic>4Oj2{E%u?HNCvY=PA=!xF^sT>htf;
ztB2I0UMZ<%(=(qh4D9Le>-Y8f<VR1Pf=1+)gFRuSKjicbr{_7{HOuP9IW;*w$?1!n
z{)N$Fzx(U8DDT+ooW9KI5~rs*eTmcOIQ=oFr#O9#QB)0f2(h`jr3vO7(zaHde4;3}
zwzcAVqPZ5HYi`C2G@t}?RtES4I%@?H{-7CrS_IBtTa~mw=gUFBLBK)4LBK)a(~JPR
z7tHe=Ug7zow{ywfBXy6~(0i=Bx7E$OM*`<J8UWoJ(S>e1defFeU+>5x2c%<`%f6ZU
z&QNrMURelvlGApm0%^j_RI()vioS{E|8MEpw=T4ynOe7_X5iOw1=m+Aim53@*-#5h
z+$Zk+`+e`mGMOkj%ZDmlNywjQ3*l`-wzW|aivl@LFG&Yxzj}Wxo*1qnMMZnJu^zn8
z7L)f1My;FNBM(0G%bIa^VIBHq)~)Melp7DYak3!jEp3<;ug2#KY!Z5R&^of2@zi|#
zmWyA%=nZsL*DpR%=#;|lcwfW?mrT5Z@f%fHDXG-E#@T{x5@l1JQw?ccDd(k#4$d|q
ze>{KW57apP9*goY35(}ghO|;N!OP^&9Yts^=dgKcBP;9u!8`xtwf_Z)rQ?%yVKmV@
zBbQ_u?aQsKz?v8}9zJv6jz|+-8ocd)A}}0lL>$0gpg}$&f&BT7Z$CerjgNpvfZb&v
zD&Olx0%<e+vGnv$caBr8b;yNxec-|y4WwnR?{vn`#4;d@nA~>O8n7s9fhtYH#kUt8
znNi!hKG$#ti{AOzVrDoy9s}*RVAVLQ-~aHmM81Og0F!Wf*9SYuU0CoXHV2;ghtEsi
z$AVGZ9^qPk@LOXfh`sG$h&<3$O%f9eU}j=a<^AT^6nP2@;smbuge-mbR|klIBS{c~
zEAq`7!ZTzH3$A0E(X019cLTW}3w(sE19InWV><aY7C^VS!j515&PlQr6#yCQfqd_S
z$L}FIESQ3uDv;l77`l%<i3KodL{@&g=O<(f8tqUMTqOyaee2Dik#C?zUm#=?kQbhQ
z^%b%O_f~_=p|e*Vd5ye;9H1B)@A>J0e;vKO=QLzL@rl}hoz4C>+uQ7Jv$@UQHe1{5
zY_qY=zBb$1>}s>A&7L+}+U#hvq0N3a+kGwn^`1|>(bjMPvG~aF^1+9*{?^-`=@ns1
zpi!2bQI&7L$)qR16s$a4e3R+aSR^y!zHT%l4F#lYq^azLbVM(!lHZ$jdq$(F(dBBl
zXC$1RnxNcre`qM+_GBjG(R3;?Vo7j&2K~XIet*#83k~)827*4He_+t>2~TC?sdOTH
z{fs+0nHhC^+(ED3>z<u$P&4TF4Pwn?8UlpS+o++c<(swKyxZe#z&;f6_9OdLGLfBe
z=hZo22hKike=soM_l0}|{h?66-$;j0zt7iLhd@n-)L1l<g&^gOJ8LNADh8xfc|tXE
zim4lJ&@!BcKxH<Xo^cn!dodUCs;n1e-O%Q=GU#FB$U(r$00*a}8d_dcpnkXC#!$^E
zMO7{+ITQ|q2!%p}Him*yg#J~fYHH;<w`ahMsR2h>H_Um;0<{d9l4=@&On(CmLtchq
zcsc=-A!<pFhBH~P5G_pB1SV7IY`l@-05fXSHJs9ps8t|fU;lV6`g(*%L^wR@`mgh@
zojsto4>@}P-;mdK4-moIrL%4hb~?GFR4T9;$+a5@dc=7Cz;0;msjw%~Jyj&z+1r|;
zUY@J%rz>~PRldM!e-=9Xv?u%e$I-qXJo}mmi=Dpx(#q`cD!@5qs%5apC4+dpWZt#y
zL1Ujtu4Gn)8@g(!<s2ME#iCZJ>R{xWanZd-<1(j0Bs*CD^=F%`nDR$9>Eg);c0P6f
z++KWtL`P;j6fb#D2+*-xwuccQI#nx%exsT*Rk+N8t5h2f1`C#lv%*Ib!QG~f=i$T$
z$H|1#J+GJxrRDZpt@gT6%uA!Ho>xsRCuI&-O|>MUA0@(NrEVbB*}bZf;AlkB$)Q$D
z8Z5dSbwgGbwd!gF?fZUq8X&ROV4X|D;epq2+`CzgSX2-6>#~s6N@~puT80X5h6r`d
z<rCqeN8tt7sQ*>ajS@Y*(DO&uL3EM8+mn3wg`YqWT%v&Wv>VG;0HXoIs{>}~?(VJ^
z_aZRe3}M^mR59v@^U-0MQ}sDRsmyCR`wX-%C>V6(gkc*i68Kvl%gA$9*EGUa^9Ay1
zHD_oQ6O2jG3<4XX8B=ajC5-1tG@c_ZRMor$m`inUTqV6ARh5!vNG1$hD(ZR#15FrO
z<(@{a?xzB{)8lm>rPp~Nu7m8&{3beO>wXI`p2%cjmtpsq&j7`H48G$;0hla`@rPkI
zEuv-h&ln_l^cTN>Vl<K($8zX|CYpyF%EV;-G{cg)vK-}LvdYh~$J+P97-VpggMfp8
zgMfp8gMfp8gMfp8gMfp8gMfp8gMfp;$BICGMC4W;5t-x~xr+E}1`NG<A*4gn(O5Jc
zO-7<JjX(@SNC1Nd5X)1ls!50yipWFEql~{vC+%3zgh~zn3&Q#}j0Hg|F!s56jRrn8
NjjTAxd)FVm<KJ=XTOj}d

literal 0
HcmV?d00001

diff --git a/DEV_inputs_test/ndvi.nc b/DEV_inputs_test/ndvi.nc
new file mode 100644
index 0000000000000000000000000000000000000000..28e2e5c50cfebff2b904645c8b5fc476bc0d1e50
GIT binary patch
literal 56198
zcmeFa2Y91NdLCHK96*o&0T3i|&N=6tL1Zw3nPASj+3aq1n4a#*d3I-K)652~R+cSU
zmMxvnC)u*jvMir;k}ZcjCppPU*4byv=b-xvDD+H^#y+iY*3bE#uBV%C{q<MXTlH1?
z>Z|(6C&_4hpl__Nt)~aI-0u_j=x;VdjO4|y&tLvnJ{c{wQGTjT+z4y$Ca+paK8t;A
z=mw6`HG+TM=xD<Lu48ZYZ1cRm`Ow9{c%y?nq8k_gxU>0YXb9u4k2as~H32n4l_Ywc
zgtjP1Qq*y>xm*7o@)ARcZpt=_KQMCRMuW2L#*G_$jq)MX!gdNS@eoL(DUcxQq~efv
zHGU@Lzit-cJEZzW!ouXD4)B6hJ}ZaTGbMzhwj&0jz<;!mtDmf_u301d`$wyrdrK$v
z;|4X(1tB5A1qAMfYaxE1Xh2?&O&1gr;zqKeav2T73oI{O#%O9hmMf?7xrK5hl#L;Z
z0Ef4u&SrvKXgt=4kA;QfV_5<G^RFepZjEm4ZmvFC40zRauQ<N@LAa8PLN<_$27#-}
zin#f%7xj`@|N0k~{z)TNFOD_Pj3xbSEX{k<e>Pr*ywLYa-C9*_A;}J^SU?v@36IS;
zFB7?QNQ4%CT>O`Yd?cA^w$7%u`T8vqi(*1h|8^*4Lue~YJN3QQrK5%%u<WoL4vzMB
z>IWyAD+}vK_4WOu&80nrrlc1UoNcS0Hd_Y@bF(6E=H11|c%!4M`H1?D50*|gmv$D8
z>TBE1X5inOoOq?-JN_59e!JNZ>uoozIH;}tlYZFOjE!2@>X)tO7LtAnv70jKL7m;$
z+pfmXzu)SLz37j!n>ku={lHBV?4jeQ_8YHdF(UMPuQd_<rA7v%$36`uFMnX5#9Af%
zfs=kV0Y%Ylq|b(b?la{99%ug_a1<F#cYnBY;8BW4<bTmopj)^4*?%$O)0ckA%U^1C
zxc5Ib@wGZq3)>mE`_O~|7#`RKk+mOLe;;J{=?Soz>F?cb4zUoKnOj4w)$6)`m=hrG
zzd!f{$i;o<4^DtMGyMOek!3r?pu*G)?bVu^KYuoClmTY04oJ`sGdZ5i<C{A>6H7a%
z_0JXJe{L#=fp-@Yf2DmPUM|PbU7Rzx8!#tU8-D+X<NxuKNsUKF^FlM)4_4}CngyPw
z?Xz>Z%G-QLM~ztl9^U^4zy4MtlF#DFps}m~H2i?g4uVi?btPr~1L3Xl*NWeof?C6&
z72e9H)%RQBt@y3(-U@H22QBTR72e9H_1wZ%{MP!^3U9@44S`m8Yg+HebqCi9Z{<U2
z20smN#m6fhcAt$;Z3aILZ}AR(CcKqTL+kmRd`NV3(8{nC@Eo5@;JE~zOW?T#o=f1l
z1fEOaxdfg|;JE~zOW?T#o=f2Wg#<n^JNz2i;!l#?D6XdMgBYMs_cl+C8-`SD2*q(%
zmu?;lk2m*L>O!;8Y}6S|@K0zox*Pux>UaIqKVTwO`J1_aVuUz)=|8<g7!lgh{QcMB
zrDUP8jc#k%B=GB&b-iUrZW&RY+TB_P2a;FI4)8SV|BMHcJx$rPeJH^OvZ-uE|Fyu+
zyo6=cGNo3=nr{ju0%fOZaUMiQY#}Pe;;~XJ7l|!oQ{@V>gPibAe|70(>9~G^5I!Tp
zrh`>5#_t_(?(d-<{JvQP;pZL(w2n|(7MQW+fAS_SQOgi_<Ed;~GHM~EacF_fT@7pa
zLM~D`+CQitoovDx1U853;XlsytbT;c;xlXa8iy4=tKOxz=vM1z^_~3#DCvoIV|C5u
zGwDrcz4<))?SBY~(8HQ6gz?7+{szJSKu|i@e4izFLhwU^zd-PxHSzjy|MMcw@A_{L
z{3`^1mEbQE{6&I4Pw-C@{F4NKvWd88&)stgJeR<82|Sm;a|t|`z!xZia3vcDt)g*6
zMNWUSKO!7L-@+W^4EHZw7)nTu^$WPlfR&y)#s4<?WgLSX-2RGT@ne+3Y8N>`Z5#J|
z8y(|-)-mbk*M0)Mg#*Zu=r1qa{w4HhZ~!P?_;6_D?KgiLvA}~;_y+|1%CG+R_s=W;
zrup>z=MO~!jrQ+-M()pahG#m$Gj0A%dq3OOpKj-8+W47X<p1kG7x|&o>3ODqcnLJ`
z8uYe(a>=T7z}mte?5;s;d-bi)Y_D3IoFD3K13a62;kONL{$S_!jWR^PsM`kr?RKTO
z+}Qm^%fcJ+WHkKrE>a`(rJr1j{XjS2TGt6$h;zXdcok`f3sa4a?9Tq)`oi84TxY|Z
zSUGu}+dtY}+S$CUuL`4^YisqRI<`zq)Yt1La8YZ2Png<UtzSId4IiH@;j3Y*pJaqL
z;g}~Eok+C<Q}YVh{x07CHkj605k6sB{e)?MO?a|V7w*-Uj)a?gYe!2*P!iw}Vq6A~
zUKkRd;madJh`PU57j7O4mxLRe>l?zePfB@kR6oWg*SEv1x<(O>@ijK#+S1C&e)A?9
z77i|`*=BosKkyUt%+qQBMLS${A=UWg_8}G;zV?8#_|zq~b~FgK3b0{(a<Q;s!@&1{
zu#KW+y1=(wTFui!+(bNbB?4bu`1CWr0rBZ)oYx=n3J%14?*99cz^CnG{~q);><-8U
z8$2U*LA(F(?PN_}gbLZdTP;8Cf7szJk6o+Z;<18V>sC_h7Sd;L9kr~9t;LO03PxAc
zwQdEqF732#{ouZi9lk}2UF+r*9)j4xXlpno!G<0S{>D{mKynidkvn*Ztd6L48y0T7
zH5`v%<Blr$L3?gEAwdyAF+mAIDM3)b8xGH<*ntM$a4j>aiae_cY6y}>MmqA04T{(q
zngBM@rh{v_4QD23A!sFNBWNe+AV}`jA}~caoSUGBAi2_peB_x}mr;N`2MLA<h6zRp
zMhU_ers3iQ69kh4$<_j;$#aHamf$$S9Kk%n0>OqM7W^ghTqal{SS2_?aFXB@!D)gs
z1ZN4>2+k3lC)n7Rfq#)aFA-cOxI%E1V4Wad;jmk80@@(BNpOqcHbEHF4Yx~hkKjJR
z1A>PHj|d(UJRx{W@QmO&!3%;n30@Lx4eTrOe4F4Mf_DjC6TC<8KEVeB9};|l;3I-B
z5`0YX3Bi{LzD)2Hg0B*Mjo|A9-yrxV!M6y$P4FFp?-G2E;QIu>MDWW5KOp!Of*%t6
zD#5Q2{5rvJ5d0>=j|hHD@LL4GP4GJezf15(2!4;?j}rVo!5<^|M+p8n!9Pmyj}iO{
zf<H;{j}!bUf`5YGPZRu;1pgGlpCS0O1pfzuf12R`Nbt`P{GSN^9KoL__-6_JIfB1H
z@Xr(cMS}k`!M{N8mk9nvg1=1gFA@9|f`6Id|3dIr3H}<v|CQjc6Z{Q=e}&**CHR{J
z{~E!+PVjFK{F?;-H-i5=!M{cDZxj4G1ph9<zen)52>yM7|A62>B>0aA{$ql_P4J%(
z{HFx}8Nq)}@Lv%89fJRo;O`RrKM4LSg8!P}|4H!Q5d60Ue~;k5Blv$2{PzU^1Hu1D
z@IMj!eS!_e6a1|&Y}#%#U$hhKAlON;i(ogw9)i6D`v~?Eq!1h+NF_K(aEKs{Af4ba
zK?Xr4K^DOgf^34L1jh(+2yzMX2=WOE2nq>`2#N_x2ucac2+9d62r3Dx2sS=E0TF7*
zvzDNapq`+Cppl@7pqZeBpp~GFpq-$Dpp&4BpqrqFpqHSJpr2rXV31&lV3=TpV3c5t
zV4PrrV3J^pV47ftVDsucLgVDQVao*h<jHe^V3A;nV3}ZrV3ptm!AXKs1g8nk5S%4g
zBREHJp5Ow(MS@ENmkF*ATqRg1xJGcD;0D1>f?EW)3GNWwCD=3z{pGeB9eA9z!T&$p
zw{|V_&%YkS={~t#-!lKuQOwhP^KZMGc*Slan#U~9-E#^2@FZ|we5?3hzLhfH)nxxo
z(n}i#6C6Y#r6NctkN-~q!>2KVaJ=Tw!nabUvdwfi&`{-f+MC~1o{;X}PO@wvd4z19
z{jO3KY`!s&2rX1h{y8o=)M4RQgevnRQ$n<i1FX%G{Ic|SMd&>oFw-31KmG5_M$}m6
z6U}{fwDzSmYE0n1W-UHq9Ouvpu12i6)Bcwo@hbWm9FS~&=k}}YUz<QTa7|LpL!G~~
z7@9?C95C6O41V)#Kfj2+i37~2@eS&`kJ9_-*Kq*Q_5nN|{`()IF<gDfqaWbo-~HoX
zMJqU<+H8lvIpF*T`cpUnDm?(O@LRG!k62jHXKipv2BEpX`*(f`y@K1)jlfcY(4YI+
z-~6kH)dRkx=53_ErMUcS=ofJeC>Rak&raF?<J~8+mO=<2|NXj8AL@|yf2Qr9Y4>N^
z{F(NCrmdf8=V#jZnf86AZJ%k^XWI0c_I#!-pJ~Tu+VGk7d#3GfuilgW_uE=ue90SM
zC$^@r8>0yQqsDA<<F~-L+gNsQT--p9@nA!L^&kH$S;O4lT9?7h3_~Ux{oPM2=SRm2
z=i5!oIsVKMkJUmcpN^Dg#gmD$&}kJKgw@Kp@Dkg?EykQ!n~3ETpS%`pqoGQ5yb*iB
z<Z@cY+HxToE9FzsRt{pV-QsYXEDo*N<<y&P4zt-}vs<*GY9*O3r7F|2;#i@a5NpK_
zqs1tmoBNzMc8l4L-xNylK|Qb~j_UQj^UakFvDWxG(VZ@%2}iHyQk7ZpY8^fzuR~0|
z(d4k&EM}M4W^%c#md`7L%Vaixei^J!%aBjUA{F?ke^y*MTG~6tH}mQX<MpGJC9wXB
zp_HK#eAr)!m1e~|aIA51x>{e@Ut56B*w;7rpd8P{Sb>v^Bltvpw|=y_3g5YlwK~f)
z5gJDe3(cbiTnU%UWq&5Z+R_g2*O!h@Huu)WTALBO&xx_Oe{{0Zhyrh$qS>vV9D!s^
zpCdz;u_?pQR0=vntd)B#RIWe^H3H!@hh&(`mnzB6YdDZuqfK#teD40+OW@h@kC)Nk
zc>U4KPmi2Fb+pttJ$pU};Cblz96%=IFLn-ShohCoLFYQmbPKyn2L~`2eR8xr44o6N
zKQJ4%=2Vyy8*^$q;=?fv1MU;Ssq)kL^s_T3`n+G7@Sp$Qzvt(G?EmlX+3}BSJOANV
zx$Q968}i)TYc5`Euuh%S_n<xg&_aB?d-5LdTSvimIQVL`D*TR>=6U`?eP?I$;CLSz
z`RN4!tc~NwZ7LP*NYv#2^xC$tbh7YY-J};EekA#s|L{8jd}!UM*)xmcAG%UlK-UlX
zAuOOqrQjFQY}Hpz>Z=PYpejAB4)dAnd`3rVhbvc&#2c{UgW_}ot2=yLaJu_RdTc^z
zh-3d~XH}S}@2}S3m)!`<H;+&1yF#o-;ZxX8Riw3c|3EdtwP{EAz5ooJ<NEF<47#oT
zqlKl5&Ep>=qu_TcQy`Mo82pq<L&O6&H;8uggJg?~2jye5Bn_WiJk^3e5mhfl!lHVr
z@;$wC0b*<DUu$sHqp^D7jRKBp$0`zNd+^B%{J{UM4p51#ttTB0|BPTZLBc<nV1D&|
zF5?#hSTNkCgrAw54vs!uo&S2Vtk?I~kCqNL;8!6%(*vI!P*Ca5ckyH@{M0Vjh%yEH
zdF-a51cOa%eUS3ke{c4;P}rYdTUOu~e#gquYBPp@?|Px>Lo55l(*oGw|DZ15iOtQ+
zKqcw*I~NY?-OzIV6z$gSMVwJzZ}0LUUJ+1lZ!Z;Z^AU>UuF!DUFwkgeM1ku@xXHnq
zfCKS=IvT)sf`;=tL)94xq@xF1vj=N~op5Z=R<B<!;kOL>@ba-7`{=`To&bALtig$W
z9Q*?-OlXvWh9Dt=VU5Qv{x3xw9uC3IZS#&#ih!9La8K>8V@rfYs$4l=#xMCC-qNWH
z`#3y~Mvpxx_B5#qKjM}|asguF-5C-$e!mtLy&b7EI0QR~*%>Pfz%nHZdMbMxe6XDu
zjPk{k75q{uRV<x1Yynb*)a}PXDB5T|>~r8BjjSyyB*aD+xCvl`pI|Hk+l~XQC{1TG
zT73on2;>~{7WX#gAViN>I(OJ;eT&|-bh(Hl>P)&+u8|^YY5@F*VNWCLpcN?_(s!W<
zgoWV$42q@jdvLv0lhvOWfbq%@^p<uvWDwQuQUFH}2-vJMFJ8{$csjE#mB$`xi-R%{
zgPFY>q&$d1Dx_$TIZ+q_a6KkaG}*M{4brV<Yal0t8QKu=R)7OKkK3W1JFE}EuGp-z
zE?&;ztaTP$x_~|07RM+j216NlHUSC?kh~!@2MR-jV`-LAG|{~5!RJ^k)?iN5;P6#<
z*D;67rk_7r8yx};t$p#bhO^dN_30w^Xrf>e!%nV@!A)x>Sf%tR!a(vyzp|qc5a525
z#~e}|Q7^Svt)X!-%-jaQZxT2lXP3n^f3%Ky&06Q;<!mFOO`j=YkE_LDvB%&t1{Zp=
zA4U6+tT8@pC<KDS<0Ffr@#duupUz^lg~z49Vf6c_fdk6rGT9c6H`LJm&06>3<qXc+
zU^iqd*prM}Fd!6!Jd#KPEt$iuFKYw>N4Vhu^B%nSV*j|-WVHq}QYekV=bPAG!yHbd
zwRX6!f*PB&&iR|uIHJy`PnWUB6($_G9x)^a4nOWp@*B_~eO45JdSK^9kr=6m)@ZT#
zQ?dq!x4hNha2PDp2Wv{OkQlZ0+Qk%3qP6K#4Gs<T(}t5PqbQ6-(ChVpNYM^bW?GSt
z3?B^&M@f5V3>LFDDT6U+@Oq0IbzBd--aN5iS2Q@Rv*%!Lgr>G?la2P^whB7RiNa+U
z%ppKQgjAhyEF!cb4{mU9JQR+Qex=cyP42k7!QsiRH#lrMQ+aPy4rvTp^YqySPNKGI
z5)BTj507MUlz|AyTL`fR@x>Z-w-A>Exp0GnqvDpKVBerdXEr%wiY7;PtuBXpShc3&
z?ux9zVVXK^aHy=Bc!L9Z?KmNJ5fHi%YKvP0X?v9*0V{Ii=Z0&wKJ;>(*6f#nEY)^f
zbaqAz<3q2|CYPsi-&HHL`OR^hLaLVsnsbBQ*vtY!Ev$$OZ4SKz+M@>sgHewS%EY-_
z)_bv3sn9wEFoVd925+&zhpMV1ir{3P2mL@MQN-t>IE_#xcH(IW^3gN^e%PS3VO9qR
z8A}?nQsIvmn(l^`bPmK9NW=;w2LvZnt89t*7)(1-zSx<K;c+hHi+#mLC+AA}CRZb|
zx;YZSuZM}81+r~?6~XV~LYbsb?`)wPI(;KPtS%frPox@wxrr;4=zLxl_yl7-W5~mT
zx^c&N)_7w)umz)PQzMO3NRk4B7!8AL8y9TH8vROeL)3cPrf4k!Ox7ryD;|b*aD>N`
zTg<~SlEzqUrQJk_x{a{dT6aS>j4`I9u_hs2WBJ25iJ*fHHM$P0kgN`n8fWYV)TZg$
zge^K6oiWVf40gk6${yh>m4mRZF=)d=je-iiblR{)FUDzx7&NZ1k(dpAq3OopG99es
z;9yM98e@=gqtOR&So;_mDbULYsY8tZPUz5d8jH`zzQF-1TgdANUg`jKOw#C`ltBs;
zIx{#rv{TbD;3%LMVl^)G4YWaxh6dZ7zHJ|dPSbEby?wp?1KsVgT2cB384T>}>**b0
z4tIAzSWho)q|vH9eLd7>@@C0#;sIC$yP!;Xz@eVTeAU<1SgLSX#{e0PsH3B^qq7TN
zaOmjj?C$HszRr%0p1$7pHV6Y>e`E21Xx$AD$rb0)g@1V_#A+<+9aux@>}+dr;;@dc
zhB|@0Ux;fEHF5$m5ciY}_4f7kQ->OI>F*n2GI}9;e_!`7Yq%RE-P_YgAED#&_xATT
zb+W;M9@Gy-21!vMW>*{d`>9xm>xUDg-nIc~Q7j}{$G{*tm!mToEUu6W*$j^`g=!@o
zlCT(5xlS<%J)cfv$P9QYMsyZ^tkEs`W3=XR!WisKSxia?BxQ|Z9gYRNgy9>Uk^Y8q
zCu!&K$>mLsKp;|^n83^x3rs;T2ZSi%vs@833wj}EjPFc1a3AD|I0`i&(ksW12oDMc
zEPWyYy$h0x6+Fc57J&-Ey`d5hAv}<5ZAx8p{XkN+O6v@9p{R1bI#ex4p;%faFFTXt
z^I=wzC^HLb2_%tfrPjuj>a$6Z7WXhK%u-r4y9bi$EJ`Ho)53jvIdbS(4dpIz!|o-Q
zIQ2HGGf|R4nT>8+X}hk6G3B(#77mx?awxRcGJig&g&cJbT{sSi#v@wf#KSNQtA@kI
z?}wzWpc$zKoRFmv#oS!HyrXFp%{JF_q_bKa(V`R*0Y`qTu802Zut;YQSCyDUW12l*
z&_QR>+O@$r9&f&guE)iKkPr%XYi)cAOd+~}1*s|a9(}(7#asdu1^G>)NV>-e2_P+)
zjl)G5L^Qc=*{wA_h}vP5P9CkP8ytqI^Cbi1t+8tYaXkO&BZeV2y9@LZ6dDxS`BYfm
z#IXKqC^njQL=PNd6oDzNio&TODd8~Lf<>$|07rUj-2g+@Zk1GyHnqTE(C8-4S4>bh
zja}`J<GF|*F|j>kFqZ{Tc!Xsa3_{zm0yd=Ww~lJWI^d9_uoCNup#-}*`)JHYtG|fX
zbF<r)*xE2aS?xAS;doo$;Lw!M>lO%8+f}|e^a3;(u?W3fs4Zwbf<A|U28HSm*paT^
zB39|Oz@b7RJ3>Vij0<u9gYub;7H?4jox|d>$F??&;IrE$<0pGY;4o^`#q)I=)KYC%
zdJ_-=tGG?+=R@H^J5(U2a2OA`fCCx&%^Ia!i{dUF3Ps@@2nA!!#Q<qc29r0hfRaKz
zfWrVBR+|(!3?NLsMq9jCvp_9XHia)GgX{)E7Nw6j)Q<cE5joUjc%ZMJ643VP`%P-4
zQ-`81b!X6ydUGfclQ-*uk7v9&1ytDNvc<Qy^teV=c@8*WSm?C+(#<s!gek3Ze+D<<
zSkSEW@`k#Qe=H(n*aZwKbA%lbF?1B8N@)jm&ZX!NDiJ*k9Fk^x=nOi8FROqGo1E4J
za6p+&7G>^a-vG;)7BsC}bpuESI0EB%Fe(B@m4{2~MSewC#<TGmRDnbqVDWSV2Bpeo
zKtY#`7UUpK8u??wX1~(wwR&Gh0hKa3tjVoy9dH=U>O61=p|$ju%2izll2lk^!Mqqa
ztbT*a&81P0-x`)Gt$gO7+F%UyDYaC+LS-=`ze_w496-_}aPXSr!=Tsd{AnfhPNUtL
z+TPJY`x;Ez!s&ql7HPG?UcFP-fFyw<T*Mq<pI+_a(g%@06qcH;eAbZ7;|+p_Kd6(d
z%_ih?33)+~asoKU8Y&Ia>ohulMm~Vr^bSjUXGaSbF@s4{Iy=w-hel_ey1S-<_E4Cm
zk)i~Kx7DZ9I5<PV5d;nkpGjAHJ^ntUhN_ckjXLCYi8!I&c6|Z`BAiA&s5+fS@6W={
zrc>)QXZCirFiRUv>dN_n7C5vT^X&Dy24(?;Ng6FlVZbYVI<1WZ6MrBWk($kX21DTY
zco};2fL^LMsFBAl9t{t78xklO;x^iYtJSFu{&5K%b*tT`@%>#b6R}NZ)#SyY1~~L8
z<HEfSEp!isQ5q}DVZt2q>vR^*zz_-qB2tT)Kg?nVJ#L;(HDHh!bqeHmOW9FgzcGQr
z!RD`lRcmx=qdzBRqCTa^m^<9lvXIbhR!-kMR0D@Wp<lkgsfFoBVUWfvxQh??b$Szr
zGK7MGsMKQO)7gDtk5j2pQjB7=R*GC687rn3G$l~P&u@$mr&^~r1qvcI8c_I*h2wqg
z2vS-s%G%|T3OG!1?b?Ga+%ihNG+9+ar$Ygq(ZK1ap-><uwVL>}F%<RK&1wb3B(iEm
z$mvxu6Q*HH62*Lw4qRLX^#~S3EU1UuTsS$<FcFkaRlC*TfO@RGu%*ZK&`T3lym1)t
zYmFw(03C(=AS)ApaI_!lVbUlDj6w_MaLYzw7KSa2qJBOVYeG(yPGbxeMGR;UhqZ8e
z1nn!b7&UWO$7&(;K&5f*#VrHul@tb9d_oEw93O1&ID^9|<O@m7CIOW_67sl>S_Kt2
zw7}t#j7Drs$2f`vcw9({)GEEk7%qwEs8eIJ7tc;K!-#J&Xy<RAXn;ejGOs_`GQpUU
z8|ASH5e#^hPh+%kXbj*8N=;_Lz=$O1^%_9OFbFL=;BZ3Ub&j|SC=%jo@XWx4>DU-4
zL*h1_*;zb4)6fyeqSr6nIfXJIwc5J=@|GFKx!k0PP6~j-;nf%&oM9#kdHqtOSx8}l
zfx@rXD~5nW4;&5wFW_c-$|xM>xD3ESRT*@~SVassRt-i^>E;EDH?~D*T)aNh!KhYg
z>>ICcS)qH#&C19W4>)4L;pQ?&P}u8}8ZDxJrpNCK>J17SaKIervh#&L-<WR#_O+a(
z2RM+{q&LUPvf*x+r97pJV;!?|*lIK`+&MN#y74$)d%SLg?xe7&!c`HY8#(+sqf5Z%
zAgeJTH(De^qXv^PU^J>&hEb(#tlwi7i;YHxG1?Oe3mWx628+QOFDcnnWU_h7=X-|H
zekz;`&)+^UNvTM0u+$&c?Jx#aHg&ip8y!Tdkip~>a|K8t4J(Zn8IvcGNdp#>hGQ7!
zbB9?zhg>F=4oWSwSWwbv4`j9*?TMm>&p;NZw{o@x9lqCTH_ct`LULp@nOClt9k9l#
z?V3<t#bY5()NFRhg%UKzj;Uc4;D|=p?10^@6B;QD`hYm#R4dr*KDKx?;a4{L6|y-^
z&SXI^fx3Hq)zi8~JksHInQFJTY;rELSWL@Ti%yt!v`%d(s}u9lK*DOVtE37vKuKy%
z7PU~;PoV@{7K6;#)z{ml3%T@a3Z;`m6Q{h|#`r)^kJ*zdm{ek9_xmSKmu#R`dcEe_
z<+@WL1?|DObTRMJz|7^=h0-PktmtWn)vi-Q5Xu<LR+!Mi7w}rm&Fvy6v1T1^Xr#%w
zO$`sSyG@=<)g<S1+Wg+B%T1?rq|5CwE?w_?pdMC>aqaG!OMxglr!G2f7Vr@@33oWu
z5+&?wQZl7QBbN6MQ2a)#NoIyUK#wlsGN_@i_YH}YUj0C0!m_zcp79y8l+kDPyJznn
zdBn5<mq)+$=-e-7B8%Cy{b<iEL8G9HCW@937UIQi7MDUG0|h$4<y%w&@yHn4C$$+x
zR@ej3#1WTH3GK%ajmG`jvBo~n>NI)s3uf^E)#7)~KfLq{`Wa5QcKemPAuwN>jHbg^
z&pbjT)7tf!iWQa(P1I_3OF0sxMlm|WqTukw3aOWEGw^I20e_SgacSgw*ngX8(Vz;;
z7IYAY$x~dl2)h_&zkTV+gAlKi<8Z4F-g+@AhcRO?pS^R(%Y~)Xrq4~;VBg~ogP}+?
zA_V(Iq_@u^V~^>rMi0fVA9YC8N^yV0rIPxfuSa^qVTIL#Im~vGr?P3~qfxWpw)*-j
zF%A;j-Kw)My`F?J81$B_4`1|+qJY+FDAgR0AQd#5{rq7*ienS2MZz3$hrAxx5-?qQ
zn?(cql|noYOL7HtI9Vu&Ig}cOAyQGX5kqNG7mf}bOaumL;rz`b7f4kG7VE1WBLkTw
zDp@F{g#>PkQfXojvY{SUNX{LiLp=;gD;?6|dLX+_z-iP2M8s>>gD>Ps%mJ`d5AtOK
zZ?&R=vnQdL9iOdgK!fIT`02T#bP&nL1YDB?W=W)x^Z7D*KONdb0U9=wG63x%K_U)C
zjN1b##UqSHdmy=r)oc$2b9h9gW>8QcgUwXBoSZ&1z#JJecwJoRfpj{<<hQf?&?t>T
z7vh}|;*ByGY)Ur;`W5Wq7}TC#=vO0%MekwZeuX%!!GT7<0#gjN*{}M*1Te(xgROQS
zrH{eo4s@cfzW$yOK9>qgNN;Z+Tfm`oA_{EShv=X(AS$J=kJ8oN1>*yh%f9Z84j3Q(
zaEj5<kH-fZ=<n`oj1MqjbT`Kbr23(|-rk;mI;*b(Yz6&23?3iq(b?12Jt7o<eWAOj
zr*A~S>FGe!?%wV}8q}<hIsoR6&bIbmDrKOnx3>%GLFw=B>+R`)|J~i)l>TlgerGoo
z^$m0lQ2JmuGz^v)sgBu;Zh&D#<?%v$^wXGwhM=D{gs3zcQy*}^pK;+*hX?s$P<@!(
zu^~FMud{2IH#$b44^w)2SZo%HNgwKA(1!-8gQFwVK@Pirkc*g{fiX4%)@dGRRAi3`
zKo=bqkNKx(R8&M86O1Gl=frHp932zJ7fPcPq-JqhIt%Pf!A!vB2#31+cp904&gBkM
zC?c6iEab8&e9q_yb3`s-jVNV9BN`-BG89rC57sxSSeu+t(^0!vFDaj2I2dR|q82P%
zU+IKk7LaJ>?(V^A@8rq(-Y6IvG)|jPs$upGDBNZ%TP7bJ8iIRNI*m*?q?C#TT!F>F
z6WC2GfeW^f9IH_Y9o(eXIA#xRBd9~;)^5D}UThTcv^K^0H@@OhB8f(8KmGbkY9@-y
zjPh(1))9L=qSV<r12j`A5at_<LI%U>cev~ZCDWo)E2N5`N2ZASxr!vRN5vs8Xz~q?
z8s_NKxO5xuy&uIKR^`PvKX55-V2-n|KLI^EfH@|lZO9l1D?z6opqrA;P=kZ%^f{b1
zJ>if9+)`!4$5kefJt7WynjBi^+>veMMwiy9-+cc|(J{o=TGcl{`oOJdlWMe%^KU%S
zGmukmRE$r;5ue7{<QO(3%|Xn;WMK}gUdgiRGzyu(@0O|}KAtL$?BS;zI_LbcjeVnA
z=QM16>B|w|&{?#XAAjJHw_}c*-+T=7l?ga<Q}Qk(Glf)IE0@YJCA0wnaIi)&hefX(
z!5ngq&n;C)ygYRb*+b%x8`s0E*SQu>Z0xokox`;K!3SU>5a=xWTi^P?EANnLbk56f
zJ~j`7!P=<IPs4eTP~GHUnG!PKFbPN4PM_0Z1`eA}t5l5mTvAQM%hyDaJ%~BXz+uq4
zmQHP>ZM}MjdFLx%33HK9Z#G_i`vaf6Q>N89uReOwI*cSrqp~og=tCSyNUgQ<21m?s
zfnNx}tdTwDYI3kKhjQ5Kl4>I!zV;aogM0bRHrn22uv>ON{IJPky8WFG{IV{&R_D6?
z@r(9h#8nzqg;^y9#$ZsbvG4}j=J=?u!NI{CCcSdhuG4|-*5i`uA|8P*g6sit$c5`+
zF&Nw{=fKf#u-W#$_SFy%iS#DRo$r1Skaf$ode@zgA35oWsWPgHH8=wtVm3L(%rQD}
zn1!QUmk%GyH#k(PevebCkGKW;Feo~h16D!{sJW{bwy_S1(PrQO`qwar!DPGpBOe5%
zJqoSfef_ORZaV12Ms;Zp7MDJ1K%+79hdAaK#VZ1iF&^eHfO*cL)2r0oZl}}`aSIJ0
zWcP_f&L)S+Q@;ru1178E;2U2L@{!nJbX<S$gOIcrbKLv(3tqfiG^)$<FdcMuH#xZG
zXqN|baQK)*uUB#~ho;@-lo}&0p)rW;UczBDn!W3nHcscD+3Gs{=o>-74T;g{y8oje
zgr$8-ox$_qJ1_V_r`8xXl|{V4g{56%7SQ;XDDsF+77<tI_PK#W#dX5AR*Qhc6mf}6
zLD)rz!%hn<$QEOx9-|#plht+f@kdY(5pX>Gv9E?DxE`KI-+SU?pgy%xJ5h&$f~bC-
z*21T9m63idD-lQF_PT9ay^8118rAAvmqTIz4s#fGhT<?BHFkBGjVAZXg`JIM<v#xQ
zx1c?QCZp@okAEX9#_i#K@<(3vvJgXU)J<*MVB+QZ^*Sq$GRBXL_=E<NXpHancx-CD
znr~N|G|FM%09{c8)pP@g+X}j(84A60?ic|MoBQ;;-wpE+-(++@`C}i4#ocm^&il$w
zyy;~lzS?M*IdlVu)@RV$xP5F!SnL<-jp9+vVO8nX0vl+i3NCPHA}+BF4j_Hvu*aqu
z9yS?`&V@54tAk>)dCtD~y$Bca%trUiKlz=os0%oJum6d6;W$pIHX3W^0Z<w|KBK|T
z>0$MS^#Os-AZGL2Zm(6L*9fh0t5zll4i#|N!6pM7Ub}K^)MPR`=1*Nr;IMlyzW-ye
zMsO`g&uc&Rqmc%O|Lvdozyq5Ewb49(I|kd-w9jO6jKLqX3cEuBtzN?BdE7p$+@KX%
zWOkj@2po!tTjI>1pkEU9JLDpvR;w{pHZAn_KCRKd^W_gfCmhvjY<EBUvX|E>mMa|(
zzw?r96zO2<svP*C&BJ!B)-=*dM;@!2qfiN$V>YW@FIFpfI+0N+)L2XciPtJHN07@Q
z@Hk9j9#^4|Ycq4Yq4r*d#<=+SiJgU5O1a_ijmI!%_+qK)^qpG<7E+4j%FL<*I<d#7
zP-qz)L&&PPvL#aP@Q6WgRPkk!F(pqU<;rzhuE3_}YQ4x};#y5w{@AEkERu(c3Tj)A
zSO(Uf3j^5YB_h?zopb$2#~5FrS-ZEbf;o{Vl7uFWu&=eK#bP<V1NMqa9g{C)4>43q
zwRB9#ACZj81fwE_oXycI*$OKdjo5mPd~Afp<8ejKxR}z`%@azKYwOA(G{omhsyl1S
z;da(2M>e%zgA+>rD3|Zds6cDef*K$igm2*#Vg-GS!=%zAVyS?|9UJDegd8SUEMkr*
z#Y~Y7$z@E1OvD-<X0zEGgPYse*2U)VJk<%%Z>XbeUSg&y8fqVAviPamIFE`(ne0)6
zUkv6C86Q#%wD%zqPdLaLp;4$j9)E0jgf%ou=dkE(E{8TO;?cM=I9#NO_?+RPq2b|S
zmQ2s;xzRDqWE-PVPCx2r439XI(XoNHK^mRyPW#|QgE2JBkeNB4knz}%0{)O78XFm-
z(CF0O{*e(jV~|c801GLdI?Q5G2ggRJEZA+Zsbg%=ekoKcb%-mWb=_#E4$-9+3j=<!
z&LEYhwp$pzH~J|9bge^A>p_DQ>L6Ff=mNVp9a41P=t9G^;ogCPzV03xjXu~<9q1kG
z9~$VR(uVr_8MMA3IB}u(4bzAE;T8|5&GZpC{Aug$qp-zdtc!x0&6SD=x^8s$^ip^-
ze4YnKZGH670Z>_~6iCr|qXU%Zfv%q3?vBp>eo9|=Z%<cWcVBOJFB~*=gI?Va2SpT6
zO8dIIKwa(V?CtM_KXwSmYZMx-s|~fmNgEs%cEO*N>gedA(g)!!9JO`!Qs5=(h9nTs
z2Ip&C@V>pRt+TVM6Z<;bJ3HDtIy>9j;Tis<TNh~e;A;m%=X3X50zVuHbbNBO3a5iV
z9M=QW5OFZE{%{-x%(lce4_;bNUm(KV(sII<=Zi+h(|ZnZ)#8Wb*!X;=>5MJU7mZxl
zdlKhtdA?{AHjRup+u#?ScPl?>!&=sOO18qjXjFb8p~4rPcPl@wzvV+yPkmpo#~6K@
z;o?)@7woa7o@Ti5)b|B@ocX629z6AZ!5&}rX@(b1eP5uby}hj?IMv>EqrDCEiSEk|
z_~i?2?HyfM{b=iGYlE;Z0MI|$J3B!e0X+sHwu2th0ZH57Uk9ijIH;`|gdaLVvjI~b
zs48$xsJ*?rtFtFM+tt|ys!TV9dezg}*3;S5GXOeSZ%0pOXZHXoRGl3?@Vc)Lt_}6}
zz-w?lJv~6)+1=UM3knqu>Ht-(tGg32?C!(n%g%04xxm!f4M!F|1IhV5um*ybHq5*;
z(9=%o?xir`KF&ZF81niVpkVfN55O%c_`{@KXkY-_6MK6pl!5-<Zn!bkGXRQV59oE>
zpmO#?rhPsAL$m?vP){GFcaSnj?E}kWKP|fq7l8))Db!KkHLbsMsE;xP*Qe?9o}vE!
zK@OKa+}}&1^i!BjxY$gm(Lh5TpbQPs1}VM$)c$@tlR5~;ECap6!}I|vD5d>_%n>?+
z*+(53V1QYeN}&x>nT0xSpqn-@2v*|zENT~HV2~k}iAPxd3@UY4B;m6M`dMI&<cz_W
z>ns*)V30Nhx51dhpykr2R2C@e)WJa-h0PurqJe@*9p-RHK)0tchDMpABlICAoi<k9
zV1QK<%$;)0!!cSndvKVo(yF+e0XB^`0vj0K5Or*rMi&a<t`lbr?4pceI-5N<LZ{N0
zG#UrYWHc~nQ+Z&)WYWM^Jt7oyc)}qjdze4U=Q4-cOa^~)n?2k&GRzv)nP22HdU<pf
z*Jw2e1%upS2FGAhi-(5<ECy391KYHK&mU%vj4-)e{upzJ!G_~=IeToF$sD1HL~#7Y
z7-cg!GNn)=A7*n$Bz%c*l*MI_h-ddXtbPt_lxuN6k&N_-n4?06$1ao7z{n=BJB<n!
zT{1et)~MlQ6RAYPVspmWLZO(?rjKw(Mx+`pp9TA#VYys3#vS2tMg$tYOsQdVg=0#w
zQpOz>a>f)3M*^_Uvcd8ZdRZ~nFB{=X{b8S4IV>F=lX!!64OsiRV;qwaF8C`I3N~Du
z<Va--F^4%O7#&ragyK=C7E_~9^8~Ol9FrI=Dy@+%kn;3!;R+5W_&m+ZshC5BlLcuo
z`Kp#Tpc)e>W9g_~$5L{6ib%q1<cw$qJigrqA9U!nS~wjN@s%n#q!{5zxLloGBIgPO
zVz$v}5J`AaF;8K4=}mTyL@BapEPAO>ArczaZp!(CGB{UC6kaz9sXCrSn=NK6Cbm`}
z)TDD^8-LU!5sJKS_?!Uj5`2+ND$?o<YSE}rE)bZ#uuB(<Ww1}TNaaGMOsMe$EOrl1
zu9Z3sPK!dKmP)K!S1J)rDUxZ%C*QP622BFFu{>RJ+Q*DyiJ@3cx<y=@Tq1>oA2!@+
zvWnoUtkh_->ZM!>Y^3cWwN5OPDS2*}Q=yh<R1!ll=J14sYNNtu@j2lLOd)gb-PK9x
zT8T<unSaM6r`g3S>&)Vm*Uhy`WtPcWE+FB%RdQt_2B#`+mrJ71YLqrOAXV~ZI+@Iq
z(3_-6rB>kg`!qV4K`Xb!)80^AqO)nj&ahXnhKoSHqx(iV!ywn1XX@|yRCKpg>s;Dc
z2nF~~h0-y<Qi;k%eyvKA&A{1Bz~_^z4SJ2s<+Eu-N|RC<$eL^ljm993M8bNL(ri>Y
zGKF9~D>u3H32!1~*4cEL(AguaifK_9>`S{}3Tub`3Zr*pe=QyndDR+EeQPGAl0=Mf
zyQl!)H$}o>mCj<)`+RUPE74drnrO-9QtR~=X)+l%SvBz4x~Di1%M_JXpD7#6#_dL@
zQ6Ibc#GxIr!}sm$Cm+NO%&^)L+&|w<CnZ6hKCpARnA6G<7K3TBii9W`kAv~XW(tL3
zK9fxE)a#RzPOr{nvMa{38Jkn@b{GN^wRB-p;|$q~v0~a~@mfsj+pl=^Y&hNX?cMq?
zZDz%__UPH2lYCYdHJTzvH`hvfW!7%C&dtJ|;A}drx4=hav3NRURhm3T<M^C6Xtdc}
zs#3Aw^q72Zb7XcYSE=beG3R7zGVgN)?bh7AH~gkipUD|KdGNKoZ6s}QC2zlQt5Q@X
zE!M=%>zzqBIds_^%ZqTQv6#=B>^_epmCnT+YHPq^DJ=)17N^svnVhJ20@jevmRMRZ
zPcIn*DNk*@RtdTzE=TEwx5HLW(CUredg+^G=V;#S%|3kf-c(hUwK+1^FCNvb+N#&-
z-l)SB$7;E3bp-?N@mwkC);hv=`@}{hY4do3`q|m(K-durIWud!Q}gSV@OWUUv@{j<
zCA{v*C-21^+^8dvy!YnEQ=YN1HBfl@oyT)C+Je)afAO`8C5K@;;PLHlgIYf`HRbR|
z!@go^GV3$A<1W|iZX)aS`=iFi#rbgDos4?&JIA$^9ecbOUYl5(O9s;c-|Q=2PPzCA
zS2*+dOW&^fIg^e^<*g52U!K=jJihX)@7`H+ndT$@;PC;RZO+$fu0SFdtWL}pf+k<e
z>s>fb=RKiN!ct#diKTqmq`!Q4v9xjIN>yXqGuz9VXf7OFc;l;CuORJ><zM^8cNasv
zS$BNu%OAhDwql&}2Pfb8@Zq-CvK$XZE>7X{`pV*>H=IsIre_u^5le795Lms)mHp9Z
z+P1a1o*EApGNH+{J8Qe=p7H6_!NS2>K2eNCR^RzXAt1^HlBIXQ`@MRUzu-;PzVZF9
zZm*kbp~&o4K7Q#SU|&x~<9BYsjrR3=JrK)f<8$+?Q*nEw5DstLDNcqGiJWtPZzo%b
zl=IQqtA|@hxBZ2=?CHwsb}3y=#y7t7aXBn5hBFgi`SI^>CIl;i?BaKS@>>Tx*2QRi
z;iDgY<23Bt&L&b1uHnk*&gNz~S<I)FmN#lCXS@=N?LMf^M$_q{`{ejAUx`na5)1cU
z+CRS^sx0L%*Dnty#;3EX-4DMz8Ix9H`I(P?`p5UuqK$B2^~Znq`=^JFdLp&*y`OyV
zGUhobq{p8;g4?u*dwa2LwLHGIzPmW?Nlhh_M^B~~61iN}ck|+`Je8WANUyy3#>wr+
z(W$la_0IL_OmRAs+Ir_}Rp5x`W<L1QkGC`8jZl8&J3s!-;|9m_H-F@{%edztpUS@Y
zu+!2>II{J%2FGMFdD!3>uX-<lV=^@}kzNLlE8wV?uea|uIo4l&uM(4Cj<>%4<;|=F
zbA0XN_m7TnJ(j-w^+%UU%#nHYpx@dD9P#W#rNNQ&q&~&rxwv^=p7?|#MmT0mP>=cx
zua;x-YIJ<&rMF*S8<%bZ$9wO;deGoleCypimub&IE|qzBf5<)v9Eow@ScQ5tI8wMC
z$$YlzIR}m^=2&XfBR*BHTyNhwtCeRmspUH_lwu0Z@!;_bt2r5PEWh^3gZ*R2S~9ix
z(ktheS<I1saF69+9PaHW$0w@|j=U!|kphmHMa<zozq}|{(@l=ktH<%_CdVA+SUkHb
z#*`DlaqH^Va$dd#9FOi@?wvS+W8uO5{mYyebKJk?xW*3m_fxs4CP%>w9BIsvDl|B5
zmMg%KS$y%0)2k<m>3S9FaXw!Gj=BAlLQFLg$<7=f9xN0TTfzMDozsJz2FLuZo9mYa
z@4<K~b$v(V7Bo2~GAs4n#U{r=gCkRQo!z=wE`P$2Y;fE@pRdklQZwt@`Is7W?9|uh
zib~*EI@?;?K6S3AQuD|Ai<c$Lk-U3b?vWnuAEa~BO^%{BRZXW4p3E+#3Yn_&?AB$u
zoPl~Qy!gi1)sy55uE*{31<WzExRQ%$CL)=c)tTB%Nwp0ed-Ju;2FKja+SFymdyq{f
zZ(nJBD&WZErzbPZP>%*j26Lo~>8kVe>atw=gd;UmpSa$>x>%e5j_TCRcub2qYPr%>
zS&cc?%B77nT#wpnt#CPkITBa5OnyCZWb-pkjuNiNQx3=J)va<7bIdpDktQ5VleJ8$
zTr6edy2(gpwi-!JR5UxlF`rDWpSw2FsoGpQb!*CZ2pqRA9YO2S!Qpseb_(jTw^a6}
zs+sJ;liG5ol%8;$-o7dqvon)f;5ff}DKlH2ytjRIb9rh$ld6=fxtPA39G{(xCMT;f
z&I`*6iS+uZZY`Zzm@lPo%^45JGwGX`&Zuhl=%`p)obk>t@2^hzla+k_;PJvrwve52
zAK$*7sEp4{<!Ud!dGq!Y;F!91aPM-ZwwO&<t5e08F+X0YO~;c{lcV)WX>~D?**Fo`
zGr8r3V&>M0^e{g@etzjr2{ulSC#q|6x|zkp_4#13QZ667xLnWW@-x1pI}fKP^E1<h
zna6Kk-F=*?)o1RVJ-St2TF#~?CTA-#bEZ_Dn@y&sr)bNua(yYD-8^QkW(sSIh3u_O
z_F<t=IJ@-Z7%OKdvy(du{E4~at<_MnJUMajVtqY7UaAET?>?^07H4KkQ%~N$bNym=
zZmo9j>gB7o`s#RQa(aF;W=%{^%+IA$Gc!F4iR$`te0=Mue<fSqT*{B%+NT_rO68MF
zU%7Ao;&gH5XnC+Sd%C+7NtS144jyf67BkiP$iek13-jff*~-LA@7%lpC_BG4fA9Vq
zw>P%dbJ^+I(oD=AsZB3{;A%BgOHFRB#B<w+XnA~McR828b&8HE)yd<_z!aLfJgd)L
z)KPBw{NNy#EHBI-JlfeVr>B+@`}f{hS+35`O;lfg@8QE2#uwHX?>~O;?$&-CIOdjT
zVvgw4<or@@yjJU-PFB}eQpN4VfyGRDeX&ruxlh^87fRd5zH<NU`OZAZikh1_+1-r9
z%d@k)_ck}mnaPFt!Szc^FsN(Q$(P@~fB%Kt{KmqAhi~28++N9#PtUK+#GJ88b#A#(
zn5)qz6P4AKY-Rh1vyd(<&6kSjo9vx*c6?>mn`6zLtV|Y{X1Mv<(bjS(TArHTzFJ?a
zWGCihhu1GG%vNi)%FN5}Twgyb%x}&=yngNO#`;nrH?y!h6LTdB<++toWxmF(Cdx~z
zg~{z><$MzARVkgWN;l!F@YywYO0;k^UC2yS)%m%DjandFu1;>8*H^3K)mrrU`tJN>
zwKiLxd-?4<cORAJw&q^Adg<=^>OwI;vsj;rxs%yqZM8DFP?MG8rG?e<%=U?WE|H#?
znJ6F68`nb7XnD~U)h_K9Qn5_ISy<SsSAD^9sk(NwvOJNi%tlVGFKZx-+0x?6Z(ZGf
zu`;(^dvx*W?%LvfsW7{|HWTxtlZD#a<m_V2P>dJn>XUOjr_owGQJk8t9!>jJe8FI5
z#u;?1>?XqCO%#`QR`Q-eIbUAgTUwkflxHGm*Qc|^iTT;m$}4Z)x_Ue@w=?_r<o<Pi
zE)xnSbH${?;Bf_0#aJ?%k`<i6bRm|Wn6uBM)457DI-4*q`@)vgl-8o1--(0+7Oyin
zx4oS71jB*w%Kp+~zBn~&JGwfb0jW$E@=spBb^FC^sg}Gxxqqj=;4~TSZl7Hx(ZbbV
zpVjViaPw-T!)tW}<I1UgE>kF4qBiN0FD6e;h~%Q_oj^DshfCPg(B$3_e08yWxU`ro
z&()NBH*ZdtQWH}-&!bmw-FXyAX6-jmA6(U!WI~}-p_Fn)cwC`WDVEBm^c+_xQ;KC8
zJ!hg^OlHT$IthExmmE!0Sfi}!jwc)d>w$1$bII+Chy<eLqot)_WqEFF>+EW}94SwX
zD;~alb^U_N@0A{&J-n>1(1(WT!*D0IA5Q(~3>thc-#gy)F-Ha}Q`J~9MPtz@3%(2`
zUV=lq!j2;x0Bg=rX=BOaj}8qE!Qa+ga!sty53C*EnX35m)y(j{Cs+3#8m(IT*7*w;
z_4?<S+p<r?N^{eZXdFMz`{09PY?#gNSR(-(l3!oa1@OD2v!!M8)XpMWIe`6yBU?$M
zyN_=_c%V>o(aOc6^R*2In20Di4H|-jRoW1Eyn};7BMj^dO(&w`<8VU_WkOk$#AZkD
zL@3_i2$aShp#;1epC8XTr`P9EZF)T(^%W<x^sD=;8yiNOZg}(b`k+`=f!Tu%Y2g-{
zNF)~v!K35zMS2zXnI>$p$qL*yL(X^xWiW>>m`gV}3_f=#l7tfnciim`%q-VXH8O6s
zL?>o)(wkSAT+S8p%MXrjZTkE^i(ZC5vw`nFOlG&$2p+ypXO8)>Pgs`37iYELlV#E<
zk2(07+3^O4Kqb$_QgA_7ZjdWtvvV^jr?p7Lnc4Y*{p`$awMNGC?z8>V6{V^ah5HGR
z7R>qK$apLWo<XlSQY&HKP=S`%TE_Rs>4g+3VU9uOMzO&$$QZ6B(qKg&9vL3a)~2RV
zf;CE`l@^yu$>W33v9av*;`r^IgBpXeULE&0P2q)dWofbq*`Q3Oa=hN~jic1bRvL#@
zlV}2Sz%PNU1{(;!7&4tnG`;1y@<j8@VS24znK{^o;EDRt(!<T|3HU3n<xDs^Hq_ro
z&6O(ai#4#w(UY0#`OX~J`i3*q@tcPeqr?5oQgUE^c8J=?WS&$b3>t+t%%aaHGpV6L
z<|v(BoGDE7C)w<wp~;n{;^grmnCU9Z`}6m=_G`oRrD7^kwu(j=ylk<uyS1(q(qy2R
z-#Xd^+ma&9E#A4gY?iY$g(PQfUBqLmHAfQ(rGzb!Yh`osOh(98o2AmiWNw00G@3-h
ziN(2G_59pqvX&NiX76nuugGQd`9yXuZ<g@2nL_3G<k)8%b4Ih1*B7U;h{%&tPdvEZ
z%eX|L@uX_+*sKu*LI)FRpG|6Yg`IQJOwOneWbF2QB{M0SNu|xk%G^w*eCu{Pl`G6`
zPTk!--*P)<$798f&4@wiN#!dy@7yT{C7JN})WbXXYBSn=(mC_^@#Sh(JDEv3&hAFs
z>dNHaWUiPn#<NqIS{RHkzRGG6X4TZBc4cWL5-iV5B#U<+tSoQjX4k5B_HG@Hr>8QJ
ziQ|VAk0qVVRjwbuv{uy42Xiw|9=>*R=vYa_7hZeg#oG(c-BdDm=jBSuwza!ERa%?&
zPA!7E8qAan<C}M8XUB7iDd*1RRi#jxEXVQ>o?P9!KR&%uzI|}_VrgzB8=X3TJYi5K
zlDX=GS6^E%NalTo+N+P>JUh^?#?#Ahzx&GFW$j@qo4Ehxlv}-hus=JoHEWF5_7<i?
z`AXhbznq%!`BQQ2+WBQ97@C@P1n#_ed2u^Gxe`9vxjI{zuH|Agw_cr=@S~~R#G|*~
z+ANICyGsjizVgn+fv6tKt$pb$Z{DkmPE&={i|^Mof}NAYh1vZ%sdaj9X(C*n$m^HR
zs?&OX#wuDmy?~om(+euY%@;1tZkDPmrtOW(llkg=K3==~R&9h4g00Jw_utzt&=#DP
z<@err|8k$T9xH8q^_%ZK*kE0zDw$Wly2NAdUL3D1951pZ6MHM=$mCRxKYvu7<_pFp
ztlIIJNGP1Ca|EXk&rVLOl~vJ7{rqUIv{+2cKX`wkzt0cbidR1TV7Jh-V4J9a@XiNU
z`+XbHiM@}$`;`~A`|qTt#^3z7PVL*fI$dAASf<d*d-Vc*+n*bp*(*#94(4h7lZPiT
zug-1{3?AG+IXavwtq#pBogUP3E9KPUiytn-gqq3CzW&kI_6umyGPUvbFMsXMKH7@T
z9)0(ZeDm=>dXS#WzxTZ@_=7ap=i8gNVM!?Lt>@wkwH%t<%E3ul3I6RL!|<Lvgx>{r
z@91E6wpd3c*donl77Ov({da5qeZllN)Z<HA`QAlawf@1IU%A@t-HJ@?fApQNKH8z&
zO-*Ir{?;0;Z~yk`di`>_pH|+l7h|(C`Jvh_)MF@5>zz0_r4G?*+Z5{I{gcDP>GJB(
z?8^CmEk8FNoxc6rB#Q}uerf8-dvCAh>5Goi^1DyIbh$me6)Em~?HgZyVT*m0DyQH0
z#<G~Xck6g%`FxQjnK)RfMkl8V!li@qv_Lp69-ckA;0YzOYoomL2WQ8pQ`J@B>gKHj
z*w?1R)r-ePksz7OPCk0`)#W^A(Uo6#^ZvURo4oBvcJs?$e)qwK^dgl@zV?*`t6=Zq
zaADzeUSyp<1oKyAvSeI4o|rZm(-!XH@uk*in^}}=uO42U-kg|RHSX@+JzSWaOax2E
z_p(}bHj$lp@bcsNymZl%t-W&j*4esZJCay?_w6_DuNhBMsrbw9&&JfdXZy3Yqnb8e
zJ6fnk3e|FO`)qn95KKnp_0!v)U}AdK;<^9A<@xRE^lI?z<iXM6Y&Gi7@7+o`ZN)^k
zdhgMLnY?bvpPG4m{QBvdX(tj{dGpoR?yY){Qqky>x2Gm;yC-{7(+4w-sm0^@#mIQ2
zS~$E}o6Q&F6UMFc>vUmyY9g3^{NnAK_bav4;@!)~$IFZ5NT76dGvTopBdO~3^@Cbo
zx8hCIUO0L6Y{R-6jV{0O`b*by;iGIk@#M8)*}Z?fJ2Aa8>#Hpv&Cf-1l|td@)@-d%
zOqNZ17k9Iz$?2(3;mPBx%lqY8y?polg_D(qVmMSiJ5K~0#bCU8ee1z|LAU0K%|AGM
z<$T+*7mF>u@#drJ$<%Q^m45OhT?iZ=?N+8Y=ORn1hx0S>LN!x3zB4mlETs$P<I6kw
z>cq@Uc=F|!Ze2Yn&8=1Mo;)~RUC4*RQ@73%5ogI8sop((uvpY@x}%HtE?&CWa~~w)
zORvBE;JT1MDP_iAdKk?{j}LdtGi&py`r5(#WU^R|7f<d@FO<vKob~+H?cyZtKO*z5
zynK24VR3$A^7`=JS$%#y7@fU)l1R8pu3+`<!Gq<Je#aGFzI*xP=ArK}nOJ=N-MiQ6
ziT}I3^LlR_N!L3j^mMm!q$rAF&N+Z22oeCnoJlZeQlbKtvs&sfy|X*Bd%VZ?vG4rn
z`URzy_B?037k<vzH(M7iGzs9xDpVCJeBb-s)r;l7eq}P5yUR17Gai+WCl`ZSw%Xu@
zyI)&7wOTRf`~2l|tyAj_;`{&ozde5atvozx|Ka@Whl4?$OAdd%N#{aUKhyYn@~8c(
z^Mc~{zI^_dPq+L{Ha+^c|MmF~@&0{VsQ&B6K(h4V;;h^q?ba`jFNT$Ty&0~2_)Bx#
zs0*3EFJC`3`_=v^dHO&8?ctX{mqsVuUrxU~91aRxtn>LYnGBVJT<z28mr>bu668mp
z?tcAn#a*P6gTMaI&tGzln|iVI+eapqy*WP-IvA)Q9iEH%sB^;gpPHk3wU{LzKYgrs
zs@-0E_h0|(hcACDqpZH1eL-0jIlg{-mf*wrFkks_{-IlNAJWm@<Ha9t&Y6=`y!YE*
zKYy&0u7q6f>!(;Id2x1JZI6dptoz`3RH{+o{A+7a6Y_b_)x*bHv(|1$+kgJ+{l~9D
z@1Xh-UC-X26yeIJC-E?oj>O8)jcw*U<50Bq;pp?lNqCZob$|KOmk;gqMIn=V{G844
zr)QH|Ykv?Qj8DaSB_zv7pIW_YDO+$J-+ioN`rF{kfBcub51%WYz4C|4pKo^ucp>NZ
zCb0mUh$V`bSErS<cNmV<?#2%%6M8Ssw?BXWay#PBGRfG@mvVuboE|oshkdTKbAt9!
zDrL(1AKINtAziZXUVp3zje3<!efjh4?PImIo4>vO`3AFAK9uSoL_J(ImCB#ppuZ-2
zOkB9>-yI$VcVdy&!`;tUhxAD@&Yk_-ECu(m*3>xa(beHG+Ig{1DDT{MT3FN(tlf)`
zLa|vZv+U!qw^#SI<|uP@`}KM>6hy`z=LyFX$=vDrUOwsWg=4~1=lbA)+KqC}hnr6~
zd-Q1{#-4t{Q2pTSu+f-wLbc&ZZxAaKvxVKecC(UCl+68$$5H_dhWqsE?ad<?>6@Fc
zSA%|8Wb85?IvE4wcqf~5^g~hMqIt1@5ZH^da2S5P+X<Z|qU`BIA?KT%Pil?BR=7Sq
z?e;NjPG`sWt!A|lFBpdxkA*xK;mDWY?(V+So1@g-_1DXOw;TmyhxSluKAAn;AEcAE
zeux*&>Sue0{sW$AK3v>?7>3TV9dUxrmpr~Wsx}XsY;$<l>BS1=WNQDR)u<L@dHwF?
zLoU~>mqXF7f4=+hOT9HpK3sjh?6ylHV@P@YSq_YyZZc^e&=KLJay&ltP55y0;q3OY
z7ruz|%*k~;<vG1P7Mcfjt~ESwcVh)19zS@*22LrO)$UzBWU^q;v0wgj|M**@HA;NE
z{Cd%8mSSKGu-lYp!RWUWDdP}~qw-{T;ydO;&4<(L&mHDUWL!j&?(^%Da_gYRcSaYj
zPOMPDzQLzvqh5|?G>2CY=}fCp4kmv6mk*!*OS3(SV+?oEu9tW^+3ERR-W-?A91okZ
zglQ0rlqaSAoda?b2{k?(Uw-U|FJpZ8_>#}LPOp!I=3#?v^)IIRh{X>-HX5}upV91I
zK4!ou1^GYz<-^C{nw>%H(=`~)QY4scb^{J?E|SblcIwf%xgU&_CdKh+-+Rb~>i3h2
zhfxTO(9u=2>^!-iRGUYwP@{X%>BkE~B6aXsZ&b^CRy(-(n9a5tMVk5L&mSItX`xO$
z-u~gD(<!mRM6C_+>s&OIIT%;@n59QWiigF~aO~M<!}Ys^v->d(27PoLuRD&f4{Oci
z4pnboqKC!kI=lZ^Z&n08r|X=5&gWW<LMZt4&mSIsZuh#;hx=bHyS)-aC4?4b@#JEu
z%-+5biCVi<q;Oam40hdnOt^NtcY1pe1S2@PNp$Q-*9VQ}Y0p<}UiL=`jMoc$54ENs
z@OgdX>`NivYUIPdpMU%C;bW)YiQYf_@v=88hp9xV88mrvnC<T#mO0*zr57-IgB|xc
z%+zkiCpX7nM5)PbYG|8W@3)%g15c%XHQY^>gj{+2P-~T;jBTi#eJvHhVBBB8_|P4;
zqjw+wbk*M#^HHk%&0;>%JEJ1YTRRk&+sn3lyPk14T)n~;**<j<<Eg{Tc*{1q+^;uI
zI=))-Ld-`go!@<^Hp*p`RsHO9F^}c`knhW{clVF&ew)9$|HDOhP-Lk@x#qRvc_dS#
zcAkw|x)hh)O*OlFUX)ekVt8;p0wWSUIO7G|<oux8IBEo&?X&J6Q79(TyAPGRP~@|k
z*4bk**RJL1;Fn)-?>@J>U|j!<W>nyUP|$VRJ^4r?IjrWGsJTaRnPI%v*(djzkZ?8_
zU-!V^=>20ZYoDAQ3XOwmxZOQz_Y(PhG`{mtsaEpQjH-KfKV?wtmtSx0er~oq{N4HI
zvt}zF2_}S!-Ra45@x&mXVPd8}#ijdvrG4l-U_<4T-p+Ls3^ufP7>?N|CzDcRuN3L^
zkD8rCF2_Z79)wCY8%-+)XE()cr<xB$zx;ZA`%9w<#_`9~dOaVd6BWVg_7>P!VvtNT
zapM3RnO&~jItrkDl#jc^>na#bcyEsmJ0{0Rg~o0^J{(ROtwb&xW_BJ0p_YlIws+4j
zi`i~9?@xXC_3HN5dcDivA3mH^D;btT?`d|3jCePh2*ykUicPk|h1x!Oz|f_GM(@1g
zKZ}M!JG<EUJvccm)b|Qpr$1@5<GBpSjqZeUB^^u22d7w-?^Lor{?jj4mtX6(cI4*x
z@wiq?g#vu0V6=&hSSuCtMa=`e6q~`Uu;V!h1&iZK8`bEPr$fULlW^>v?iXqY#ZbLF
zX>?=xbc`Qe2_+#FOKo+Iuk-niknx2cf4(~X*r+zStMgAs^=2aA;}SW8#ft$-tPbk~
z-ZJpB@rFND95@ex!NRc6m<+rp5sK~)*sNoGHZD|;N`Z2F((FYG>12F(hCTIUJiXO8
zy2%&X<rL|A_<S|FYgQZV#nqQdqZ1=ZhR^6sA|qN!bB>60;A3MoZ#>tr9|Wj;FH_y$
z^&E2ny4#CjhIzhIsGO9&x#pxf;ESnja(E;Zi-~wzRyn-Q72BmG>AZiu+PiF3YRuW)
z*TZI?cY8xjN^2kskyx~pVXT~OK(eu_i%&PL2R<s-P6)dP?jz1mw>#0AV{|bp2<H`N
zvT@WJ@ugHTGd!pi3yDNVRyw%HHc~M`TJAnvk51c_Ds%Gq%VBH8dB8}hEo31aiRAg9
znX`^OG+(fVl4#U^U#7-qM}6{`_xfs85@mHhVj?YGp8}&tauiGV_sgYXAV9BiqpMW0
zSIPR-XE*1&mmQ2yjz*vM8$HGgAN{z>;VV)sm!*Bih;{7pN780ryyG|yc+zDiHmC(o
zV=h-IZ{(c^=Oc$J>Xu_~q1`e%Ee>mUuT(CZ&DKR{_bgr*Hp<lI!NtkmZMR+J57T$!
zdOPF;BdYNDN@R#ld)<bJ<-qA>V|rJt>pZ2LseCZp%ZJVr4r?wUqnt<QI|`LewUTWf
zb!u>eR&KY)#Y$ysYkNl7IO6l;PCc|VIGG$i_IlOWzUyMAS`YetA%@@f`^&B%leF1&
zyyXZgHmt_N_uc0~TQcqQbz<yg%4$fm^G?_4`Of0f`qKM+>!@3XUfSY9bGJ~f&Cf5q
zSrqo!+`*vBylWorAARYMs)@1mXjrLGlt08q<a9u=!xlns)<!L-Hd{ET)JKNiE83ih
zdmOD$^g3f!M=2;;T%2Prk5p_Q^_~`K>pS^67Tb{WJE7e1?oNa#@9m!caky7Y?-F|h
zp#<M0r`s!W*djXEzig}Iq`IMnpo(<cn>TeFs(|05OFGFet&)*8LSbyRiqreCsaobq
zktSXU$DfQ^!v$x!RxEglNV(XU!0sSeB+9u6r&g=9dc%@N?b?FPlSI01+>FR%I?am2
z5LCo8QnN|1MQV+4t8~?<CM4U@C_QC3d|cqkeA!Fey-!BQc3NxaXxc!y!gQn<4126J
z5vE*DxGHZ*H)l5_sukF|Ex>|r$-2BOmB3ihux8&_+S*!M(Qe4>+e`B@39-0pvppGK
zD`B|>H*+M7L7~x1j}(Rtg~e#tgcq{WN*YXR7;3}hO9_{27%;wprW$-@=HJ4i^4<Eo
zcXP8a>|C9Z%)W!i&HLrqdC9`NpXO(Yx9=tEtJC`}F0Vu5M@-ABi_6PTpRL0DadmYL
z!h5UhvemVv1!#Z6(R%vj@BEMGT_d{C{0Qq(c+ZGlHLyg5E6tBTz5DS!Y*Jy23hUGv
zob>kHkMDkftthNSVR{NLADE!Rqerw!75{z@e;?SYiUzYN9GIoTdKG1aKQM2dofW@Z
zgnJM?(%;U$U4k>dXqUS5ZgmDO?K5+*99=_^E-t@cdW(IGg{AeynZ<?Gr4_V~`MH&~
z#YOy`|8ZeuV}5=a2D9%Z8*6K@;ar-5a^mU|Og-meUi$w1{07{dRu<=Hmu5E?7gtth
z7FMKl+l!0qu<_j5no}$<E=x8Sq%(4bTq0F&EN-lCZ*D=hZf#{t2^-S&)wQ1_TiUf%
znPhc-PNPvMG+R>XqEfAwOE+b!D>!9tZdnO;t<ANSMaiOi4F;JDtDCB2<Jy{hYh9|>
zFPk^kBnr){YQbREZ>ubdHRZNZt&=IYx1>6=LZR4}ZN1;tIkse)ZQ05a98{qPq*kq2
zY#sv~;H4Y#E&1Z2#0a-DrEFtWzG{)ll*$#U(!B1L$+Vg+g~PG#-ImGq7O8331)&m?
z&#-CKxh-}TJmeL2QV-=7wPsdtr!^`#yGz!npck#iY?Ao`P|LCy<y&gCW_3m8&}eje
z<(5Rf;ZZ8}kZ{qvH)*BPVpiz_flWrM)Y!dp*E$@JZ7$BCa9TqiKRhf=IzI%M;IVC<
zxA~JMBWW{kN#hBIN%-9^B^S*F-B9FF84X6sy25WVS!_C$%(zACbT;dz#uik>bvn08
zZ{c}G8p<nVNKHv&Y2Hg^ecFIGL9=$h*J<I>BrGf)t|c!gIPGDtT`4b>(<y<4$8fgL
ziqRp`ptITRn^G+#dR!i}PGM8V%w~^UZuIb)qS;IZtgdWUQ*oN@;kY3p&sTDxNR2UZ
zG(v1c^2F~-SHs~jO$Js%>As&#gh-vL+p7q@bTneFwe~Q|qpi51SH5LTlBAz<m^B`4
z&f%Z}DjStCv>XnOb@{7RL(k{%#`4y*veBzW)5C;4#kX?;o#Yw6(2vD2J7L$N!VyCk
zqqIqPaM<fkglyU|9A7oEdESK^22~1M5iYD8=`c{n8cA}j&dpV<yCj)P_?dp+dKe;u
zxw@;Y9UTs`!g0Y{%8Y9rzL<$K-AOK=%S00!+3sDOZ)SPB`TFMY;I>mOV}|<ks9KM>
z907k&sk9Y>!7z)>sVRen1QA=Y?0UV)gfHChSq^BD%GI1D?eMUl5snJpVrIA6=8NfA
zxO<q(Wz&3IlI`BcB8?1hGhJVcj7ph2y7)4z)Selt{}UsT@P+zakrBw%97Wv=BfDE|
zePd+7h>46S7*Q~;4iD~L865KA3o6vM0)>2d9J<$hRLx_r_iq_hdqFpx^wSmOqmTn5
zT1rPl-Gkf<gZ~eV&%;XfnW1@M1WkAzOjze+MTRGy@X`H_8Ts(%s`k8o%9s>Lk<pBn
zQqfTFAeTWtqGCQGV8pDZ>&wH#`;H)>tUeEgZw#&dIUm}fss7BUnf9I-1Cika!<I7)
zkNW8<82LOHv2u!kW#qcI92jw{>E?2Bh<$FscXaUy3>IZY2Df$gLLeAs+>TcUGiAgQ
zzTluOGJLs;HES3g^-@ztVW-*<8FX(ym(Hdlyd>AXVIy@gj5in0jLF3(teHMD^xqg7
z+KhaJLpldjGwyla(S(m0v`k=-U}TI_#$myi1EU@*CnFdq=2BopL<R#!!eYERpG+P)
z<%(}|@o`YDFvy2D2sh<|ANg=O{)SOE?*Ds+U`ZPX$Gv1tWE8=O3&}{ZzxT=r|2>BA
z%rJgqs9*C@H-Nz>e1Snz%!fB8m{X>2jBIJAS_6X<8A;@W6Y~)QBZ+*RO(u^|`FI?Z
z!ElQCFgc#`;Z%n#4bm45=^aeNaPY(!iVVNVz@pGeFIhh<`eAueMW`w;23S{{GV)K1
zl*I`1mZOi|a@9Y%oH9Jnm!P&Rwwy0WGj@k2Y;0nf9WprB=5tvEh86Q+*zYY0W=w#N
z`ic5}fhtr+m1;srhWx$5Tr!{JnDtWcDvZ7)WwG5|9v*#am#e;`)B91m#E>p9WOi$o
zgsHdHu4VKsFG+`tj!4_IPr73Xe_&LwOgKMTsaP26;Bc70d@v-`hC(G#P6oW)<7~2&
zWJ7D!!37&_=Ta8;{msGgXV`;~CzG4qa)E_bm^ZNLv}8z{V$3!jtL=Eabjal3yXHO8
zlMoqc+Y#^gH7bx%L2#{jeK*h48>lc~c-`%@OiD;H^jc$d#>Kn&l-2*?e*ff49ZL&m
zdsllxj-}jSNImAXmkuyStDaMLJ>Fo*?1=R(dp=Jx5eV$W?8kAxzf*;H)NnkGV-y`}
zxB4X^QA#=;t;=+(mJCy?o!t{Y(J!Q|u<G19{kdAMc`ioh2SSFUJ$^E<;WMS6{mJO9
z242|*BW!gfhSoj5H=U#cyR7ph6$lJ!cDHM^vmbAa^09sg43v@0*1Sojo5@gMb+C6F
zOASkD8~^xaeD+lk>Yl6q>7<fk18$!;utAzrWH=BuT6Gas*Nc2u9MO(-?DOLN5g10D
zlej<7794)ZVD~8>jdo9{Cdw(dvwfaU*W+Pot%ZCf`zR~s{$YIbsamdj&PV6_SRWx>
zP?g>AT2elS3R`R%PTe8NAZ>O;8rCuCO~y&8mvrt&eSvz(8M61sI|=c5lv`b~4MLo(
zdz4MrVhp`r9i2rI-E7hlyty48eW;fN&&l{|TrV@wVili<E$wH6VY^MuY1=S}r%eu4
z5E+=`P@SxE%=rRB!4<Lg_C|@?I8SGroq97?$%Ooay<EB)WkZs}@SKgeQ%Q^O@?!7s
zwpl8BCkHp9W-aV=dcknp(*ceSyX;C<*YYD8iO~@*Sn>Xd$6Uc?(K%wsKq2Q&Si1Xz
zWNj}Wh&S5JcC?y}(4%1?UEw)qBi+AbV$DRt>^?i$+rMlVi{AaC`$4-Mwu3=!x}2F{
zBoy{K6s(~c2n2%$C!Kj>IO%%TIS7$~OvaltcMkfg+J1p#>#cS-TFWKbakrEaB0MLJ
z^{&IQMl4~rpB(P*opp-^a(sqOw}Hs;1U4OxDZ^t|gu(C!DZPWvTJb)QM;%nX>KNc1
zmC3ksrp{hJRoyFiBDGelgM1{|onA3hj_|Bx%BaU;2J7+R-ri|DpK$D)KJ>bM2C{#i
zfXwd5(2!d6*c72>hMmq?cVPFzI|B8peSr5zCgaK(JL7(;vR8CRt0JQU#&A%~mUy07
zkM(ZC@megZGapa(_D@?m#ymXvFlD&F*s|HPBExM}&?18l=xlV}I>IyO9sYXNHlU%R
znsw$3ot-Je6|b~KMkc`yM#XFqj5RQrcs0grj7O8bgVRRFr|Tcz_qsjSDKg|%do~nd
zz)%KX7zOLl=S@WHzQ#9(Bd70<1}R}&awaR{GZ!*(wzpHv6{382CDyxT5|tRI(jOg-
z4^Qf8hq8MlGC0hfJ%MeDJs0Acu*;&PjIDq_6wq1eqR2=_>|~>A9nb-9CTGv<x&x6>
zc4UN3vmFx|-ElErz*6~gynoAr!78+q<MHIAnli}RlRJ^&a(mtWZL=*G=Gm~ztn`E7
z3;DGcx?~;ry-CjQZB#7-n(}0F*1WdcAEe6TGOSTVMmZg4+c2ZdV=%fD@87XW0h&?j
zDWj55tThjBd);ni%21e|7*4aw_rfS!`(Us(PxBkYlGpUQ!&GTkuolXlMk`iI$JpjU
zv5<@N^um-;j)i2(!;{^k<8pjsv3@XRxIM^+(rC>yJQsGFRHVoV`!r@+u=WDp6l-%g
ztJo6>x-)rmUft_V8J1G1D>70swlOIdvS7?(r;<yRVzhJ{jN{`{baAG-f79!B_$fnW
zcw#t=D(^GHL|3d`Fd$jlteOTi?aJg$c~!3k=YrjeStyE(A{e#fVlk88gL6|xF&30;
z9h~l-92fYvZ-l*TFrpqW>Gq30UTieNIE-rO$x=SXr!j@9wk}1c8Jn|J10&?h7K{Z|
zA3Isa-Kwcp>@`|(Fu3Yzv6M;j!P&&%18g1R)SB$zZ1?myAAw=T-gO@gFX?mpkdGY8
zM;N<N<1uxEBtvS9OwHD($P8n5_ll~1IOHi5b(zgUdz3Ek)a-SXRWngab8O?ZSk5J*
z)O<!{l#*eo`uKeJ;-U~+o`s?3pwo)^eAp;b8f_UcSesGjG4*MZ#WD<6w+(}2mT{0f
z8EFrSmF1E;u{!MTq${Jkr-6FaNR~5UrgK)R6p}G&sVFjpG{0>=zu3LGEvB@R?&N0F
zX~Y8%`tmBw_7sag2V!3y(=hDiyjo+lVcQLoImQ{-k1zH3u)iuOqI0AEZl*SD23kU=
zRZj_7iW!{2>L(SaR;oiV(7UMJ*H^puA4^%7=N#Yebn1y96(l{|7DoaM&SErqOgoG>
z;?+Uk)wV~IdB#N_bMLxwhH6x{*q?TW<4j}F3b!lmb}e1W`NF&B<%W<-Q0py`QOjoy
z!Tamok6+8hr0wwZZrrIRLv+aJmfM_hXedK~(&9CbIWO<k8&fUYVaQitJ<J77bW?1o
zUE2&1<6Ufx3_8(nt=*|*YI#ri;6i9s(<w^YAA!**maP28+cA7GDlO^^!QMKR1RbQo
zusWg)$MaT`1^HmT5v)EXo7R2ggY&Xyq3_!XmTuK#q3?Hg#+fE|Puj3RuBEGaPiXJF
z3?G~%CG8IGBdL0!Xk{O6c0d0L530b)+1*}SNYE7R7t4wVBMKodk8ucw-)KqTeU9hB
z`?%BakF6L>H|z4yo8f3T)9AM&t!lGXOI7k-diSi{s3em?X}f<9yNg1>8os~T{q!f1
zaeQ_+ZV7Rk3i`d<7$SfXb6dd}FkTKK%JI4t&paCq@F(osW|R#!8;a1|J{a|WE7Ghq
z@ysiEGB`S&GH7YDch4v4`GO4|M!O&X2!{XY>~6Oy#DjrgfK-?q9K-T)(q{1(dJM_Y
zW_zq=9Z_U9PDPKn_w@)HY&2A%neK2WQ|q_5dZp2XQAQrKlaq2?NF+njdgqRh*YkNB
zb$z|__)EFc@J~)}Ve=?5C?6Os7>R%bj4l{qt20`)4g+K+LB}TiT#W;xsR_@)2tHHo
zH`!XH(P)6-r+UX=#L@I>?Yn5Kp3B+&m)AQFe*nXG2u4GQ1^pD|QyT4IFp|(}_2}DS
zl0`4YS1f(Vl&8Xp!{|blVX1mc7hY@*M(GL|RiWN!B+B_fpmS8Jfx$>Ct-B~j5jnf}
z{CelZSCMgWdONBMG0KMrqBPh;V8E3Yj20N2-4hWkJwKUBGfA{TA<R;>wjsP+?+?>L
zzrkRCsNRg1@|3?dDOJm{Br7d9Z)1E7K8)_utDXCwOO*z>e|j^l2~o<22cp8-6T?K(
zoCj8ujgXJw9bT?%>0m6H<x+cz)iTWnYHd?^t=b>JfV>edmuodJ3Y4!QGNK7?vslME
zLoJ(eIFGMJcV9}CI=Odp1DlN~MFwa}rLzaY$V9zh)WL{3eQe3n#=ADhXU3`ZQjiT)
zJC?9S=nc}PUOilb2`l>cLXfN-mMUOGHuJUHI9JPLoc77p==L)hUhI+ft3otD20}r#
z&Q5`mP53MxZ551!JHQkzO^izl(d<rIS_pvAwK1EeUO!#z)kB4HwJI`bvT{(Wl=uX{
znXBH!*;*#;v>sfJZax*k*gb(+pTGy;HXKxIZGJFvDL)t$gepsVsBppDz!<I=&yBLO
zyr1<~dJbl*(Cww*q)F#Xl}bBWD1^M_{i0BcB%+&{%1wenui`T8U5>6lg5lXYx$a@V
zIN<SzX$@iyGvP=+8?bmZW!fLcj3<;g*T5(z3&Wf|N3wpQ=VBCa3`*s}$d!ajo6i@*
zp5k7y3`T4-RlZ4tMTTMba(ML!TV2oS__`|y5x)nsSdH3-HIhi75Co$X^oO!!h|Za-
z9&b!Y75fEc#>@JIzK2m|I^9&RR}E&0g3tj2e!AmgxyU8rn~BnOGE_|^-MZ1m@baNh
zskw*8SDi8#Zsb6#vU$NMmT8MeT>vBRV}e<;;PytVnNqK)PP<uOc|bDSRJ)VP_Nr96
zST1)WxdQ9T!V9IyCK8*m;&qa)rV?(=;9_|3L1YY$uiB+411n-Gn+J?iHEi{$^T9x<
z7~rU^xeP`<Cv?lYlnacZpEV@G$bgXoqZ`Q;Iam6Lk(5OX*Qp>F9##Khcz%!F9e01q
zz^l*VcF8v(Ch78WQJ2Of!%WACQI##EiWp5~JRC_!<dnf}iR7(rgf&T1G}gf=HWEpD
zvEpS27FmvU2dFiY=~DJ66)<gPTlv<2iN=iCS|`Ot!m4ehUbns~Q|NRCr(3tFShZ`^
zS_cGHZBo~^!U%zKziZth*DGB?m0n`;`7PFfTWxW;Tz-Syp_f}}htr`qTNaIuxJBc(
zm=&vmXfT>|xINneF6M&><T|v@q|j+tlS((0YU%RYyg~x00k|)0FR517WEzD^r&^Fh
z06}F~moCa7JGZS<OXNy8Lu@LPo2zo8QVET6`HW0qmrK<0t+hG5&0upX)aq5e#jaJV
zf-~4#nS1*alqcrq7MGyq`_rPhM=&!tyYLSCF7wNaODjun7v9g!EUnJWy<LEm+5(i!
z78c<QF}Db<2`Jj)5BAmO7bOd`OW0`mX+;8k1)T7H1=0<ROE7}}-IM2cCqNqS{{iD^
zm+bciiIUU5cY2kK`u#|+(pSG9W9j$1^~;@KrN~}B`cFT2^Phgl|35c`jLyvL4CH&>
zzLwhS5p*A))4Oj6&+)Uw|B6mcPwDR+eftjL=<xEAtSv*p@wJA%9zE?bi)X*qjVDhh
z@iWn>=_$RofY;ByeFtrQxm>AMU{_W0+8SSvUay@B{7)}~5cYF=?T6sP^qj2U4*uq^
zS#L1gOd8qx_Ujd|M~0U(U%9TA&!#~!rvH2Gy`QWEUQU@FkWLp#A;!s;{`EXjn*S-i
z-tYC7GT*)=vM)zdpSJ1$UZu2Nz9g_g35G>2Q&SLidiVB61d8*c2E2v9Jbb@?M07ap
zrz56?7$*)Am&8vDcGSfh|K_LcRLK&em532A-Gj#leC7@D<M#yi*WW|=@`tC$vmU$1
zFso(T6a<%+A>Z@^;h)3yIE3=SK=j?EWe7$Q7DzKK&A&y=7iMmD9v>iDZ*Ho{DGJ)B
zh*?O(i-rjqL~?sWKz0c_%{U#xs&IR=xyT?Go0}BGkYwv~Gj9n0%EAWryHO!s3kwiP
zUBr=PX<>bH4YhiWS%r`iu3TQpO_fDwAOoKwAoVmmzqY!NStbw=ZhdKe3!&g<r?H<A
z<<sR0dQIw02_=!rwelrM_4}pR9^YPqNUdvSRj$^oLlW1rj-xHWwwO)H7W9M(i8L<}
z6;eg<|1_}y`T6B7*;)qLr89^`i7-OjvW01MRzwAZMllJ-3ME%)ZMro?hVm<=M(l1w
z(ATvg)tap`GzsgL6hb?|YFMRQgMc}NTv3p3qt>sem0)adudGV28=sLA(s@F=rF9T`
zM5~!rM=0ZfVN^>OD7DVwgY3}^;n#0FZMID`4wqbxSS3m{EURkU;xVIjnpuNp>RhPN
z6|m*5kgOvRTZVuP(}v7!+|FnSM24}bEdhj^lRQPBL?|IYO&ZkFSqK2RSr2sa3BMKF
z<KAr*7%H`!4(qjG=rwNIg_h{T!i`Ihd>Ac7qsgL`Z6X+1hA^8Ez0zgZW{d=K>o&L}
zsIr?73PG}&-F|GMt2f?af1gSPRq_SG?>2FvkQz5~=?zpgZZUyjHBzwvG8ACF_JG)$
z?XIHT;W8=YJ}+8|6EgbR0BOqD3FA6JVP~8O`L@N76u8zTGW6<=A7Jmyl_NUD+3^Q#
zsc2M>X6Le4IatS|;agocF2|yHSvKfp#kT47m%NCYqEUr|VEBA$t&yc|883lQP#ku~
z5kN;xOuGWB$rA|spbUl}KuoHgHX)RhpY@coS&I`4rz=%$_%H~w`Q0ffcOV822g^)9
z7(u2Sq{Ci=E*2qjL<my-wgm6V1_}2TkqjhzL^h@ugU*U&;0a(g7iH-Ghh=4ZDJ$-!
ztO_q>m3}F!(BCO5n*)rXPCHNdUg}jqulV=%3I<B{mwLs$)T``Ey-I$oR|?zTsaFVb
z6QDyn1l%TH+937P2Hh`h(DBj+BQI?*{n7^GFKrM}dj3Nj6yp!lj869>;eF|^f-n7*
z`=!5fzVug-m;Ng8(qBbi`YS};5&J8c^1Sv}*x3%yoDL83b#(pu4tV|SwU_6`$2d)~
zZwEx=<>>Dni20ycjqV4+jX0yr^Ka1G(zEZ;@BW0&)hjCKK_L${Ow7?NzMSLdXJ_Au
zF-GQ+v(WNGsF?*K0&A!r-@jX!nSos1%$r37Efr<-rh|a_B|Dx+fKh3GB3#l9rDkhp
z<vW7jSiviJ9vvSE<vI2C3i?DB;-9M|=&v?;C}%Grnx%X*BiWX&AQ-lsh-^sK-mNTe
zudm5AHrHkp2=J<uD&=DLOvqHISCpSoXl4kPPU9rai(6o5wwxBrnhFstWOA1W0^nfC
z<yMy(eW5;T)M_?o)(DFsueRyLXjK-1*Q%BCTboc{Gr>G=(T-rai0EdS4lj%@h~rbf
zuv|TlI9LuY5m-aq3fiILC~6BGMnJ2N#6o&A7zRBZLHIPn6D6H?_52p1UzThkmmWWW
zB6z3GysXoNEGF2~nAXAw9L%~{3cnI%5kXY2AQ_gamJ#2XE;SP}1hJ%PM=6=opfLw+
zPN9}{qSt}=TSlPJ#lZwO6tpaBh-9QfWf`Z|h&NO;6d)zoyY_ivbi^ymBj|N8Tuh;#
zgBW|+Y;gEolyn<`h$uhy;G~xC5v$RX;iF1aDxcY!gGLJ~6>KHrITBqM%`yR>aZyR|
zp*)$OtcWofC76KMwPrCT+z!?sa!Ql9S(-{m#qOD4vIecyZ6_so!<k4=JQtD}zb8B%
zBOHyk@uS^X=1oTsloWlpiA%aL>+;jk`qwQW%2XibijXFSnqUxV+-YAm>LYe5?FrZ<
zJZ=^vW6acSF__S*Om?$tZF3$)`0aT(UY>{j<#|Bx@A-N7UY-a4@;u0Y_&ih!o!MYp
z-IyaZHcOE6tg639R5f+bOG@#h%_>cZ@oO-Y(`q%okQsfu!|k$~6?0O&g2HCM4RU~Z
zNw};g!;)I*H^73;p<fN)W<gt!oDNJ~^!64!^i9j_vxL%QAZf?4;(NkiP>>F%1eM5u
z_$`!Mi9D-RN|#@c&d~yEPJ<j#QXLk*-fdED$OuG4f??c(Lhjb8Z3c^G*@K&rCekSr
zpNGvSlPWYy-J;l!y|fR_OZ(8gv=8}9`!Kz<5BE#^Fn((vioe%BHr8c|E!8|mBb#s~
zcj)G&-w_Ih)S}m~V9cP{*tFV=n`nBQo0|r!7_(5V)hZQ|_bY^6W8F5YHqi;_2`!8W
zXQdLOT&9I#;=BnrGpo!9o`$5gKyzLumCs>RwJzV%nN)M@-w|6|>pGQs5d#F-+B&ot
z#h}YDbA)^$F4|Tqwq&bs7YWt2L8g(fV;H6;5PnBGv%aQ5ceJ&=F{j1NbaEZWK1d4v
z@)D*WQbbY`V@N;uR~s+=6~r)~`zr{0KKEC$m;Ori(qG}Df6!mSloMmicNl>}s6e*8
z@ea}TR#%~jzkx7R1VlZdFS-cR*`=jL>6U0syDov!|H7MDVq<-41;fPIc|@z0AnXw|
zOr^_9Sn*kYw~3p{&^=B^yAs*#%>3f~5^P+FImicWto$?weFVg$URimMfhWXC*CnfS
zI2S8VtLqE+a%p*aadGy?_fRukhh6F%I{p;`Dy8#p=4N4h3cJ#I1d5njT3cG1#+xFR
zVGoSy;2dP@r*W#E-XF6s?+=KqKEFR8=KB2pn16YHU={26{eh4E!TSSK2rQ8OfH?z9
zZ(($cmk$AJTueoO#6$pMvG@!V42a{4hPpTq#bn>UK~!hNeioJbA!z=dz(+qpl@|K3
z@5EWk3~naM$-WVv`MVj&+rmKjyYIh)K)<MI`yHk+_yijL--}06om)|LUpxz@#MAKD
zkiY#-RLPwx^?#2e0s^={KwTGqpxXb#cL<RUQQSB9;Yfd@sxJQbpXr|^@Xr$XX9@gY
zmH;{y(We1Di|EOKo()HEZNRb%LY0Vpk{Cb)_ZIt~X%~lnl6X2G5Jp5?fI&iM`Q|&^
z0t+?hHm3<|GU!yN=`B>ZXJ&Cn0zMAV?G}5mx!HG!#rHiX0_a<35%>gWzMqAH_VfzO
zVz6u?rUk6~;l(2!%@MPJ%^+6a^u2_SFvt1k_#tj&z#-z@cLajStgf$!@Kwm(tStYC
zvlpOuy}B&ka()T#F!6T8`ZDq4QUc}T6#_45@v{{IQND>q0&!}@|08V?;$p}&EkbpB
zZ3Dj=A4z3$sW_k6lq^Z4Vn`y;);44t;=NZSP`935fjFzr7jHwmcnP`1+KH5qAVcfK
zx@<b8!bez>!IWZ}Ff7@UZ7(kTK*+Wg(Crolo|N)+_*T8eoum@kwrX4aX;|yodRkM|
zZhu>>)I;-ETyDbnMnP;LLo#AZg`wm$*$KpLeM;iSl2W^ZoGLXsqeZuXx41#G3BB-l
z7*%b{lp3Q(TqoUvi1U+=k3~yppRTpx$6}EO$)Y3FC}bs}GJd<U9|6IiPQy$}r`MU+
zSAQh*CX2&uSwdOcjN3-5>AeaINeZ3GZW6PvG$6Fn^cuIBz??@+Ud*7L=3{OPfjD&-
z%<1jlZp>n_<LNYvu*_DoQ;O2WZX8LuR#E71rm(vmvl#ChA@)o<#NV&8X{^&J3qd(C
zz@4TLv{_ZtaT?_!Y<Msh!b)O3G)-yDRGtqtxV6K#wfTl{lKv3qUq`_Ny?SqeoWrlR
z866}|YH+Rzk^Ah^Yd9a_d%75*;0WEcA|U3y8<oXLxajFB489i#KKXPS)T>B88B)l9
zB7Af>nqs7Qu#uo8$kGeQwcAG0JRG;MmJxtr^Yofz=xH?vVYi4dVV>5pRG0{$viJ}W
zF#3Gs4v{H7ADUs~8jPqFKM^z+O&4P_Jn1CsVB_4P2N${E2g`}Y?+!ad)byG{gor#{
zEMS&qpH{YFlnCSD(?lqhe)^6W6Xp8zp&c8qY*L3>!6nkAdPa_UNj~n0XA?^STtu>o
zVnQq<f7Bacrq|SxSk9fM7S<85kc0$vtmp#fqeP@wdistyC7S)`LpKvk#_~o)PlpIp
zy;oAI2%?&W*>`RQCCde4xmwPRbLmtd`BV(uJduC8ctjKlthFIQf+vNR8YhyqZ$JKy
z80<ftM)>%2E?u!0-w|1%)*IK=sG#i<U1|uc5nRNig{~kzVzv-|7VX-t5S6Ek4~ZHP
zREfvq8d1QaLYByNdr#jHUx<UV=fhyKC=^?E%X^~K?2M1Q;_A$(&ei+P^#lUqMvKjz
sW)SDbDiPu7=O4AcMc23>I>gi5^|(V+P>D-~urq$X@n49uyXV9I4TH4|0RR91

literal 0
HcmV?d00001

diff --git a/DEV_inputs_test/soil.nc b/DEV_inputs_test/soil.nc
new file mode 100644
index 0000000000000000000000000000000000000000..005a8d4325786221abf7ea9b1890480682e2d1e8
GIT binary patch
literal 14251
zcmeHOeQZ<L6+gBi1Tc`c6b1|QVL**K6CB(5@Kg%MvE!uqYU}`8pss7bB(JsKbDv)z
zLF-^rTTyA!x~d;*%O=XIeTc2n+G(`aq*2y&s`z7Ug@Rbeidt#fM@3OLsYq<>ocrPU
zWu_ukw$jNrdOq*nd+vSb{?5JkopbKZ!}0dmf_cm5IcsY{rpqE7YUMdTQLuiiTkDIr
z4>;XFb$VEyt40*bDpTn^Cpd_uu8w{u7SS)jHcsHidYFZ`E(c7lC!cw}w_LnJO4O@a
zOgD~MPVwxRcvReS{<jzK0<S3N#PgT;W-;*>Jn3+-7Mu=;W1irV)U}u-zE0T)9FR>A
zoDadvD8cN7RN}g25~z>%(EDb?TmjTU(U$DjoGbRr9VS*6ECx4Tkm^;fOm;XJF$^=S
z>1j(XLLI(`&;d|@1jTo=Bq=?(zko1cGc_nI@<zJDi3Hq$1q3%i(%#t<?M-y{^`;V$
zaCa2k^_1^InOTM2a8H!+6E9uktVtJgEJVq$Nhw*ome;a&nUI~5yPf28!F=opbtR_d
za$C~5F|{IySLzACM|%L?{QQSOuA=%GRS`h5n@jr)C#9kHUsVyQYQ|MGfsOwC`41Z^
z>4RRl`t4iS^hM&gyYHZ`k|LR^Q@xRX(<rE>rQzVBG17<^>3pYZ7BxfP;P>im6mAGM
zh8hFwvg%GXXB6_PZh1$u!_6C-8k_u${?v~TAI3`X!1lFaqCXP!grJuMU9*FicMEC>
zx?j-O1pOPQd;j?U9n{|5GlHHJ^thnM1bto5R|WmKpoay0fm0F<&h38ew(525VDR))
z%WkV)*A513l2m`H2viZMA}|dE$h_c|uP?x>9Tc_6+PP}$k(oy`Xg&Tbb7wub9>L*2
z+m)L)E;pI&WKBzhzTA)}0|?VzX|i9>ZD-7Sl({SbF9}+UMc6sVN!GG#4YMAh{A(>8
zId^#usi|fqN#>X(JSdo7ETk<holBYO@Nr=iH~izkUt@_x6wZpSI+g+aiPdnt3$S7i
zYhv61`<RooWygKb#^RluOGvT4`&KXlKCDa0cO$3LlRNa>_kUAT&KfFnqAwu5t}9XL
zjon1P2qp!V*1}Ts6kzr|QmgA0o_~aSVgMJL7#p<G*uo^$`Y~pHl=!CQiHR}8jAyD$
ziQdPZPv0t*0dH^uAe>S8i1c^&o<6~pBm7ud<QP9MZfh8Nj#*-?fio7#FXUyiXtrQ=
zZtYJ7`TeO|GZ#ee)somiyq`@BC-WIaHL0{ag?}#xMu|s^9P|-cVpHQPUWwh2Jz|If
zG$3siJ#)aHzxmYB&B=H>N`xR^hb?9sd_*u^kAFP9^6S;zjB74&5xmGn9u1~RuBYqb
zCt?X~3oM~9IR_-=HAtj>bbK4|k&l{MA<ufjAnD^Tj3+iHyJINb5-RHE<+uM}n*w*!
zJ*t6PIJV}(DtLl2zRZ`PpI!9(-~-C&prM9t`Sx#h@iDrABXV$Yu@@|=jKh&M`@yW(
zAiPW&aX|MT;P~f$y9FH7lTHAdk?-6U9)d2)*vc29Q{Q>@E_jwQ`~Y*2JpNFZ0dG<U
zzLsX#z8l`#56ejaY-2u>?_PNE+mNA*LG)CS{C+{(0eFcru+xNW>fN=!f+eK1olbNm
z0d|~w=MDHeN%Td4MM!@6%Bi<t39YRrA46wv-hCQgCk_;hly`Jn^FQ}Iy!II0e&!R?
z|6H5?ZMwJV-KKM!zHPd;>Di`Zn|^J&wdvKSQ=2|*y0q!hrbC<lY`VKY`?a;7d8N%@
z2V(K|&66KK8ubr8^os@;E(t8kk}->^UH4inzA}h2k0xLLKwnoRG343Wkx<%#%6erm
z*`w?>bX5uXdOhBbXkW)<vB%pUP7d}k?o?B#E$Hzk`s2}ozRq^pg2&quXl-i>w0ixa
zwnl$*tKT1JZV7n9gUNW`KxgtxL!M}VqQm3$wE6-*&yF1xN?HQ`7Aol<zyKk>t;1AR
zeV3LQ^>}?1*xN$BCSo7#?Mx1Nvg!!3BWI(psWsRf@Q3`(O`%XQP)UYRliy!ihG0pC
zzOHB_i9yOCPtr{5MG8o%sUFqLq%Fhrpp@YO3{)ng14Eu1S})dER!tehDZ|u8G#%xz
zabz%HWuk$TS4}OerLlZnz{XI_q;qO&IGv%!u@j+CsKv%GoX+9<)pXI)^bwD@*+;1Y
zhi;hGC}Y7=4ozOQOhl%s0)`<U$1uFD69+?7wjK>9lBgk;*;{I`w{IXBucSD{jOla+
zhjg{6X&_*q|MXq-uXP-9;m4EB|1#^kx(1Ziq3RkS7UY>-16*kBvQamJD_ts|E);Mv
zO6lpm%Elg9#Pt2a)ljafxF)hS)deg0ZcWb6M@sAI)Rl9Z8#uK`@MZt+$v*$7w`W?<
zo^~PNZ1+pk`NwM*zf+c~qdu-0#C>_|i51V8n_RGrYZX0o)l~Hieu{EAtxz;j$tC5I
zc}>dYMuiJj@%NX1+fr#O^@&wFzW?CrmoJ^)Nc%^uXSR>xRd<R2+WWxDFagAdYQZ$_
zQ8SjBO_8Zo>JGWPliitOBe~GrW|c?ri;sqrg<tn<+8WDG)?Y5`8)h!6bf`vFwX}?q
zm?&CmULiZmh0aP@N94DAT0NoOxS&ozt>!hHboUr$Dm|_hr_<1;?<Kb(BytXxnKS|(
z+>R60&2*Y2@mM57Ag$%qk`*)w6?ujPb;;y&p`*uYAZ*P38aB*4TS3V3C)O4+k<i*x
z{PaUtuoX?BpnO_K`BQ+=5aH<obNtEr#jkB7VCp%-w#BJnmUriq!!n{8BWAiVs%7jW
z&_1EC)0Gq3S=~f_e&rqS{|r=B-tj)N$BhApvaR^iR3ImtK#yrE4&_L?kj`k<gquYU
zXs4eI;xS#bim(h11yIP2E)lw@cvOJG1s7N04<_O}(B5^XBhuH+>d2TME9IdDy3HSb
zaXns46nZ9ChthdYu=_ks&6n{VnJ-UfP*8&s$wEoCPq<T({IY$?H#J{AODfyvtzg`d
zViMI<MWBj66@e-ORRpRCR1v5mP(`4MKox;10#yX62rvXFkRgvW{XsbGkMJM@xbNpj
z`o*S3RxM^stzcnPiiL%?V*?bSNMH=k#Nd^oD+NQ-EyWmCtWj0LNM+8e<S{CiAImHE
zjiq(mvzwsOUHG33W{F41q;+M~z$nMvD#o0!VyVGKx0lXAl=f9JW&apKacw%hKq8dw
zx42QP8mzGh!y$OCX7WEDx;qm|jH&W4$x@^gI-sZwH9$#H>b(q;P<W&KnGou^?0LT&
zwg&}l{a8lH{S_&Jk0cm~#-an!-bggXBxu3cP^<+bm+7Kv0Z37@7}(P31U}SKId4>{
zD4{$9B=iLG;E#Nc{KWGOlB6V2j-$8;W=bOYCP}AC^w`d4p6AEBe&}?pAz5bjn3o)~
zc@Z~~L-v1q(96APMgo90)6gsZFHR`{NXKqyNN&|~v;nR3)7+Y|_#egL+&e#Ze&bmn
zCP0s-UttDo`s5@Y)3Z3IcVR$vR4L?et{2(7)(wT?w^*c(;|^5D;*<YVpV$;vq`@mE
eIxPlu>5~oR6TA6So7gA}OLAh%=dTto`Th-R26mhP

literal 0
HcmV?d00001

diff --git a/DEV_inputs_test/weather.nc b/DEV_inputs_test/weather.nc
new file mode 100644
index 0000000000000000000000000000000000000000..0d7f57fe048eeab4f18a23da0b727ed72d6e8ceb
GIT binary patch
literal 167106
zcmeI52VfM{)`0KKY|EBx8UYk!jer;_A%T!!1430PLO=vHhGZd;G?UP*fT)PrQPgMm
zSx_v11?;{1)F(<&vG?v%Y$*RZ_nx~uJJ|q%EZLcV=YyQ=H|LhUQ+IZ<d;Z9w!+q|3
z+^WmP<g}fpr<bZ+<<N}3uReF)*pWlWtG>fjHH%A4mXc+R(XBECGgDIV=Q$ewWNd^8
z%zGgIqA@RXF&kq=wKcPYHI0?EH4Q8sW;ozKcA;o^ufDwpcqcHZE0kPW)6iJiI196R
zaxNtgn^j+17c8v_6hT$>!NBy|`oQR5dF8C?UgJw^W&{Jr21^?&g7tynwe`WW(uT$c
zZ-D~MXwV1Vi;BvI#@hOMywZK4^VoyjoZP%#IRkq2%_;7iH?Ut`V8r-QLjsMpfuf<K
z4=x%Wm|a;KD4RZ`zH(ZxoZf{%-#)>zqT2FcL$BQ4xdVFV1m;vWRs>GVYN(x6Ulz<7
z7|7~VGrgftS#5RQtj1ukoSdAzKINs2r8#{{%S-E^{od83^$oo%YZ`;~HKkR(`sNoD
z6c*<9E9?u^_3xY4H@8=AUVgv+Fr!z$z6Axj1M&)c<>u$)_Ucy}%<olDkdp_$O`n!C
z0A@}rAJ7~7men(mHNCP5sw|l{uViL$UITRGv@^g4Ys#RrC1uk=OwFn*tuL((LK+rm
ztgH@VoDr;@QD0hDF%KWi4%RnTmX%gzop6SCe_VnZbBfUuPmm9uT`rNL(Z>$w3zM;7
z_<PwQFBOXwi7843WHD7ylp4OGFeans#Z@N-)-bV>4Py3WG+H)37ON(kFIKEqX~p^$
zD}$+BCtRNxte(adQa<($#ypK>z{&?*C>axMEGwVhuc)@RzPz%gv@zJgo5org3-As?
z+kBD{A9xQIE2Tp``=X-<Pnf`;Jg)E9$DyM}4;wRK)Yvg46N&~OJ&Z9o77I^zZOxuF
zwrJ#0vd8vjdn}&mvpA|S9U%g?j?)`G|9%gv!6XNBoFMGYl;rJ#oza+G^m5GL(ZeE^
z5H4f>RB%Bi!-&AsmLoE>vbwT7ax7%6g$xokXpe6<vuY|E#df<NB;a!|FP+yAXsE0y
z3kF~V%ITHU7ybouat88$-o6Yvq_HGX9M+~_O?hcO?-*|VSVLWXZB?)iR<)8D^}!jn
z^_8VHusisE#UD*qA33fvY8f~*sevP)dYlHaeje6uJFF=B)rikxzn>)rSHnFHOhZHn
ztwzpf7fLIvM!2Y2z8c}^z%*^udG)y59qsXvB34)U&4!C~9wmt<H^+%7j2$Wh|G9SJ
zoE~zuVQkevG3yx0@OFy6+F-T#e;)k*B;~I`ftuRJKxv>R*jO|a_JyUen5I<)dk1Lh
zo5NOQI5UQiA32UM3zhMMfoMc-kE_HQ%MPGk^6eK>LQhHLCyB7rpHM7RQdHdo)`zS`
zG$ib$&7&s|DIPf#RzFZasIR<7)<Q8i38a{QxZUyk@ClH`WdEp}wAqEbVPoUiT-F1G
z$k<NLV_!;<@z$k3&Hte?Jm&j~wUdt}raJikP_faP$dtJ@zGy`?*-tpqc&}a9>r{}?
z>(h<BKJ>9SZ)oz|&NODd({19_!zWA_78-H|p5s`aweZq2dSrzAn9@$PK3oLmg%t!h
z{A97^mW(MHhsRtvB*0!1H{IUwA1i{xL<1b356kTd>jvyz59<eq#=KtT!P&v8+PZ2u
zK=-UDpWgqlzP<axVG?Wq<l?^yFybVKT67$D)?$*&mG*>mWMOaw5)+MRXd%iyf`{gN
zDdIxMDOm*OrzqD>9A2E#MBvX1f8|REB8>{mr1^3=b)d%`Degm2t^%-C@u@RK++E@X
zJ83=zDDF=&i{b$kyHe~%u{*^BDdI|q(?KGzY>IHG;!{tGy(r?+!l{o4EQexWihRMq
z?|JmNAI1I@^C=FXC>J0c9{5Ci6?O<cmOG%s=<(qc<soelJswPP2*o0bur=}t4yJsP
ztK|rKJd)xlibqmBisI1}M^hX_aV*7g6k#pr(|C$-{>i6eiW4bLqIfLD$rO*HIECU=
zipNttf#Qi2ODLX1v6SL8ie(hbDF!J{r#ORR1;t8=CsUkBv5I0f#TtsW6zeFSLb0A=
z1I0#)vnbA{IEUg~ig<Xy=~NNeX%tVVcm~BYDV{}fKE(wT7g9W%;v$OYP&}97c@!5@
zJfGqciWg8^O7TL97g4;J;w2O>rFa>|%PB6Scm>5PDK4jY6~(J5UPJL;6tAUt9mVS@
z-azq2iZ@Zbnc^)J|4s2$inmd`o#Gu7S5Ul@;$0M1QoNhuJrwVycpt@86jxJRL-Br!
z4^VuN;zJbIQd~#zVTzAXe3as26d$Mf1jQ#QK1K0qiqBAdmf~|1pQrc&#TO~QMDb;c
zuTXrI;%gNDL-BQrZ%}-b;(sZ=MR7gFw<*3uaRbGT6gN@aOz~Zc?@@f8;ueY@Q2dbM
zM-)G%_zA^NDSk%rbBbS3{F35V6u+kU4aILMen;_pia${Nk>XDjf2R09ioa0&mEu;4
zzft_1;vW?Mr1%%bZ502e_z%VH6y<?mQG`$_Y7}i0?GzmpofKUZlPJ0=dMJ7+`Y0w-
z^ixcsm`X8?VmiePitQ+7QfyDL1I3OM_o0Xv6L9K8u`|W}D0ZP3ptwK9EQ$wE>`Jj4
z#qJajq}YSvK@_to9!#+(#a<M9Q|v=AhhkrfxfJs#_M_OJVm`$I6yb6Tp9(1sq<9F$
zLn$6c@o<VqP#i>YFvTGhizp7IIE><OiX$kFq&SM=kra=jcr?Y)6vt2;OK}{<V<?WN
zIDuj@#fcOrQ9PF7WQxa8oI-Ib#p5Z`i*$A(Juabm62(%A(<qivET<TxIGy4QiWL+q
zDV|JmCdDd>)f7XorER#_RQSJZmtI_Hn6D`-vN6ibG~tUcX>49cF~_ryh<zzCd59G{
zPJYgud&<i2i=G!UMYj&4-Z}QeUlf~0k39`_rHz%PRVDSo>9!1#5Ze85!TbDW)v@$s
zl~l%l;H4DU(`Pv`RcTb%QaAo~_f}u0;Uh;435^v>wxT-^Y(*AD-Z%?ghPdv%x*tM|
ztsmy7TKEbYuNa(+HMF4RV;N~NV`;1cUY6TkxbxMr@$vTUasj`62%pV<q&=3rR+Yqx
zXwDk`tQgxwkxdsRTIr{_fFWN9ze`UCy9(29OGY<@7u%b(wH;1PJKEtdijA+^@SECL
z_Xt;Fe$s(!{bC@O>n=sg0g|!h*OKhaPn-MAqM^NWEeV%fi?1c^>@v^XmR%a;eHbdP
z5P0ClGHe?q!z-(*CY4sr3P!ykhgW|imnn@IT$p2QWydpa;TzN)w8K<hQ&=9h!t(zf
z@VZ>|cDBke;@(#5#2ZKYVB=*0ej%8*4Uc4$E0lZ=&(WS=KA!j3W`3?a^ZP^PxhdQ6
z)v(bHe0q$@eMsnvnOAtYM9U|>2PU<n0BgV80BEbht6}HRhM=&0pZ|S`@w(Cq+JB{m
zE`an7-f)ufy3&eHvL~`Hbw}QEkGN)5Mgt^wnKH_r8n2m6qW<)u1R1CP{Ih>Z@zKNK
za??Io$I;?i&{20x4zTGcU=N5-@0G4U?a%6P=%z4!HRr)Ua&y>6b`T@PHNK(8Uv(sV
z6mwy>o&@3MjD|664vJyo0QRG8_(awN>li7n4876e%}H!AK9f=6lGKMM4W7#Ku+n20
z^FY{e{Zl8gOspeUyk~XAxkuKrm$8n?;v(JR^iQv4Z&Do?Zux1=b?h4K`b4}4&e*%Y
zf*aXLtOK4j?nFyI>GmAUbir@4#r4t?KKk$__6KGR6npk(9^UYdxb{?!RWi2az*9G}
zCDbuoL7tx6|L?+?AMoP?&oJ>1zWkFSo^Qlbi5T99)g)pSi5T>Vg%&aT5vxSR3KOvi
zB36Tlp^sSWBG#vf6*OWAN31gcdG}gac@mQxf#D;E>dujHyZG?JYrYA+8>~AS<JV~?
z!X|-R#klcfk1U!nHS5?B69NVO0(}A#i$@1et*r?L@^Z#x^&By5?1=E}te!&$7f&3`
zvrGCG7WB*NIbqz$VdKY+8Y*kZ>Y1N6prCKwfS$R91-*0o56I2U>z|+3bMVCCkz>b?
zDn4#%*06CCMr8HO8jzEhlXb!gQRd|5<>up@apU2=L1=A6eK1%vr?RXft7lG>^n${i
zzL-97%&6k2S>-`^Ck)<@&FP)fcR;`XdAWtT{reUc_REXbhr+(OxzYR3FVu&zM-D40
zhW96@W);_$)->QH`(VlFV0~FBybGNLy&OCq-YzU2HhyYW6`b%k&Vu)eYp0jMdvP-=
zYoH$y8D;pMJDd|&2kR@#;mx|Np1txS5*qMrzLM#sWmq{pi9&c&J0f9vX%(~|g!knu
zYi4Bi><{mZ=R~$4GHPn;8!LDU%oRfeZ&}yF$n=dehJ`s|3<pmh1uMfaS^KcT6N+I8
zaXBVb;h3@Gi$_LVa4=?knaXuC5%=B*L~j3h7`!9@_2u~fJ%97Lt{&d+46Sx0HE{m|
z-w)t19`^&GZ7ftm-p>H<Etl0dgiGLI8%u~ih6YC*aaRbe@UY)#gq0ia)0kBqE`X;=
zyg+SzRe4}U5N@P^H>?8_<~1}1s{^<j(V#2T%i5EJWsSk|$R_c$C)}iAC#AvZ7cUHV
zaFmx+m)6z6S|4sboyx*yCpTW$`1r=FvHis9w|={Ychs>9@U9$f!Yw=_$5Ss^ZEekr
z(Eh-<<+rd}*OT~mB7d<Eu?J{THS-tkH&}xcKV3C)H_@VcU`XLv?l0=m@3Bfs8%uWk
ziO#+D<^$Gk`{4);-fiX!6jl|!%Y{3Xs`<*APs=jDFfX@XVeg!R!u;I+{Ri~x2e(ha
zEkyXg-+-Ll0=NmNAUCHU9DEA%x$TAXvDem}BDa5hp40aE_3+S5Teu0&nc0{)_J__t
z632cz5o_IJzs6E%zYwghgx&GUwe=;Xb1NHKvV$3PTk&MrQAJc#345xRRrJR4P8rW0
zffc2T9#z8KlX&tr9rpK)!5TPB>~g0TS2teS`QP<}G}ck<o8Y2mupaITgj2Dqs>-^C
zS~&KE4ywjOWGB&nmHv1Tc|hLfj_h;?Hlk>(xPG4^Qsp%kd_{nhyyppzr#KnCyaplt
zT0E72i9(Znd^wGiF(1bOC+TOBkQYFum;I5CWt1LIP;io7_D5cNlwRiJ1qYm@m-+I-
zko2-YK9L^Em-QqYysXD>@G_r<J5-O%hl|^MlI_ZTc{NCSnUCiyI7u%r>cUk*K1na@
zp-V!cd|3}&oez0gPy1%->0t0OpWdMe)g#-bi(erx^Y=B>BlGDpSSVliv$G*z)<YNc
zL;13vE{1$rPr%@1J*3%>as`kg5fdYj7=gqHBt{@H0*Mhwj6h-p5+jfpfy4-yd<5j}
z9`Noqy`(#to)O;nvde$JCN3K+|LIWq`gQk)Wkb#B>vw!72i&W|CwcQq=t6_M;6M$_
z>i|3J2NVmH6jk%TYvMO3VA*lqKKu91>PXDryCYzSTRg0Q4OTrXfPd(11@`V*oM`R;
z%LufCKW(%ENlwH?IQe8o_y1+ROFYlS2<*-XC~z;c6~I?7R^S%rEzW;CZ*|`0yxn;R
zu)=vKn6N_|3wLvrHd^jPud#p6+pNI~@15SeyeqwTd++hy3*6^j<=vUXOk>{eJJE7v
zeJyvpvu}vUGBsX{(6JV%>v&4XdZ3|WBTy4FhpiafXdU-${*U>;<o}lcSN`^Vb%0}l
zdw_2My8kun{)v(&ToNz)1C?#+{>*RA*#}_RZ}LT1k*9ZC^P6+Z0L=3a=mH-hf&_G!
z(QRU3HZnojY$CH!m(}`C)~Ii{r7q8r%jxS7<T!F1W^`{T*5FVZa#)DNZHIGXHuntE
z^E5;B<{56{-(9R(gG1pn)vUl2#}vnLz+~W9$0T4PnB8|%qm<nnE*{VMcBfBPcukAX
zp~xca&N57tNsPeGMxZ0yb#DcB_B<2w_wEQ3Gi%_5kGrq}_o(-%_X78+tAN$&8o=n#
z?uonG&<57m+D?VPR<#0qcWX?vwzCm%!zWf*ft@{1v&^^3>kqSxccPz(5on7Mu*!CR
zE)$%?!kEz|YXYAZ#*FR*^+N+6scq+sEn{DH&7IAx!S9OX4@Dx)?;-O?D66&HZ=o7x
ztA8ka;uM==Q#qPVU<0!^of@Uxmb>c;ZYsE~;I4vu3+^viTku%H(*@5LH1F;V<vj&?
z*YdLH*6^}Ba<>-9yhjUU9=dyY+2-7H(8i4gGVh@R%#&^fFDu<~nNv7ZGX=+I9-rA#
zHzi8hz3sZ&j2zg+bhpL7OYx`r(|~k;hQA$<>2L4XJG3ozSR>Y_FB@xjvHxQKMZksr
zrN9Ni5-@wyVavvLTkg^{YtRQioYV^ZPu26Yis;S%xQTx^gIR;Nc+C@wHVc@ZLWbx~
z0XOmQ9!VBFoGe*K=8=$jI60boI7(jIxmYaYwy5*eBxD|EZW1zw8?(8EC1;hKRx+og
zv81k~s-&_cSW;Tj)J=<U)jUtSxxB1&XY=~PZhFbIlCV3a1oNajy#(`G;!chzJBybU
z?g;xJ`(XPJdy##peVBc?eT03a-Hh%DHO!N$AZ(r>^Q0~->>k%yS+=d572Y$z=U|oV
zaQKY@R^W3D`GWXTL(J%IaNg*=$%%gn6aE{V9Xb8q%*edfb)q-&Z-hh}J9jQe65s|r
z9IwOcXdNeNjI1wq?y=OzQ=dqEGWDs{r&FIveKz&E)aO&p=>Amy0=B7ttN*Cm;rpc(
zP1T}1w2^s#sPW@`@QKV;;4EehtZ=UnK8L@Fe^?kZx{otI0X_vj%lsVp0{9r%n+{tx
zw%bzIy<OLKS?%_3*QH(OcKf#L*sguMcJ0s|5aBYLaowWiMY;~{uxwLzV1(PJ9p>$h
z>(ow^eZcjQYn|&+*W<1yT~E87bG_ht$<@4jgz}!%<+a>BA8O+f*JHfhXI<#-=Y827
z_oC}yCg{n+#L68UrH$5dZE;N?(Pr}_@D!X~J(l`N>cgpPQy)ydKXpy&Dmd3d_vjAX
z15xr??j8@d5iTp;L#Zuuuee_ZUUI((yx@Kwc+S0d9JXR?w`C4X%RUcQBtsO5G=oEC
zNGK~_-LFirm4z{*JBEcB#{^+>44HAdtk!q3Mt!?2br0Ddv^`+E-?qlK+IFAqUfVsk
zyKSxIRz<0A*WG<h`mzdoXzZVKy;+!EOb|A`$@J1?wZ4-z>f3Fpvnh7Mp-9@^*`xGj
z$8KxtZ>hhh{sH`%`WLV*_3u>S(C!_Ftr**Fsms^$wE;i@r%)@@27+nciOk5ncyO1f
zmjagomvfe>%hVk^kr|oS_U>Gl;5=8D9l3L(<i&$KGwqDD(}B~{P6g(r%>~TpR;8~=
zzd!xK^oP>dr9YhhX!>L6Po%e<TOFk@=pKOnt%d$QLj5zk$I~Y$6O>|LqB03M7MKiX
zM-DTMd9t*=taQif>WiKGzVr8;e*k{${1fmq@INrE;jl(yeX(<YFyv3-FNQSlWTrkZ
zKHQC}HCU@i)+rKc)`rZwP!`T>n%DO3c*pUM6NnQX0yDbx-nYH)csF=AdN+ADd*AiG
z=iTDnRrgktvhRC$-K{sYi)DA;9o)HB=iZ(Bbk6DAw{u?Sex3Vw&hNZ4mmA9K87kZQ
zE|2=6ZzJq_QJ-_6*VyNawTplYfu+C&z!LammRP$B3~h*(MZ1x`!B2@X_9mK}#jJr}
zO#xDYG$0+wP}|0(M`?FgoapsK?zO-=&cp79-H!lfbW@o%X#R?+x#pxYYv5qkzzUyP
zi_gI-?>k?_j9l2nT*Ttv&G#W^MPNpEQN?8yS5{n8aec+j6}MH~S#fv8s)|UrqKUh;
z!su?QFuLn1jPB|Rqr0NQ=q{}=x(h3e?!F46TUl{mMTMu5bF$}T&rF~SXibMT8tW5I
zSlRBLx`FmWqQEXOqq{i!qU;N^mu6p(eSY@h?DMkE%|0hPc5X?ijY~r9wvJ0qvj&H_
z4|flC4|R`lAL$<L9_OCmp5!*8>*VR|+0WC(6Y%Ws$?_cF>FVj`F{3lrS#9&ZwaR_r
zRyj|!##y_=68Gri-x*kC>dZQ`{S27VS>^LyM<yScd{px2<gv-eB#%!nPM(xJIeE|A
zD5`Iqu0EqX7TVZ-H-gHJNtSIKoBXX}4MwMoO*tlId`fZ3q?E}iQ&NskIWeWFJAu1#
z5qZ1jrlw%tm=t~4U315!VBXOYWt(>wC`%1osF=|OnKkII^jGo??5Q)>i1o$V)!J)_
zYCBQ|?3t6L_4UQx{VVgv%ztOD$h<pqb>>5vk7PcX`D~`%Js5H~gz_TY9lR{MCwN)R
zi*(myO1COg=0&=Pd0BL~@UobP?wL&KR%S+b4^ZBR3i6STkGY8t_feF-M0Y*A_3YNG
zTW_FGw;Z4^usaS*8_Twy!yfkGJq))I>8iA9twyWUPSNVM2CY$>rOnp%+|@>@Z`WO2
zlfKkK4~_kcbX}P>Xz<VS&+(t?Kiz+(f4={0|2h8i{LQ<0l($fq*K)Tg)W$sjX}sNq
zesuM`FT3N;^*7iXfmy(8V2*t*5HAjUW9*A?Xv2L*w}dTW7Z})6XRHzHi?y2!_rDz%
z!4xp&bJ@%qY;|pQ{pR`|_yhRU^%u~Z4r>(kVJo7IShzf!o_;n&Z}PZ_f44x5WTCpy
zU>1PceTStZTfj%H@3t$*KM};nMYbZ^Xg$~T`@6(9ctw7vSIb;meEZia&QqNA&IV_r
zbCz?qbB=Saa~_zcu5OpyTxZj~-FMg<V|~pzJXTnN%UzedmbtC~u5>L2u5w)s#(eH7
z`_;fT_J7&0wO?ny-hPAqMtkeJ8=yX{QQt=EI_&izb&z_50W&(QTyM9^{XSM%9zB&_
zL~jM@1MHcTrS<j2-o2q(gLSIlVKs~yosC(8f2o3N)i7psKPf*G|5GG8c0WhS+uN>b
z59hh&xMsT=UG=V0T(z!hSC#8zSL?Z&i253M8#{7yd0!&)cFoysh?wLj>hF%T(Wm0;
z%kE8gxON0ENE-|c(Taef+Rhzj8uMc7u!lo<55w*1ofW=Y4WC0|imwr{%IBoU*QoA(
zZ~V*}%(PpB&uw1-Ujkp*z6QPl^bT!2IqZ$Gf9NJNYhZUc98UOzSNx1u`Dv`J=debv
zBeK5Ma@cG6Nwr?bRoYeBa_vgs3SgObInbI8Yc$r^a)&(>{R@|kboa7**?q(+CNQHL
z&8)%hA9)fb6C=?42w3B4iH)6T5f-sUp>MR5-)Ct?w}M%NToxve3Bo3qOr9<)K3(&B
zfW(}{2<-j{Smp1dzE{3izT?QhpZErh`CL~;&`m+Qb9U7o5T$HeIr&+&<Ji~`W^}v1
zWyM!XJF^DyHL8i;Bt{@H0*Mi5egwL(FaahAn=WJmx~%wg3zdb+*}x*sISLXhCt8VY
zqqW@Ynvpk}nAf%VcVDtFUok=0d`adjT~>TLm(68M;<#;y2QZ^cSKFzXYJ0V#x{tcA
z+F9LC?E;?#xhF0oN`1n0h-hPHt`qe|pVz#*OOdQpB+}d!GAlz_IIn45+q=!$CSs!|
z*r4qlCtK0CyZcUjgYfL6SxJpa4N3J$rzF+Uw+ipv&Ewxb9GTZv?g_;j3}?gHhzLf4
zF`u)_=M*)6=9QR}7=hg%fyDN)`(u|VnHT}HjKCT8Gwt*33+xN+i|psv&$XXtKi|G%
zcXpJ#9l6E2`ea`7?kwJy$h;l7MZC|8?O$tOYhM{KqdQeORXHt!)4`a}nPvS?^fNI6
zi4jPQfRkB+#Cdyst;^@ClJiuFH0OrQd7-R$bv>9hNN3g{zSi<YZ>%x`$7t5zU}g>E
zU)`m!FlKZU^2X*JnKwLdaNc2g1$lXSee$yNB3(c3iX-yS9h`^ms8HERH-y*M)D`CS
z%gfP~#k{Z^2lb7}i|!7Ol83H8ukXU5ON*{3y1MAPqMM3tEm~1@chRb%_;XhkU0ZZx
z(Z7rCC|X%`U(wQ{ONy2iMY=bcHJHJy!BNZ_sEVp+giR5klh!Ow-0g#R{DklNF_#l}
z4AuCp!K-fM72;(#Vn)}VS%VJD8d&AL-74381~6-I@2EAS9vt=XsK-V<HR{<>FOGU;
z)a#>~y8js53y`OGPeNI}TL<;&-D+rK$L?ibUpVg>s88=6h5GdF0jN*!zR36@<I9Y%
zfUkjXGQI^`(_xLq`daR=hoXPsvXSo64i|U0u)_r%&hKzuhjTid-C;q8vpTeuyDUmy
zF6(efhl>pT6K-ilpKoH;z#5<WF`3<A5v+0??8~fyk6DB7nKjtR;+qXo`WMwH8<cl|
zw}JJXx4@XkiC({@yh^;Hgo%}VHA)+;<t8v|V3p%X;&~)SAQmGqp~HCK7+@SQw!;`8
zUL5wu*cai@hWm_e-^|XL`(*|)voa6J?3Q_8W{=G5%w2b#qLgjk?H|#G$P=zRwYzI>
zpUkFtpQIz70-td{Pyd`7qx(2YUOc&h8HZ#X3LKVkc*YUHpp3y8dWW{94r|2v^krl1
zVzJG&MV(o`12|oMNH>ld-3fMUuz)RK3k{gj$<OXGdu)qsG8S!G<$8N8#wJnc?vH>q
z?m~95<aB>}3Xqzf2BfEF0DH?}%f@zF=CHKr^OwGFeLwhq_WkPn-M7uR-KQqolfU$R
zvqSfnF7JCP8{PdAs_%Q>PrQx4e7o)x>I;_r1^S0&(R~5^i*&Z+6BTPP%~j?qcLiP3
zT{B!2u1eQT_~sF#J9&@X%qaEkx+^vG?_}s<G&gkMkb#2+9x?E+frkzpIIv*gfPwu7
zq8r>Qt}se{J9mfkzC`9lx+2~`z1#ggU3`_io*~!ixR#sPy6d8}Azft4qJKMgWvZZD
z6$I5FnKClvU}EpE71<Yk%hN(O+{0+j8h6RGkFS+H(Hr)QVhyfPtU*U+4HB;h<7?-V
z`2MK)S{{ksBt{@H0?m(rRnB*dnKgLdF4$rhd|>~8%==`vfQh}sR%BoFEx#AC;T}eF
zKQU|2{FalL(-tFOjjg<x?QaoIV%DI4YX8*y)B!+2YGLX?;1Dov<N8Iki@iA{bywY=
zsx|nR4Y`JMwM}3~_i5KpyMD&`yeskrfKJ$Wa-T)%U)W`P590Lj_V6AEbO-j#iPDkv
z#oqnOtie&)*5LDO$!FOTX+95`&q7)8>Kq*qrw$i4@!_0N`V!q)<$B`5ww?xh+05t`
zr=6d6LE6%^i_$JmyEN^xv}I{mr0uFZFG^W-OAPhx$X%#wSLQY6E}{PE^Q3Ew*9wU?
ztv><}Y7c4;7%-!&W${f-l>S9?%a}FLm^B#2tU(d81}8IXki@LPAVo4*kw`NrWCn+_
z;?<4rJi7B3U@R~WI0i7IOS7lh(}@hbz>MxfRj^bIV@9{sX${7X96Db0F-7Tz=ftX}
zC^d>WJ9g7avW(4n&C8%f3YO+2QkdO^M$5*>+qb*YP-{x=iK5~kJB!$^6SLurBsSB5
z8I1jTf|&C+J!aE|Z>66p*hl7{DjNAs1cqIO>F@4O3>{lE@+e0pK2mIU3@~ZKigTVj
zPUJW;&?q)12AGt2Yw;Pt^mEDlC$Y4Au31aPSd0_1nfG+Y9B}63XE_Wi89!|Ju<^sj
z6b&mmdenqsmJi9;xAM})(uQCogO4O*%!Q(|I#|<CSzE(gIA1gtV5*`hqJBB%j6Kp0
z5@3>E_E|pudULmczK8W~)<Ze&vWN2VutleSM^B({PEL-urx-FXj08OKvfA4E^2(ah
z#$W?{Vl-nON^PhsZLBP<Dya`n_f|3Z(0Cl=6ddHCqec%KGhx)&F=CKALIE5^Ikq^6
z8a&2HM(oO}Y_1rt_WH-zp#HHJ9^Y`v_B#EbF2jWeQynJ5go|DDDAtQzw8xknRQdSM
ztsRdXK4Iig#$Kawra-H`%)^yy9GeRtDlVQSW4Ct@^Y*0_8Rws6zkNtC)ab%>EfwaL
z_sE(rRx>>H5!0C$*F7p)k}Kc?cTKzPYDI9362^?~YsDIjQLVvK>XQbZP|fJpDEAw9
zKry4cQPFdgg6PeS+{C}@pz7%eU#8nJ!gSy!{@pREH8@<=bC`<g&Eed{zxzzp^QnsH
z&1c-izk5ru2A3<=;3vfze5`(~eq_Lmu0SnN2N*D;TdXeDeKnH!+9Wf&)ry`~3Zgfw
zxru+5t?D^QMf4_{oA`GRtCDr9M4E>~W?d)?=QYi1duPVi8qxavr=sUC1<{*7xru+*
zNwo&+6v0|0j2YcY>Pc#efrKl8QE8iT>8@IXZOS&<BZ++#{xqYjRINe2sz;pp>P<eK
z3&*?5Q?0>l)f$|x>N!J2^yYMK;@@@A1OY9K8QrVetJ-S@%;-+DS%d3UYw&>WegkW4
zW^~(BJ%6c)-fZJ0{+%_Bf_dr{#4<IE8Qsg8HE7iIG-!z4G;$OFPMnFD|K~i_8rVXU
zW|M?e$Y@~&zctNkd-q@MzuKDy%;?Uy=~-+;^yYkS;@_R2S%a@_lCNwMX}%7buR>Wk
zuW4S}yAC!<dz(a>4k6P%l!fz}=C!@MO0x!^YoBYM88D;kV-xhYg)yUh(Wd7G8=^Na
zaufe9%dTgCJEAvP+{C}zXtM@uY}TMevj&f8)&Oq6*S4#Xn9;qcy{Nrlz>Lm#@0D2}
z-=^NC-fF;%?m|s)kru{`Zl$W{E)~(6mE6R?>#pfJP($>lJ2&y~-mtx4d((g!U4>oG
z3_GGX72G7;-Wr!@9fD^Zf~Or%lX;fRGhpK1eQnqCl^xNWuephT*T-%RMr%hKI7&04
zTVO}#6KB~GGrBc)WHqtMj+oJHP_4l*%^Emt*5GNoo~P`H-aO4s{JTROf`N`Ov2uq*
zX`{8=NA?d5d|)@DE4S$>vmts@&Q1KgQk&o;TNpFC`!qfGYKY$4$4&e@^WVuJp24&B
zXYJ1ra(`<ECwndXY{vI&qV@SVn_#Ogj2Ybmr=IyvL~j;w6aUU*cLPZrm)&J|0%mk+
znx1qG(VH}G;@=I`^oVnKy%|d9_3`d*u-#z0o^zc|a4i_~IrywD$9B8mAG^Q|PWD>%
z`4*?(W@ngKxm%*N(OPbzYoe=|Gr=Vo55{~h)ou;$u}SW>Nu;?aWbO`S;k>4KZSRhA
zOg3<=!;Ef<Q*fM9Fxfeo%oH-mfr)=N*D09e6wG$cCNr1J95C_kRytNX?lNFTx6HN7
zbvbdFOJD{kdoBCC)gk!RA^64d3z@BCegzZ%uE=f;7C9tmJ0#L93YoJ*Svaq0Ufa8$
z9g?3M5@~)8nV&*gIIn45+q<z&J!71R-i+lY{@q&+Yp}p(4bFBf<ScL?^Eqa8@2Y~$
zY8W%RDUK;{ek+h?zGiT;*Rs#Kc585%Q*x<OBF$wXb7?3G=QYi1d-sClc>~Wm%;=Ul
z^qlWN^kxY+@$c}?5^+~a_@^1&OAc%BsYA~v4n%K0<tF~!OougC;kd)V?G7`#rAc}&
zNJ8{xDL3)&Mkfi5P6}g2cZL@^9XO41s`pfG;>DdArGL>}dzT>76~>J26_=iuU5MVi
z!cF|Uvpowr3p~htjv3t#-XFZ*1K$DPa=ro6I!@FWSzoN&LGFXxJ%9s&?woF5TE~eR
zBkPNm`(3pL-`geM*(K6^A2Q#CvT$D0yta1}9fD#<7&E$?9XC5}GGIn`r`;M1cIX-8
zK=fuXH}UVjaed?Zn)8)Q@Ff`YxsC3P?hV8{Zh;w`?6vH32Tyw<(}T1F_RPuB`ubw;
zws<7(dnD3q37PjpSvaq0Ufa80-GX1-VPfTejnYPIxg2+nyN>}gy3KCMCbvYI%^|ZX
zl!fz}=C!>W<CcteOQabSGNVITIIn45+q;XLdM<P#dUFvs@$W`D^^9;LdNY!n_;**j
zBv-g3(p(ubSA?=~Uemm`cl)^{om~=X_6wQLp)8!&G_UR5T=yJewi}rR?3t6L_4UQx
zjdde(4<+_QW^{#~LQer_fJcxI#(eG}_e1UnIS;r6_k%H?lXoA>dyvPw#?$@Cat@mK
zcYnJC+gxGH=*BqojCLS;GlrY^cST8pAxUA(=#Fy>rntkH(S7dv-1V6OGrA94f-SBv
zW_0V^g0=22v2yF8w9#5_iF=9rd;?~5uO>-eNs>tOYRJ42%EEa~^V;4WV;79GhcTnu
z>XiKIlt{BRWPS~0;k>4KZSS@@_59^T^ky43@$c3<^}OXo^kzLb@$WVz>Dics=*=c>
z;@@Sv4kCKEkOP4|bF#F)zSz5fThPTV*w4KmnE;tCVB+7I@o(m%_4x?r;RX(Kn$Z;|
z=_yD;^rn!T_;=rUzvg`9MZV;i(S7Rq)bk1FW6#H)kAM%sw7wHHMz+zKuGWXt5Y;}!
zjP4sRg6mdz&59@YZIu2+bL+foIS+Y}2RUYRy?woXy*NF6J$(lQ*<f1Vi5er@Xic}T
z+ZrrM5}cnTSe&$&%n~x^gNc9F#e?hzbmnyObmAsnTp&vSqB;Hb3US5cA-IAPzII|3
zC*LiQ?-|_hzMtMdkaN()zkAXB0`a^Xc@Ee!Crj(=i@iI|cbe~1U>-1+GY3rTI8kF{
zeX(-#==Xs4eqar-I&?&g?$8?ZWTw8XbZd0=#m?1v>%6s`8gGra8mI!(`cBjs*+y%+
z_Wt(%Oinw$AOno~T#2v5ccSkE;CSCuV2bZJFm2;biqdZL4tpK$^W{FlWxg=6a+gPG
zqqSV6M^A+Z(VI$c;@_3H1t+@0#LAUKX`{8=WcOtEvA`r?B2WyN(JgoDxzdg3&2nzy
z-(`B+0U6$OZyJ#5P4W7@@#U~L*gx!{?4Qx4cyT+P0ZjJ<J>@`|XBuqhyY4X4ST<G;
zdno%S+mNo4Z(m>^Uq_$=&>o1D!&Z!Kw3eIbo9HX{O#sIGjseE`#)4@ZHz`WH%{%OM
zxX=Id{p|Y*_|f-+?|a}o-?zS&I4o@}8!LxBl>L)!Nauyg!*P4v-Xy?<s@0vSF|v);
zba*y$Dd!U3CBBP+i}=|?Ydg_OWE-vNI{P~T`vUv;I|3a5y+a#M4trzlAG#0y9{^kY
z@B7~a-t}+xZ;CUAJ;eTD4`u(1ZnJ+{aw%{UPy(C?oRHjh4qG<%W!K!3-X}PZdmr~c
z#(C6hCMSCRsP{VmwZOmp*8o@huktSkqC2$4yx2PI&2_vt;dYyIeor#s<KWH~Co{S$
zJ$jaV5WTsQoA`HqQhTTNO6{3?aB6nyK|qhx15;z^uoYt)t>H3!f_6SZrZ1CB2AOtX
z;@{nwtY<|sqBnPP6aOwd{h;(7=?A8FPw$rAHT{6}tn~fUcg<x*$=iL0EgRcyiJR@8
z?VklS0u4YtU`FTjdHGJ)?Mw2x0H@F4Yl*|s#<H<;*hASr*@kp)```Aj_rC@F*Z(H)
zhW~XiZR6gF(r)t(dmZlcyD9IcYz8)^Yy>s{?||8n!%SnIEUhmq-3DEKv2*`Pel7Xc
z<X4hkPJSu*#pD-~pHF@+*^KTSw?*)u<6q=I8(8RH0L=HF1*UD>xl!6}-eIr9ect5T
z<l6{r@Vx`P4Xg*VBZryBJXu;_R=T%!^~KKJ;k})>&5PU$?3t6L_4UQxiC2o7{2Ton
z{Bj@HoWs|=#=NFZl#a}6>f~!wecqn9RsQ>cd;RwScl%fR?*gJbw8p&HI_%9V-kWf{
z%{jRg<WU8WgRysSL}^1haSX;I+#7)C%g&wZRXEJcyxTp<KOQqV(d&OaM<pRg5~Gq3
zGrBsrpw=y@ao3QkBU1||{@rTNYR@XpeV+R~_X78TX?-VZjBKMdUAiaDK&r=#F7mH|
zBmX8?JTiyV_tnh5hu?X^_jbn7`%&^8s@-$qoubIHJ9hGYuwwq+*m&LhOxVA0+2-@`
zZCX#7J1q%tr8xlyzH!^E!%Sn|-gDS;xPPKu;ePc9e({7cquW4d?BZNpoR6E)E%Kf1
zTL>)h%?Hi`&h*93VJpTqTF0I5S?oCvI2Sm_vj{lbv!@PgG}gCchb@bC!+9-pX8aqw
zXnp=L>BFQCfGxoLz<Yoh-E@9;h`hH1@7|aO2JaX#mrH`lmE;5*{L5c$NgAp)aqcMX
zHt(?4;XdQ*rhhztbGCW5dHw>-=!!iPi18le7+}wwEUm9E_U^l+?~=X+z5%}Gd<CX;
zoTxFfzF0Z8$L&c1TplOj@W7WMgK6qKaqbG;1;7AzKG5IY&z%QEcW8}yv31y+Lf)Hj
zyUn??lJuOJgy_v#+{C}@<LwRf0(yE62C@OYLmN*Hdt>Y$I{Z|h<3jiFiaU6B<nV4^
zd8e<;+j|aMF6N&^gnbt63b)#~+P4b0&v!3ykMC|Enp+bkFP>bz?-XAh{|j6E%i1d6
zOy9}A9lM$+d9io+H^9dJ=^cM5i!WzQ&NI%ZIZru}Cpl(x@?8Dhq|E@H*>4P;zejgy
zjd?OtUsk$Jy82@0@Xdd?)zt+DbeEfxgtr5<iW4<P))y-`!6_(qhKZG%5T%XQa(Jf$
z-s@loY=8#D%3&+UHd@PFntVy}#lS_$7bY)Fz94x?^7+ZlIm|TXHFctNWL{Hu8MUI%
z+Y`4uX?fC>z!gc$fXji)!0gCjrZG>J)|Zv;a$S9~b62^R6XMf@uYgYp-ZLjl>+6fX
zd)51@_Z8q}?@Pdo-WPyq?zJd+@#LOxA&(P}xeznD%WT%bjNf4yt<U0ngylC0i|-ST
zmD9H(+h{HKk6q8-c0_Of;U@lFgI!Q>4-+fb5T%XQaxXhxZt_i|W^^ao^_*Zw^yWlv
z;@{cr4g*fR8Qu9#WHE7`6EUM(>XKaGl1Q^OWG)D0;k>4KZSSn{t>L-$)xatPW^~p#
zySmKQUK3<$Va(_jYkJPp5WQK<P5isVHEXcmrsr)NqBrZgiGSD0*2&h{fEnFyc59Gj
z*R#JJ(VHx8;@|N%XyHv-w+=Hp{N55P5Oi7thtuJ-8!)50-EIw@wCQ=uhUm?c+{C||
zWw!<^Z7Xee5qH`IW^l6Cvd>#JJ-=#*-fZP2{@u%(HTcl+0kOq_ybtV|lcn|b#ooQB
zNnX$-(!3ZlFNCsiUemm`cT;VGDYh`Na#N$U(OT|3?LF;X17>vVRBPbXtif06mj=F2
z&FHRD^(<Ery}62;_;>H9dN!zt-n_$2{5xxG9bPTlB7CM<g9B`WuC_2{bPF_VaFk{Z
z4py@b9Hg4j{iIlfzcs-&EsPo6db{8)yWqd}|B_ix<}EPs@BFGYU>ei5tC5(|ovP}Y
zry_cDDmU@(Zcqi+t6|LOHmV!x%v_wa<4-fXEvlaPRYY&La1;N|jAI|I&--h70ve(>
z`*RcjZh&eH(lu*vscH?Dsn*~FP0tn$(VGvriGTN*?K9h_oKI|mkHMJFjZm#YPgPGZ
z715iX+{C~8Q_=I6g6Pek+{C|oR~2kl!<f<Crb=#ACDPm$GPj1Za9-2Aws%F^5CelX
zGrE)1lhhIeW^`XD*5DV#8vLeMgCjJ7{;p~>x_z~M40P1Y=;kYW&QcJ)na@r9yD^HM
z(F&qBW4MWb_l{x>HY<8IDTv-|<|h8#RK*&sSJo?U88D+0uQ<&AbD5&&as|<wW!%KS
z`?q2ZCaIIui3ZH*)+*K@pz7(OB6<_xCjMPeu?DZGFB^DCHKP+}BIf^jh$<MUhB2et
zpa}G5y=HVTDc0aQ#Tu+qtifH1HL%86$y&vX$Gm8L{!;sb_*_HquH-#)vb4Uw*t_M5
zH7HOe;rlbBDbU@s5q4o4Pi~8{MS0(V8Qt@W9&t?6o9F2`8SgG%u?9nmj~>oktc1ao
z&zJ(?mnB`gvehVZ7~DF_*o$-fjAt$temj^!E@MIO;g_?I;TKlHm<_^feNMWC9f^gz
zGiHaNjy~j8b_M*xI>Q5umDjIgb(qnUF$aWKO3rwRU5lcMF(-tw%P;+kWx_9P6#fA5
z+8bY$9kYIfc$D}VlM&#{KPh6|BUXmUg{WUEB36Kig&(n;BbIT*TJ*pBkG>@)&nZ#q
zu8+VDxBnw_+kaQ%_Wu($GE2ApPv7bGztp(>uie@9zvPq={#}1EiOd#{z;0~+v{o?v
z3c&OWpXrw~(=TJDUyFJi`7ItQ5>s1_!0?enhp?G=K$g?gT<_}|%YZMy&0wr#Ot7)6
ze0sm4+S>Z^%9_%~U;|5mQvf{ZGZtW~q9}NZfRl{)J8svw!_b(ADjHr{RW+%!YF3b?
z&^!jeVR_8ONyhElH8?|XVH>k*DjOSEI?Tk%@PBo+#*)R%hQ`vG^3wY9lA6-$AW!N<
zN#kl8gEftnrB#98?9#g0#`@BlhPul7(#Fc#8Vw5by31;-t09XQoV*JK1EpnUv#MuR
z!5EYWYG()Q1C15IKzZrBo`JH`sxq1}r?Rnv&lo;<Y_I<Lf$CsmMQyppeAH)14F;x{
zmNnMa^OkWS@p<*l$;r`RsHkh<(Z(n5p`?a7cnYOeCH2AS8uP>O!4eueYV@!%6Gn|4
zQ*!jE3B{}<q~ZX#hyQr07Ht9HK(ES%A296)|Mic%7#;`UF^;$w7l&LMhE`TrmPf7!
zkd~NCI0Ah8FFHim^N*&V`@KoGe;+&z+?dAI2g@q!DjS>M;E$rrNwrn8s)K>r>4C=D
zM%Z7JmsVEI3)Jb$2ZD3qaYF-a^`#AgIk5TH2P%S<Gb$Pbl{JCt%Bo7(#)I_@@w<)h
z3r|+RjpJ@o{}^7PCAM)Gga5D7qT{;_#^=70;(ChPD0V$T%s-l9BgJbeK1=ai5$AvK
f-7#3-{LK{KruYWMS13MD@hOUrQhbQwYKZ>_IU|H)

literal 0
HcmV?d00001

diff --git a/code/modspa_samir.py b/code/modspa_samir.py
index cce94e4..3957598 100644
--- a/code/modspa_samir.py
+++ b/code/modspa_samir.py
@@ -10,7 +10,7 @@ Usage of the SAMIR model in the Modspa framework.
 import os  # for path exploration
 import csv  # open csv files
 from fnmatch import fnmatch  # for character string comparison
-from typing import List, Tuple  # to declare variables
+from typing import List, Tuple, Union  # to declare variables
 import xarray as xr  # to manage dataset
 import pandas as pd  # to manage dataframes
 import numpy as np  # for math and array operations
@@ -18,7 +18,42 @@ import rasterio as rio  # to open geotiff files
 import geopandas as gpd  # to manage shapefile crs projections
 from parameters.params_samir_class import samir_parameters
 
-def rasterize_samir_parameters(csv_param_file: str, parameter_dataset: xr.Dataset, land_cover_raster: str) -> xr.Dataset:
+
+def xr_maximum(ds: xr.DataArray, value: Union[xr.DataArray, float, int]) -> xr.DataArray:
+    """
+    Equivalent of `numpy.maximum(ds, value)` for xarray DataArrays
+
+    ## Arguments
+    1. ds: `xr.DataArray`
+        datarray to compare
+    2. value: `Union[xr.DataArray, float, int]`
+        value (scalar or dataarray) to compare
+
+    ## Returns
+    1. output: `xr.DataArray`
+        resulting dataarray with maximum value element-wise
+    """
+    return xr.where(ds <= value, value, ds)
+
+
+def xr_minimum(ds: xr.DataArray, value: Union[xr.DataArray, float, int]) -> xr.DataArray:
+    """
+    Equivalent of `numpy.minimum(ds, value)` for xarray DataArrays
+
+    ## Arguments
+    1. ds: `xr.DataArray`
+        datarray to compare
+    2. value: `Union[xr.DataArray, float, int]`
+        value (scalar or dataarray) to compare
+
+    ## Returns
+    1. output: `xr.DataArray`
+        resulting dataarray with minimum value element-wise
+    """
+    return xr.where(ds >= value, value, ds)
+
+
+def rasterize_samir_parameters(csv_param_file: str, empty_dataset: xr.Dataset, land_cover_raster: str, chunk_size: dict) -> Tuple[xr.Dataset, dict]:
     """
     Creates a raster `xarray` dataset from the csv parameter file, the land cover raster and an empty dataset
     that contains the right structure (emptied ndvi dataset for example). For each parameter, the function loops
@@ -27,14 +62,18 @@ def rasterize_samir_parameters(csv_param_file: str, parameter_dataset: xr.Datase
     ## Arguments
     1. csv_param_file: `str`
         path to csv paramter file
-    2. parameter_dataset: `xr.Dataset`
-        empty dataset that contains the right structure (emptied ndvi dataset for example)
+    2. empty_dataset: `xr.Dataset`
+        empty dataset that contains the right structure (emptied ndvi dataset for example).
     3. land_cover_raster: `str`
         path to land cover netcdf raster
+    4. chunk_size: `dict`
+        chunk_size for dask computation
 
     ## Returns
     1. parameter_dataset: `xr.Dataset`
         the dataset containing all the rasterized Parameters
+    2. scale_factor: `dict`
+        dictionnary containing the scale factors for each parameter
     """
     
     # Load samir params into an object
@@ -44,7 +83,13 @@ def rasterize_samir_parameters(csv_param_file: str, parameter_dataset: xr.Datase
     class_count = table_param.table.shape[1] - 2  # remove dtype and default columns
     
     # Open land cover raster
-    land_cover = xr.open_dataarray(land_cover_raster)
+    land_cover = xr.open_dataarray(land_cover_raster, chunks = chunk_size)
+    
+    # Create dataset
+    parameter_dataset = empty_dataset.copy(deep = True)
+    
+    # Create dictionnary containing the scale factors
+    scale_factor = {}
     
     # Loop on samir parameters and create 
     for parameter in table_param.table.index[1:]:
@@ -53,6 +98,10 @@ def rasterize_samir_parameters(csv_param_file: str, parameter_dataset: xr.Datase
         parameter_dataset[parameter] = land_cover.astype('i2')
         parameter_dataset[parameter].attrs['name'] = parameter
         parameter_dataset[parameter].attrs['description'] = 'cf SAMIR Doc for detail'
+        parameter_dataset[parameter].attrs['scale factor'] = str(table_param.table.loc[table_param.table.index == parameter]['scale_factor'].values[0])
+        
+        # Assigne value in dictionnary
+        scale_factor[parameter] = 1/int(table_param.table.loc[table_param.table.index == parameter]['scale_factor'].values[0])
         
         # Loop on classes to set parameter values for each class
         for class_val, class_name in zip(range(1, class_count + 1), table_param.table.columns[2:]):
@@ -62,23 +111,24 @@ def rasterize_samir_parameters(csv_param_file: str, parameter_dataset: xr.Datase
             parameter_dataset[parameter].values = np.where(parameter_dataset[parameter].values == class_val, round(table_param.table.loc[table_param.table.index == parameter][class_name].values[0]*table_param.table.loc[table_param.table.index == parameter]['scale_factor'].values[0]), parameter_dataset[parameter].values).astype('i2')
     
     # Return dataset converted to 'int16' data type to reduce memory usage
-    # Scale factor for calculation is stored in the samir_parameters object
-    return parameter_dataset
+    # and scale_factor dictionnary for later conversion
+    return parameter_dataset, scale_factor
 
 
-def setup_time_loop(calculation_variables: List[str], calculation_constant_values: List[str], empty_dataset: xr.Dataset) -> Tuple[xr.Dataset, xr.Dataset, xr.Dataset]:
+def setup_time_loop(calculation_variables_t1: List[str], calculation_variables_t2: List[str], empty_dataset: xr.Dataset) -> Tuple[xr.Dataset, xr.Dataset]:
     """
-    Creates three temporary `xarray Datasets` that will be used in the SAMIR time loop.
+    Creates two temporary `xarray Datasets` that will be used in the SAMIR time loop.
     `variables_t1` corresponds to the variables for the previous day and `variables_t2`
     corresponds to the variables for the current day. After each loop, `variables_t1`
-    takes the value of `variables_t2`. `constant_values` corresponds to model values
-    that are not time dependant, and therefore stay constant during the SAMIR time loop.
+    takes the value of `variables_t2` for the corresponding variables.
 
     ## Arguments
-    1. calculation_variables: `List[str]`
+    1. calculation_variables_t1: `List[str]`
         list of strings containing the variable names
-    2. calculation_constant_values: `List[str]`
-        list of strings containing the constant value names
+        for the previous day dataset
+    2. calculation_variables_t2: `List[str]`
+        list of strings containing the variable names
+        for the current day dataset
     3. empty_dataset: `xr.Dataset`
         empty dataset that contains the right structure
 
@@ -87,35 +137,157 @@ def setup_time_loop(calculation_variables: List[str], calculation_constant_value
         output dataset for previous day
     2. variables_t2: `xr.Dataset`
         output dataset for current day
-    3. constant_values: `xr.Dataset`
-        output dataset for constant values
     """
     
     # Create new dataset
     variables_t1 = empty_dataset.copy(deep = True)
     
     # Create empty DataArray for each variable
-    for variable in calculation_variables:
+    for variable in calculation_variables_t1:
         
         # Assign new empty DataArray
         variables_t1[variable] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'float32'))
         variables_t1[variable].attrs['name'] = variable  # set name in attributes
     
-    # Create new copy of the variables_t1 dataset
-    variables_t2 = variables_t1.copy(deep = True)
-    
-    # Create dataset for constant values
-    constant_values = empty_dataset.copy(deep = True)
+    # Create new dataset
+    variables_t2 = empty_dataset.copy(deep = True)
     
-    # Create empty DataArray for each value
-    for value in calculation_constant_values:
+    # Create empty DataArray for each variable
+    for variable in calculation_variables_t2:
         
         # Assign new empty DataArray
-        constant_values[value] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
-        constant_values[value].attrs['name'] = value  # set name in attributes
-        constant_values[value].attrs['description'] = 'Values which stays constant during the SAMIR time loop'  # set description in attributes
+        variables_t2[variable] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'float32'))
+        variables_t2[variable].attrs['name'] = variable  # set name in attributes
+    
+    return variables_t1, variables_t2
+
+
+def prepare_outputs(empty_dataset: xr.Dataset, additional_outputs: List[str] = None) -> xr.Dataset:
+    """
+    Creates the `xarray Dataset` containing the outputs of the SAMIR model that will be saved.
+    Additional variables can be saved by adding their names to the `additional_outputs` list.
+
+    ## Arguments
+    1. empty_dataset: `xr.Dataset`
+        empty dataset that contains the right structure
+    2. additional_outputs: `List[str]`
+        list of additional variable names to be saved
+
+    ## Returns
+    1. model_outputs: `xr.Dataset`
+        model outputs to be saved
+    """
+    
+    # Evaporation and Transpiraion
+    model_outputs = empty_dataset.copy(deep = True)
+    
+    model_outputs['E'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['E'].attrs['units'] = 'mm'
+    model_outputs['E'].attrs['standard_name'] = 'Evaporation'
+    model_outputs['E'].attrs['description'] = 'Accumulated daily evaporation in milimeters'
+    model_outputs['E'].attrs['scale factor'] = '1'
+    
+    model_outputs['Tr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['Tr'].attrs['units'] = 'mm'
+    model_outputs['Tr'].attrs['standard_name'] = 'Transpiration'
+    model_outputs['Tr'].attrs['description'] = 'Accumulated daily plant transpiration in milimeters'
+    model_outputs['Tr'].attrs['scale factor'] = '1'
+    
+    # Soil Water Content
+    model_outputs['SWCe'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['SWCe'].attrs['units'] = 'mm'
+    model_outputs['SWCe'].attrs['standard_name'] = 'Soil Water Content of the evaporative zone'
+    model_outputs['SWCe'].attrs['description'] = 'Soil water content of the evaporative zone in milimeters'
+    model_outputs['SWCe'].attrs['scale factor'] = '1'
+    
+    model_outputs['SWCr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['SWCr'].attrs['units'] = 'mm'
+    model_outputs['SWCr'].attrs['standard_name'] = 'Soil Water Content of the root zone'
+    model_outputs['SWCr'].attrs['description'] = 'Soil water content of the root zone in milimeters'
+    model_outputs['SWCr'].attrs['scale factor'] = '1'
+    
+    # Irrigation
+    model_outputs['Irr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['Irr'].attrs['units'] = 'mm'
+    model_outputs['Irr'].attrs['standard_name'] = 'Irrigation'
+    model_outputs['Irr'].attrs['description'] = 'Simulated daily irrigation in milimeters'
+    model_outputs['Irr'].attrs['scale factor'] = '1'
+    
+    # Deep Percolation
+    model_outputs['DP'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    model_outputs['DP'].attrs['units'] = 'mm'
+    model_outputs['DP'].attrs['standard_name'] = 'Deep Percolation'
+    model_outputs['DP'].attrs['description'] = 'Simulated daily Deep Percolation in milimeters'
+    model_outputs['DP'].attrs['scale factor'] = '1'
+    
+    if additional_outputs:
+        for var in additional_outputs:
+            model_outputs[var] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))
+    
+    return model_outputs
+
+
+def calculate_diff_re(variable_ds: xr.Dataset, param_ds: xr.Dataset, scale_dict: dict, var: str) -> xr.DataArray:
+    """
+    Calculates the diffusion between the top soil layer and the root layer.
+
+    ## Arguments
+    1. variable_ds: `xr.Dataset`
+        dataset containing calculation variables
+    2. param_ds: `xr.Dataset`
+        dataset containing the rasterized parameters
+    3. scale_dict: `dict`
+        dictionnary containing the scale factors for
+        the rasterized parameters
+    4. var: `str`
+        name of depletion variable to use (Dei or Dep or De)
+
+    ## Returns
+    1. diff_re: `xr.Dataset`
+        the diffusion between the top soil layer and
+        the root layer
+    """
+    
+    # Temporary variables to make calculation easier to read
+    # tmp1 = (((variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr'] - (variable_ds['RUE'] - variable_ds[var]) / (scale_dict['Ze'] * param_ds['Ze'])) / variable_ds['FCov']) * (scale_dict['DiffE'] * param_ds['DiffE'])
+    # tmp2 = ((variable_ds['TAW'] * scale_dict['Ze'] * param_ds['Ze']) - (variable_ds['RUE'] - variable_ds[var] - variable_ds['Dr']) * variable_ds['Zr']) / (variable_ds['Zr'] + scale_dict['Ze'] * param_ds['Ze']) - variable_ds['Dr']
+    
+    # Calculate diffusion according to SAMIR equation
+    # diff_re = xr.where(tmp1 < 0, xr_maximum(tmp1, tmp2), xr_minimum(tmp1, tmp2))
+
+    # Return zero values where the 'DiffE' parameter is equal to 0
+    return xr.where(param_ds['DiffE'] == 0, 0, xr.where((((variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr'] - (variable_ds['RUE'] - variable_ds[var]) / (scale_dict['Ze'] * param_ds['Ze'])) / variable_ds['FCov']) * (scale_dict['DiffE'] * param_ds['DiffE']) < 0, xr_maximum((((variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr'] - (variable_ds['RUE'] - variable_ds[var]) / (scale_dict['Ze'] * param_ds['Ze'])) / variable_ds['FCov']) * (scale_dict['DiffE'] * param_ds['DiffE']), ((variable_ds['TAW'] * scale_dict['Ze'] * param_ds['Ze']) - (variable_ds['RUE'] - variable_ds[var] - variable_ds['Dr']) * variable_ds['Zr']) / (variable_ds['Zr'] + scale_dict['Ze'] * param_ds['Ze']) - variable_ds['Dr']), xr_minimum((((variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr'] - (variable_ds['RUE'] - variable_ds[var]) / (scale_dict['Ze'] * param_ds['Ze'])) / variable_ds['FCov']) * (scale_dict['DiffE'] * param_ds['DiffE']), ((variable_ds['TAW'] * scale_dict['Ze'] * param_ds['Ze']) - (variable_ds['RUE'] - variable_ds[var] - variable_ds['Dr']) * variable_ds['Zr']) / (variable_ds['Zr'] + scale_dict['Ze'] * param_ds['Ze']) - variable_ds['Dr'])))
+
+
+def calculate_diff_dr(variable_ds: xr.Dataset, param_ds: xr.Dataset, scale_dict: dict) -> xr.DataArray:
+    """
+    Calculates the diffusion between the root layer and the deep layer.
+
+    ## Arguments
+    1. variable_ds: `xr.Dataset`
+        dataset containing calculation variables
+    2. param_ds: `xr.Dataset`
+        dataset containing the rasterized parameters
+    3. scale_dict: `dict`
+        dictionnary containing the scale factors for
+        the rasterized parameters
+
+    ## Returns
+    1. diff_dr: `xr.Dataset`
+        the diffusion between the root layer and the
+        deep layer
+    """
+    
+    # Temporary variables to make calculation easier to read
+    # tmp1 = (((variable_ds['TDW'] - variable_ds['Dd']) / (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr']) - (variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr']) / variable_ds['FCov']) * scale_dict['DiffR'] * param_ds['DiffR']
+    # tmp2 = (variable_ds['TDW'] * variable_ds['Zr'] - (variable_ds['TAW'] - variable_ds['Dr'] - variable_ds['Dd']) * (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr'])) / (scale_dict['Zsoil'] * param_ds['Zsoil']) - variable_ds['Dd']
+    
+    # Calculate diffusion according to SAMIR equation
+    # diff_dr = xr.where(tmp1 < 0, xr_maximum(tmp1, tmp2), xr_minimum(tmp1, tmp2))
     
-    return variables_t1, variables_t2, constant_values
+    # Return zero values where the 'DiffR' parameter is equal to 0
+    # return xr.where(param_ds['DiffR'] == 0, 0, diff_dr)
+    return xr.where(param_ds['DiffR'] == 0, 0, xr.where((((variable_ds['TDW'] - variable_ds['Dd']) / (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr']) - (variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr']) / variable_ds['FCov']) * scale_dict['DiffR'] * param_ds['DiffR'] < 0, xr_maximum((((variable_ds['TDW'] - variable_ds['Dd']) / (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr']) - (variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr']) / variable_ds['FCov']) * scale_dict['DiffR'] * param_ds['DiffR'], (variable_ds['TDW'] * variable_ds['Zr'] - (variable_ds['TAW'] - variable_ds['Dr'] - variable_ds['Dd']) * (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr'])) / (scale_dict['Zsoil'] * param_ds['Zsoil']) - variable_ds['Dd']), xr_minimum((((variable_ds['TDW'] - variable_ds['Dd']) / (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr']) - (variable_ds['TAW'] - variable_ds['Dr']) / variable_ds['Zr']) / variable_ds['FCov']) * scale_dict['DiffR'] * param_ds['DiffR'], (variable_ds['TDW'] * variable_ds['Zr'] - (variable_ds['TAW'] - variable_ds['Dr'] - variable_ds['Dd']) * (scale_dict['Zsoil'] * param_ds['Zsoil'] - variable_ds['Zr'])) / (scale_dict['Zsoil'] * param_ds['Zsoil']) - variable_ds['Dd'])))
 
 
 def run_samir():
diff --git a/dev_samir_xarray.ipynb b/dev_samir_xarray.ipynb
new file mode 100644
index 0000000..99173ab
--- /dev/null
+++ b/dev_samir_xarray.ipynb
@@ -0,0 +1,1257 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import xarray as xr\n",
+    "from dask.distributed import Client\n",
+    "import os, sys\n",
+    "import pandas as pd\n",
+    "import numpy as np\n",
+    "from typing import List, Tuple, Union\n",
+    "import warnings\n",
+    "import gc\n",
+    "\n",
+    "from parameters.params_samir_class import samir_parameters\n",
+    "from config.config import config\n",
+    "ndvi_path = '/mnt/e/DATA/DEV_inputs_test/ndvi.nc'\n",
+    "weather_path = '/mnt/e/DATA/DEV_inputs_test/weather.nc'\n",
+    "land_cover_path = '/mnt/e/DATA/DEV_inputs_test/land_cover.nc'\n",
+    "\n",
+    "param_file = '/home/auclairj/GIT/modspa-pixel/parameters/csv_files/params_samir_test.csv'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "\n",
+    "def rasterize_samir_parameters(csv_param_file: str, empty_dataset: xr.Dataset, land_cover_raster: str, chunk_size: dict) -> Tuple[xr.Dataset, dict]:\n",
+    "    \"\"\"\n",
+    "    Creates a raster `xarray` dataset from the csv parameter file, the land cover raster and an empty dataset\n",
+    "    that contains the right structure (emptied ndvi dataset for example). For each parameter, the function loops\n",
+    "    on land cover classes to fill the raster.\n",
+    "\n",
+    "    ## Arguments\n",
+    "    1. csv_param_file: `str`\n",
+    "        path to csv paramter file\n",
+    "    2. empty_dataset: `xr.Dataset`\n",
+    "        empty dataset that contains the right structure (emptied ndvi dataset for example).\n",
+    "    3. land_cover_raster: `str`\n",
+    "        path to land cover netcdf raster\n",
+    "    4. chunk_size: `dict`\n",
+    "        chunk_size for dask computation\n",
+    "\n",
+    "    ## Returns\n",
+    "    1. parameter_dataset: `xr.Dataset`\n",
+    "        the dataset containing all the rasterized Parameters\n",
+    "    2. scale_factor: `dict`\n",
+    "        dictionnary containing the scale factors for each parameter\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    # Load samir params into an object\n",
+    "    table_param = samir_parameters(csv_param_file)\n",
+    "    \n",
+    "    # Set general variables\n",
+    "    class_count = table_param.table.shape[1] - 2  # remove dtype and default columns\n",
+    "    \n",
+    "    # Open land cover raster\n",
+    "    land_cover = xr.open_dataarray(land_cover_raster, chunks = chunk_size)\n",
+    "    \n",
+    "    # Create dataset\n",
+    "    parameter_dataset = empty_dataset.copy(deep = True)\n",
+    "    \n",
+    "    # Create dictionnary containing the scale factors\n",
+    "    scale_factor = {}\n",
+    "    \n",
+    "    # Loop on samir parameters and create \n",
+    "    for parameter in table_param.table.index[1:]:\n",
+    "        \n",
+    "        # Create new variable and set attributes\n",
+    "        parameter_dataset[parameter] = land_cover.copy(deep = True).astype('f4')\n",
+    "        parameter_dataset[parameter].attrs['name'] = parameter\n",
+    "        parameter_dataset[parameter].attrs['description'] = 'cf SAMIR Doc for detail'\n",
+    "        parameter_dataset[parameter].attrs['scale factor'] = str(table_param.table.loc[table_param.table.index == parameter]['scale_factor'].values[0])\n",
+    "        \n",
+    "        # Assigne value in dictionnary\n",
+    "        scale_factor[parameter] = 1/int(table_param.table.loc[table_param.table.index == parameter]['scale_factor'].values[0])\n",
+    "        \n",
+    "        # Loop on classes to set parameter values for each class\n",
+    "        for class_val, class_name in zip(range(1, class_count + 1), table_param.table.columns[2:]):\n",
+    "            \n",
+    "            # Parameter values are multiplied by the scale factor in order to store all values as int16 types\n",
+    "            # These values are then rounded to make sure there isn't any decimal point issues when casting the values to int16\n",
+    "            parameter_dataset[parameter].values = np.where(parameter_dataset[parameter].values == class_val, round(table_param.table.loc[table_param.table.index == parameter][class_name].values[0]*table_param.table.loc[table_param.table.index == parameter]['scale_factor'].values[0]), parameter_dataset[parameter].values).astype('f4')\n",
+    "    \n",
+    "    # Return dataset converted to 'int16' data type to reduce memory usage\n",
+    "    # and scale_factor dictionnary for later conversion\n",
+    "    return parameter_dataset, scale_factor\n",
+    "\n",
+    "\n",
+    "def setup_time_loop(calculation_variables_t1: List[str], calculation_variables_t2: List[str], empty_dataset: xr.Dataset) -> Tuple[xr.Dataset, xr.Dataset]:\n",
+    "    \"\"\"\n",
+    "    Creates two temporary `xarray Datasets` that will be used in the SAMIR time loop.\n",
+    "    `variables_t1` corresponds to the variables for the previous day and `variables_t2`\n",
+    "    corresponds to the variables for the current day. After each loop, `variables_t1`\n",
+    "    takes the value of `variables_t2` for the corresponding variables.\n",
+    "\n",
+    "    ## Arguments\n",
+    "    1. calculation_variables_t1: `List[str]`\n",
+    "        list of strings containing the variable names\n",
+    "        for the previous day dataset\n",
+    "    2. calculation_variables_t2: `List[str]`\n",
+    "        list of strings containing the variable names\n",
+    "        for the current day dataset\n",
+    "    3. empty_dataset: `xr.Dataset`\n",
+    "        empty dataset that contains the right structure\n",
+    "\n",
+    "    ## Returns\n",
+    "    1. variables_t1: `xr.Dataset`\n",
+    "        output dataset for previous day\n",
+    "    2. variables_t2: `xr.Dataset`\n",
+    "        output dataset for current day\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    # Create new dataset\n",
+    "    variables_t1 = empty_dataset.copy(deep = True)\n",
+    "    \n",
+    "    # Create empty DataArray for each variable\n",
+    "    for variable in calculation_variables_t1:\n",
+    "        \n",
+    "        # Assign new empty DataArray\n",
+    "        variables_t1[variable] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'float32'))\n",
+    "        variables_t1[variable].attrs['name'] = variable  # set name in attributes\n",
+    "    \n",
+    "    # Create new dataset\n",
+    "    variables_t2 = empty_dataset.copy(deep = True)\n",
+    "    \n",
+    "    # Create empty DataArray for each variable\n",
+    "    for variable in calculation_variables_t2:\n",
+    "        \n",
+    "        # Assign new empty DataArray\n",
+    "        variables_t2[variable] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'float32'))\n",
+    "        variables_t2[variable].attrs['name'] = variable  # set name in attributes\n",
+    "    \n",
+    "    return variables_t1, variables_t2\n",
+    "\n",
+    "\n",
+    "def prepare_outputs(empty_dataset: xr.Dataset, additional_outputs: List[str] = None) -> xr.Dataset:\n",
+    "    \"\"\"\n",
+    "    Creates the `xarray Dataset` containing the outputs of the SAMIR model that will be saved.\n",
+    "    Additional variables can be saved by adding their names to the `additional_outputs` list.\n",
+    "\n",
+    "    ## Arguments\n",
+    "    1. empty_dataset: `xr.Dataset`\n",
+    "        empty dataset that contains the right structure\n",
+    "    2. additional_outputs: `List[str]`\n",
+    "        list of additional variable names to be saved\n",
+    "\n",
+    "    ## Returns\n",
+    "    1. model_outputs: `xr.Dataset`\n",
+    "        model outputs to be saved\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    # Evaporation and Transpiraion\n",
+    "    model_outputs = empty_dataset.copy(deep = True)\n",
+    "    \n",
+    "    model_outputs['E'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+    "    model_outputs['E'].attrs['units'] = 'mm'\n",
+    "    model_outputs['E'].attrs['standard_name'] = 'Evaporation'\n",
+    "    model_outputs['E'].attrs['description'] = 'Accumulated daily evaporation in milimeters'\n",
+    "    model_outputs['E'].attrs['scale factor'] = '1'\n",
+    "    \n",
+    "    model_outputs['Tr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+    "    model_outputs['Tr'].attrs['units'] = 'mm'\n",
+    "    model_outputs['Tr'].attrs['standard_name'] = 'Transpiration'\n",
+    "    model_outputs['Tr'].attrs['description'] = 'Accumulated daily plant transpiration in milimeters'\n",
+    "    model_outputs['Tr'].attrs['scale factor'] = '1'\n",
+    "    \n",
+    "    # Soil Water Content\n",
+    "    model_outputs['SWCe'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+    "    model_outputs['SWCe'].attrs['units'] = 'mm'\n",
+    "    model_outputs['SWCe'].attrs['standard_name'] = 'Soil Water Content of the evaporative zone'\n",
+    "    model_outputs['SWCe'].attrs['description'] = 'Soil water content of the evaporative zone in milimeters'\n",
+    "    model_outputs['SWCe'].attrs['scale factor'] = '1'\n",
+    "    \n",
+    "    model_outputs['SWCr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+    "    model_outputs['SWCr'].attrs['units'] = 'mm'\n",
+    "    model_outputs['SWCr'].attrs['standard_name'] = 'Soil Water Content of the root zone'\n",
+    "    model_outputs['SWCr'].attrs['description'] = 'Soil water content of the root zone in milimeters'\n",
+    "    model_outputs['SWCr'].attrs['scale factor'] = '1'\n",
+    "    \n",
+    "    # Irrigation\n",
+    "    model_outputs['Irr'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+    "    model_outputs['Irr'].attrs['units'] = 'mm'\n",
+    "    model_outputs['Irr'].attrs['standard_name'] = 'Irrigation'\n",
+    "    model_outputs['Irr'].attrs['description'] = 'Simulated daily irrigation in milimeters'\n",
+    "    model_outputs['Irr'].attrs['scale factor'] = '1'\n",
+    "    \n",
+    "    # Deep Percolation\n",
+    "    model_outputs['DP'] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+    "    model_outputs['DP'].attrs['units'] = 'mm'\n",
+    "    model_outputs['DP'].attrs['standard_name'] = 'Deep Percolation'\n",
+    "    model_outputs['DP'].attrs['description'] = 'Simulated daily Deep Percolation in milimeters'\n",
+    "    model_outputs['DP'].attrs['scale factor'] = '1'\n",
+    "    \n",
+    "    if additional_outputs:\n",
+    "        for var in additional_outputs:\n",
+    "            model_outputs[var] = (empty_dataset.dims, np.zeros(tuple(empty_dataset.dims[d] for d in list(empty_dataset.dims)), dtype = 'int16'))\n",
+    "    \n",
+    "    return model_outputs\n",
+    "\n",
+    "\n",
+    "def xr_maximum(ds: xr.DataArray, value: Union[xr.DataArray, float, int]) -> xr.DataArray:\n",
+    "    \"\"\"\n",
+    "    Equivalent of `numpy.maximum(ds, value)` for xarray DataArrays\n",
+    "\n",
+    "    ## Arguments\n",
+    "    1. ds: `xr.DataArray`\n",
+    "        datarray to compare\n",
+    "    2. value: `Union[xr.DataArray, float, int]`\n",
+    "        value (scalar or dataarray) to compare\n",
+    "\n",
+    "    ## Returns\n",
+    "    1. output: `xr.DataArray`\n",
+    "        resulting dataarray with maximum value element-wise\n",
+    "    \"\"\"\n",
+    "    return xr.where(ds <= value, value, ds, keep_attrs = True)\n",
+    "\n",
+    "\n",
+    "def xr_minimum(ds: xr.DataArray, value: Union[xr.DataArray, float, int]) -> xr.DataArray:\n",
+    "    \"\"\"\n",
+    "    Equivalent of `numpy.minimum(ds, value)` for xarray DataArrays\n",
+    "\n",
+    "    ## Arguments\n",
+    "    1. ds: `xr.DataArray`\n",
+    "        datarray to compare\n",
+    "    2. value: `Union[xr.DataArray, float, int]`\n",
+    "        value (scalar or dataarray) to compare\n",
+    "\n",
+    "    ## Returns\n",
+    "    1. output: `xr.DataArray`\n",
+    "        resulting dataarray with minimum value element-wise\n",
+    "    \"\"\"\n",
+    "    return xr.where(ds >= value, value, ds, keep_attrs = True)\n",
+    "\n",
+    "\n",
+    "def calculate_diff_re(TAW: xr.DataArray, Dr: xr.DataArray, Zr: xr.DataArray, RUE: xr.DataArray, De: xr.DataArray, FCov: xr.DataArray, Ze_: xr.DataArray, DiffE_: xr.DataArray, scale_dict: dict) -> xr.DataArray:\n",
+    "    \"\"\"\n",
+    "    Calculates the diffusion between the top soil layer and the root layer.\n",
+    "\n",
+    "    ## Arguments\n",
+    "    1. TAW: `xr.DataArray`\n",
+    "        water capacity of root layer\n",
+    "    2. Dr: `xr.DataArray`\n",
+    "        depletion of root layer\n",
+    "    3. Zr: `xr.DataArray`\n",
+    "        height of root layer\n",
+    "    4. RUE: `xr.DataArray`\n",
+    "        total available surface water\n",
+    "    5. De: `xr.DataArray`\n",
+    "        depletion of the evaporative layer\n",
+    "    6. FCov: `xr.DataArray`\n",
+    "        fraction cover of plants\n",
+    "    7. Ze_: `xr.DataArray`\n",
+    "        height of evaporative layer (paramter)\n",
+    "    8. DiffE_: `xr.DataArray`\n",
+    "        diffusion coefficient between evaporative\n",
+    "        and root layers (unitless, parameter)\n",
+    "    9. scale_dict: `dict`\n",
+    "        dictionnary containing the scale factors for\n",
+    "        the rasterized parameters\n",
+    "\n",
+    "    ## Returns\n",
+    "    1. diff_re: `xr.Dataset`\n",
+    "        the diffusion between the top soil layer and\n",
+    "        the root layer\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    # Temporary variables to make calculation easier to read\n",
+    "    tmp1 = (((TAW - Dr) / Zr - (RUE - De) / (scale_dict['Ze'] * Ze_)) / FCov) * (scale_dict['DiffE'] * DiffE_)\n",
+    "    tmp2 = ((TAW * scale_dict['Ze'] * Ze_) - (RUE - De - Dr) * Zr) / (Zr + scale_dict['Ze'] * Ze_) - Dr\n",
+    "    \n",
+    "    # Calculate diffusion according to SAMIR equation\n",
+    "    diff_re = xr.where(tmp1 < 0, xr_maximum(tmp1, tmp2), xr_minimum(tmp1, tmp2))\n",
+    "\n",
+    "    # Return zero values where the 'DiffE' parameter is equal to 0\n",
+    "    return xr.where(DiffE_ == 0, 0, diff_re)\n",
+    "\n",
+    "\n",
+    "def calculate_diff_dr(TAW: xr.DataArray, TDW: xr.DataArray, Dr: xr.DataArray, Zr: xr.DataArray, Dd: xr.DataArray, FCov: xr.DataArray, Zsoil_: xr.DataArray, DiffR_: xr.DataArray, scale_dict: dict) -> xr.DataArray:\n",
+    "    \"\"\"\n",
+    "    Calculates the diffusion between the root layer and the deep layer.\n",
+    "\n",
+    "    ## Arguments\n",
+    "    1. TAW: `xr.DataArray`\n",
+    "        water capacity of root layer\n",
+    "    2. TDW: `xr.DataArray`\n",
+    "        water capacity of deep layer\n",
+    "    3. Dr: `xr.DataArray`\n",
+    "        depletion of root layer\n",
+    "    4. Zr: `xr.DataArray`\n",
+    "        height of root layer\n",
+    "    5. Dd: `xr.DataArray`\n",
+    "        depletion of deep layer\n",
+    "    6. FCov: `xr.DataArray`\n",
+    "        fraction cover of plants\n",
+    "    7. Zsoil_: `xr.DataArray`\n",
+    "        total height of soil (paramter)\n",
+    "    8. DiffR_: `xr.DataArray`\n",
+    "        Diffusion coefficient between root\n",
+    "        and deep layers (unitless, parameter)\n",
+    "    9. scale_dict: `dict`\n",
+    "        dictionnary containing the scale factors for\n",
+    "        the rasterized parameters\n",
+    "\n",
+    "    ## Returns\n",
+    "    1. diff_dr: `xr.Dataset`\n",
+    "        the diffusion between the root layer and the\n",
+    "        deep layer\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    # Temporary variables to make calculation easier to read\n",
+    "    tmp1 = (((TDW - Dd) / (scale_dict['Zsoil'] * Zsoil_ - Zr) - (TAW - Dr) / Zr) / FCov) * scale_dict['DiffR'] * DiffR_\n",
+    "    tmp2 = (TDW *Zr - (TAW - Dr - Dd) * (scale_dict['Zsoil'] * Zsoil_ - Zr)) / (scale_dict['Zsoil'] * Zsoil_) - Dd\n",
+    "    \n",
+    "    # Calculate diffusion according to SAMIR equation\n",
+    "    diff_dr = xr.where(tmp1 < 0, xr_maximum(tmp1, tmp2), xr_minimum(tmp1, tmp2))\n",
+    "    \n",
+    "    # Return zero values where the 'DiffR' parameter is equal to 0\n",
+    "    return xr.where(DiffR_ == 0, 0, diff_dr)\n",
+    "\n",
+    "\n",
+    "def calculate_W(TEW: xr.DataArray, Dei: xr.DataArray, Dep: xr.DataArray, fewi: xr.DataArray, fewp: xr.DataArray) -> xr.DataArray:\n",
+    "    \"\"\"\n",
+    "    Calculate W, the weighting factor to split the energy available\n",
+    "    for evaporation depending on the difference in water availability\n",
+    "    in the two evaporation components, ensuring that the larger and\n",
+    "    the wetter, the more the evaporation occurs from that component\n",
+    "\n",
+    "    ## Arguments\n",
+    "    1. TEW: `xr.DataArray`\n",
+    "        water capacity of evaporative layer\n",
+    "    2. Dei: `xr.DataArray`\n",
+    "        depletion of the evaporative layer\n",
+    "        (irrigation part)\n",
+    "    3. Dep: `xr.DataArray`\n",
+    "        depletion of the evaporative layer\n",
+    "        (precipitation part)\n",
+    "    4. fewi: `xr.DataArray`\n",
+    "        soil fraction which is wetted by irrigation\n",
+    "        and exposed to evaporation\n",
+    "    5. fewp: `xr.DataArray`\n",
+    "        soil fraction which is wetted by precipitation\n",
+    "        and exposed to evaporation\n",
+    "\n",
+    "    ## Returns\n",
+    "    1. W: `xr.DataArray`\n",
+    "        weighting factor W\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    # Temporary variables to make calculation easier to read\n",
+    "    tmp = fewi * (TEW - Dei)\n",
+    "    \n",
+    "    # Calculate the weighting factor to split the energy available for evaporation\n",
+    "    W = 1 / (1 + (fewp * (TEW - Dep) / tmp ))\n",
+    "\n",
+    "    # Return W \n",
+    "    return xr.where(tmp > 0, W, 0)\n",
+    "\n",
+    "\n",
+    "def calculate_Kr(TEW: xr.DataArray, De: xr.DataArray, REW_: xr.DataArray, scale_dict: dict) -> xr.DataArray:\n",
+    "    \"\"\"\n",
+    "    calculates of the reduction coefficient for evaporation dependent \n",
+    "    on the amount of water in the soil using the FAO-56 method\n",
+    "\n",
+    "    ## Arguments\n",
+    "    1. TEW: `xr.DataArray`\n",
+    "        water capacity of evaporative layer\n",
+    "    2. De: `xr.DataArray`\n",
+    "        depletion of evaporative layer\n",
+    "    3. REW_: `xr.DataArray`\n",
+    "        readily evaporable water\n",
+    "    4. scale_dict: `dict`\n",
+    "        dictionnary containing the scale factors for\n",
+    "        the rasterized parameters\n",
+    "\n",
+    "    ## Returns\n",
+    "    1. Kr: `xr.DataArray`\n",
+    "        Kr coefficient\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    # Formula for calculating Kr\n",
+    "    Kr = (TEW - De) / (TEW - scale_dict['REW'] * REW_)\n",
+    "    \n",
+    "    # Return Kr\n",
+    "    return xr_maximum(0, xr_minimum(Kr, 1))\n",
+    "\n",
+    "\n",
+    "def update_Dr(TAW: xr.DataArray, TDW: xr.DataArray, Zr: xr.DataArray, TAW0: xr.DataArray, TDW0: xr.DataArray, Dr0: xr.DataArray, Dd0: xr.DataArray, Zr0: xr.DataArray) -> xr.DataArray:\n",
+    "    \"\"\"\n",
+    "    Return the updated depletion for the root layer\n",
+    "\n",
+    "    ## Arguments\n",
+    "    1. TAW: `xr.DataArray`\n",
+    "        water capacity of root layer for current day\n",
+    "    2. TDW: `xr.DataArray`\n",
+    "        water capacity of deep layer for current day\n",
+    "    3. Zr: `xr.DataArray`\n",
+    "        root layer height for current day\n",
+    "    4. TAW0: `xr.DataArray`\n",
+    "        water capacity of root layer for previous day\n",
+    "    5. TDW0: `xr.DataArray`\n",
+    "        water capacity of deep layer for previous day\n",
+    "    6. Dr0: `xr.DataArray`\n",
+    "        depletion of the root layer for previous day\n",
+    "    7. Dd0: `xr.DataArray`\n",
+    "        depletion of the deep laye for previous day\n",
+    "    8. Zr0: `xr.DataArray`\n",
+    "        root layer height for previous day\n",
+    "\n",
+    "    ## Returns\n",
+    "    1. output: `xr.DataArray`\n",
+    "        updated depletion for the root layer\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    # Temporary variables to make calculation easier to read\n",
+    "    tmp1 = xr_maximum(Dr0 + Dd0 * (TAW - TAW0) / TDW0, 0)\n",
+    "    tmp2 = xr_minimum(Dr0 + Dd0 * (TAW - TAW0) / TDW0, TDW)\n",
+    "\n",
+    "    # Return updated Dr\n",
+    "    return xr.where(Zr > Zr0, tmp1, tmp2)\n",
+    "\n",
+    "\n",
+    "def update_Dd(TAW: xr.DataArray, TDW: xr.DataArray, Zr: xr.DataArray, TAW0: xr.DataArray, TDW0: xr.DataArray, Dd0: xr.DataArray, Zr0: xr.DataArray) -> xr.DataArray:\n",
+    "    \"\"\"\n",
+    "    Return the updated depletion for the deep layer\n",
+    "\n",
+    "    ## Arguments\n",
+    "    1. TAW: `xr.DataArray`\n",
+    "        water capacity of root layer for current day\n",
+    "    2. TDW: `xr.DataArray`\n",
+    "        water capacity of deep layer for current day\n",
+    "    3. TAW0: `xr.DataArray`\n",
+    "        water capacity of root layer for previous day\n",
+    "    5. TDW0: `xr.DataArray`\n",
+    "        water capacity of deep layer for previous day\n",
+    "    6. Dd0: `xr.DataArray`\n",
+    "        depletion of the deep laye for previous day\n",
+    "    7. Zr0: `xr.DataArray`\n",
+    "        root layer height for previous day\n",
+    "\n",
+    "    ## Returns\n",
+    "    1. output: `xr.DataArray`\n",
+    "        updated depletion for the deep layer\n",
+    "    \"\"\"\n",
+    "    \n",
+    "    # Temporary variables to make calculation easier to read\n",
+    "    tmp1 = xr_maximum(Dd0 - Dd0 * (TAW - TAW0) / TDW0, 0)\n",
+    "    tmp2 = xr_minimum(Dd0 - Dd0 * (TAW - TAW0) / TDW0, TDW)\n",
+    "    \n",
+    "    # Return updated Dd\n",
+    "    return xr.where(Zr > Zr0, tmp1, tmp2)\n",
+    "\n",
+    "\n",
+    "def run_samir(json_config_file: str, csv_param_file: str, ndvi_cube_path: str, weather_cube_path: str, soil_params_path: str, land_cover_path: str, chunk_size: dict, save_path: str) -> None:\n",
+    "    \n",
+    "    # warnings.simplefilter(\"error\", category = RuntimeWarning())\n",
+    "    warnings.filterwarnings(\"ignore\", message=\"invalid value encountered in cast\")\n",
+    "    warnings.filterwarnings(\"ignore\", message=\"invalid value encountered in divide\")\n",
+    "    np.errstate(all = 'raise')\n",
+    "    gc.disable()\n",
+    "    \n",
+    "    #============ General parameters ============#\n",
+    "    config_params = config(json_config_file)\n",
+    "    calculation_variables_t2 = ['diff_rei', 'diff_rep', 'diff_dr' , 'Dd', 'De', 'Dei', 'Dep', 'Dr', 'FCov', 'Irrig', 'Kcb', 'Kei', 'Kep', 'Ks', 'Kti', 'Ktp', 'RUE', 'TAW', 'TDW', 'TEW', 'Tei', 'Tep', 'W', 'Zr', 'fewi', 'fewp']\n",
+    "    calculation_variables_t1 = ['Dr', 'Dd', 'TAW', 'TDW', 'Zr']\n",
+    "    \n",
+    "    #============ Manage inputs ============#\n",
+    "    # NDVI\n",
+    "    ndvi_cube = xr.open_dataset(ndvi_cube_path, chunks = chunk_size).astype('f4')\n",
+    "    \n",
+    "    # # Create a daily DateTimeIndex for the desired date range\n",
+    "    # daily_index = pd.date_range(start = config_params.start_date, end = config_params.end_date, freq = 'D')\n",
+    "\n",
+    "    # # Resample the dataset to a daily frequency and reindex with the new DateTimeIndex\n",
+    "    # ndvi_cube = ndvi_cube.resample(time = '1D').asfreq().reindex(time = daily_index)\n",
+    "\n",
+    "    # # Interpolate the dataset along the time dimension to fill nan values\n",
+    "    # ndvi_cube = ndvi_cube.interpolate_na(dim = 'time', method = 'linear', fill_value = 'extrapolate').astype('u1')\n",
+    "    \n",
+    "    # Weather\n",
+    "    weather_cube = xr.open_dataset(weather_cube_path, chunks = chunk_size).astype('f4')\n",
+    "    \n",
+    "    # Soil\n",
+    "    soil_params = xr.open_dataset(soil_params_path, chunks = chunk_size).astype('f4')\n",
+    "    \n",
+    "    # SAMIR Parameters\n",
+    "    param_dataset, scale_factor = rasterize_samir_parameters(csv_param_file, ndvi_cube.drop_vars(['ndvi', 'time']), land_cover_path, chunk_size = chunk_size)\n",
+    "    \n",
+    "    # SAMIR Variables\n",
+    "    variables_t1, variables_t2 = setup_time_loop(calculation_variables_t1, calculation_variables_t2, ndvi_cube.drop_vars(['ndvi', 'time']))\n",
+    "\n",
+    "    #============ Prepare outputs ============#\n",
+    "    model_outputs = prepare_outputs(ndvi_cube.drop_vars(['ndvi']))\n",
+    "    \n",
+    "    #============ Prepare time iterations ============#\n",
+    "    dates = ndvi_cube.time.values\n",
+    "    \n",
+    "    #============ Create aliases for better readability ============#\n",
+    "    \n",
+    "    # Variables for current day\n",
+    "    diff_rei = variables_t2['diff_rei']\n",
+    "    diff_rep = variables_t2['diff_rep']\n",
+    "    diff_dr = variables_t2['diff_dr']\n",
+    "    Dd = variables_t2['Dd']\n",
+    "    De = variables_t2['De']\n",
+    "    Dei = variables_t2['Dei']\n",
+    "    Dep = variables_t2['Dep']\n",
+    "    Dr = variables_t2['Dr']\n",
+    "    FCov = variables_t2['FCov']\n",
+    "    Irrig = variables_t2['Irrig']\n",
+    "    Kcb = variables_t2['Kcb']\n",
+    "    Kei = variables_t2['Kei']\n",
+    "    Kep = variables_t2['Kep']\n",
+    "    Ks = variables_t2['Ks']\n",
+    "    Kti = variables_t2['Kti']\n",
+    "    Ktp = variables_t2['Ktp']\n",
+    "    RUE = variables_t2['RUE']\n",
+    "    TAW = variables_t2['TAW']\n",
+    "    TDW = variables_t2['TDW']\n",
+    "    TEW = variables_t2['TEW']\n",
+    "    Tei = variables_t2['Tei']\n",
+    "    Tep = variables_t2['Tep']\n",
+    "    Zr = variables_t2['Zr']\n",
+    "    W = variables_t2['W']\n",
+    "    fewi = variables_t2['fewi']\n",
+    "    fewp = variables_t2['fewp']\n",
+    "    \n",
+    "    # Variables for previous day\n",
+    "    TAW0 = variables_t1['TAW']\n",
+    "    TDW0 = variables_t1['TDW']\n",
+    "    Dr0 = variables_t1['Dr']\n",
+    "    Dd0 = variables_t1['Dd']\n",
+    "    Zr0 = variables_t1['Zr']\n",
+    "    \n",
+    "    # Parameters\n",
+    "    # Parameters have an underscore (_) behind their name for recognition \n",
+    "    DiffE_ = param_dataset['DiffE']\n",
+    "    DiffR_ = param_dataset['DiffR']\n",
+    "    FW_ = param_dataset['FW']\n",
+    "    Fc_stop_ = param_dataset['Fc_stop']\n",
+    "    FmaxFC_ = param_dataset['FmaxFC']\n",
+    "    Foffset_ = param_dataset['Foffset']\n",
+    "    Fslope_ = param_dataset['Fslope']\n",
+    "    Init_RU_ = param_dataset['Init_RU']\n",
+    "    Irrig_auto_ = param_dataset['Irrig_auto']\n",
+    "    Kcmax_ = param_dataset['Kcmax']\n",
+    "    KmaxKcb_ = param_dataset['KmaxKcb']\n",
+    "    Koffset_ = param_dataset['Koffset']\n",
+    "    Kslope_ = param_dataset['Kslope']\n",
+    "    Lame_max_ = param_dataset['Lame_max']\n",
+    "    REW_ = param_dataset['REW']\n",
+    "    Ze_ = param_dataset['Ze']\n",
+    "    Zsoil_ = param_dataset['Zsoil']\n",
+    "    maxZr_ = param_dataset['maxZr']\n",
+    "    minZr_ = param_dataset['minZr']\n",
+    "    p_ = param_dataset['p']\n",
+    "    \n",
+    "    # scale factors\n",
+    "    # Scale factors have the following name scheme : s_ + parameter_name\n",
+    "    s_DiffE = scale_factor['DiffE']\n",
+    "    s_DiffR = scale_factor['DiffR']\n",
+    "    s_FW = scale_factor['FW']\n",
+    "    s_Fc_stop = scale_factor['Fc_stop']\n",
+    "    s_FmaxFC = scale_factor['FmaxFC']\n",
+    "    s_Foffset = scale_factor['Foffset']\n",
+    "    s_Fslope = scale_factor['Fslope']\n",
+    "    s_Init_RU = scale_factor['Init_RU']\n",
+    "    # s_Irrig_auto = scale_factor['Irrig_auto']\n",
+    "    s_Kcmax = scale_factor['Kcmax']\n",
+    "    s_KmaxKcb = scale_factor['KmaxKcb']\n",
+    "    s_Koffset = scale_factor['Koffset']\n",
+    "    s_Kslope = scale_factor['Kslope']\n",
+    "    s_Lame_max = scale_factor['Lame_max']\n",
+    "    s_REW = scale_factor['REW']\n",
+    "    s_Ze = scale_factor['Ze']\n",
+    "    s_Zsoil = scale_factor['Zsoil']\n",
+    "    s_maxZr = scale_factor['maxZr']\n",
+    "    s_minZr = scale_factor['minZr']\n",
+    "    s_p = scale_factor['p']\n",
+    "    \n",
+    "    #============ First day initialization ============#\n",
+    "    # Fraction cover\n",
+    "    FCov = s_Fslope * Fslope_ * (ndvi_cube['ndvi'].sel({'time': dates[0]})/255) + s_Foffset * Foffset_\n",
+    "    FCov = xr_minimum(xr_maximum(FCov, 0), s_Fc_stop * Fc_stop_)\n",
+    "    \n",
+    "    # Root depth upate\n",
+    "    Zr = s_minZr * minZr_ + (FCov / (s_FmaxFC * FmaxFC_)) * s_maxZr * (maxZr_ - minZr_)\n",
+    "    \n",
+    "    # Water capacities\n",
+    "    TEW = (soil_params['FC'] - soil_params['WP']/2) * s_Ze * Ze_\n",
+    "    RUE = (soil_params['FC'] - soil_params['WP']) * s_Ze * Ze_\n",
+    "    TAW = (soil_params['FC'] - soil_params['WP']) * Zr\n",
+    "    TDW = (soil_params['FC'] - soil_params['WP']) * (s_Zsoil * Zsoil_ - Zr)  # Zd = Zsoil - Zr\n",
+    "    \n",
+    "    # Depletions\n",
+    "    Dei = RUE * (1 - s_Init_RU * Init_RU_)\n",
+    "    Dep = RUE * (1 - s_Init_RU * Init_RU_)\n",
+    "    Dr = TAW * (1 - s_Init_RU * Init_RU_)\n",
+    "    Dd = TDW * (1 - s_Init_RU * Init_RU_)\n",
+    "    \n",
+    "    # Irrigation   ==============!!!!!\n",
+    "    Irrig = xr_minimum(xr_maximum(Dr - weather_cube['tp'].sel({'time': dates[0]}) / 1000, 0), s_Lame_max * Lame_max_) * Irrig_auto_\n",
+    "    Irrig = xr.where(Dr > TAW * s_p * p_, Irrig, 0)\n",
+    "    \n",
+    "    # Kcb\n",
+    "    Kcb = xr_minimum(s_Kslope * Kslope_ * (ndvi_cube['ndvi'].sel({'time': dates[0]}) / 255) + s_Koffset * Koffset_, s_KmaxKcb * KmaxKcb_)\n",
+    "    \n",
+    "    # Update depletions with rainfall and/or irrigation    \n",
+    "    ## DP\n",
+    "    model_outputs['DP'].loc[{'time': dates[0]}] = -xr_minimum(Dd + xr_minimum(Dr - (weather_cube['tp'].sel({'time': dates[0]}) / 1000) - Irrig, 0), 0)\n",
+    "    \n",
+    "    ## De\n",
+    "    Dei = xr_minimum(xr_maximum(Dei - (weather_cube['tp'].sel({'time': dates[0]}) / 1000) - Irrig / (s_FW * FW_ / 100), 0), TEW)\n",
+    "    Dep = xr_minimum(xr_maximum(Dep - (weather_cube['tp'].sel({'time': dates[0]}) / 1000), 0), TEW)\n",
+    "    \n",
+    "    fewi = xr_minimum(1 - FCov, (s_FW * FW_ / 100))\n",
+    "    fewp = 1 - FCov - fewi\n",
+    "    \n",
+    "    De = (Dei * fewi + Dep * fewp) / (fewi + fewp)\n",
+    "    # variables_t1['De'] = xr.where(variables_t1['De'].isfinite(), variables_t1['De'], variables_t1['Dei'] * (s_FW * FW_ / 100) + variables_t1['Dep'] * (1 - (s_FW * FW_ / 100)))  #================= find replacement for .isfinite() method !!!!!!!!!\n",
+    "\n",
+    "    ## Dr\n",
+    "    Dr = xr_minimum(xr_maximum(Dr - (weather_cube['tp'].sel({'time': dates[0]}) / 1000) - Irrig, 0), TAW)\n",
+    "    \n",
+    "    ## Dd\n",
+    "    Dd = xr_minimum(xr_maximum(Dd + xr_minimum(Dr - (weather_cube['tp'].sel({'time': dates[0]}) / 1000) - Irrig, 0), 0), TDW)\n",
+    "    \n",
+    "    # Diffusion coefficients\n",
+    "    diff_rei = calculate_diff_re(TAW, Dr, Zr, RUE, Dei, FCov, Ze_, DiffE_, scale_factor)\n",
+    "    diff_rep = calculate_diff_re(TAW, Dr, Zr, RUE, Dep, FCov, Ze_, DiffE_, scale_factor)\n",
+    "    diff_dr = calculate_diff_dr(TAW, TDW, Dr, Zr, Dd, FCov, Zsoil_, DiffR_, scale_factor)    \n",
+    "    \n",
+    "    # Weighing factor W\n",
+    "    W = calculate_W(TEW, Dei, Dep, fewi, fewp)\n",
+    "    \n",
+    "    # Soil water contents\n",
+    "    model_outputs['SWCe'].loc[{'time': dates[0]}] = 1 - De/TEW\n",
+    "    model_outputs['SWCr'].loc[{'time': dates[0]}] = 1 - Dr/TAW\n",
+    "    \n",
+    "    # Water Stress coefficient\n",
+    "    Ks = xr_minimum((TAW - Dr) / (TAW * (1 - s_p * p_)), 1)\n",
+    "    \n",
+    "    # Reduction coefficient for evaporation\n",
+    "    Kei = xr_minimum(W * calculate_Kr(TEW, Dei, REW_, scale_factor) * (s_Kcmax * Kcmax_ - Kcb), fewi * s_Kcmax * Kcmax_)\n",
+    "    Kep = xr_minimum((1 - W) * calculate_Kr(TEW, Dep, REW_, scale_factor) * (s_Kcmax * Kcmax_ - Kcb), fewp * s_Kcmax * Kcmax_)\n",
+    "    \n",
+    "    # Prepare coefficients for evapotranspiration\n",
+    "    Kti = xr_minimum(((s_Ze * Ze_ / Zr)**6) * (1 - Dei / TEW) / xr_maximum(1 - Dr / TAW, 0.001), 1)\n",
+    "    Ktp = xr_minimum(((s_Ze * Ze_ / Zr)**6) * (1 - Dep / TEW) / xr_maximum(1 - Dr / TAW, 0.001), 1)\n",
+    "    Tei = Kti * Ks * Kcb * weather_cube['ET0'].sel({'time': dates[0]}) / 1000\n",
+    "    Tep = Ktp * Ks * Kcb * weather_cube['ET0'].sel({'time': dates[0]}) / 1000\n",
+    "    \n",
+    "    # Update depletions\n",
+    "    Dei = xr.where(fewi > 0, xr_minimum(xr_maximum(Dei + (weather_cube['ET0'].sel({'time': dates[0]}) / 1000) * Kei / fewi + Tei - diff_rei, 0), TEW), xr_minimum(xr_maximum(Dei + Tei - diff_rei, 0), TEW))\n",
+    "    Dep = xr.where(fewp > 0, xr_minimum(xr_maximum(Dep + (weather_cube['ET0'].sel({'time': dates[0]}) / 1000) * Kep / fewp + Tep - diff_rep, 0), TEW), xr_minimum(xr_maximum(Dep + Tep - diff_rep, 0), TEW))\n",
+    "    \n",
+    "    De = (Dei * fewi + Dep * fewp) / (fewi + fewp)\n",
+    "    # De = xr.where(De.isfinite(), De, Dei * (s_FW * FW_ / 100) + Dep * (1 - (s_FW * FW_ / 100)))  #================= find replacement for .isfinite() method !!!!!!!!!\n",
+    "    \n",
+    "    # Evaporation\n",
+    "    model_outputs['E'].loc[{'time': dates[0]}]  = xr_maximum((Kei + Kep) * weather_cube['ET0'].sel({'time': dates[0]}) / 1000, 0)\n",
+    "    \n",
+    "    # Transpiration\n",
+    "    model_outputs['Tr'].loc[{'time': dates[0]}]  = Kcb * Ks * weather_cube['ET0'].sel({'time': dates[0]}) / 1000\n",
+    "    \n",
+    "    # Potential evapotranspiration and evaporative fraction ??\n",
+    "    \n",
+    "    # Update depletions (root and deep zones) at the end of the day\n",
+    "    Dr = xr_minimum(xr_maximum(Dr + model_outputs['E'].loc[{'time': dates[0]}] + model_outputs['Tr'].loc[{'time': dates[0]}] - diff_dr, 0), TAW)\n",
+    "    Dd = xr_minimum(xr_maximum(Dd + diff_dr, 0), TDW)\n",
+    "    \n",
+    "    # Write outputs\n",
+    "    model_outputs['Irr'].loc[{'time': dates[0]}] = Irrig\n",
+    "    \n",
+    "    # Update variable_t1 values\n",
+    "    for variable in calculation_variables_t1:\n",
+    "        variables_t1[variable] = variables_t2[variable].copy(deep = True)\n",
+    "    \n",
+    "    #============ Time loop ============#\n",
+    "    for i in range(1, len(dates)):\n",
+    "        \n",
+    "        # Update variables\n",
+    "        ## Fraction cover\n",
+    "        FCov = s_Fslope * Fslope_ * (ndvi_cube['ndvi'].sel({'time': dates[0]})/255) + s_Foffset * Foffset_\n",
+    "        FCov = xr_minimum(xr_maximum(FCov, 0), s_Fc_stop * Fc_stop_)\n",
+    "        \n",
+    "        ## Root depth upate\n",
+    "        Zr = s_minZr * minZr_ + (FCov / (s_FmaxFC * FmaxFC_)) * s_maxZr * (maxZr_ - minZr_)\n",
+    "        \n",
+    "        # Water capacities\n",
+    "        TAW = (soil_params['FC'] - soil_params['WP']) * Zr\n",
+    "        TDW = (soil_params['FC'] - soil_params['WP']) * (s_Zsoil * Zsoil_ - Zr)\n",
+    "        \n",
+    "        # Update depletions\n",
+    "        Dr = update_Dr(TAW, TDW, Zr, TAW0, TDW0, Dr0, Dd0, Zr0)\n",
+    "        Dd = update_Dd(TAW, TDW, Zr, TAW0, TDW0, Dd0, Zr0)\n",
+    "        \n",
+    "        # Update param p\n",
+    "        p_ = (xr_minimum(xr_maximum(s_p * p_ + 0.04 * (5 - weather_cube['ET0'].sel({'time': dates[i-1]}) / 1000), 0.1), 0.8) * (1 / s_p)).round(0).astype('i2')\n",
+    "        \n",
+    "        # Irrigation   ==============!!!!!\n",
+    "        Irrig = xr_minimum(xr_maximum(Dr - weather_cube['tp'].sel({'time': dates[i]}) / 1000, 0), s_Lame_max * Lame_max_) * Irrig_auto_\n",
+    "        Irrig = xr.where(Dr > TAW * s_p * p_, Irrig, 0)\n",
+    "    \n",
+    "        # Kcb\n",
+    "        Kcb = xr_minimum(s_Kslope * Kslope_ * (ndvi_cube['ndvi'].sel({'time': dates[i]}) / 255) + s_Koffset * Koffset_, s_KmaxKcb * KmaxKcb_)\n",
+    "        \n",
+    "        # Update depletions with rainfall and/or irrigation    \n",
+    "        ## DP\n",
+    "        model_outputs['DP'].loc[{'time': dates[i]}] = -xr_minimum(Dd + xr_minimum(Dr - (weather_cube['tp'].sel({'time': dates[i]}) / 1000) - Irrig, 0), 0)\n",
+    "        \n",
+    "        ## De\n",
+    "        Dei = xr_minimum(xr_maximum(Dei - (weather_cube['tp'].sel({'time': dates[i]}) / 1000) - Irrig / (s_FW * FW_ / 100), 0), TEW)\n",
+    "        Dep = xr_minimum(xr_maximum(Dep - (weather_cube['tp'].sel({'time': dates[i]}) / 1000), 0), TEW)\n",
+    "        \n",
+    "        fewi = xr_minimum(1 - FCov, (s_FW * FW_ / 100))\n",
+    "        fewp = 1 - FCov - fewi\n",
+    "        \n",
+    "        De = (Dei * fewi + Dep * fewp) / (fewi + fewp)\n",
+    "        # variables_t1['De'] = xr.where(variables_t1['De'].isfinite(), variables_t1['De'], variables_t1['Dei'] * (s_FW * FW_ / 100) + variables_t1['Dep'] * (1 - (s_FW * FW_ / 100)))  #================= find replacement for .isfinite() method !!!!!!!!!\n",
+    "\n",
+    "        ## Dr\n",
+    "        Dr = xr_minimum(xr_maximum(Dr - (weather_cube['tp'].sel({'time': dates[i]}) / 1000) - Irrig, 0), TAW)\n",
+    "        \n",
+    "        ## Dd\n",
+    "        Dd = xr_minimum(xr_maximum(Dd + xr_minimum(Dr - (weather_cube['tp'].sel({'time': dates[i]}) / 1000) - Irrig, 0), 0), TDW)\n",
+    "        \n",
+    "        # Diffusion coefficients\n",
+    "        diff_rei = calculate_diff_re(TAW, Dr, Zr, RUE, Dei, FCov, Ze_, DiffE_, scale_factor)\n",
+    "        diff_rep = calculate_diff_re(TAW, Dr, Zr, RUE, Dep, FCov, Ze_, DiffE_, scale_factor)\n",
+    "        diff_dr = calculate_diff_dr(TAW, TDW, Dr, Zr, Dd, FCov, Zsoil_, DiffR_, scale_factor) \n",
+    "        \n",
+    "        # Weighing factor W\n",
+    "        W = calculate_W(TEW, Dei, Dep, fewi, fewp)\n",
+    "        \n",
+    "        # Soil water contents\n",
+    "        model_outputs['SWCe'].loc[{'time': dates[i]}] = 1 - De/TEW\n",
+    "        model_outputs['SWCr'].loc[{'time': dates[i]}] = 1 - Dr/TAW\n",
+    "        \n",
+    "        # Water Stress coefficient\n",
+    "        Ks = xr_minimum((TAW - Dr) / (TAW * (1 - s_p * p_)), 1)\n",
+    "        \n",
+    "        # Reduction coefficient for evaporation\n",
+    "        Kei = xr_minimum(W * calculate_Kr(TEW, Dei, REW_, scale_factor) * (s_Kcmax * Kcmax_ - Kcb), fewi * s_Kcmax * Kcmax_)\n",
+    "        Kep = xr_minimum((1 - W) * calculate_Kr(TEW, Dei, REW_, scale_factor) * (s_Kcmax * Kcmax_ - Kcb), fewp * s_Kcmax * Kcmax_)\n",
+    "        \n",
+    "        # Prepare coefficients for evapotranspiration\n",
+    "        Kti = xr_minimum(((s_Ze * Ze_ / Zr)**6) * (1 - Dei / TEW) / xr_maximum(1 - Dr / TAW, 0.001), 1)\n",
+    "        Ktp = xr_minimum(((s_Ze * Ze_ / Zr)**6) * (1 - Dep / TEW) / xr_maximum(1 - Dr / TAW, 0.001), 1)\n",
+    "        Tei = Kti * Ks * Kcb * weather_cube['ET0'].sel({'time': dates[i]}) / 1000\n",
+    "        Tep = Ktp * Ks * Kcb * weather_cube['ET0'].sel({'time': dates[i]}) / 1000\n",
+    "        \n",
+    "        # Update depletions\n",
+    "        Dei = xr.where(fewi > 0, xr_minimum(xr_maximum(Dei + (weather_cube['ET0'].sel({'time': dates[i]}) / 1000) * Kei / fewi + Tei - diff_rei, 0), TEW), xr_minimum(xr_maximum(Dei + Tei - diff_rei, 0), TEW))\n",
+    "        Dep = xr.where(fewp > 0, xr_minimum(xr_maximum(Dep + (weather_cube['ET0'].sel({'time': dates[i]}) / 1000) * Kep / fewp + Tep - diff_rep, 0), TEW), xr_minimum(xr_maximum(Dep + Tep - diff_rep, 0), TEW))\n",
+    "        \n",
+    "        De = (Dei * fewi + Dep * fewp) / (fewi + fewp)\n",
+    "        # De = xr.where(De.isfinite(), De, Dei * (s_FW * FW_ / 100) + Dep * (1 - (s_FW * FW_ / 100)))  #================= find replacement for .isfinite() method !!!!!!!!!\n",
+    "        \n",
+    "        # Evaporation\n",
+    "        model_outputs['E'].loc[{'time': dates[i]}]  = xr_maximum((Kei + Kep) * weather_cube['ET0'].sel({'time': dates[i]}) / 1000, 0)\n",
+    "        \n",
+    "        # Transpiration\n",
+    "        model_outputs['Tr'].loc[{'time': dates[i]}]  = Kcb * Ks * weather_cube['ET0'].sel({'time': dates[i]}) / 1000\n",
+    "        \n",
+    "        # Potential evapotranspiration and evaporative fraction ??\n",
+    "        \n",
+    "        # Update depletions (root and deep zones) at the end of the day\n",
+    "        Dr = xr_minimum(xr_maximum(Dr + model_outputs['E'].loc[{'time': dates[i]}] + model_outputs['Tr'].loc[{'time': dates[i]}] - diff_dr, 0), TAW)\n",
+    "        Dd = xr_minimum(xr_maximum(Dd + diff_dr, 0), TDW)\n",
+    "        \n",
+    "        # Write outputs\n",
+    "        model_outputs['Irr'].loc[{'time': dates[i]}] = Irrig\n",
+    "        \n",
+    "        # Update variable_t1 values\n",
+    "        for variable in calculation_variables_t1:\n",
+    "            variables_t1[variable] = variables_t2[variable].copy(deep = True)\n",
+    "        \n",
+    "        print('day ', i+1, '/', len(dates), '   ', end = '\\r')\n",
+    "    \n",
+    "    # Scale the model_outputs variable to save in int16 format\n",
+    "    model_outputs = model_outputs * 1000\n",
+    "    \n",
+    "    # Write encoding dict\n",
+    "    encoding_dict = {}\n",
+    "    for variable in list(model_outputs.keys()):\n",
+    "        encod = {}\n",
+    "        encod['dtype'] = 'i2'\n",
+    "        encoding_dict[variable] = encod\n",
+    "    \n",
+    "    # Save model outputs to netcdf\n",
+    "    model_outputs.to_netcdf(save_path, encoding = encoding_dict)\n",
+    "    \n",
+    "    return None"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/auclairj/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/distributed/node.py:182: UserWarning: Port 8787 is already in use.\n",
+      "Perhaps you already have a cluster running?\n",
+      "Hosting the HTTP server on port 38039 instead\n",
+      "  warnings.warn(\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "    <div style=\"width: 24px; height: 24px; background-color: #e1e1e1; border: 3px solid #9D9D9D; border-radius: 5px; position: absolute;\"> </div>\n",
+       "    <div style=\"margin-left: 48px;\">\n",
+       "        <h3 style=\"margin-bottom: 0px;\">Client</h3>\n",
+       "        <p style=\"color: #9D9D9D; margin-bottom: 0px;\">Client-03e36644-264a-11ee-9c09-00155d33b451</p>\n",
+       "        <table style=\"width: 100%; text-align: left;\">\n",
+       "\n",
+       "        <tr>\n",
+       "        \n",
+       "            <td style=\"text-align: left;\"><strong>Connection method:</strong> Cluster object</td>\n",
+       "            <td style=\"text-align: left;\"><strong>Cluster type:</strong> distributed.LocalCluster</td>\n",
+       "        \n",
+       "        </tr>\n",
+       "\n",
+       "        \n",
+       "            <tr>\n",
+       "                <td style=\"text-align: left;\">\n",
+       "                    <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:38039/status\" target=\"_blank\">http://127.0.0.1:38039/status</a>\n",
+       "                </td>\n",
+       "                <td style=\"text-align: left;\"></td>\n",
+       "            </tr>\n",
+       "        \n",
+       "\n",
+       "        </table>\n",
+       "\n",
+       "        \n",
+       "\n",
+       "        \n",
+       "            <details>\n",
+       "            <summary style=\"margin-bottom: 20px;\"><h3 style=\"display: inline;\">Cluster Info</h3></summary>\n",
+       "            <div class=\"jp-RenderedHTMLCommon jp-RenderedHTML jp-mod-trusted jp-OutputArea-output\">\n",
+       "    <div style=\"width: 24px; height: 24px; background-color: #e1e1e1; border: 3px solid #9D9D9D; border-radius: 5px; position: absolute;\">\n",
+       "    </div>\n",
+       "    <div style=\"margin-left: 48px;\">\n",
+       "        <h3 style=\"margin-bottom: 0px; margin-top: 0px;\">LocalCluster</h3>\n",
+       "        <p style=\"color: #9D9D9D; margin-bottom: 0px;\">a5d645bb</p>\n",
+       "        <table style=\"width: 100%; text-align: left;\">\n",
+       "            <tr>\n",
+       "                <td style=\"text-align: left;\">\n",
+       "                    <strong>Dashboard:</strong> <a href=\"http://127.0.0.1:38039/status\" target=\"_blank\">http://127.0.0.1:38039/status</a>\n",
+       "                </td>\n",
+       "                <td style=\"text-align: left;\">\n",
+       "                    <strong>Workers:</strong> 4\n",
+       "                </td>\n",
+       "            </tr>\n",
+       "            <tr>\n",
+       "                <td style=\"text-align: left;\">\n",
+       "                    <strong>Total threads:</strong> 8\n",
+       "                </td>\n",
+       "                <td style=\"text-align: left;\">\n",
+       "                    <strong>Total memory:</strong> 23.47 GiB\n",
+       "                </td>\n",
+       "            </tr>\n",
+       "            \n",
+       "            <tr>\n",
+       "    <td style=\"text-align: left;\"><strong>Status:</strong> running</td>\n",
+       "    <td style=\"text-align: left;\"><strong>Using processes:</strong> True</td>\n",
+       "</tr>\n",
+       "\n",
+       "            \n",
+       "        </table>\n",
+       "\n",
+       "        <details>\n",
+       "            <summary style=\"margin-bottom: 20px;\">\n",
+       "                <h3 style=\"display: inline;\">Scheduler Info</h3>\n",
+       "            </summary>\n",
+       "\n",
+       "            <div style=\"\">\n",
+       "    <div>\n",
+       "        <div style=\"width: 24px; height: 24px; background-color: #FFF7E5; border: 3px solid #FF6132; border-radius: 5px; position: absolute;\"> </div>\n",
+       "        <div style=\"margin-left: 48px;\">\n",
+       "            <h3 style=\"margin-bottom: 0px;\">Scheduler</h3>\n",
+       "            <p style=\"color: #9D9D9D; margin-bottom: 0px;\">Scheduler-2d9db510-69be-4ff5-86b0-f61af0227954</p>\n",
+       "            <table style=\"width: 100%; text-align: left;\">\n",
+       "                <tr>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Comm:</strong> tcp://127.0.0.1:37671\n",
+       "                    </td>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Workers:</strong> 4\n",
+       "                    </td>\n",
+       "                </tr>\n",
+       "                <tr>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Dashboard:</strong> <a href=\"http://127.0.0.1:38039/status\" target=\"_blank\">http://127.0.0.1:38039/status</a>\n",
+       "                    </td>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Total threads:</strong> 8\n",
+       "                    </td>\n",
+       "                </tr>\n",
+       "                <tr>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Started:</strong> Just now\n",
+       "                    </td>\n",
+       "                    <td style=\"text-align: left;\">\n",
+       "                        <strong>Total memory:</strong> 23.47 GiB\n",
+       "                    </td>\n",
+       "                </tr>\n",
+       "            </table>\n",
+       "        </div>\n",
+       "    </div>\n",
+       "\n",
+       "    <details style=\"margin-left: 48px;\">\n",
+       "        <summary style=\"margin-bottom: 20px;\">\n",
+       "            <h3 style=\"display: inline;\">Workers</h3>\n",
+       "        </summary>\n",
+       "\n",
+       "        \n",
+       "        <div style=\"margin-bottom: 20px;\">\n",
+       "            <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
+       "            <div style=\"margin-left: 48px;\">\n",
+       "            <details>\n",
+       "                <summary>\n",
+       "                    <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 0</h4>\n",
+       "                </summary>\n",
+       "                <table style=\"width: 100%; text-align: left;\">\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Comm: </strong> tcp://127.0.0.1:45285\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Total threads: </strong> 2\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:41621/status\" target=\"_blank\">http://127.0.0.1:41621/status</a>\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Memory: </strong> 5.87 GiB\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Nanny: </strong> tcp://127.0.0.1:36255\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\"></td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td colspan=\"2\" style=\"text-align: left;\">\n",
+       "                            <strong>Local directory: </strong> /tmp/dask-scratch-space/worker-jebkacuo\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                </table>\n",
+       "            </details>\n",
+       "            </div>\n",
+       "        </div>\n",
+       "        \n",
+       "        <div style=\"margin-bottom: 20px;\">\n",
+       "            <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
+       "            <div style=\"margin-left: 48px;\">\n",
+       "            <details>\n",
+       "                <summary>\n",
+       "                    <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 1</h4>\n",
+       "                </summary>\n",
+       "                <table style=\"width: 100%; text-align: left;\">\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Comm: </strong> tcp://127.0.0.1:42337\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Total threads: </strong> 2\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:37155/status\" target=\"_blank\">http://127.0.0.1:37155/status</a>\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Memory: </strong> 5.87 GiB\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Nanny: </strong> tcp://127.0.0.1:36249\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\"></td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td colspan=\"2\" style=\"text-align: left;\">\n",
+       "                            <strong>Local directory: </strong> /tmp/dask-scratch-space/worker-iprr4c5j\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                </table>\n",
+       "            </details>\n",
+       "            </div>\n",
+       "        </div>\n",
+       "        \n",
+       "        <div style=\"margin-bottom: 20px;\">\n",
+       "            <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
+       "            <div style=\"margin-left: 48px;\">\n",
+       "            <details>\n",
+       "                <summary>\n",
+       "                    <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 2</h4>\n",
+       "                </summary>\n",
+       "                <table style=\"width: 100%; text-align: left;\">\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Comm: </strong> tcp://127.0.0.1:38697\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Total threads: </strong> 2\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:34703/status\" target=\"_blank\">http://127.0.0.1:34703/status</a>\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Memory: </strong> 5.87 GiB\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Nanny: </strong> tcp://127.0.0.1:40667\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\"></td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td colspan=\"2\" style=\"text-align: left;\">\n",
+       "                            <strong>Local directory: </strong> /tmp/dask-scratch-space/worker-92iu95t2\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                </table>\n",
+       "            </details>\n",
+       "            </div>\n",
+       "        </div>\n",
+       "        \n",
+       "        <div style=\"margin-bottom: 20px;\">\n",
+       "            <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
+       "            <div style=\"margin-left: 48px;\">\n",
+       "            <details>\n",
+       "                <summary>\n",
+       "                    <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 3</h4>\n",
+       "                </summary>\n",
+       "                <table style=\"width: 100%; text-align: left;\">\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Comm: </strong> tcp://127.0.0.1:32795\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Total threads: </strong> 2\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:41767/status\" target=\"_blank\">http://127.0.0.1:41767/status</a>\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Memory: </strong> 5.87 GiB\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td style=\"text-align: left;\">\n",
+       "                            <strong>Nanny: </strong> tcp://127.0.0.1:35569\n",
+       "                        </td>\n",
+       "                        <td style=\"text-align: left;\"></td>\n",
+       "                    </tr>\n",
+       "                    <tr>\n",
+       "                        <td colspan=\"2\" style=\"text-align: left;\">\n",
+       "                            <strong>Local directory: </strong> /tmp/dask-scratch-space/worker-bx7cwqe3\n",
+       "                        </td>\n",
+       "                    </tr>\n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                    \n",
+       "\n",
+       "                </table>\n",
+       "            </details>\n",
+       "            </div>\n",
+       "        </div>\n",
+       "        \n",
+       "\n",
+       "    </details>\n",
+       "</div>\n",
+       "\n",
+       "        </details>\n",
+       "    </div>\n",
+       "</div>\n",
+       "            </details>\n",
+       "        \n",
+       "\n",
+       "    </div>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "<Client: 'tcp://127.0.0.1:37671' processes=4 threads=8, memory=23.47 GiB>"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "client = Client()\n",
+    "client"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "2023-07-19 17:36:42,921 - distributed.scheduler - ERROR - Couldn't gather keys {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": []} state: ['waiting'] workers: []\n",
+      "NoneType: None\n",
+      "2023-07-19 17:36:42,922 - distributed.scheduler - ERROR - Shut down workers that don't have promised key: [], ('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\n",
+      "NoneType: None\n",
+      "2023-07-19 17:36:42,924 - distributed.client - WARNING - Couldn't gather 1 keys, rescheduling {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": ()}\n",
+      "2023-07-19 17:36:43,297 - distributed.scheduler - ERROR - Couldn't gather keys {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": []} state: [None] workers: []\n",
+      "NoneType: None\n",
+      "2023-07-19 17:36:43,298 - distributed.scheduler - ERROR - Shut down workers that don't have promised key: [], ('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\n",
+      "NoneType: None\n",
+      "2023-07-19 17:36:43,300 - distributed.client - WARNING - Couldn't gather 1 keys, rescheduling {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": ()}\n",
+      "2023-07-19 17:36:43,454 - distributed.scheduler - ERROR - Couldn't gather keys {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": []} state: [None] workers: []\n",
+      "NoneType: None\n",
+      "2023-07-19 17:36:43,455 - distributed.scheduler - ERROR - Shut down workers that don't have promised key: [], ('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\n",
+      "NoneType: None\n",
+      "2023-07-19 17:36:43,456 - distributed.client - WARNING - Couldn't gather 1 keys, rescheduling {\"('astype-1dead4f4f28400d17d384d6a2b513c87', 0, 0)\": ()}\n",
+      "/home/auclairj/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in divide\n",
+      "  return func(*(_execute_task(a, cache) for a in args))\n",
+      "/home/auclairj/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in divide\n",
+      "  return func(*(_execute_task(a, cache) for a in args))\n",
+      "/home/auclairj/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in divide\n",
+      "  return func(*(_execute_task(a, cache) for a in args))\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "day  2 / 366    \r"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/auclairj/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in divide\n",
+      "  return func(*(_execute_task(a, cache) for a in args))\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "day  42 / 366    \r"
+     ]
+    },
+    {
+     "ename": "KeyboardInterrupt",
+     "evalue": "",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[4], line 14\u001b[0m\n\u001b[1;32m      9\u001b[0m save_path \u001b[39m=\u001b[39m data_path \u001b[39m+\u001b[39m os\u001b[39m.\u001b[39msep \u001b[39m+\u001b[39m \u001b[39m'\u001b[39m\u001b[39moutputs.nc\u001b[39m\u001b[39m'\u001b[39m\n\u001b[1;32m     11\u001b[0m chunk_size \u001b[39m=\u001b[39m {\u001b[39m'\u001b[39m\u001b[39mx\u001b[39m\u001b[39m'\u001b[39m: \u001b[39m-\u001b[39m\u001b[39m1\u001b[39m, \u001b[39m'\u001b[39m\u001b[39my\u001b[39m\u001b[39m'\u001b[39m: \u001b[39m-\u001b[39m\u001b[39m1\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mtime\u001b[39m\u001b[39m'\u001b[39m: \u001b[39m-\u001b[39m\u001b[39m1\u001b[39m}\n\u001b[0;32m---> 14\u001b[0m run_samir(json_config_file, param_file, ndvi_path, weather_path, soil_path, land_cover_path, chunk_size, save_path)\n",
+      "Cell \u001b[0;32mIn[2], line 734\u001b[0m, in \u001b[0;36mrun_samir\u001b[0;34m(json_config_file, csv_param_file, ndvi_cube_path, weather_cube_path, soil_params_path, land_cover_path, chunk_size, save_path)\u001b[0m\n\u001b[1;32m    730\u001b[0m De \u001b[39m=\u001b[39m (Dei \u001b[39m*\u001b[39m fewi \u001b[39m+\u001b[39m Dep \u001b[39m*\u001b[39m fewp) \u001b[39m/\u001b[39m (fewi \u001b[39m+\u001b[39m fewp)\n\u001b[1;32m    731\u001b[0m \u001b[39m# De = xr.where(De.isfinite(), De, Dei * (s_FW * FW_ / 100) + Dep * (1 - (s_FW * FW_ / 100)))\u001b[39;00m\n\u001b[1;32m    732\u001b[0m \n\u001b[1;32m    733\u001b[0m \u001b[39m# Evaporation\u001b[39;00m\n\u001b[0;32m--> 734\u001b[0m model_outputs[\u001b[39m'\u001b[39m\u001b[39mE\u001b[39m\u001b[39m'\u001b[39m]\u001b[39m.\u001b[39mloc[{\u001b[39m'\u001b[39m\u001b[39mtime\u001b[39m\u001b[39m'\u001b[39m: dates[i]}]  \u001b[39m=\u001b[39m xr_maximum((Kei \u001b[39m+\u001b[39m Kep) \u001b[39m*\u001b[39m weather_cube[\u001b[39m'\u001b[39m\u001b[39mET0\u001b[39m\u001b[39m'\u001b[39m]\u001b[39m.\u001b[39msel({\u001b[39m'\u001b[39m\u001b[39mtime\u001b[39m\u001b[39m'\u001b[39m: dates[i]}) \u001b[39m/\u001b[39m \u001b[39m1000\u001b[39m, \u001b[39m0\u001b[39m)\n\u001b[1;32m    736\u001b[0m \u001b[39m# Transpiration\u001b[39;00m\n\u001b[1;32m    737\u001b[0m model_outputs[\u001b[39m'\u001b[39m\u001b[39mTr\u001b[39m\u001b[39m'\u001b[39m]\u001b[39m.\u001b[39mloc[{\u001b[39m'\u001b[39m\u001b[39mtime\u001b[39m\u001b[39m'\u001b[39m: dates[i]}]  \u001b[39m=\u001b[39m Kcb \u001b[39m*\u001b[39m Ks \u001b[39m*\u001b[39m weather_cube[\u001b[39m'\u001b[39m\u001b[39mET0\u001b[39m\u001b[39m'\u001b[39m]\u001b[39m.\u001b[39msel({\u001b[39m'\u001b[39m\u001b[39mtime\u001b[39m\u001b[39m'\u001b[39m: dates[i]}) \u001b[39m/\u001b[39m \u001b[39m1000\u001b[39m\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/xarray/core/dataarray.py:223\u001b[0m, in \u001b[0;36m_LocIndexer.__setitem__\u001b[0;34m(self, key, value)\u001b[0m\n\u001b[1;32m    220\u001b[0m     key \u001b[39m=\u001b[39m \u001b[39mdict\u001b[39m(\u001b[39mzip\u001b[39m(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdata_array\u001b[39m.\u001b[39mdims, labels))\n\u001b[1;32m    222\u001b[0m dim_indexers \u001b[39m=\u001b[39m map_index_queries(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdata_array, key)\u001b[39m.\u001b[39mdim_indexers\n\u001b[0;32m--> 223\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mdata_array[dim_indexers] \u001b[39m=\u001b[39m value\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/xarray/core/dataarray.py:840\u001b[0m, in \u001b[0;36mDataArray.__setitem__\u001b[0;34m(self, key, value)\u001b[0m\n\u001b[1;32m    835\u001b[0m \u001b[39m# DataArray key -> Variable key\u001b[39;00m\n\u001b[1;32m    836\u001b[0m key \u001b[39m=\u001b[39m {\n\u001b[1;32m    837\u001b[0m     k: v\u001b[39m.\u001b[39mvariable \u001b[39mif\u001b[39;00m \u001b[39misinstance\u001b[39m(v, DataArray) \u001b[39melse\u001b[39;00m v\n\u001b[1;32m    838\u001b[0m     \u001b[39mfor\u001b[39;00m k, v \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_item_key_to_dict(key)\u001b[39m.\u001b[39mitems()\n\u001b[1;32m    839\u001b[0m }\n\u001b[0;32m--> 840\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mvariable[key] \u001b[39m=\u001b[39m value\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/xarray/core/variable.py:977\u001b[0m, in \u001b[0;36mVariable.__setitem__\u001b[0;34m(self, key, value)\u001b[0m\n\u001b[1;32m    974\u001b[0m     value \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mmoveaxis(value, new_order, \u001b[39mrange\u001b[39m(\u001b[39mlen\u001b[39m(new_order)))\n\u001b[1;32m    976\u001b[0m indexable \u001b[39m=\u001b[39m as_indexable(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_data)\n\u001b[0;32m--> 977\u001b[0m indexable[index_tuple] \u001b[39m=\u001b[39m value\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/xarray/core/indexing.py:1338\u001b[0m, in \u001b[0;36mNumpyIndexingAdapter.__setitem__\u001b[0;34m(self, key, value)\u001b[0m\n\u001b[1;32m   1336\u001b[0m array, key \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_indexing_array_and_key(key)\n\u001b[1;32m   1337\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m-> 1338\u001b[0m     array[key] \u001b[39m=\u001b[39m value\n\u001b[1;32m   1339\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mValueError\u001b[39;00m:\n\u001b[1;32m   1340\u001b[0m     \u001b[39m# More informative exception if read-only view\u001b[39;00m\n\u001b[1;32m   1341\u001b[0m     \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m array\u001b[39m.\u001b[39mflags\u001b[39m.\u001b[39mwriteable \u001b[39mand\u001b[39;00m \u001b[39mnot\u001b[39;00m array\u001b[39m.\u001b[39mflags\u001b[39m.\u001b[39mowndata:\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/array/core.py:1699\u001b[0m, in \u001b[0;36mArray.__array__\u001b[0;34m(self, dtype, **kwargs)\u001b[0m\n\u001b[1;32m   1698\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m__array__\u001b[39m(\u001b[39mself\u001b[39m, dtype\u001b[39m=\u001b[39m\u001b[39mNone\u001b[39;00m, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs):\n\u001b[0;32m-> 1699\u001b[0m     x \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mcompute()\n\u001b[1;32m   1700\u001b[0m     \u001b[39mif\u001b[39;00m dtype \u001b[39mand\u001b[39;00m x\u001b[39m.\u001b[39mdtype \u001b[39m!=\u001b[39m dtype:\n\u001b[1;32m   1701\u001b[0m         x \u001b[39m=\u001b[39m x\u001b[39m.\u001b[39mastype(dtype)\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/base.py:381\u001b[0m, in \u001b[0;36mDaskMethodsMixin.compute\u001b[0;34m(self, **kwargs)\u001b[0m\n\u001b[1;32m    357\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mcompute\u001b[39m(\u001b[39mself\u001b[39m, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs):\n\u001b[1;32m    358\u001b[0m     \u001b[39m\"\"\"Compute this dask collection\u001b[39;00m\n\u001b[1;32m    359\u001b[0m \n\u001b[1;32m    360\u001b[0m \u001b[39m    This turns a lazy Dask collection into its in-memory equivalent.\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    379\u001b[0m \u001b[39m    dask.compute\u001b[39;00m\n\u001b[1;32m    380\u001b[0m \u001b[39m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 381\u001b[0m     (result,) \u001b[39m=\u001b[39m compute(\u001b[39mself\u001b[39;49m, traverse\u001b[39m=\u001b[39;49m\u001b[39mFalse\u001b[39;49;00m, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m    382\u001b[0m     \u001b[39mreturn\u001b[39;00m result\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/base.py:660\u001b[0m, in \u001b[0;36mcompute\u001b[0;34m(traverse, optimize_graph, scheduler, get, *args, **kwargs)\u001b[0m\n\u001b[1;32m    652\u001b[0m     \u001b[39mreturn\u001b[39;00m args\n\u001b[1;32m    654\u001b[0m schedule \u001b[39m=\u001b[39m get_scheduler(\n\u001b[1;32m    655\u001b[0m     scheduler\u001b[39m=\u001b[39mscheduler,\n\u001b[1;32m    656\u001b[0m     collections\u001b[39m=\u001b[39mcollections,\n\u001b[1;32m    657\u001b[0m     get\u001b[39m=\u001b[39mget,\n\u001b[1;32m    658\u001b[0m )\n\u001b[0;32m--> 660\u001b[0m dsk \u001b[39m=\u001b[39m collections_to_dsk(collections, optimize_graph, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m    661\u001b[0m keys, postcomputes \u001b[39m=\u001b[39m [], []\n\u001b[1;32m    662\u001b[0m \u001b[39mfor\u001b[39;00m x \u001b[39min\u001b[39;00m collections:\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/base.py:433\u001b[0m, in \u001b[0;36mcollections_to_dsk\u001b[0;34m(collections, optimize_graph, optimizations, **kwargs)\u001b[0m\n\u001b[1;32m    431\u001b[0m \u001b[39mfor\u001b[39;00m opt, val \u001b[39min\u001b[39;00m groups\u001b[39m.\u001b[39mitems():\n\u001b[1;32m    432\u001b[0m     dsk, keys \u001b[39m=\u001b[39m _extract_graph_and_keys(val)\n\u001b[0;32m--> 433\u001b[0m     dsk \u001b[39m=\u001b[39m opt(dsk, keys, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m    435\u001b[0m     \u001b[39mfor\u001b[39;00m opt_inner \u001b[39min\u001b[39;00m optimizations:\n\u001b[1;32m    436\u001b[0m         dsk \u001b[39m=\u001b[39m opt_inner(dsk, keys, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs)\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/array/optimization.py:49\u001b[0m, in \u001b[0;36moptimize\u001b[0;34m(dsk, keys, fuse_keys, fast_functions, inline_functions_fast_functions, rename_fused_keys, **kwargs)\u001b[0m\n\u001b[1;32m     46\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39misinstance\u001b[39m(dsk, HighLevelGraph):\n\u001b[1;32m     47\u001b[0m     dsk \u001b[39m=\u001b[39m HighLevelGraph\u001b[39m.\u001b[39mfrom_collections(\u001b[39mid\u001b[39m(dsk), dsk, dependencies\u001b[39m=\u001b[39m())\n\u001b[0;32m---> 49\u001b[0m dsk \u001b[39m=\u001b[39m optimize_blockwise(dsk, keys\u001b[39m=\u001b[39;49mkeys)\n\u001b[1;32m     50\u001b[0m dsk \u001b[39m=\u001b[39m fuse_roots(dsk, keys\u001b[39m=\u001b[39mkeys)\n\u001b[1;32m     51\u001b[0m dsk \u001b[39m=\u001b[39m dsk\u001b[39m.\u001b[39mcull(\u001b[39mset\u001b[39m(keys))\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/blockwise.py:1080\u001b[0m, in \u001b[0;36moptimize_blockwise\u001b[0;34m(graph, keys)\u001b[0m\n\u001b[1;32m   1078\u001b[0m \u001b[39mwhile\u001b[39;00m out\u001b[39m.\u001b[39mdependencies \u001b[39m!=\u001b[39m graph\u001b[39m.\u001b[39mdependencies:\n\u001b[1;32m   1079\u001b[0m     graph \u001b[39m=\u001b[39m out\n\u001b[0;32m-> 1080\u001b[0m     out \u001b[39m=\u001b[39m _optimize_blockwise(graph, keys\u001b[39m=\u001b[39;49mkeys)\n\u001b[1;32m   1081\u001b[0m \u001b[39mreturn\u001b[39;00m out\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/blockwise.py:1154\u001b[0m, in \u001b[0;36m_optimize_blockwise\u001b[0;34m(full_graph, keys)\u001b[0m\n\u001b[1;32m   1151\u001b[0m             stack\u001b[39m.\u001b[39mappend(d)\n\u001b[1;32m   1153\u001b[0m \u001b[39m# Merge these Blockwise layers into one\u001b[39;00m\n\u001b[0;32m-> 1154\u001b[0m new_layer \u001b[39m=\u001b[39m rewrite_blockwise([layers[l] \u001b[39mfor\u001b[39;49;00m l \u001b[39min\u001b[39;49;00m blockwise_layers])\n\u001b[1;32m   1155\u001b[0m out[layer] \u001b[39m=\u001b[39m new_layer\n\u001b[1;32m   1157\u001b[0m \u001b[39m# Get the new (external) dependencies for this layer.\u001b[39;00m\n\u001b[1;32m   1158\u001b[0m \u001b[39m# This corresponds to the dependencies defined in\u001b[39;00m\n\u001b[1;32m   1159\u001b[0m \u001b[39m# full_graph.dependencies and are not in blockwise_layers\u001b[39;00m\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/blockwise.py:1341\u001b[0m, in \u001b[0;36mrewrite_blockwise\u001b[0;34m(inputs)\u001b[0m\n\u001b[1;32m   1339\u001b[0m sub \u001b[39m=\u001b[39m {}\n\u001b[1;32m   1340\u001b[0m \u001b[39m# Map from (id(key), inds or None) -> index in indices. Used to deduplicate indices.\u001b[39;00m\n\u001b[0;32m-> 1341\u001b[0m index_map \u001b[39m=\u001b[39m {(\u001b[39mid\u001b[39m(k), inds): n \u001b[39mfor\u001b[39;00m n, (k, inds) \u001b[39min\u001b[39;00m \u001b[39menumerate\u001b[39m(indices)}\n\u001b[1;32m   1342\u001b[0m \u001b[39mfor\u001b[39;00m ii, index \u001b[39min\u001b[39;00m \u001b[39menumerate\u001b[39m(new_indices):\n\u001b[1;32m   1343\u001b[0m     id_key \u001b[39m=\u001b[39m (\u001b[39mid\u001b[39m(index[\u001b[39m0\u001b[39m]), index[\u001b[39m1\u001b[39m])\n",
+      "File \u001b[0;32m~/anaconda3/envs/modspa_pixel/lib/python3.10/site-packages/dask/blockwise.py:1341\u001b[0m, in \u001b[0;36m<dictcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m   1339\u001b[0m sub \u001b[39m=\u001b[39m {}\n\u001b[1;32m   1340\u001b[0m \u001b[39m# Map from (id(key), inds or None) -> index in indices. Used to deduplicate indices.\u001b[39;00m\n\u001b[0;32m-> 1341\u001b[0m index_map \u001b[39m=\u001b[39m {(\u001b[39mid\u001b[39;49m(k), inds): n \u001b[39mfor\u001b[39;00m n, (k, inds) \u001b[39min\u001b[39;00m \u001b[39menumerate\u001b[39m(indices)}\n\u001b[1;32m   1342\u001b[0m \u001b[39mfor\u001b[39;00m ii, index \u001b[39min\u001b[39;00m \u001b[39menumerate\u001b[39m(new_indices):\n\u001b[1;32m   1343\u001b[0m     id_key \u001b[39m=\u001b[39m (\u001b[39mid\u001b[39m(index[\u001b[39m0\u001b[39m]), index[\u001b[39m1\u001b[39m])\n",
+      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
+     ]
+    }
+   ],
+   "source": [
+    "data_path = '/mnt/e/DATA/DEV_inputs_test'\n",
+    "\n",
+    "ndvi_path = data_path + os.sep + 'ndvi.nc'\n",
+    "weather_path = data_path + os.sep + 'weather.nc'\n",
+    "land_cover_path = data_path + os.sep + 'land_cover.nc'\n",
+    "json_config_file = '/home/auclairj/GIT/modspa-pixel/config/config_modspa.json'\n",
+    "param_file = '/home/auclairj/GIT/modspa-pixel/parameters/csv_files/params_samir_test.csv'\n",
+    "soil_path = data_path + os.sep + 'soil.nc'\n",
+    "save_path = data_path + os.sep + 'outputs.nc'\n",
+    "\n",
+    "chunk_size = {'x': 5, 'y': 5, 'time': -1}\n",
+    "\n",
+    "\n",
+    "run_samir(json_config_file, param_file, ndvi_path, weather_path, soil_path, land_cover_path, chunk_size, save_path)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "modspa_pixel",
+   "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.10.9"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/input/calculate_ndvi.py b/input/calculate_ndvi.py
index b24a581..0796bec 100644
--- a/input/calculate_ndvi.py
+++ b/input/calculate_ndvi.py
@@ -97,10 +97,9 @@ def calculate_ndvi(extracted_paths: Union[List[str], str], save_dir: str, bounda
     ndvi['ndvi'] = (ndvi.ndvi*255).sortby('time')
     
     # Write attributes
-    ndvi['ndvi'].attrs['long_name'] = 'Normalized Difference Vegetation Index'
     ndvi['ndvi'].attrs['units'] = 'None'
     ndvi['ndvi'].attrs['standard_name'] = 'NDVI'
-    ndvi['ndvi'].attrs['comment'] = 'Normalized difference of the near infrared and red band. A value of one is a high vegetation presence.'
+    ndvi['ndvi'].attrs['description'] = 'Normalized difference Vegetation Index (of the near infrared and red band). A value of one is a high vegetation presence.'
     ndvi['ndvi'].attrs['scale factor'] = '255'
     
     # Create save path
diff --git a/input/lib_era5_land_pixel.py b/input/lib_era5_land_pixel.py
index 0a8f7d3..93b12c3 100644
--- a/input/lib_era5_land_pixel.py
+++ b/input/lib_era5_land_pixel.py
@@ -429,7 +429,7 @@ def calculate_ET0_pixel(pixel_dataset: xr.Dataset, lat: float, lon: float, h: fl
     return ET0_values
 
 
-def era5Land_nc_daily_to_ET0(list_era5land_files: List[str], output_nc_file: str, ndvi_path: str, h: float = 10) -> None:
+def era5Land_nc_daily_to_ET0(list_era5land_files: List[str], output_nc_file: str, h: float = 10) -> None:
     """
     Calculate ET0 values from the ERA5 netcdf weather variables.
     Output netcdf contains the ET0 values for each day in the selected
@@ -491,4 +491,5 @@ def era5Land_nc_daily_to_ET0(list_era5land_files: List[str], output_nc_file: str
 
     # Save dataset to netcdf, still in wgs84 (lat, lon) coordinates
     final_weather_ds.to_netcdf(path = output_nc_file, encoding = {"ET0": {"dtype": "u2"}, "tp": {"dtype": "u2"}})
+    
     return None
\ No newline at end of file
diff --git a/parameters/csv_files/params_samir_test.csv b/parameters/csv_files/params_samir_test.csv
new file mode 100644
index 0000000..3ffa399
--- /dev/null
+++ b/parameters/csv_files/params_samir_test.csv
@@ -0,0 +1,24 @@
+ClassName,scale_factor,Default,class1,class2,class3,class4,class5,class6,class7
+ClassNumber,1,0,1,2,3,4,5,6,7
+FmaxFC,1000,1,0.9,1,1,1,1,,1
+Fslope,1000,1.4,1.5,1.5,1.5,1.2,1.5,0,1.5
+Foffset,1000,-0.1,-0.1,-0.1,-0.1,-0.17,-0.1,0,-0.1
+Plateau,1,70,,,,,,,
+KmaxKcb,1000,0.98,1,1.1,1.1,0.65,1.13,,1.13
+Kslope,1000,1.6,1.6,1.6,1.6,1.2,1.6,0,1.6
+Koffset,1000,-0.1,-0.1,-0.1,-0.1,-0.18,-0.1,0,-0.1
+Zsoil,1,2000,1600,1550,1550,2000,1550,2000,1550
+Ze,1,125,100,125,125,100,100,125,100
+Init_RU,1,0,,,,,,,
+DiffE,1,5,5,5,5,5,,,
+DiffR,1,10,10,10,10,10,,,
+REW,1,0,,,,,,,
+minZr,1,150,1550,125,125,1550,125,125,125
+maxZr,1,1000,1550,800,800,1550,1000,126,1000
+p,1000,0,0.5,0.55,0.55,0.65,0.7,0.45,0.7
+FW,1,100,30,100,100,30,100,30,100
+Irrig_auto,1,0,,,,,,0,
+Irrig_man,1,0,,,,,,,
+Lame_max,1,0,,,,,,,
+Kcmax,1000,1.15,,,,,,,
+Fc_stop,1000,0.95,,,,,,,
\ No newline at end of file
-- 
GitLab