From 206c8a172b0a3832e425e590d66fa4939a065b93 Mon Sep 17 00:00:00 2001 From: "Seb (MBP)" Date: Wed, 4 May 2022 10:27:19 +0200 Subject: [PATCH] v1.7.2 Added OpenStreetMap plotting --- CODE/VAR/prefs.mat | Bin 669 -> 685 bytes CODE/VAR/tephraProb.mat | Bin 2344 -> 2376 bytes CODE/conf_grid.m | 4 +- CODE/dependencies/plot_openstreetmap.m | 208 +++++++++++++++++++++++++ CODE/dwind/analyze_wind.m | 10 ++ CODE/dwind/dwind.m | 2 +- CODE/dwind/process_wind.m | 2 +- CODE/get_prefs.m | 25 ++- CODE/plotMap.m | 3 + CODE/runProb.m | 30 ++-- MODEL/forward_src/new_tephra.o | Bin 29855 -> 28264 bytes MODEL/forward_src/tephra2_calc.o | Bin 27107 -> 27064 bytes README.md | 11 ++ tephraProb.m | 3 +- 14 files changed, 281 insertions(+), 17 deletions(-) create mode 100644 CODE/dependencies/plot_openstreetmap.m diff --git a/CODE/VAR/prefs.mat b/CODE/VAR/prefs.mat index 663efed0b2526aae643d64199bbccd307e345967..388012fac97f6baec7e634c692f8ed6d9e96b363 100644 GIT binary patch delta 594 zcmV-Y01+4{;Hxx}_c_1J%ATcsJF)=zaG9WTAGBS}-BavVRk#rD$B?15dc%1E3 zzl#(x6wdC>F7Dwi7FTHRA0XOjH4CSlB4S~qHe+VU?vBnRA<1~B7Fr0F)>`=&SZil1 zXk}#~D1snZ_$QbjFT2@wVdo+Uf`RbzgIggb;eQh0s9TBYJj`8-2i?V46KH z$VA$8VmW62xZP(tLN(@p+BqpJB2l3MLL==v+i4dcLENK3+!{-nshJ>hW^?$piTeUA z;ttlt9R+cRMx&k&c3$p%x&JtQ{>+Tmsq=mY`xETW*x8d{&w;%J_J*@#=Ik2*_@T4> z3E-z=d;SeN`+?3b;H?vgn&_*|+H#J3avZIPVL$|23W1@eTan)A#Sp z_>r3Yw>IP-tmHp`736ojzVCd$L4VErs$Cz2eK3T5-|w+}uX6;t+w!}5<-fuIynmC| zi{vGLS?8bTdlc~6ys!Q}F_z}%%W0`3xTAD*J+HUmbyj%WE4y$t&ilQq|^8<`pTxRiHPUpMm05!=kM!3 gu+M+kx=jVX7sJaLC9!%7#(At>UeY_i0LsM7d2Ko{6#xJL delta 578 zcmV-I0=@mM1)T+uHxx>BZXh5xATc;PH846bG$1lCGBJ@+BavVRk#rD$6#@VNc%1E3 zy^a$x5S}EPh@80SiYuV!0T47$ZKMMoh$nz1cvnucuxnekmph<<;v`CHsCWW;I$9)< zkf`WDNJvN&JOch1Cz}u{d!{{0*5}WYZ)QB6bw&Unb^rzj&zac;Kl;c!jcNCUq*LYA ziPPBq3vQoh02L9&I4w$lsvtK6U}S99I_>-e`aR$9TM?xudM>G&yBJ|@{l0;v-@yjI zqlVw1&C$%K`)?1wKRHQXzp~?9>aE`p{~+#GBez_#}DaMH@+bMy66MH|0A8c@e}R8(oY}l_?_DL_qW6!ti`|F zi0^kppN0Oyep>W@QMuLx0EWL!KkdyRE)!R^gjF)1B$s=#u}8_q?HN^Dc&W z!6M6U6thArG)Lvp4Z7Y^*ICo;uIbjl%Va~G=KHaSTo7583b{s9z38rN;g-J z;&i-Syktd}K+bs?TOTE4!DjLGRdrQS8Lua`t}7O=KQjLvoFV>c>o$|*af~kGjK=1h Qb0K2$6$N|v7a%jp8|gnAumAu6 diff --git a/CODE/VAR/tephraProb.mat b/CODE/VAR/tephraProb.mat index 0ae031d474b9e2ad29bb26fb602b4f349267264e..99174481f0f3480cfa41a1f44a358399ec3740f1 100644 GIT binary patch delta 2298 zcmVvGdZ6jmU#{n5z)Z7H;!ve22<28`vHT>GOyQf^by3?WVb!A4q3JC$}d+LcNA1bv)7 zN*|+-(Q{{GaX>rT{Cdk3|y)_&dESVK1+#zCxxaEAra*NLniWu z3%@W0`P@_Z%B#v(>iMKb&(Dl?DDM>hL|uk5@~Dtwy!$nH7Y#hTe`7)uMwLbbTXB0X z93j&xPF)=O22Q`k*?;+*wBB$XLe`LoX~@GnM6iHrK4gRovPVJ|bCrZFW$niPykhWO z8QUMRAz(`ZnnUz|qqZ|wTy6Htb4al`i%Ak?5a5vUAk-6W=D?*#wYYVDyJX;As^}Y2 zrAX_Kpk1-w3a`-92JfnYceSFQpQS})A(>Luw#F;;n}ugnZ?>dRI)%bQ6pfPzcTJ}7 zXYPPQnF*h!oT;Bk>!Fkpqk_2IUi-b~<+sKXh>n{BFv!Uf~_(B4Fmti+wuK$ z3Vnp1)ZlM_8u%zz@?FzZgU`Wd=&ncni1{Q>hdqx8NP0QLZH?^u-rz-fMxL=WOY$u| z(Zn8qVppD5K2Mo9st?mPG|8{Ffsb-o$4vCyR1U#g0>-x6q_^%?>!lS>QQywOQ@0O^ z!GrRXQ3onQp;HO_n?fZid(!)yr!E`oE?4Y+`FL-Co3xUUx97ww;x!Ae&R^FIylWM_ znE4!%$%M7XE&J^$xbyT~#eJT`q{JeP1d7{Hv2jk#oEm59!>uWb5>yGc|C#{m+J&K5tC z-X@8EBvj$zK?T2VepoT^QQj-3Ca5G+vu%d}btX9Umbf4dpsDCE=3c?~R}J5R!FOQr z;fMyK!)Z)+$*&?h#u3RUZ$c)>O|NVeHwjnJQ#D4)ZFFv1GW0=tGSA@P(JW3s7&s^& zfWiz9xdR@QqwQ!iIyraR(kf1} z{EdJHf`j>;6ZtNOi}F}NEgvGU?=Eo^YIbtt~u0G9YV^7GV*g4FQg*K-@QzX zZK-yn_cp!^HeN$x4)GpGyt7?1qjxNncQn?coOc6hVqjx7mIMS#7U1w~m@nkA~q zy!Ew#hw|=8pDND9Le@;_S5UO(S}!wg{a*XU)9tRG#qseg*dCL=@LhCMyDk{GC?7oC+3q^? z)h`==Y)K6Vup@4(ci-JuyVsEaV`DwaE9>t1YHy|1Zd`YMYW;bq)9*K0w>Y(bZqZrZ z*jRg?exBv0o(tC;>z6C*&uVYaWgGT-vtF+u@50pj1?QgSx>d_{P4q0UcNN#=-bSy~ zx3158ZmdVy?RB0}(Sh*xi0f|V-=#hu%HK`%%G`fDa%i>gL6Xj;=l|{Jrv7?Ye(tZj z_H}2aS=}hFS^J_+-*Z#z&N-cbRXhKz+wv`c*f>_}%zpp?0RR8hnB8vMAQXVzG)>!W zt+L+iEh=|SlV0X>Z-J`J_ET`80W_M5)fxw4&1gt+Iu>G6E1ZATrNHn z4mvL9(%vF~YYSvuPkY=niXY=e2zS}zE~D4gGN?V5TP{|hIKMCB7{WcpVVs+nUSe22 z^Kml|*PZu?d9&|hxSOvB?!w{L4rhP=I?%7KXTGb;njwZO0_AzFNAbC&DH6JX^1mqC z^U(dZ9%%QA!`b>zUM84-@uu&iJcMaPu+o&}p`JU(bJuzpm%7-BLOi#Q=eG6y;f!>~ zXJ}ro9nW>^DfFnkH`MOT;cUI4%r*<>@$X>Xw;pZ#6;=w&oAyI_-?-cM(^`V(&9>k2 z>jB=keiF`k%5Cw(v zW;|OS(*GtuleBSvj2Xu!PoKhL;!*lE@wgu7_~JYtHNTvWU_$?}6;QgqP-AAObor!rt_(pffH*aa=hb UTgRm+a0L5*)*I3IFF`hIn_L;W;Q#;t delta 2266 zcmV<02qpK(5~vc8Hxx!?VjwaxATcpIH8wgjG9WTAGBA-*BavVRk#rD$pa=i}c$~yp zUvtzp5O*N7eduI>77Ei&5p<@#1H$=oBqZrWfu#KDK$;=Y^bHhSXKQLp9!oxwKJ~F5 zsUM>ssk{Dj=TnZ@(E)YhBcb`!&q!^S?Wj&o9p4@6U3I+5pC&uLF!WH~ zDdU6wjAZ0dA;)+&@%fDjO&C?0i-2_E?p!!RrgfaK0cT*~3@V&gFQ1dnTdqUM0unI| zd3Boz7EsMcjBr8rNXTNYl8~jW+d7_?4Zh1`|06a6Y$?EULRHkAZ(VKj|2b^xIEzUV zWw5{@<3Xq=*35x_OAi}yo8orCz`anjH>OIF&L2Uy=D!tQX{Qa|6$9@|%|1U%%aw&> zN>$eyue5I#UQ?dmGw@K}l0xZ}LM|&BClT(NOySQw00%M?K214OKatKuDI-P&aeMvl z2hGcG4Ly{1viN|dSaGH*AC|&Qu>#j-pKnvm2t)u~F3Z_}0SmYgwFd>GOtL5sd22te z8GP4H@FRhTK0V5D-!$OfGH_AeV+nLV@FmTJB1342jCvtaP*O2Trq-^qzp4LzXYio> za2%Xi3TQ=Th%keWO19Sc*A4vZ@5cAjDf9wAX~5q!@KLViyQZmzm_y7kT#xu6^GTi# zdma;z^m0aj*joAZy}^s}j67p$mgHM_Vu?Ne#IHQBe4a9|s1MUVw8^iofsb-k$4vCy zRE{880>-dgq`&Ug>!lS>S>H~>)3gtX!GrRXq63vnp;HO_n?fb2deVoSr!E?L7i)gM ze7v_!I!Va8bK;fjH4CpPURMpgt2Mls`5cnTgtf0F|{DPQR*xw4DMo}^A9mjuv(||v#P89Kh%qb&-6l!XwJU!1%gXiWu zc;4a4cs-hUZW}x(zk*hd)<+SGAg$sg%ijoSAUKeR$C+&3$4v5nYxmXUI5JnAYw~y2 z-1M{WA3^aQ4&n5MTL+gf#=t@lH23mkk5_u?H1nBj@Swasrl|K0q>E@8uw!(#`S=Kb zZPy&?sSY9KLmBxwix*N6Lm9R@fp+weS}B;UOceJ+s&Zr zGCRCkbH|p1rXs*!&4RKlD$NqrX5RYRz(aX=ygtVu$9ZA%UiLwhgf!oAMf)!Y>xWEJ zSSdXC9ArpkTmN#PX87Fxf0s`89(uih!JazK`BK(Q>DN%S=F*p$wti3j;^}tJ&*J## zHEfT`Uw1McxAp6sfs69~)1B>}Ghh9(iN}`IZ~!~xwssHPjeB=n@_%aRqrAHAuCMi1 z8~sMTvs3zKo#nxxmEPi%-lDUzv2pK1#(7qp`Yv2^=&#iDPwQ{rWgGT(vwpvSCGWzN z{(^JYQg6*tuZ^A6^}g11rN7ay>}~2ZUmE%-yZz-~s91*O?GxAC%)cvpK32Tj=9Rhs zcI3!P??IBzW#|7L=eF^BUvVC+x%PEuwO!q)u37t{$=)+ldS{&FHM{t%+w!e^+;z0A zhxmVw{{R30|NqsP-EP|;6o5&8oBnNUl=WtBQMqfH^fH%w3q);fBiR@M-j+T6z=;zZ z(U}})m14G&u6FH-c`tT9?K3y6$7h6Q5N;L5O?un} zy`W*PeGc6Af%NxuCfqZSi_V0D=H*=aTLf_J16kM89ydkNW3&k2E_>X6C3;;ggZgv1 zxZAG>ZsKrjhqJ%`Ixw!TXTEF8 znqy8%3gxR>kK!}Q5-fBM<$p20=b`&+J<#tLhqLXURS^?jHDe#;pnE3H5w-JhrV&ZKyaYnBD}(Gj}}m-ng)Tot;gHXXbcjeNQ4c z+wXo1c)mFv+upAeOl2)JtIV<`nu6N9ohV<9eX^C0RCV zena{9&cB3_@}uv6b4A57TpV`~^PPXW#)_1p=-!`(JpTDqs%si?9-1q@f}e+yX^Kr% zOb_D$-P0RSrl9qG%cm^Eii*Iz{@X)}Vzj5IP7{yy9O60l?{7ii_dUbnSk$A=fuV8x z))h^X4Q>wGah`;wsqxS>4fpu=ZK|=(B~|$m= 0) && (x <=1); +validScalarPos = @(x) isnumeric(x) && isscalar(x); +validMaxZoomLevel = @(x) isnumeric(x) && isscalar(x) && (x >= 0); +addParameter(p, 'BaseUrl','http://a.tile.openstreetmap.org', @isstring); +addParameter(p, 'Alpha', 1, validScalar0to1); +addParameter(p, 'Scale', 1, validScalarPos); +addParameter(p, 'MaxZoomLevel', 16, validMaxZoomLevel); +parse(p,varargin{:}); + +ax = gca(); +curAxis = axis(ax); +verbose = false; +baseurl = p.Results.BaseUrl; +alpha = p.Results.Alpha; +scale = p.Results.Scale; +maxZoomLevel = p.Results.MaxZoomLevel; + +%% Convertion from lat lon to tile x and y, and back. +% See: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames +% Example tile x=272, y=154: https://a.tile.openstreetmap.org/9/272/154.png + +lon2x = @(lon, zoomlevel) floor((lon + 180) / 360 * 2 .^ zoomlevel); +lat2y = @(lat, zoomlevel) floor((1 - log(tan(deg2rad(lat)) + (1 ./ cos(deg2rad(lat)))) / pi) / 2 .* 2 .^ zoomlevel); +x2lon = @(x, zoomlevel) x ./ 2.^zoomlevel * 360 - 180; +y2lat = @(y, zoomlevel) atan(sinh(pi * (1 - 2 * y ./ (2.^zoomlevel)))) * 180 / pi; + + +%% +oldHold = ishold(); +hold on; + +%% Adjust aspect ratio. +adjust_axis(ax, curAxis); +ax.PlotBoxAspectRatioMode = 'manual'; + +%% Compute zoom level. +[width, height] = ax_width_pixels(ax); +width = width * scale; +height = height * scale; +zoomlevel = get_zoomlevel(curAxis, width, height, maxZoomLevel); + +%% Memoize downloaded tiles, to save bandwidth. +memoizedImread = memoize(@imread); +memoizedImread.CacheSize = 200; + +%% Get tiles and display them. +minmaxX = lon2x(curAxis(1:2), zoomlevel); +minmaxY = lat2y(curAxis(3:4), zoomlevel); + +hgrp = hggroup; + +for x = min(minmaxX):max(minmaxX) + for y = min(minmaxY):max(minmaxY) + + url = sprintf("%s/%d/%d/%d.png", baseurl, zoomlevel, x, y); + if verbose + disp(url) + end + + [indices, cmap, imAlpha] = memoizedImread(url); + + % Some files, t.ex. from openseamap with transparency, come without + % colormap. They have three dimensions in indices already. + if size(indices, 3) > 1 + imagedata = indices; + else + imagedata = ind2rgb(indices, cmap); + end + + if numel(imAlpha) == 0 + imAlpha = 1; + end + + im = image(ax, ... + x2lon([x, x+1], zoomlevel), ... + y2lat([y, y+1], zoomlevel), ... + imagedata, ... + 'AlphaData', alpha*imAlpha... + ); + set(im,'tag','osm_map_tile') + set(im,'Parent',hgrp) + uistack(im, 'bottom') % move map to bottom (so it doesn't hide previously drawn annotations) + end +end +set(hgrp,'tag','osm_map') + + +%% + + function [width, height] = ax_width_pixels(axHandle) + orig_units = get(axHandle,'Units'); + set(axHandle,'Units','Pixels') + ax_position = get(axHandle,'position'); + set(axHandle,'Units',orig_units) + width = ax_position(3); + height = ax_position(4); + end + + function adjust_axis(axHandle, curAxis) + % adjust current axis limit to avoid strectched maps + [xExtent,yExtent] = latLonToMeters(curAxis(3:4), curAxis(1:2) ); + xExtent = diff(xExtent); % just the size of the span + yExtent = diff(yExtent); + % get axes aspect ratio + drawnow + orig_units = get(axHandle,'Units'); + set(axHandle,'Units','Pixels') + ax_position = get(axHandle,'position'); + set(axHandle,'Units',orig_units) + aspect_ratio = ax_position(4) / ax_position(3); + + if xExtent*aspect_ratio > yExtent + centerX = mean(curAxis(1:2)); + centerY = mean(curAxis(3:4)); + spanX = (curAxis(2)-curAxis(1))/2; + spanY = (curAxis(4)-curAxis(3))/2; + + % enlarge the Y extent + spanY = spanY*xExtent*aspect_ratio/yExtent; % new span + if spanY > 85 + spanX = spanX * 85 / spanY; + spanY = spanY * 85 / spanY; + end + curAxis(1) = centerX-spanX; + curAxis(2) = centerX+spanX; + curAxis(3) = centerY-spanY; + curAxis(4) = centerY+spanY; + elseif yExtent > xExtent*aspect_ratio + + centerX = mean(curAxis(1:2)); + centerY = mean(curAxis(3:4)); + spanX = (curAxis(2)-curAxis(1))/2; + spanY = (curAxis(4)-curAxis(3))/2; + % enlarge the X extent + spanX = spanX*yExtent/(xExtent*aspect_ratio); % new span + if spanX > 180 + spanY = spanY * 180 / spanX; + spanX = spanX * 180 / spanX; + end + + curAxis(1) = centerX-spanX; + curAxis(2) = centerX+spanX; + curAxis(3) = centerY-spanY; + curAxis(4) = centerY+spanY; + end + % Enforce Latitude constraints of EPSG:900913 + if curAxis(3) < -85 + curAxis(3:4) = curAxis(3:4) + (-85 - curAxis(3)); + end + if curAxis(4) > 85 + curAxis(3:4) = curAxis(3:4) + (85 - curAxis(4)); + end + axis(axHandle, curAxis); % update axis as quickly as possible, before downloading new image + drawnow + end + + function zoomlevel = get_zoomlevel(curAxis, width, height, maxZoomLevel) + [xExtent,yExtent] = latLonToMeters(curAxis(3:4), curAxis(1:2) ); + minResX = diff(xExtent) / width; + minResY = diff(yExtent) / height; + minRes = max([minResX minResY]); + tileSize = 256; + initialResolution = 2 * pi * 6378137 / tileSize; % 156543.03392804062 for tileSize 256 pixels + zoomlevel = floor(log2(initialResolution/minRes)); + + % Enforce valid zoom levels: 1 <= zoom <= maxZoomLevel (Default: 16) + zoomlevel = min(max(zoomlevel, 1), maxZoomLevel); + end + + function [x,y] = latLonToMeters(lat, lon ) + % Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:900913" + originShift = 2 * pi * 6378137 / 2.0; % 20037508.342789244 + x = lon * originShift / 180; + y = log(tan((90 + lat) * pi / 360 )) / (pi / 180); + y = y * originShift / 180; + end + +end diff --git a/CODE/dwind/analyze_wind.m b/CODE/dwind/analyze_wind.m index 6b41f42..544c40a 100644 --- a/CODE/dwind/analyze_wind.m +++ b/CODE/dwind/analyze_wind.m @@ -53,6 +53,11 @@ load([FilePath, filesep, FileName]); +% Make sure time is as datevec +if size(stor_time, 2) == 1 + stor_time = datevec(stor_time); +end + scr = get(0,'ScreenSize'); w = 650; h = 600; @@ -719,6 +724,11 @@ function CHG_DIR(hObject, ~) % Typ 1 -> Years % Typ 2 -> Months +% Make sure time is as datevec +if size(stor_time, 2) == 1 + stor_time = datevec(stor_time) +end + sel_val = get(w2.subtime_table, 'Value'); sel_str = get(w2.subtime_table, 'String'); subset = str2double(sel_str(sel_val)); diff --git a/CODE/dwind/dwind.m b/CODE/dwind/dwind.m index e4ec421..a2b7b2f 100644 --- a/CODE/dwind/dwind.m +++ b/CODE/dwind/dwind.m @@ -566,7 +566,7 @@ function download(wind) % Offline mode elseif strcmp(wind.db, 'InterimOff') || strcmp(wind.db, 'ERA5Off') - + copyfile(fullfile(wind.ncDir, '*.nc'), fullfile(wind.folder, 'nc')) %% NOAA else % Work on input coordinates diff --git a/CODE/dwind/process_wind.m b/CODE/dwind/process_wind.m index c0d04e5..61c93cc 100644 --- a/CODE/dwind/process_wind.m +++ b/CODE/dwind/process_wind.m @@ -156,7 +156,7 @@ function process_wind(varargin) stor(i).HGT = ncread(fullfile(in_path, fl(i).name), 'z')/9.80665; stor(i).UWND = ncread(fullfile(in_path, fl(i).name), 'u'); stor(i).VWND = ncread(fullfile(in_path, fl(i).name), 'v'); - stor(i).time = double(ncread(fullfile(in_path, fl(i).name), 'time'))/24 + datenum(1900,1,1,0,0,0); + stor(i).time = double(ncread(fullfile(in_path, fl(i).name), 'time'))/24 + datenum(1900,1,1,0,0,0); end % Do some tests to see if the timestamp is correct diff --git a/CODE/get_prefs.m b/CODE/get_prefs.m index cd2c7e6..401354d 100644 --- a/CODE/get_prefs.m +++ b/CODE/get_prefs.m @@ -355,7 +355,7 @@ 'ForegroundColor', [1 1 1],... 'BackgroundColor', [.25 .25 .25],... 'Tooltip', 'Choose the basemap',... - 'String', {'None', 'Google Map'},... + 'String', {'None', 'Google Map', 'OSM'},... 'Value', 2); p.GE_export_txt = uicontrol(... @@ -397,8 +397,28 @@ 'BackgroundColor', [.25 .25 .25],... 'Value', 1,... 'Tooltip', 'Attempt to open Google Earth to display the results'); + + p.OSM_scale_txt = uicontrol(... + 'parent', p.pan_map,... + 'style', 'text',... + 'units', 'normalized',... + 'position', [.55 .66 .3 .05],... + 'HorizontalAlignment', 'left',... + 'BackgroundColor', [.25 .25 .25],... + 'ForegroundColor', [1 1 1],... + 'String', 'OSM scale:'); - + p.OSM_scale = uicontrol(... + 'parent', p.pan_map,... + 'style', 'popup',... + 'unit', 'normalized',... + 'position', [.75 .63 .22 .08],... + 'ForegroundColor', [1 1 1],... + 'BackgroundColor', [.25 .25 .25],... + 'Tooltip', 'OSM scale',... + 'String', {'1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12'},... + 'Value', 4); + p.prb_title = uicontrol(... 'parent', p.pan_map,... 'style', 'text',... @@ -623,6 +643,7 @@ function but_ok(~, ~, p) prefs.maps.GE_show = p.GE_show.Value; prefs.maps.basemap = p.basemap.Value; +prefs.maps.OSM_scale = p.OSM_scale.Value; prefs.files.nbDigits = str2double(get(p.prob_digits, 'String')); diff --git a/CODE/plotMap.m b/CODE/plotMap.m index 4393867..efc5274 100644 --- a/CODE/plotMap.m +++ b/CODE/plotMap.m @@ -198,6 +198,9 @@ function plotMap(mapType) % Plot basemap if prefs.maps.basemap == 2 plot_google_map('maptype', 'terrain', 'MapScale', 1); + elseif prefs.maps.basemap == 3 + a = plot_openstreetmap('scale', prefs.maps.OSM_scale); + uistack(a, 'bottom') end % Plot vent diff --git a/CODE/runProb.m b/CODE/runProb.m index 9e091b2..dbf5854 100644 --- a/CODE/runProb.m +++ b/CODE/runProb.m @@ -149,17 +149,10 @@ mkdir(fullfile(out_pth, 'CONF', seas_str{seas})); mkdir(fullfile(out_pth, 'GS', seas_str{seas})); mkdir(fullfile(out_pth, 'LOG', seas_str{seas})); -% % Should probably move these out of the loop -% mkdir(fullfile(out_pth, 'FIG')); -% mkdir(fullfile(out_pth, 'KML')); -% % -------- + mkdir(fullfile(out_pth, 'FIG/ESP', seas_str{seas})); mkdir(fullfile(out_pth, 'FIG/MAPS', seas_str{seas})); mkdir(fullfile(out_pth, 'SUM', seas_str{seas})); - %mkdir(fullfile(out_pth, 'CURVES', seas_str{seas})); -% if ~exist(fullfile(out_pth, 'PROB'), 'dir') -% mkdir(fullfile(out_pth, 'PROB')); -% end % Get the wind vector for the season considered if strcmp(seas_str{seas}, 'all') @@ -261,7 +254,7 @@ % MER and mass calculation for j = 1:nb_sim wind_prof = load(fullfile(data.wind_pth, [num2str(wind_vec(j), '%05i'), '.gen'])); % Loads wind profile - + % Get the closest value to tropopause (assumed 15 km asl over Iceland) level = find(abs(wind_prof(:,1)-data.trop_height) == min(abs(wind_prof(:,1)-data.trop_height))); if length(level) > 1 @@ -280,7 +273,7 @@ speed_tmp(j)= wind_prof(level,2); % Maximum wind speed at the tropopause mer_tmp(j) = get_mer((ht_tmp(j) - data.vent_ht), speed_tmp(j)); % Calculate MER from the method of Degruyter and Bonadonna (2012) - dir_tmp(j) = median(wind_prof(level3:level2,3)); % Median wind direction below the plume + dir_tmp(j) = median(wind_prof(level3:level2,3)); % Median wind direction below the plume % Here, either sample the mass or calculate it from MER % and duration @@ -558,6 +551,17 @@ function load_data(~, ~, ~) prepare_data(1); +function choose_wind(~,~,f) +pth = uigetdir('WIND', 'Select the ascii/ wind folder containing the wind files'); +id = strfind(pth, 'WIND/'); +f.tab.Data{4,2} = pth(id:end); + +function choose_grid(~,~,f) +[fl, pth] = uigetfile('GRID/*.utm', 'Select the .utm grid file'); +id = strfind(pth, 'GRID/'); +f.tab.Data{3,2} = [pth(id:end),fl]; + + % Get run number function run_nb = get_run_nb(run_path) l = dir(run_path); @@ -593,7 +597,9 @@ function data_gui(tab_data) % Menu t.menu = uimenu(t.fig, 'Label', 'File'); t.m11 = uimenu(t.menu, 'Label', 'Load', 'Accelerator', 'O'); - + t.m12 = uimenu(t.menu, 'Label', 'Choose wind file'); + t.m13 = uimenu(t.menu, 'Label', 'Choose grid file'); + % Main panel t.main = uipanel(... 'parent', t.fig,... @@ -636,6 +642,8 @@ function data_gui(tab_data) % Callback for ok button set(t.ok, 'callback', {@test_param, t}) set(t.m11, 'callback', {@load_data, t}) +set(t.m12, 'callback', {@choose_wind, t}) +set(t.m13, 'callback', {@choose_grid, t}) uiwait(t.fig); % Test input parameters diff --git a/MODEL/forward_src/new_tephra.o b/MODEL/forward_src/new_tephra.o index b86068be0023fa0925c315354b0519d22e3a4c9a..7a6be60e11699c379f00c2e28cbef1193c23e539 100644 GIT binary patch literal 28264 zcmch93w%_?z5m%gdp5gyZ3wUMh%Qe-Bne6skpPhdHV`0%072AcNp_Phd70htlJLkw ziE(SFaF&|Aw$40A6o-afPZI%hj~^M@fqEsP<%qru|tOwk@8|>#Yrj{Y}2w zX)-ra-YF#~OI6x&k*2YY=9SnaueY+HuCAfpEcbeg3d;*kO4G_zyW;CvW5vREaU|MT z845QA>#J2}qJ56hvfvMN60SY4FaC;`&3vzSO(+z{kZ50pvTuJJPkejn!uZSUt@5vF zuJ-!YtZDMEqx_;3g-hl}S>B&0lPgul*4${3ndoJrDKA(bXdpwPyg8F(wYQxTU(b{1 zsnImmNBzo!69r?Uys`PR{`D&E_D`2|+;@UJUt?pf-&^mi^V5%s@>*5>u&OnDcAO-c z>{njHny*%}I8ok5Gi0T)611;ftr9WlSDvu1(N~!$Z+Jm(c|Z7^^1>S%`)%(`CGR0+ z#qobC590%Ng}u%7>x1=GUT%2Noq6hXV6}{}c6vd27S(b_f zKDfxl{4&%dPh&0SM=-%2tqdPtcJ*PXbN(BH&lkOxc{<*(EdL|;RxQHjgp+z~^gdR>a;KsMsRs6|y&kOl&f6nM^`NY#D z#*n9L+DR5CkX%w-HYpx5?^P>LtKkoq@t)O{2i_aV>{wbL{TaLo?gKaHG zZLY}Zc&=%lNdEXynl|O2kozpn6R9a5c_J;x_VJG6iKJ}+OKHRi1NO>tPh{D#s4hym zKEABw7;ntLsV%yqq*9Xp)S@obPCHeAb-eb*jSQi<2G4wb_VR-kzg;%mhF%>A9hDA ziPzbyU(at<7OdI^2|_*vs(>&(w=BV`jz+|sGAzr}*|7m+Z`%$?j$JH~7MA_*EJXYm zu8z0j4VWf1m$_l}xR_MsVdg%5VSjBoVsNZjC`!abpH{7owaSfF)%0zZ4l$&Ka9?lA{zo~LUts$r;znOqbe zDHgZ?uHyYyiLQ4=I7fTB#J$o}S-js9JlKfY$z%CRRBT46k|^bPA`9=5eObJ}?VzQz zcz@4h!Y^hlVFD_?3ttQOz|K*$)7smxG+G4)V}ihI4Y$CH9j9H_A3E(CQ2el`>&mVS z%6S+fhEs&)Jm&92RnyYSg)sjGY5q~vK5u~s6Q#5^u}WFS`ARA#z>kky08?zTHIah7 zq~w@@LzOQW0j?DXd$*W9)*j=YA_s%13Uf`@!ZV&LkGeY2(VVA51a3dP2kGmc7lK?V zfo?UQmZ3DRBIZ+1=Ydma&z z4vLB$l!f3v3J|)Z)87Hfu!|)!xNT6V9B`!L6J*o!o<|VIo=EXA2ziUFkpr5Dk6l#! z%u%pN=ZNX12F$K!I>^y>W{j&H5@U3-t34y0cDp**Z9U(nyvP%xrBUhBB?WhO9HZW@ z0|FxXI}t$pI21%N)poRlnsTu7h$*FOSG4++7dn3q)8j;lYOVf_ zBl$-dfN~u5yyuM9{R~R>4gB(S?}+IZ$w+=X?Uy*j2_BqRCRak}4e z#J4-^OV(Wx(=C#b{2yn>d|uT{_d%2-^fx!V=jrc$0j z_Pmc}T`b6|sTA&ybM}{|bV!`i#$HNOsg$?2SPy4nqWA(W6?1i40Y*J1uB(KXMV_Zi zG~K+>bq+wpn>W4uJ(&8iOneMrdSe%=N;#I#FQ-Dv>5P-(iOR9MIxhd5_GCyoRdI4q zQ|l35(|FU3&R9+0NpM*W*v{j!^UPSB1>Nq!uXv~P2fyD?I7fn(#!*}4=m?2IY zLRZK7?tf=UBVnPc01K0nrAsZm$GYnxd44DvT z$k(7Y=8<7EggdL4A*(=rv95`mr#)0C=V_6Q$}$v#H>m!vx@rduQ< z`D^2J--_LItlj=ES+_f;TO=d-1LAb2^wQ0eFSFefV!A~#l79euO4H}JU$PuDKyC+-NY1P68mLxa|_! zdL9%q4^XCPqNh??5j?{!q0Ex3p6D&Hh~}7gA{_a2ILz9yZ4a|NSZzh}U%=5|biZ>6 zvU&Afj4^wpV9>*Jr*z7UN|5{cO5w@qP95t&xe7|i zKjdjUNl}BaEWhHXmDU ztaD~hIMR`YeBmtVjV_UKmBkgA^%>PYl@ZpHQreG};kF|yiZ-yQ=TUrh36G0Bs*h_X zJG;jbKXj&b?Gh$N@@N;bhxSg|LXYhrB=UAPbWre zrwIL6B*t7KhaM_I504@u7}Uzm(UnMs;X(On*R{rMfzi_#xf8BDAR2Y(J^gbIB%Tr^ z^f|NSK)N7iz!BjOxa3^V01VTZ`D!2u&$}Xz3**Ir?$Xc^p1VA2&txfEqC26=2>Zj+ zIp}CCUPMmYjHB2&ERPA{6xZ4D6fBZ@Iy(;cHCPnyVu3vGxE=R}zTtTdtP;aCF)XID zQBWDzC>{mOMh;@O<5kEt=>&7TY&dKffO_W(>f7$Y!ca!7IM$&&9^4Ssp-Ub;(PxV7 zJcN^LoHVsRj9Dvg29gIWKQ=8kXPij>S*&$EkzHaEj6^tWOCw@%D@T06Y2yK5m7JJG z7d#A!agN(7NX|i>9gicO5JG}`mtcypSWX)+V0YQKU!KEY7Vd|W#k`L58sT(T>pF~u zKDix!bd2j3ujrl((Z;^rq9zQ_+844s7IxBXo>O+vFQ13m(eq-Qc#JmHYZ(z44cA6a zLegp1N}PjcMCf$zs63vHR>E=4f&xFU5)ne5Ly3cJr_j&QcCCkp z%&Z`rH9!Tf^mHtMr8CVy?z!~S&m%npgFzjl-zs{Mi@Gmj9C^P`t2x>}N2~F>jCqs9 zje*12q}5#2Pv^xQT+8{aKc2WV?DI3>I6IQ>#Hr$5Ic=kr#(oUGe{RD{#x;T)YsTn; ziat6fBB?@6UUCvkykjo@N34}NkJ-f_ZF|nrC2qAa=w@ehxSEED+tQv`zPu}-52N?X zZ^ro>6G!jmsVLljnD@rS`=vYne#sy8`nyqb1wv&<#HeS=oB+>1fpS?9H-WYfiPe&~ zrM5FxR~u)#p6@}ti~=qp#N^j=cYL0w@{H%ud+LXU#5HRC`UNfT6q8NDyeQ@$u}1B{ zO|;Jsr?I1X{oHGHsdbS~E zWLspfG4Mqz%!RfHY0sfQ*gTeBqLzgW;tqKphC=+{qy~c75He5UeHD*boG> zN%83&Gff6}L%@x?!FqQnxXB-ynVhVR3u*NYRerCwxY1wlht>uuYx+1yX=rkn;!0`P~KXan+&Zg}x2vFVxBkmy~;# zmKT>TEf`njUK6bMt`Am)YYN6yqnC2>0`9o4YPXiG%?;KELpA;?xVo_xQF3Vx+NrIN zwy`m(Uwb^{9(X#?xJ9oo3Yd6_f@&! zsLSGvQ*~F7Zyf3elC_0}72bJE3QHC(ExD$cF04Xc$)Y}aDmMDsSGUnEBdfPzpVzp` zhK8o9V7)Kw4=Dk~6(!}Xs{LWF6nowDC2R>CC29Dtcu{%r62ZbKy$(6&wN8X_WdnMy z-jx5@yiE=3Lr~VoQ^HwJ(-a}L%C|A(ULOqC&^VD(>tE-u4Y?7iZWPi|vn{5f$xgRv ztQg#KVrg#j$$&raFexgTJ9pXAlEsU>a~2oRom(=e1WwTU%`ISfL3L5lN-!^4S{|1~ z*`+9(TUc7^EiYME>@D{!DPHPXjAsw~$y1U_OUgYZbG${xh>`M@T5nSElr39WGN(Ao z)0;GTqGK;^7^;w)v$#|aOD&pGC1s_{VCTZcMa89BTpF`U)peX~memI_h1MbN<7N)g zE8&_Zf2gLRwo1c%JGmZAzS`bw&0-c?;}83^tMS#lP{=4OE%Ow*d8fF#xV$hnl<*9w zG-gR*dCB7VG|d)H7F9JVfKyYlMiJtRJfNinPjSgSPkApw$Pk^lWLX*HFY+!dT)I>% zj;EB2UgNj#Xs&ioXsiiNOc>LBM|g~{xTvHsZun~p6EcKeQ!1vk*i51IqI&TyN3?p2 z3zwFcESjgq6VaB6%tecrAk<72lj_Y=Tw1)GK3}Bdm}Hd8GQct}eW0PIw{A}V+t&1q zp(D}bky zGz6fhvjQHw9ca()$ZCCg~zKrkwmgL2TcO#AFCY0yQ(*cMG2Y-m=MMQuY8xyUng zYA>FwUOa{=X;@T}K}o}-l?18%9*ZrkT5&nBcGR>~Ef3#Z4co^Oz%euvBhe;89tXnnbP#|ksw>g}{&GhzEXExH`AWv(C^I|jImg&qxdJps4 zvz+sQqHf0oCx5>=#QYt3PChsN9MhfCoL=Csn7(DU^Q%BUuWaq|IClf5!&vLBWzGkY z7BX+EaDE%<>zUtO;4CB<*aMW~0^^ zZRd+2!snXp10;1nX^ZVsa#69UWfv(5(UP)}*@#l>SCS?mxgJT|yrj#Jya`D?Pw!?` zhiw*e;`Ee*RC-n-nL=f@N~BWvB9+FbrSruBU2{E;v?D2&WzlT6q9o}~e4R#_k&)48 zH*_cafN^6Onv)Nei?rS0)m&QVC8#m8{EotK*V$rGGcGY|1-q~r}) zY^GE3ONV|5D4Q*f!h~^Enw6b89?bS}aUnd3waH~a-iXYpt1iMJzLO zJu*8X(K5rSoA0C~-GPjaD8RV9%NXTsRTYgfn)7Q&i}sFnIPWpj<1(DzHq+y?oKKnQ zi`>qGW_rQ|=TS4Am*@PYnVvGu`MR0Dbhh(-GkuxIN&g95`DM;zm?YAdS2%~6=_^({ zCz$COHO{GKy0FnX-ArG-!8yxJ&)Mc&W~PgFID=-oxLeCy58~^Xek=t}{nS$U4wvmt zWY~DPpY$z!Z3K(4cwm;Xz`4b&TBe1=MPL#*J#eC&GagVW_=z+J49LL$NJIwGc$8R zT;^}g%+k2bH_Xh1vCPc(kSPM@YGj(G{>{vs^77> zy!GaI@~*8rNmjPPlE9MYbSn`u`Hj;lKECD8K&<%)pAhMbjfgi z=gF3zk}kEjqkK379VfZa`a?)U@mnpwl{8HTy=}{zq!(F_%JSc}oFKi}`m&^twY*Qd z%qrfs%j(Xg9dJ21#;#Y7@%S)o_QWm|Gp!y|CyHu@I{VYLMlvb*uboUl3dpb%h(|>N7 zzGDl=q2(7A_MPxrD*2=InHG-C7=5Y*iPl!mClbd&}h%J=9dSV@qK#S~We(YKKjut<-Xcouk>9mlQ#Can3WAV+va8 zycbF48|9kwJ7#*BLrdcgRotY8#icDG=YqJLqj70eFXS72?#fypj-uAd(-!xX%X0E^ z^4t?d?&K+#<$=n%^vaye@+NBa{`GhZmFF`#m0B)q2P?%!Hk4aYiSw=;yds!8r=dQ~ zv$fnJZjQ0(ZOSdH^@X`-3*|2KRW4pyk<+L$;&_%8Er=F}s@7_XlgZp=AsCk%@~_D) zYN(VacTiL159O8%Ys#7$*5ryEV=mOJ_cc{{Lrs;MIOdSAq2f?R3k1VjRYNn5RJ6Gz zrNwyTFeI^2YpxFkt9ejYS>tQcy!czMc>}e+>X3#FWm6agcIujUO@Ou5Gy~T7Ldf8W zpXRMy6OaiUZmbO=qq4Svr96_+q{xupzg7!0HiWzo=B@HI;ef!~5C{ko)?GdosoDl~ zg!rPpUY-&{O7j}ci&uY}e5S<}$b*T^!VYiYk|l*Jy*$nFmai-;_IfpMy(p1ekZdI_!U&#|9tud@sdmHf% zuLx!ht(GwdVH{-|Jt0n~nW7Y)eaL7YM8z0!Pw7nZuZf4@54QqV0AuV=lCz3)x^Q4wWo`rfbBE4ug^l#8MYB+kYZQ_;MbXA$qXa>Ij39v}P$E%?a8rz`o} z4VMsHyYZJ(%TsCV-@Rw@)~0QX4^O$^)S%ouuKCJe#(a41wF^eRocu$p))Cl|a(w-! zv-jPzBjt$)9{uU7KfeF<7Z=yPIdS7k?P$}El(*m9?E31~f8YGb{tRz9$QIva_L3t&{yzJePaFR`oyeyV8hlUP2Vj={zA1k!v{m={;o zq-Q98vBm?v7soxk#NrefR`tbV7xWRuFP1}~pHTL2Lq~Z76u($AfZprC7cMUCS$>hK z&n+8iuj1#9g7h)mO?mN*gY+U2xVX6>{f;UZOJvZal)YT$v;0OSpL+q)0Y!7gPr6&t zV&MRK0oDb)xJhC;SB$*Om}4HV#F)b`hv&Ws4^FN`{F!?QsDnN`;?}@MgQZeG+=3kJ z-0?AEnYyT-CuKNY(kk7%HfN-+d;F12Md4MUTYpz}=%Y_Be|t*W6K{?GTmHGm|33cq zmQ8!U_0q=&Hs}8%Z`+QP@Tk=EyWhUS^_x@IWGuew*J-z(J~m>7@wLtFRRteE_=EP% zBkug+y-QCu{@I4c@&t${JO%jh6sHm&KE>3D&pr5j4Q?qL-40@Z8a8%i>DjR;(>?C;<*R2cpSo!ii<}W%nHZ28dVmLQWO_= zmgG8zu2WguyD@7oT&c2n!on;byr?XmfHCV1eCStPJQiXW4^os}+yF9*dqS1Pb2w(r z!$)OtZOJU|rj&RddoYWeeZ@5vG_!b?rm}bx!>n6lc5xHRtOsM(bA!$-?%0)h9_ulS z%UC6gCwt85#7A-QutQ`)fXbQ*npu1BrJ4;Ydaafza`d`ZD~j7#Tz$&-h-|rq?XUb4 z!x|U&l(_oTkKgI_S029|>#sb1&(d%CDCFYe#sOC!{pQHDY+9mOHe^^bQj;7EQ~DNf zCuwJ%;W#{ZZZ*>4DooF~$vP)ZO+bRL2nBi13bu2~Pnf#<^3{URqxiVx7gs-gjLgQv zQl9c#zw-8j@8mJb$DOjc`jK}61@4z6AGfvV<4Z-mC%i8CxCcBRA45rdQ}QXl_OUM+ zw_uxA$~d%G{*j(WU^RR!|JtojZgZ$9H&|a;+g#=U7w@v-RpkaRn|AsCrM}R{kd~91 zlanj|agf-u{oi>u7wU7p>c6}6+SsC#v0$wVHssW>wfJWD!zb_G;j6!>~K@ z?D`#BE5N3IfW<|MbsT;$Llu|ZG0I4@J1UF}dx5kBb%wnHGK^|Ja@5%MPSi`b>mO^6 z8dg*pY<~^*wYKW*TLWOxKh-LXL92~)d#jNOIRm5|{eQLuAXvZI5;pA3MwWe_;iAM3 zUbXB0y=AwNX~$0mVQ8Guiri`=^C2S@5|Saor6jy7CEQI3_d&u5DZ#Mo>1D|Y=CkUzk#|EH9{&#>7YgAD5)yZ*wK3M0j???VCnr~hQjaX3UjAW36{-D`}3 zEIka5J%g5njR7+ZJzxw$6SF{#JZ@aP+dv}Q{s9DZ>rX)*-1Hj=SOsevDhz`@bnG>T z!Q%0)#yA-n_SMFvD4&^QSS#TESJ1?;kv_w)24Gta8Xh)Mjw87smR#9yvRnt>a6Nn+ zo!QNJNJoR-VPow7<#-G+q(`;O<%GmH@dV+^uq z7&!ss^048YVay8Xdv2^S#@O|@ZO7pjw5bA-(k$vb9W(UZP>s?*+8&igbsG~1Cx~5hCcQ9S8ziVTv6@2=8u;5BKFNvY{iZOJhsGJUYcS?Eh$K<&oZy4kaMZDVE zjGTU}eO*@jTddl6RJ#z>vSqa^`mGj`)jo_>8-;4aQ7ubWV-NQ=a)+$;@fWPthHB&N z`dPDDg>ivuTc%?^8kcKaLUlQbJ~_Kd)N3=8k1G(X9KZ-=Rc}(i^*%;D^v#Y~J+^Kn z>J6@dOKZ?KVg2u5!r;8sfT)ZD5z>wTqRO$rxL_qG9sL8zf#^FQ$6qAJ9g1TrtZ@YN zKW#z8XQVlsw824`7}KU2o*LdJuTF|$7zuSppX$7Cdg zUH|d(vaI%4tXe|YP3br6ej%&e< z)Hc0eBM-`Ihho)6qS|0o8^~(yM*kympR9HyR_#JmV~h{^tdaPHtoGAbwS?*7vVM)+ zE35r7RxP0e7{Y4YC$GeN9zhSd_Eye)cczOTj{pNbZa-Utlf0K2LMSsJ<83EGW zxHDUEV;hs~c+al4x9sMeyItIX!p2xN=T`Tgb1~Ik!g09L(Eng0Ukl0TRJ)#BfjfpH z2g>9`e`>Rs=yyWc0^HQD#w>pV^QW!HNWKrX1Gus2C&4oU;|n~ztvPY0?lrwjRY~@1 z?D~n#ymRThm2D2YW1wi)5Jx8R7ev#8>~k^cAIHdCYa}%onRQ}99vLu3?1qVslS0~fNKUhN+Vx|baXZ)dfO+oxpIs#<0I{==$mtcC<|1=~2m6VCeqe3Fkk%na)3iTe_w{yBW@Z3FS#x zaSeh;b{j)>L+&?uH^nXAvBH}RRRM!|ez0O#CxIO#OBOaDG3{rfs3EJEF4<=>Qd z!3le@x(ev`gJ+0cf0J6+H>vc?0$2*H(D$GeOS898dNryn+D9Q)$9%ZJc-Od4M<5yY z6A)5?rGfsQj>Xj;)NjSKtUFtg*J@-r+{jsI*Sp{&m;QH%qUIbpuG&a{%t&j)oW8jg zcS8Lv$~Z;_pWw|p9pytnp-V0ZfI<(S!_9kyz8$yQVb(F3m}!GX`VGd=uObyTMg@$E zDvW_Mj41*AA7IR~>o-lbUygD8x!qd9AsYw zn_7*2xdwylB65ud7nil8h3EywAl=EhYsLK12BTXcrcKaJyI#?%=Zb|Jd^6v0V5J8e z?cf>iI3X5g^z;d=5x6RAHL@Ki^kZ-W?uGveCm0170RbZe9NF`Yae%xEV;aC_Utt)I zRgeZE1%t~5Po{}=e=7(2fvu-S4|mC6eHyYBAofZFMwwW5rtr?nwcV>CdJo;m=zZ44 z=sk2JqW2Cf7T3v+ckYMZckAzA_P7YEYGiT{r=Kvan3!FkfH6$JAI0M< zxKo+LN#r`is&71D4Aj#&ym1pwMR>D~DX5r^Wi#gm-Tnv1lxq#=eOPsd^`k8n)U)q= zdJYRc`{VW47+=q!-T?g>%9|Rxfhv&eEnJNXI zP46~F?lT5*CF;1>vC6JLjsJ@`)UF@3?t^XUQP`b@eqwrvXv6hJQUSOF#z2thl<@%e zhVwOeB&>fICSu+GGE~6F*!mzkxt@Y+b&S4L?9xzoI5f{NQgif&CHrstvg;4F97h%X zEv)OU`ZWk*eG=oM6`|KA&}zTguzc69-`;|Ko&FbVD?0D?7IfZ6mTtqD{8KqGa8^jN zSLhvb*YpWC6s~pX2hPh+Ve7LH)6NWi`0ZBgzoCLfdNd9fW$@p6i z#Q)ZrJQFwpI01+UJttoP#NXth_#(d@NcrdBG|D>#q&;iHwQ>R4+qlkAK`G3^lL!s{}qt>j{>Rhkjmeu^8Z8S*8^$)EFiAz zX+ZiZ4@fx^fs|uYxm3R8KyIpk{m@<4=9i^8yg<%KnkUO+e}i z0a@=_AoVN*QqOcC^^5|tzDuD2WIz2Kg_Qq1kaG3|Dd!O&<$M!JIlF-PE4Ue7xazh^ z%-#s3oRvV@^%aFBK+3%uNV!)5Dc1n9f6wBJ{(c|Ge*X!O`QOGD`R@mk|LZ{V%LlIkQqL+N^^^jsXQnK#yFy}i9+3VS3#6QEAmbDdcMDoq_;U!M zzTH4HvEW7^Oq#qx(bIvPCl8|WwBtuW+W!L}{r)Y5_W^0gUBF!MZ3pH7eZU;x3Lxth z0V$_I(U$`$XOhA(K*|{jq#OfCIq#q`l=C`}a*hHi=a8cLPn#&`dkVh=q?|oK%IO5= z0fRuw@d7EQ6i7M6ik=CioXZqW0#eRcAmyY1^MHSGNjZN6QqD0T<-DNi=YW*+q{8n3 zDd%BeF8Fo=Ddz?t?WqP*PC1Zr7Akr!kaDh4co~p#E&=9(Z!nN@{+=%7{0T@ouL3FO zmx_K7NIB0cd=f}Gdx4bm4It%g1m*x6fRuABkaCtOdJ&Lv<|w=hNI6r1l#>mloKNwJ zY7X$ffRytVkaAvA^sj)Fb41~@K+1UvNICxjq@3+Q%2^MjoGKvYT&L*eK+0JFWc>n3 z*Il9VbAgnX1*E()sZ!owft2?;kn&zu^v{5l_Z*P&9+!08cUAsFK+3xfNO=uF+UEyS z-Vz|?El_k3kn*ksQr;v<*Ns#8*+9y(11axa{0L2XzXMX(2zT{#2l_3&{Gz6`cxXeFKOoq~H_$BFy?dK-PZ&$okI$S^rTW>p!UI z-9Y;Hc7q#QgwDda2wQqD9WL~|O&r~4w zTnwb#!76`0e$eH8=@uaMwgFkb7I+EhDj;r6lji}ce-e;OyzxNt@r5bkSu6&5Upoop zeeK7Jeg=ppPks_e`Fuu!{0{*~ga1BYI`aAN87=~1mmSfPXI9_3z zLM}#m|K&5%)YAi`o|o`N|2?npX&~G2IFRxB7?AuA0r5AOuMV)BPsU84ejw}D0l9v! z1menG1jJuf=*9)VC5yeZ@f5<2zH-cR7&yCIKm*&$^Hv z4kW(|h^60Tz5>U5?uCheLoiV9J3!WZ1IYZBfvopkAjeSzDC`4rUDyH~h5U7jt_Qk7 zmm)|Y+WvhwkN883YvNT3zmFgwJz62hDe0pM?*-z|zC&RMhxM6-R5aTL+1<(b&05}G?3ie`amOcxJEcyCCG|ieDW=xHWFyQyBJ^oWJ3K%6X@{?{OJ2$ z_2(y)PfsXE%#&Y77>ny8Nh=CtGNJzKvOLOh;~9U^@;^yBs*m#u z{-W2f6X^T|n)63I|7z5Xr<)UKK5-adep^EMrwR0P33O9J{reNj_av0_T>$kwHXiDh z6kkt3(tOU_8~0eY7tbTB2gJPr5mi`N9a{-)+eK7Q`S0Id^f=xP_Jo3pP2XJZHsux)y;U7L?{h60Gs^c0e9hEnS#MGF>emuOyGa=z#~H9 xl@u@El+e5@G;i!bkZAb#4Ca45kZ;k;ms7m*O$690cH!mU>Ou#a2*zqQxQK68>_ z#d7C^ncrS}?X}ikd+o>UnK@~_Yl@&Fy;F!y5~YcfcxXplVm{Ru3vomj;@}2O|5Ovg zC+fGBvvKBBg#_E5sEyGa1F0Z_$-0eCs_zWYh_< z(&j^RFyEC|rSKIYwYf2c#C(gFq!4{b?HR|1`d+m(g)bMWVekdqdUL@yzbwT-;(G*q z>&Njegtw9h@r{6Q^EkfLikv_NkAttx=A&}uqI`8zp#RlCe@$+nzo{@V;Lm+4C)dQa z5LeNbQ}WcdNrzA4IC-g^DeO$s&LqE`^c*IGwmc(+^`7#Pyg@h^)f^Ck{-<7-a2#cz zCRgU{ZQ9U4X54;i7>qfWCtk_9eAk=(H9ig2D)#~8_!eEyCuUx@lxjB1OAgUUG0c2 zvH-Ppy2xH=;M^2PfddZ7KPE4-AEW6gcqx#wd?3q0NB!hEadKT?z{9-x)aL^O8Px-x z5ukylmjf+74upT)lP*NH_H^~2ek!RnBeC(`f_^dCcWb zl=B3j4A-o3dENY!?brbjuYm3wpL#22iXTIQY(ZDdkCVd^Y7voXle-t}hS_GgEYyx~Sn9u0xsQyqc9Lq&T) zp$4Ijw~`z^m{sn`K~+u!dUu@=iKWu63CCa_!Ro1~Jwi*~z(6_XIsXYbDQmyfk?Yxu ztJCf@`CYV?Gx9rdwO#xa;*cVA3tZu%~q;F^UH>UV^ zF#o;@&06X|5znNuj>UTez0Y>s0JmUny!})e1+?v&K>r|ho}+z-aHivE-$4dI58LSE zoSbj^N0!fw&kV?(XwSjD?`YfglQ{#PpMbA-MC+G;SgHs5Nb2LqaB^%2U?4lJcoET0 zy^^N#`t5;$FW|gNI@#4AYWn06WiWr(RrTa-maIr%K-Nsk8C?!$u0zh9WmIr%Jvv9H zgWUV%Ff2Zu$hNijzLS>FPu7pM<+%aruEC9(%t0s=4(y|nf&O@rucH4)+^>(jEJQVKXm`815nj#v8oNG98)=k!_VCIKEX2li2O(HQyEuhMu- z=;I)qf>A>{aUvE{ia;vIKyW6Wcgr)AeUn1=yF7FrM%9`QW8ffAbC~tYaP?uF`#FsS z&sWGLQl_hJy-&)3O9akr3QV0l(l5i&<>rH^1B)wki930n9a*nyOJ&<6*8Wq$6x}Yt zlw`qkU4nb1H0@o2I;Za;tOzjnOG{C4PdOw}nj}&>fGdq}7@w%?$?-{c#N5=5_z3up z_T7my9VcfXU>`xJn=Y>z2t;(X9-;tiwen%x!`(4jC8|7Nep8TYPldjM(Q`ckC%o>55fw7s#st2h@oFNk( zm6{+x!?Lq9$^WcNdw(BU$F(zZ`toHD<~ew*CM*ix7@^KN`)J=iC`8A}#mHxEbiT*PeWNFJb+ryQ*lV)1$K|x$83eFllkziVhfP!Hj)8 zpp9Vk#Y{%J(T$`<`^W)2HXDJ#3;DJ*&~HqbkJ7lifr+bR-CYx2@-A>dO04RCiJQfj z2JvpKAf4t~XjX5Oh^!Aq#I{{(J^q+@dKcpkx;O1ifn%nn0h_ZfX>4#=<{ z7atm}_54`9Pr+Y4hMKOpvA{Et}@8oK3h7$*sZiaq}Rx z72TM;FTET4SU&6~t;Za#)X*UOpJ`?JyjU6EuP994)mye*Q^d&tnN7=R+cDeWD)wQp{ z@_!+a1Mjl~$hGPF2nE-E8?9S^oLskVan`L(Shvyuj9OI}&Bocjk;O_L>x2Ai*O2@E zeExpx?D3go6?JAga12s-%BdpR_cr#CQ@|kOqF0VT?ll0h=X|*GucuG*>@XN`axwRC zz2!X%M$i1M?wY1f*0~VlOx8nSlU9&sS{Kktz`pUm{Zt1M&axoq@^8q_)jvoJfTR8+ z8Iyg;@FBY9B20slb^+CQQ@u@QLeEtkoowb<@fb)wA{Fl1vg zb*!Q~-Am+usZNhT0_?bwu4ifObeA!s! z%C`tE^6f?!!@03`k8Leo&(hjjb-c-k~1RxVCS3Pv<{eI!33sjhKuYp@EOM82+z}c+RC`^gbw7o-Z1h{D|V8kK1>-OdGu)a+K}xI<@*ITnuFU z$EEmxWiHwa2~NYo-!p9EiF}sE54&s|eTfo|<1rhN@rPUt=f-$yI+yXSP(Yq1@>v@H zRhMm}7rKqdd_=|{aWRnb|N4yC0l{P~UC+|m=Uv8(cGF1UxEgakS^G;D!@03`iES-i z&(hjcE@MVVXlXLWT0xO+k0j^D+H%`kx}K%A*)C&7r%qrkZ5y3iYXi2mbUjOJ7rTrZ zZKW+P$G5osea+b(b-t*UQy_eIH}{_Lz&|+*mu?wwA7EY3++HV@405 zyov;S6@ALZK-S(ZdH$7|ZYw0vcuUu_wDxZ&czpa$QvCPy7>{pu zon2h;<0&@ZDBu1p-9F=C?cTk*A;m*)lI@sHz!EiU|sE(GV zd(sT8UNa4kxquL8`5B)49h1eB4pU@l(?q&ql<8Z+lbzO=ZKOY~mjO>r@0)65XqO2) z&rV-DCSBsB3n%pzcI1mG#$;`tAPMA5A!_&{uYJHK{?qkd%s12Q>u1n%bAjo41r~rb zebMzKIz<4OuEu#n(C|Rl+F~I(fNR3BSh&00hlUAuxB5E6-63Bv+MehN!HaVc(5cEY@%Ch2^uPqYwZI6cIWKIc=46%KVs0?;@N8-K+*nQP?>!jcnHx$RVETdD2 z&%Z+!0JAM|LP60KIkvs;B1@$3s~SzLxQQK@h{Ypac1`J>&t<;iSVoE*_5%52VnO5= z7Pc;!K(hq9eXZeGJQ{9J#KVzp>7aP1rz099dQT`C166tx?`+nDyM6V+Xgu5km#+vN`kG~5X^nTRC~nUX#U*8JKG&(;C)S3$!?BJ~ zD>abl=|re3ptp6&`lR0y@t#E7*Aqn0gyOK&7mrAP+s#%VMK7g>2DzbeO3`xNm;s%OfRAX*56vt(YZ%w$> zM-dba1zUY6vCw6dU8cAbRqSjNHxPd@a9f6`sn}?(TVGLC+fa4ApS;)#T2<{hnjN3x z?d#g<<2bXGuzTK%eJznlv^CrvjE7>W21+AsIsT1RjW@K1;`VsAW!ZYFAKb=D;Z1*S zqkp|*q4?bb3LR&(R|a`Y1kKs)sQu^s(a81~lsdjw{*x}0`X)PeN%^h8oiX3`aJ+*o zm6XoVwoqrxhhX+0lgxH1=C)edEyn>KeG<_i_lh(2r!`e;*ETg&)zw;+b^f($t17GD z6EPum1;rkeuUWGR%(V@TE{MDZ$Xi=cU2Qd1)%dN(zr@pDCs?x7`Qj-pkn(JKM)Q(3rK?;TJocih}*jZDz#$PR40pm({7~^Xf zx`S&#ku4GLh=yVvkPdd(r{CH>VvsT@*xMRk3k z!bda9YJX!zvg?Inwzl;Zja7ASOf``HmPG{%piwMCkO+Aq9;DgcfWK;8pfM#OWytni z-&7CzwN_0}+q1hPQ=EZz*k*b6R~g(p)M-wiJ=NJs54Z zV$l}Xi!&`)oI^J~oe@kPm_vj?s|qVF^b-#tTQ?`koJ4mFtIZIN_{h>HA;6DS~`MJBu!PE3A%@+N)>4yzD!a~(m1QSuJT%|Q5cykuo6x) z4J*>tws;ZjO}p{Cu-CBKdLl7OMrlRINzUJ~3MC8+>i}85v>{^fU{s7#ZA03{$WBj7 z``p)Xon%?tTRKGA!QCp+EYiL}S;G7aPEP|(CI=ZJ!%TbVB$8gsYKbJ^p0v*v5FgUp zf}M$wNc-Z4kxql|uI5-A(-4{WA0z?hVHUob(qAb-I@8R<-69+%9g;-}!?K#H(Kf;} zoea#jR)^ck6FyJrN09Bcy6O8L&wNTCb+hm+Anm5t>I&}I-Ws#hnO=*WZ>Ms-Rx)dv z*V-XG3->`lAuK4x)%4-?S^Sq(B&}AwOB8voYQRPT>J{&@Vxbn{S?W|2gsw28FC#)4 zb96)!QQ^7TwuZ`dVMY?3xS%LI@ofjl6jtHawQ(vW|}0^JP$ea-xi)}UWh~<{2psgw253#q6%rR z)f|k$>TQ%jiXs^iELA`d#pSdJ^mQRSU!KV3{47);e+Zf2*@3Jq$~w=AhHy(FwodYF z)eCUrQ(%Rn(e4N=(@7xGu{c&u!qcvQ1~*8iOh*!MGZS(p0c9LQVF9X)pa{$IKFI3| zMKOA{pv`a6TX64GRA8Yu_k>!)ZQ)R>h0qUnhES*tTDxSa@&cPGJl*JX?$Y?H-03ap+h(F)roLr;O# z1`l*%`WBuY`uz~Hqrhs}7Viw9&L5H`_DBLXfag}7ym(Imy4!Z)c|czZ{=J1x2zVaS zsnOrj97>t+Jgig03cZjzp4FBh*BsP8k8*{cLR6}qh6c}<^eu3%Pg0=fEBe>Lu*8Fw zs)In;npdN7T&*dSc57Y$o3et;X6}!P@29Uknuy$ z31kc+;3=NI=%nQ!5L;VS?O(h|9QM3%11JUPdCu7U6k|AUb(lXOQ~d z+Zjlri_^WD@ce;NQxVOm@gFZjJ<+i{o7*^?|3oEcdoB5ypXVKl)iN(e*2v9a6#W}T z^b#*xJBZ%(zZAITP>8!+WS(A&BqAmj3Zd0?ja-3XY3YnmuQ4>LHX;kiSnISy2;spFBGtohGii1y*~ot1BpHJs{z^%NaXJmFGQghUiI^W39<9wu%j z*JBWX4H)2sr(Y-cBRIO6Ke7qqWn$ zz9h*9aF!k7831XCE__zFt0xf+Sv_&IVTnOQmJdBI(u~%2XO$W;Wa(*XSoL9lM>o^# zrSS~QYHRHbH+NavcoRx^(r6S)ze8`wBZEs-b~qFYiQO|*I<_s`Dnv59gYM^rmLy;? z9Nrwm&DO!YiOPmx4!dvR3uwuZ7>@9QBB-}fnjoER{0XJ_Oa2rhU22&QPbwtRiXC_JRHz{O4sC zaYQ0rSOkj94S$r@{?snm1_LduWo{Oka;ZqHkOWtuBH>7sx1$ggA~PiIME(a8Vc4Dd zK_?9vcou-08u+y}MBC0HbDNWSQ*c`l)2`%@t{0iRDDSe}cwRC~8t;^kqE2yInngwC zod=MOVRhxS}?j=(RLv7gvaIP8)Ww-C4k47bASSWhR+r-uM3 z#ezm!F6sI{%8_g;5srm1|D>?BO19$KrW&hmt%dF7ikb$IHB%aI6F2zRH)8V(|Ur}4>7g-BP=j@*JY~w>yus951^2UxdRt!WJUhI0sJ0u~U z2&OlLy-E(r>`6q)^79=2)%Mc z-s-~l{hcCPx?Fg-Q9CBHcz-L=O=xSKdM@60qwc_jyJ3l|-xe z_9VK()UiU?;fe1Q*^7z1(c6hl<#0=qd=-(mqAh4c+mic~*~^L6;SF~7bOc3q1*Lkt zoVA(~JG`yowl?f+qMx@!LTzo~mN4COW>-?~9&a0sbh76q_e`_bP~JV>P;2W>6zWDt zrV@Ul?iGf!hE}>90VcAmiAh|JX-Z_*5H`KE>5cA3zf}S&M0PEa zbA=(}KD%xu%4XzdH++Lyn%aa+B8xtvNt~s_Gv_$ipP!l@3Ov~rITtPj zeKOrm$^&=$PWjDLo&-&F5?A0Mu*fO?;F!`&k)HEq#JILNDRj;OB3$Jn+)IQdNy6Na z(3Vq|m~8lLOXglDv=zst_`^0S(|i%>Op{t=E~Z;;B~@y&G2#?ym&+t=)g-(sIN3O+ zijfQFb|<=!9TVDW^2cQ3xKoNUH1Khlx}9U8lB1+Lwk3ZVi7r?2(+rvhssOo`J2#c^-jq!va3{1ghH#+Z1+ zp&fXz0M)~;koJPc6+o?o_I=5%4{a{JDVf+lKSXwyMjN^ zxPr%bl}SxM)F`Ztp<}yBojxoyhd8pUj5@hWAKq1pK3?V!^r8CFj<_osfYV$_IF+7~ z3GKM8$T@L za@4JEM>z@GQBJZRArRvY4#nTdR@jj_Rhmt2CeYDu4-DFYmmMRfDkGpXU+I+X;Uw!z zr8FL~z(Ve(cw{OKT|!ZvC@i!q^X9Hnny|35qfNJ^Z!D$Z(@kj<^Ld4Yrw{Tj6Wo8I z%@jJQ57VI?xOBP@3$RQ=CYvDU(k`R+4#(FOmB$#+EvY6spR>(G6hNIy&ySGIKZpzk`> zCoZRycHnY63tF%P*#xK}I~6uDq>`mVo1zt)d0xZ3z?`B@os^6Is22&ZcO^(L={#YW z%e+^6t6+0l4(RZeP4f^fy+HGNi9thRQkIL`)U>KRqE4P6s`5%vFeA^!>rL`5L0%@w zCV^R(fgAH}o@SUy$vH@7dp(Y#TvuU_5K|69}g4b@)hC097|E}adY3o&rclmM0$%|N>q+eN|2MK8qwiI5_6Nel$=+ht(q<)z(dZX2q`bTf3-Be#G{NSR{TI3 zGgjqMvJ6RT1XPWZ^OYd)tm!FY7T97?5n^2EW*nn&QPMa>{zQx)W7FiDRy>y4GYyXkJgtJUkx7JW?tj ztzsf>=}feSgygjorP#GdA&8O^#Ry1m+F>0lAEBl^d4*m)DZ$Twm2{8+^gs^?wXx(Q zCC2h$3&XBpyaNcMHa%d8FvPA!H$fBw73MQJfj5KVtxU({jQtCyx3}Pmy`<#(KuE(U zj*;{y zjqUo(X%D?^{>*&!wd?efX~X7!o9&;|XXU+i-4nN)fB3AqDVg{T62tmsdB}a>HGNqg z@U(C1GxN<^`pokn-58>iUO7c-x+dMoI;3wsPru@P{i;05`8J)*MmpahVH2INmv94} zH%WLMoei&!_aOI2MA2*->)&_lc{(}!3A&q+a<97aLT($JQM0-4a{&o7w zd}HS^BST-N=`-{6qUmO6lNr+M^8ZaQ&ZEjs!`1w=aedsFW4`tFpQ&`Uk$K4Md#IFT zu--e;OMhj4%gBJ%C-k{_55Z{izpeF-bYt>HD(I*!)hqKpZr)a7?yffX7+LxvO~>2` z`>FK|y_7u~PdBICVr-2Ym*OGrrG+G5laX=Aywm)FUUZ&5Z<_f4DI0!OUzKP6Ouv$x zu785q(&Nx*`i70kZ|Rrizxu=zBjjo={f_h-<<+VE2CDszPaE0#3e7AvpU_w2>8sDv z7f;a_=G}gSx~|GICYL<_ys`0^QEI-WH{?N~UM|(?NX477MFVK9IiVx|RPfH@ZeDwoZUC*2LC3$`o=aPI!vefosx=Vqi{N~*$ zF7~s4ul;jg|1_h2x6k+s)(zTk$&yVb#!tD(?g)k3JK};c9@FCS4s9vQ1rC-}F|n?) za+z;7-cq@oUr`x5DxYl{cR&KauW(uu+74glBy*U?a2!Q^QDnMC&GvbESI(5>_ zl=5S5GscyFgvvjn8HaSc{3R*n$KJk;E5Fx3`MtVfn0EO&DdorB-8$uajdj_3sOTPj z-N(|6vK;anYR{$df;yy7Tcn9LoXaKZ!g;AgJ8&+SXb;X65?^J1H?31Y5OH@2GBBz<|woTXsx7G0Iid#3Wz1s z(^Zzz2E>x-WA+Nkw;hP(-2=pYdx4nmt3VZ!?-`Z*0ub9mzsk!V`cH*U0kJKA1Y%q0 z;}N!neu9j3%m88?3xHV47L^+UVjU46_V`YP?gC;R^n8nTJP5Q}TKaV$*6{=o>-ZiJ zOF0d+TJmKX4qrYH^UYGUQXsZqDbUrj#A=`mB&t!g>w(yo79h4|8xZTb8)&91@ev@_ zF$~1=jsr2@D9}vF_dB3UiT()0EkkRhfTT?VVhgerx>%u^Kvj}&p`tBQXpKTQDHH*^ zMwZyAXoEm()r&xE)h~h8NWRyA*n;0F^cSGTGB?NLp!0xO%2hxtWfc%hq2Ci`DP0QP zqtHVN<)=Gk&j-3z>L>!LmZ$`XrK|&DUs*t`jXsC0lzjAQCeuEz&>Q;y*0&UA6MvEAolb3f!Mo01!DVt4aD}H z1Y-Lx#HJ!^s|4b*A5iFeAeM475KFlmh^2f4h^5dEQLrsfD)e)OJXubOEFf;TSwLK( zNudEC9vQx%Xh(p!-F^tfdVdDQdW*B2+!a8ayBerc#>O>3tYf`GTY$K~^#ZXc?^U@& z3VjlYz40&*``{1|``~dP*7h_I+wvY9N+e4aAZ=fLL-55KE>9AWVB0h+F1&KrDIA6bC(| z(DbPeZIePj2C9`F|COTQCd-x*SLjiN{-DrYEQncKt3say;x>N?i0%6?Ah!8+#pla+ z$~FQqUl^!f+VT;gD<%385X*Z8i1}UsVtGGSxxZDpZva)wvUBmMfJ-a^;yIueh+FVF zAojBb#D0zgaS3|tz&cIWTQ&o+;qSV##elEP1!0-3`Q&KL^B;j{r5v zmi`XVVu^mB&`*I_@~?qd@<||;oCZUfmI=g?=K!(f3Lutz4bYWRa)UzG1F_@|AePbt z#FG0JZ3u`ZKMKT>zYD~Yj|1T!)*CM>^j|=WB<*z|x*tpvzgM|`0%E<|g^pgYLU}-J z!NoveMcFHWSl)b~6_U10<*ou^dDj4OJsK6-48-!$u=zN5@Pc*0EgC)&Q}Nb|BVq2N3J{Fc9nbs6wAr z=!-zC;{*`b;}sy*@fStYaFfnDihx+hDj?QzEfDKyP-r6%OTG(;C4U%*C4XMgz5>LO ze+0f;4AKrFdkp-v!{{0SiZgN_GrVaZP@+A}~b`86PxtQ9$uvw&D~ zzCz~%vE%>{3PEx`&>D%tiuOSumi$p5mi!eUw)qJlmi&xD$ADP!pMY5MKY-Rq$rlzo zE2u?4<&uVzz_g5qCTR-cMoA!~LBrzO&c$-gM$=R-CR3Y+p~*%*l?$~tZMH&lRc@K0 zEmCNS%3ZE#T2nA^0>A2BI2hVh<3_gPIruqPbQR4+2peY2q6S9aG3zV`}2GqFtQt zlvt_|tvks#n&<*rD$(6QG?QrJGeB2M^erH&xh6h^RUe}~T)0L{6#51b=PpANGI~Ry z60A#@_C1CEsu1}P!2|y}5upkCp>HmM46yrhs0RmL<3Pc{&8A}c5UqANFhR; zEqu5lMDU+#NylOQX9`}plxc?R#={8bwna}-q z*kpTYCz1keY~Qzq~g7 zSY+3QGQop&x^Em)xy^%ZDDco#r%gC@YhV0CISIDU=lDk_W%9B0sy9f+YpI=xb2>kTzqui8%x%d2Wj2 zXo#MlI!w~$SWB?msi-5FzHf8Y0=l}BBdIJ!-@v8F3C1H~M{crSiU~_Xcw9U!E7cJfG_@=_nly5WI7`6%?4~lNNcAu_&%#@`DQ3#}@KEasCnm{n)PDICH5YV_Z%(7(-fsE({m(*WUQx$@wr%{x(5O r5G52%;7chy>bGh9qY@T6hgzwNmVZs@41d@4e4t5(3)q z``vqfoU`^?d+oK?UVH7e_u2bQ`0=}MeUM-n4%5H`cRcPkfk!rOTvKqLCGk9$>t>a4 zkY#XvdYyuhY4CO{)CdM6;kJlI0D{3%U!ZgZw^M7fREqLgZT1Gkc$K7hqm#hCg2B4l zNVqvvw@6XO$_r>Y*=!W9(-#>A+h|{7n-mOIg(4w4m-;PVTDo)uZ=0dgE=}^x9xZQ# zQ>B8zik6lb;#hgRwVgU+cw+U_YDaT}!OE6Mb8W+BO&QzXPOablIzJQo#a?l-2-0Y8 zxT1A)FjP^|9Nr?^d-KxO*GIV~%v8y5YegOzBM-TwMB4Id8>$;w673x&??>6B8pf|Y z3RcclaG=64wCHi=q0y3Yth}%0tMYpi6Uz`s6Nbo87ptWHLBzz zx)WaT#z<1hapg5u+*--(vGRVkSQWC(wq?1(W1i#6llnDBiFt}n=3DzfSEheK)A_F%0XZsH`#V4c7## zo6#F8bF4q_UZ&`0wY==Y(c>-IFlrQsU0(TfU2TJuIM#lhKE;<{DVTYqg5D|>_kfMk zZvhCIxK;INfkU-|Or-l!E*EeS24p9Q(8nmj&5Kp=W(kI&blm^wwbsAqxBeY(4xyZX z&-SzaJ*}txd)A)x?rsrGyj~XTxAQ$G*%mM;wVyZ-{qS=>ITQHzrjmcJ@6~G@?tNxF&G0+~ z?qdT{_O!pU_*LL{50+psLjhlBwL9mCdly|wX#iD0FNTJgWGhs!@xo#*~EAA^&a9m-rnN?9Dsl4vgso4|}XEN>m-;VGoWFRj-UvX^h}S z6s*c_igs>|w8e#r!2*FENRFw30#9m4mI@qUB>`ld;ux}P8pEwPm`M?|n}VKwswHI_ zQrb-el&HMkBMk+Gq=nqObmMIENh+l7Q)9&=apgqk;glumT;S^5O@noX@YL?s*(;l1 z_2fw{RqYpYmw-1L)u&ct`6cb*o<|?j#mi_KU*9PPU{J*(aH)#QIy)6jDt^HV{PSQf zC5STKL{J04*GJI^93E+Uzh5E){K|5j)2 zCNiuuwOvw=DLj0ZwnQ7PwI}h|IIN;IZQx-{wR}TN4ORQy?F0TN*(O|_{`^5u8U|o9 z**$|~8nj2pL%Na!Dw-pPxux0+@ny$Z%`w^;Eh%I0m5)u@o#h1KQDn*W1nPK(iuAD# zV5;k*90>T#=#!q`m5R9#sBoy}*$yIpg7|=<_Q=weX&s=5M}}+YQjAJfn+7Nr>RCDh zczmQtKD(g@wBPMx2yo9sQ59++=UFOxm?EHF8B;N>0H#}^)xHTvrSMConx2Uoit8k< z)Pxa~rF^JI#tBYET0T-r>pqfzLB8}rqNkA=ooGE62s@6_r{HEe?wLE;p~i^6b-!88 z(r57~MnT&zkP4$-tA2X$} z^}y3(4#r}wTn2RUl~jr%*YK%53KgqS7!1C%BOFrU;ZqfwV_uYFSCJy6X2}|D(Yvlt z^udg7w#=rwv+ z!sx@@iyZDAz4U-d_ii87{w^bh)EZKkQ;3au$(T;HJsGZ-Nn^M_x2YJcMk`Yw&PyoGT-1H--{?JEW87q zabH>tM0$d5o886$T?BKW^!vTOLFpcYr&`fvr{?I|W2)nbp;Q1syBt*vGku5h(SS19X>_SO?b@hDUCrb5r#^CO*pnB@*f z#q>+55EE@!A1NSj@*$;X;6W8zN9Gh6AMk!SpsaB~tzS-}sJ#H5_$71;1aMNZis6{D z2d#}z)CwQg9@76@PXzW>Hiae0b&&NBz#wvEu$vglD6Es%VL9r1`VorUcYO)Z5h1n@ zNvU+BcC}{U{1Q;pUI5Xiab$6*1rKk`eW$8wF$HuUk8 z`+*~_!PyS?;{)Ai5ha2w@fkc0JE8)(gTgjIKy+mwr{Hr+&tmBkJbvlx|Hbyz9k5$_M$|I22 zi#EA_;l>MI>HK@-;EWBw_W1L{zx{g-cwjpJ9!wICsZ>uU2!(r4`nya;^t{|{~bW*J% zda`9V%fO;Sf(l@;k={1yLrQt{a<}zxo?`ErC)?3_mbF%U{d-lap~OyY@$1xvE#5Oe zciA8C2>ce@)*RVdpZlg?5aWyr6R6^wYCtNJq)~Z3uIiApsYY`mHHv|m>lUGkO)rqH z?J9+qUBAdKbp0Y`%FdG{ciH2;((;@Bs&Gr=*8HeeZcjtG%G}TpfF)!OplXHIAGq#A zUG;MUucb1Ms%gJ_AK&>U_d4|6nUAWBAC68 zaj9jOU=Z-{mFus<5-`f9Q#3QcqzhX67D?QaTRi@~bThr%|hxf znK$+3Re+0l)B&k-EpX!00aZeA4nU}S;X{hMhuzmA`$I&QX&i+1uH7R4k>8oJPCaZ# z_zI++${qIfxKLQ~?OmkU_kNco|DI>#tmU!oRoR zCw?02^#7~DhHZnrd9I!^?*KED!R}7@XoFo(gGm#8w99`5eMx`hYZZrwA2CBFN&?bi zcDJ&#<;JqbJ#p}#FV|&{vkMHM8Yrk=1^YMz4)W1!Yv2V|JcNfF=Fj>&x zIvm|2xBQu4!=|}3SehSj1HT9B@wOm+wHlm2LZ0hWU zp)OdEry(O_7u`8ZYmYY*lJG)4+b|bz9y=YWs0(;W2oumbQ$KTOP$x zUabYt!)5wJ#XgtW#8bcXdxhooLscIlzy)})mozLF2*anc)Xx8T7Z5kRi%#Bgl&kcO z&r&tx{jdJ<<>HExB_+PqrSHv3^h%=P(J)@aOz+m-WwjHoerNH_p5Nc}#G#A5zj^rc zN3#$ANu8R<4=6rwhxerydAYo{TqV!Ya;9r|`umDL+s2n1cx~><~@mDf*Iy+-!H?zpQ8021sfWh>#v=&IlCni zZj!LIKG@XUSi!W0Nyf5JLqlUEyFOG`*I1cdAFgj~-kM$A*qmM15D9N?4n=Ak8?sA7 z6?NfH-7=?T<3Cd1mN``$bs>l4n3I+v%(Sj;X>WI=c_vRwvmEJ8oZROiE8XeMG}Atn zd`W?Eqe32`{wiWf$l&M(A1rMPkw>YNSc#mUJy#!na}Ok6UhFoWV; zQE`Ie5~Jb-#U(|>Ni4a1Y+Fw>HgSp?q? zq(-%01T>v{Y#DQ9p$QqtDu9rQ4oEYMNwYw{8cdVfh~(tONM>d@l7*5qB}NiP)l}r@ zI;M?kDa@u?Iz0z^u&g3vLtK`CWRh6lj4Vff2JIb4E)zE+cP4F-kDTmL0%1BD&iFYt z(6I2p<>FErpGD;1S&aJ;qPV`TfPC@;us9oE4}~4;%2C>Jqxm+d!;s+EkZ9VUW0~*9 z5bbZO;^XdV%V)_OZ-Uu?H7$gbIBq-!IlyWbNvw!!#obe0EwOHe1@4yEVPGQep3*N9 zyZJXb+Xq(CAB*|_jaYyRAu;W-P z+S~CTvFPrOH)7Gwj<;jcv+aM5MNhZ?7c`GbKG|+!73^5AD;jQ}9E%=nzbqCVYQHWP z9c(WpdW$YQ(7ry38qRmn2l1|iZ0tAiN@E|bdz_sXAIAc&4?tzYkr0JXV37?oPL^L;8=T@65EyZ2q2kRt3`Os5NePo^>nI z8BFh4Wc?G;7c+gA*V=}3A=7*P)?uXgGJSW6#V_}y%>4oDDWtzndT)*OBc!En_cU2Q zM*2z8?~hogkbaTr1GigmBQ47wY_s06`9HtgdLQXZ@;umUq#TFB&dH#-9L{+FBHnow zL7eA$08iWos~Aw4WTRoNL0VS3$z`OcWI{-)G%A~Cq|6(Ir(FwF+T2LF4oT2m zJMDB?)9iGD$C_iOU76MuNUuV$5{cQ?a+{ts&niWF9qGw=Q$lt+1z$8d7XyorpU6t2 za#coh6G*p%**Pc1g|~sATv+7%o0yzv_S?kI+QgJZ^Po+9 z$|j~Jnte9$&o(hF(HyXeAKFBBqWOeP%z%Am?Vd#QuuYt06D_pmW;Gnm`8F{<(L8N) zdTnAxqIuFLuC$3063ro-SZNa{CYpmbaVv`IMZ+#5P>9U1*t}1NabGYO7QFwpHoulx+jt=5^Hvhqn{lxQaytm_v#1}~1 zNEiabi+ekcb$7DKZ^$e*Sx$Ptv(Y9w>7BIEnr~6gDH|<+Y7`ofv+IA^xg~NM|Gh*x z<;W(#WuvOe8fBB8p{yv5CBi7Yzk@z7lHaG60gB(4Hh}qyR}i;j-5#?3k)6nu51Xqn ztV^0+A$Ep&#UHwjSkvq9k|!ZXwY8DA5Er$MD4f@|m$i`_XubUZ6)){oPTW!;n$kkyK^~=>^6!o5o zVP<>Ij=~Rj$n~LOPcrZ6?~vxGv)kL(!9}T~-o_7h$c|{R@!pPa(q2(~cgN$S@Xn4W zN8xANpB{z(sa<+yqg-;%cm+=RZ8$!!Nw^*N#U;*+yJy3H&{=2v9?w7NTzp$+=eGS- z=9*Ym&*bOixadf_Y`I%5x3fES7SE#`>yl}NhvMj1-={tg#c}E|vimgGHzCro4zF{D zqhE5#_oRL-o8%18^ehYY$MM|=GHT=y=#SH5wn{eXsJ)6Vm+wPSG%(1rKTrmQH%fa( zQU43f{h^)feJuuM@BL;Jez@b0qwxNY4@TjA9ZtkD$8&aX9qgDe3h(WhHVWU}anUHe zvqKKwX!*14(i;YCd8ga+$REW|wtsRIKHR>H_zO1wvG!8pKeF+m_RYj!wDG}q8Czbm z@qu>f+y8FkPqep^KZ+l2?;t*8=l8eo9Kns0b?7a32#>4Lk)m+wK{3nnBFXV3E;g(c zNORnY*)EH(&@w&8V{Nw6b2F_icKVWRtKCk|n`iB^(|LK;emkAN$oia}Uf{L*?erCX z>uEc^u*Bj-A?vy-U~vkO>8s1FH|%sljm2Rr=}Vg|4q2JLuFd+ponF>yWnk$g>C1Op zm)L0^UO865FWu|$C;(NgnwlwYw3gZoD-5gHP8YkZo9*;UkCC#Q1@FV7js@3E^@ufA zv&~RySRHnHt;?iXEFqc^eY=-rQb;wTN0;gQcpA}_S)`c22XKx8Nw{WK` z|B9TvoV@ILk~%;CiacOB3qF~1McyZjoSanoo3g_T=OS#J$v5(deX~|t4-V&~BYzjB$mehqJ*kP-0$t@05u38hw zY0`uko;AyFjApk~-DYS`CUe)eK-t`ua7FI&#>&=uDSa90YsoE@I+b9LFSmHra^H$v zRC9Z%xrz!K^|cK~Rbwl5@{IaWo53AGqq;U?Twl1thaW+;DA;7QHnh}kZU|RpSJs4@ zjbO7;+Ym8=w;RFgy3pnpqot;?8F+OqccDGL+bsaQHT(x~Z`x2x0gIUo(E16l|=nmMChx>Pn>Q8qqcKp>$}*4+o84 zYlRWSFA|$Ww$=gAplPerBe-<+>ZN7DHHGVa!P2r4UodC{8)Qx$%2m|erZhuUE#Zh! zvZ`=psSysfMEC2+|rrciUFwh}+>tBDw)mYS;CP`#w#2afQ; z$X25oyVAi(ZGG5is)M5kt7@yOTRBeQehuNEUH+Wh%EnN0OH6*GF%qf^j!2U~gKNQd zc$*PTZZ)c~`Q4BYr`xQ&Tee31m{j_p9z67)Sf|yvWQmR)`m2f%GC-(ijVOTMMusB6 z5TZ_KGa6S@yFj&ETElMr=5Qp~T+>Lu31TO{7B<1)QkDQMZzP$Ikx<*fAJ7GBTN*+Q z9F=WjXjNnIx{=}w3=D{j#j)aU*T~kfSkjW8Os%#8bSJzh7q6H%k7mU<2 zhY{-Ps*G?|)z%;eeoG_-j;4rF720atT)1+%@|~J+?PmBd{ivlW9IjH-s@mpoC2eP* zSfs^RQ@XTt4MVW>Lg~_M7Ck~Do0~cI3SrcmS`|c5BPR}{*0|MZ2;UxTYr$_g!!o&5 zB`Xm;S*_6mu{r@0DnDMJz!e-=vKpfg)#brcPw&~Do~ ztVUaqVz-P%wo2NuUuKCE{~c!7b5 zG@p+JYN*PhYP{Lyf-O?0l7q1gqg6#}WW&&G9=4&c!%dh(8p3(<3EBF{P0qHVidy*Q zRs-!-vRZ<5;p&LdW@yRuaYGIlSsb%X9YZ7i6~CFU$6Vc(Z-nxVihQGO0ppbPHDjwP z6D)z3)!U@?4>IAUKJ4s}%AYSbV`RWxo1qc`A*R7b8FqzPmS3@ysf zwKbG=Mv=%fVlJ7K#%jW7Fq*g&@vm9x=Qao5bS}oNPa_P{90rzbxNv07jYxy)-J z_bDu+{N$~^x$d_2URpWyo&Eif0?)mF{Y}$9o&MxEkxt!Tu&bwMw&%X5ei1lW z|KPSa9-QCtLHeW}OLu&`yR(0DANT$5#{D4fM{qxl`$gQxaleWCecUN%>~!3h<1WCx z0`~^ob-1_Tz8m+0xF5m&H0~F1AIJSB?)P!0;Dus3?#po(;9h}y1MWK9+i>5F`$61t zn*}zi*ZV=-9zgX*oZk815exAzKqfEl8x#LI7X7@q5lnoF>L#8yO1xW_=Z-Y<4{QFn zktTivZeHB~CSIZQc`8U8e>_57+&(70Ov{&>fWSMjfafK5t^j{dr@$rmMFIKQ7%y&< zv;60^{BqqsegMsj+Z)VZto5zb^>yj`xYf-3&uUyAjsgE2iMY7A&g2iYe7QFVe7Dv| zzSaZZsO8H&B;fLA=K$Hd{xw=(IyvRV>-r4HB|c3L3T~1Uuh;o<2NC#cT_1PWnIF*k z^6eP719EwBPn!9Qb$#6SA)c-E;}HPyw{?5Cr%yab%d^KSX}rdcRgS-Zgoow%t9*`O zUZ@WBIRS+c4GJbk`z2~vvxowj|JEIzYs z+GihsD6eyOvL|KBxB9;Fqra`X=9SdS>;v0wX}SCOn$t~hW}@Mj;a-fJd#6>ndB$}Y z?uT&m9Pe4&ui@so9LKU7HoXYXlV^>y0desd zfFzy{a0?0-PsvE);h84!jh7_80&5bF5=d$sDabQ5l6cgnx%gg95?{16iN_Hnai?!K zATFL1lQfKnCh^pcBp&AJf;@F0iN`RS#IshCc=)PGJkTJCXB=8qiss_!EIvoz#ltdk zt%4_NF76#j(uk_uM7pj_)RZJ1ylU-%>5JPqxR^icD%T^7=Q+5}$&Wq!7%xBeBw@V#*fWFi z@?%c`#>CTaLT;h~m@zG%jBq_>MiR_)f{vc#I?O zevtc)Dn6b7$!i?Gr;%|?^ZiP*jVJE}@_K)uQi7`27>KfGs`+uG8cTiZ}s*IE_+U#u#(_KHPU{cqOSvbE)Z zxw_n-{@X}uhb6iPiyp?SMpbQNP7OaqQ(Q>H@uKN?%vX97H`3f%8EM69K3Gr`FEn`f z;vZd7xO`a*Mt~qt*M}lCNXby0r~)=z%6!3ah@>Pl#r2xFs8uX*o)Omtgyk#~S5`D_G?M}% z9XuBY#H=EbsCZad5v~k~1w|s2JY`}E$dH0cQmh)K?yMe3)W z<~!{Xk#6pDM8t#w;mm=E8jjE>XCem#uh4U_;V9N}M>@JkrA;hk+ zl$2c}*?Gb|*8Ux_(Cl|~i4+JR<#kEacO`ZI~{#2kTazS3VlH)r*(;o zoaW;lu596S_R_arHeguuU`GH*F9evUK`^2JWH$LmJL!Nd!+h2_AtooqMW%|xZWtv1 zbKU}jK~+~ZEx;Nxfli9MC|$Vd7UfdU%cNVt=98QzTx$kcW=kz66oJ1?c6JKj1Z2R0 zoz9>nIg8A%ZLbz?r}^_Z*9nY}E^&E4EDng(UE&&4UW6jkoP*+`fSBARCW0lIN|(V( z0Wk-}B4pFg10vD6LBS&P3lN#c zQte3`7@J^~fH0iq-wnvXS737|!4h)Ra3t$w8cIzB6+Sf=BUZ{!L-tr&D){L57RCp+ReW&&bawOASuR|d=vw@Z2#=nI?y^8=V{ zp-F5BwUKy&6W@+;Laz^QX@6q9y{ba?LHAp-&F({YEi9_F77grubB1;^r@ zFcX?uAd(71`XMTu1cg0L1PzrOolUdr>;Srw?SL8Lir_1?vThIDrpSDITQ!H%eJ0!i zEHHX%5o#~nAd>Hq1I%3@k_$w}A(5yoWT8TxDtLqCj*^_?nF4xK%cuoL8cp-v*Fno z@W2pN@zn0M@Qg@Fay=~)cecsVMiV4rl)9ZyIn6V7Mns1BnfQp9SRfK}V1&(5q=liL z0nrvm+h&*nR=Kkp>Y4X~Hy@#QPJ#J{I}6~>=6&&97_;c=B80r1a)2c|`^3fO?Hn<^ zXdCB~RIMvw8lu@0XP0^6PBzLs6d$z*CL~Xi>#s5rvAG;`7~=40*bI?@KAnJ)--9Yy z<};|5?WR#+UKU71G6kxWWWi1_5~nwMB~%~>8->VgF#;m#mc~M!Hr6EcNJ>)Nk_jTQ zr!Q(Q&QIE*Xm6QZKVrnY?mPq|{)I+Fl{8|FG-3ucv7qj!p)O*@yYXk_h@1*FI5MZh z;j;lTDP1}!a=JLcXP==NQ=Mf_vm55YsQLh{l+m!r?7FksX?_SP#4}WoVtYP9&g{BR}0~eg=YQ zQAR}uVf0O`povB1AHkRi#}GNr!FU&#&F^hPL~#b<9i%^Fy3(DmolrSPwjCPD`M&AO zCfPi?%~b$$fhzE{dB|B{zI6wt2d8(xrx(m+Ti#Z0r5kYG(c~*<7i^TQSqPR#@ zR*O&(;~A_p5hi!uBV)!?z*%6Ok8bUS8Sa6jqTnHsQzSB~(dc`{b7tV!oC8 zsi=<2&ot-1U`%$w-Oa6s&`rJYL=$6{@!o@(kl_OnORYQ;F=T$dwVHj>hYFIN-@#~V z#h_p~-Yjl>U0erps`EFlJcLI~>(C{`{4!`*6ei{Xr;*LCNTw4a{SmYP^9j~%7N$N( zD0*G2g9^n^-sviYl*IMoB6F)9t;(^U0fOoe%j#V@7q0&WS^wkF`d?7>e_Ph?8vF+p zFke93aPuQDiMBZ67DfnWY`D0WwSx*RzWqft`0{V=B z{;Zsk;4N33;6#4I3G|%P^&3#JY|O%x#ih?ocs3#%G?~YC1h6zN6AQ7*Ny@Jaujw%w@257q?yE3fz+VHC4 zH1}aZ&x3csH85ybAuwRYh!!LEqZiozS-|FSwU1R722GmNJdU+J2Q%IQKHH&2vqP@< z;}WswVL0R;=P>EQc{&c9E?B0x0x+7v+S{mtj!c4YaV0b);S{xceX3J{yL^X~_o1^#lt zxxm)}egbe6Al~>VE&|K}o)35t;B-K2GG{pfDfh2xi#Y2oz?sPZnZ{oRq#gr+l>aCo z<@WyG<0hC0rnSJ{x2H-LBlsR{1qVE_c|c$It0jZwI7i7+6YLyt^=f9 z7Xe<3{4OX<|G!N`{`E4}lc@1CD9rqy0h0cN#t&&2(eyeE*J^mB&Y!B`pI{u8e^bLJ zX)M%#oyO0=h$IaIUV{350!X`l4v$NL-v!9>_;o5z|FV(+Y1iXGh`*}gvw-CLCW(;$ zARzf#HGUQd`;FgFv!3IC?9T&$Fv4VyJw zsNshwbUEa`3rPMKHGCM5^(F(($G`U>bPiwzAo=(^Ny60{z6)oh9R5&~e)$p}^wXyR z>9=14B>(3%TnR`&wE!I^P2!E z=LSu`T;s`r)X%B$H&8jteH9S=S^G8Y0Hpj{4SB*l8|B^u6Xm{x2k{?hoPYg*^c(P? zJdcLoh7jTpQyBVbKOlz4{Oy`vtKkg5S;%(+($0TCp*es%;Wr>T`4=?_muWah!vqbd z!(W*1fIkp6z#j;If&L}r|06)SM#E0@C-HsgH^OqjY`FA|fG~~oU(p|gJTN7^L&LC! zB^q9%;am;bkK}(J{X+OV4PVyqyBhKnYvwlsUWR)4g~mKUo&uA85=v38TG)YjDImJQ zc?BTpb2Tj2d|nN+H6#=MIfpghfQEeCD(Ay2i5A6`evhLNpt`?rE$)$;8h6y|B%{wd%}Q>&4yqLv&`+W6TjhDked3_Z(^X>E5X6P^Z;m5p6G;W`d@~5LRA3>MbCvg+E&o@_U z+&&L|Q{(pea5{L&Z(sZ(AqF+xLq;slaTtqpa08EnN#qq6lYiTo@>|E?|1t)@at!{7 zF?eJQo~v*Z!&Y98DO_Xlt1^ZSYjz1qb<|U$7eAI5tG^VK2WT=OgBmVA-IU&zo5SYS-ag`OzO5 z2X9!0k8t^yJ`5)DIlfBA)DPM6Q(JyuoFl*S3)&x(U+9DMV5y-$^Dgz3tTBSiR;>uG zQOR+>9G4%JpZi_7B#(Si9vt~rJ^KCoIUlrt>8 zGXZ<0&+Yeo=P~E3zrEJlYp?x#&Pen#K0%im5u%x)0#Tib4Q3dYx7tgEc*qoD@p6g( z5|4mrYO9V$GAO1m(zI{9Mi~f-roGke>2xl2eJ-Nl^&}Tf#X2$rF^3%l@Y-4WYWG$T z<{N$#a^;sfrQ|axYJ=6ue7S94?S|avZj&+-@nE@BzT8d&^a)bNxuYJznwbtJpO-`C zFqE9SUE+)I2nyfzHPdG03EtomqA;9Cl5#U@^N55khLXd!B%64!zL{6$$vsve#2U!G zr{v(Sl)Gwnp4_Xz#x~`Vq};5!JYoRY@w4Pm-`R8W4wmHBd2iW8Q>lLtZ@XVPr4x9uW-V*XF zJXlG{sS;n%0^kD>woaAY{sR2i*sVlC6b=&-lG!GiR3aEuI>MUUBDn$9=U4!@^#w;! z@Kni4vJ%<_N+7%?h=^gJ2J4-5mfxrF3Yk8>~o`p-+?hS7a9a4#K>~Zw3a~rU>gpP(g4u&@$48L&LB!6Bw zTo5)WNWog%8V{W+Y2iq10v6hXxP?F}GTAnKum_i$=5*k1&X-8<0)Y{wDVpB34kFOt%l8ZvQl3NT@QIfUURM7ijF3x z^s5eS%Ha_WndpP?rcHpk5S35euE@K$vB>!bHvk{Wysv0C%Coa@{Ur+F;p452$@S8t zdk;Tb_n&_wRw_HPV;+XVM!Km8OX^M#9;2ZG&+RmxkkoqkuoplG;-w2FLXQmBf{B6Q zy`#J^+dCdUe;32l{^UdB!w3O)i)}DYtsZ3e6;>*T!R(*_0U8imLM)I=FD*E@A#~J$ zF9@92BIbmSmfY9$mosNJgihpAPn6uZPwG}k${6^$HnP=*#bagY zkPM<45_fAu1f3lRP?aG+_0Q9Uo}@A-!y8VlciQ_zEzxmqV|?4Y4KkTKHj} z__j8kMD6{U2T>l;hN3KU^N)t{j%Imz1%noqJe=j`mXypx$G*z(M~{=~9w3JgA!WG< zivk=T*)XS%ZG!9u_?p$2XYcRoDEUQwIdrmQ!(ut&@-R@V1KY%b#Z3(Zv;+-rRpZ%) z(?xfny$vW9K8fO42A9{mTl$8HfF7nEgBM6km;Q{4CIyLp0Kf<9$>B?zhPCo+@2L z3??IqUPy&EeDXcfUxKNJ=vzcp=;Y?md#6gSpc_UdnfDkmc-(dm%EHdG5y6wSKD{Ay zGIxGFx#r&EsM2tAyyBBe)EGf69?&hGJXaLRTMo}jf^ci4AgttZlN-%(FLc_8!|4=T zs%7_K2^o1B8M%iP?8$2iiV20Nm2Mw3y@a{iGP|AWD@sYA_|S{MhD3PSXio( zh5dx?N3D~p)}K1D4Mk|PK*r=uB^ybmE^Ds@OEsOOV*b$S@@2#% zi*9@Xs8W3_#o=OSRmSkHlSZ3EsU(!r?BsLrWH)5H?;(P6fWli!I}obh7TV|XLS`2c z^gN78`a*lb-hDP35tHlKW4}^0C3{>GPS##3%OvT7tlZznh2tw=OxFxXv2iE!R5L{7 zuhr}v71uQ*6-TK<9wv=PX%^9aZ%*j7lKa+RB>wx@hFLYe9sNWr_W^8T{UF{!w-t1Y z`EYBb8}|mac&~-4SiE;Ol&qya?rQ|>4{sfVHWJ~@YQKa8;~H^Uwdke1R!L4rBG?Uu zH@mw*+y?4$oS4w6BfRw+6W4=QGUYq;2yfoHn;+Zv;IVk~#=Uqrogm-O3m&% zD9g%+yV)um;BkNlhAck_ehd16_5rw*qroPJUPpSuY=iht$6o0SdOV$G<%5vV(+gbc zZOPNSi5iDgH|3~qVR$8WBFe4O^^RDZs-=~uCD^7#ec={0f4Y52b+;QIA7m=@+dKs4Y0Z0i`#1I7Cu^vcr?Rr_NP zfsnOy)}h&EWHH;@vSJOL1=9K6)^T#>7~bqUo5_cNs6+qU;a-uZZgO>e)BK zor-0B8y>Wh`)O53ADNI7!o*{aE*gmWqsq<+@zaOz$T4N|in zixtWRA!i|9-MOoS4xf9*t^0>_{!M4x=L(>;B{XZ24H$XfX8^DQmEFJ1qWHYsYW=Sb9Is8W6eX43(5rt)U>@Fp-^K zVyp~Vlu&;|&*?@voOBN6E(8Zl?%VvI{O04Ytrkn}i<~}l25&#)rtVYw{1{koyW_=e zbfG!K!Og!1K_RhVmGoHEZ${Bg2WVCZQY2mW^h9&lb`sVt}uAFX^ z6uL{gi-pIyL=?D;Cc|>MtrS5*Fg!X3K9vJ6%Yh&G2zYu9{Z$3CJqz|mOM3;Tn>dDl z#72L&Pn5bw8WV(BOZMrSQkZ?N2kZk%&+D*)83&l2)gy(d1WeD*juL`C-Js`99{R>W zKCks+M+}&r^Kp4;9L2&x-$A0hto`SD{dX@(_nLY2n|JYUY z6n*aszJ4|D|3vOR+5NlZyiYuF<@|9agTY6>NY{aR^7-p|Q|FI6{@Yo1&bw@a5CvS= zG8Qgb5tu)1T3~kdjOyCJ)ZR>gU+mhhSYIp^?*JM}b#z}lST}3xteJtSUBX*kT|-?c z(PyVq9W`jy?o?#D-4W^OsO}Io9jm*dsrW#w=AN2nW;LO)YZ{j>2raHbd#;M4qEt@$ zjXNcbbD5UAC(+jxi1%fLD;>WlW@m(X=xMM`YXo!p`_kyau_%j}=7@B6B9#f`Q*HM%iGaIi2iP9FkG_p_Pa_1_UBW$`tXm$tH!`>?npU{iV^i=I=JVKDR=LL>6z&;M zK!G0~%)`^F#qB%z<0_I>JJTyBx@UK|xL%oFJ00r~?m13HLFkP``&=RvpynWSC;C&u zeYNt2D813hYT>?yDDa0V2^6}P2ylrAt7AwF_jN>pUqo5Oh2{~(uVk1KB+58RiKnA! z1fiZt0a1eCR#>L{#_zGWoqds17vivpP~o$q>4e=K>5C#I-HVxK`0UQU%<5zeQcI}n z!t~kkwB0cn5$>CaP+$s+67B$0a7->0r8Zh}FrE>9m-`MMy0p_2D+vJ>?~4kzO)jA^ zI{TnIasfEL(d3@HokHm|5OC6IT@DWYDdk@4N+4y3d%C^9Q~2Edt-w8Ydn6672M7U1 znT!l(6Hr8XS-uXpPl7#$Z{+isi&23wUq%VYYyzu@SQpr-SPXg8Hp1OzPR3+aVaHOb zz63loNg?obCQ5DEWjt`ZiJ?#@toCxLV1=Rzcb|Dbii9f=f@mC$Akyw63pIh((G!dG z^(TcpW$uJnO9iTv?(cz2ze!!8%^^^a7~o>v6?XC}w8APDNLJXLh(J%A5@*nS5=sUu z?2ds*P4oZxs8UqCtsXMW4J-R*6k!i@I!erC{3#s&Y!OX;2cE4Q$9WG;5-PVZ)jK? znl)4GalbPcZhLGrz)6*i$b&y}CSq(6?o$M!lXvz+sPDadAy9Pj&N%hD_Xrz}Y(|a$ zG!^wk$L?wGlh+q+>C(x@k}-kSrZN3F2DYALWArgODKm%7*q#j1Pxf0K@K1n`nlmf z#v>@T1`W@s*PsZkLz+TK2`fRCwD;U4ogdAGw;{Kv-qVGr#ujZJU}%xGy%S z<;L2aWq$Vri+Ypuk6tewJ;IkQd_ zyXP5yVqJ(bnjpv+boFxYt7o*Z!j442yx!PHqmHoX;l|2WrPhUq-0RJ6VJ(R$_w}bS z7Laz~{}mY2v!F!tRRn z_C|Qt10>v^amEfZr2_F;^G9$3DUm|W+?&krqGD|ndJF>a0Ry~nZ#F6XNRHn2Ok!1s zL;f6-ktVA;GKt=HhrET!>*3X`nGxyeedb#*vI&*K1RUv+#I5ED5RoWdu?#XjMuV?# z{{#ETji1VK;eJ;74RQft$@;L3>E9-Lni*-wjUnbdPu_u%b?*VD?;|>uKsl$xgUtD9 zh28DU0}rt?uL)k0P@K5&#_A>E+3P}Neb8G^b#MlbML_|0tt` z%N>YtKW5H>z)4|IpCpi#T`Vk04H|Y8J`T#(^0acXn#$bUP0F#0RZa=_<0cuGZftN$ z`o$%-GxcPU-Nd+*<`C{yAZ!QhwP+H#ltyc(dwrkD-@u(a#4iD`C#DG4@!n*ADrP4$ zXv1n34Os#7yhJ-%JC{}C+5nbWP=Hlm1Zz&KKrM~O+jeKPC*Iy`55!Vw?6`!xfJULh zPn)Z-7racf3=RBr$y%SppWr@j*d8N()@o^GmZ4P7S6YunwJUhB${prTa`#fg?BEk z_$LAF>cIY5cyIV;m?k_KNnryW>k-}t(pLTi=&NnIVnyVcRl zl-OXBt$p61Vl20@s*IyEbjKn|;eEv6boQfGSaIDgyu3Xkt$lD9k!)X`!OGkFv=hv} zgdObx;XAr47HW1n8R>}G(fB|-oq&Viccih+LI2OB;-uphhl;_w+@iC}qRJ+T3qlEov4;`3nn5lanV zInW!yFSLre5`IrdqOUU^#YQ3~ibt9d40!0!qrV69rXxCvZD$oEyW^tRLw?lya)g=e zC32HX#^_6cc^vN$uBN38OPWQo?=V6&s;Rha)o7gU^}+YzF|rDgbayn4#YFK1T!Sh1 z5U)b|-bvk-(nu8hiQ@Mpd-{5RR3?tyizpsTQh~FiIuel-dQn#H zLXxYT@-JlTB7$o@IJCfF4i5ANMe#Ty)n!SmMe)T%3QlRI5|#;5@zH}Gr zJL=)$@df2*66s-YB8ri#nA1sk1}I0fZ8ed~Tav0hdO`+0 z*(U5r#5z0U9dX(@7tbf|dXH{+JJXFRr906R6~zmPx5*QWMpt9+(ud)QN`#2IO}IkK zRx~v?EL|ds8}3Dpjq(>SYNqz!cI}SEySgE~n1scYN#u7U!ItDlF}h+MjH$1936Xum z#l5(AX*qO^@)a+CkaSTikZUp}n^_qLFiTv)T7_$QbA9u2QG9coEZ_PHl}B^P1BT*T zWm!}s!}Y#R72`8{QGEMOh?f)XB!Dq|qoHBU`#>?a-4Bk}qJDM-&i+NoxKeysHj?(n ze0i-mobq$Q;GTA_o45k62}H@HrNc^B0IzCVyrMB=H;0#nmWP)vULXt}53)w8>Km5X zH#aPwzjS#+bE`dQPpcD#hkPFC!peuQs>Kb>;fDG4g3uBcLIC^IcGHSR466=N7+zTj zBZMsE@EhtIL(QROSs&;?4CNiAXYqwmBvsMZQ*>1=T|!l8SO)itB`d%Tu!g%oefm|x zD3L;yju6`8W_$Tf%NC%1Glk)k>{^EnPP5&7%Tl|kK^UcyS?4gX4>i}LSf0#iDHC*L zs4UBo4Plf?cGzJzEx5tvlvv(yD+R1^>5^%};KUf|YI3AV$skQ(m{=MXHrC^QYe*Q} z;v-%3#)Gb^W?@_~r-GWNjc{BznFQ{&Xi%=|N#IC@&Jv6Y%D0iOZYdB`KVf(QZ#1PB}vT0%LGU=TOU4O>hPM7iM!LfKd!;-`&%p4Nc)FlajopGzWD@E8_ zjQMKf!q~^#By(tS&7D~!T(epjKQ<{s=1{O^&_P^7d)C#+VY3P2LXR+BVj#)|96+R@@=)p=DFDuKt*46O111BBQPDMf`9MZ$< zYPvIaq&+bZ!yW>IS(XE9vEF2sux)J(y3j8~MNJy1+VhOvL5*!yK!ycn@n!JjJ5?zOcTC4})ah3{Ry%z~FgFz`lX-e)p zCCnkKTFz!lg*oJ_<_S&X{u%VbOQV6>mk8})gXU0+YQJ8EH&FujDI{Doh;A)0*AS8? zl7k|<+jCH4=pZWM2qUe~VG^ptEXb)X%3RI*>2)Vv@5={j`ypyu4f~SMb?u2`DX0ln zH3uzA@KtjU%M$$6*lFaH2yi-Ljy0}TmFkuOkSKRjikB@n1X`+%6CJ)@r;4>Op*Bls zs%?einjfxh%SD4?rfy#pDqt=2vX`zxOqp~UkNy(xjAi>uiER5ixWyIHDyV&_hVaah zPz~&$s(}lg>Z%4F!T3nmpfoLN`Y_U@DZFWHN;?rCOEp%MT5hC5+9)4JIWZ8^aatF_ z*nqQ(NEgEVZw8Hf^kSQ?6Y_z2y>v0^wF9W}dzL`GW-KK1s6wskLEhyFvY4RwA8E~Mp5nP2!w09$f$z` zk-5=M8)q7EycWef7~^R+M(^|KdPP1^-yAt|nB2HWxnW5+KB?WH+93tkI)-bAfbQ=QT4=h0z&a(+L^${k9;O^y&u7u* z5RDEr*5^!ES8_S}n1Zg8+TFGQ@`|=;hLj*jdm^iGjG+k_5Y^UAlYBF1_oh*Q4Tq#! z`4+C}zDMAdu$;klKX(XCPEtvNYX~D9UL6UeMwJKxMQtR?!{k5)-O0f*5~gCS`><8O zNm+kyJNa76MS1X{n$srH#vu@$2}1ESrMR1?KxXEIX1e0GiK}#3?3mN>45skEP$lRM z4_!Zy59H?S0U>7bDqoJLEQZ+4rk#>3ghZT8dngCt@V5-YSft8nMaH9-PJjZnHd~~YX0vdZ zXDp&k2Hc!Q?}Uv-a+T(rMQ8Ugu!zc;L;hKmw~J)&2~?q7ZlMi|Idpp!)@%xHsluy0 z+Zljb$Tio{p|wnP1|dZ~xVEYT9~sD{p`z_tSCz(&E04Bw-CWwe8A{|~62VSa7_GA6 zO>^6b&6sOw+q$_t$|!Wxx+>aPQzL7U_3NrsVk@EjYpZk?;FSw%qg03I@-V5yPOq&} zqcb~*SFkuZL38lV4hz(Ca1l-BQa;Y*L0h4_*Ueo>UMuwE+NwcmnYV?mp`&Z7R_758 zudU)$)a_L1?83Z?S5k_Yghsx+7HwX!KQSg*WgeGxkyUDp9^u2v-<&8sp6ftDyIde# z)?Clko`y2KnlAwzZfEIh)jv@RQ(s2uKa78H#Ywl=pi>1@VEgK^yV-rz6#+W=- zdNjt(HZIV`+~jA4Qi^(!NB$fXTE=RJSolI!7e+1uLe?1>hlpt8;;a&UNLB`~uYoI< zk}G5~FjLk8Cg+VM&iDdU`#KMW!qInGC1l27u{MDLcs3DzTuf3$At{SqKE@!GlWCJ% zsT|@~nUWI@l)fUjpb6=!EP>D~XNrYn+(XFJYGdA54z+voSjvx}|MS(y@UItk)` zkR-(#$_T9mN=<-@5jyQ`HD8cwrsr}R%JNLifQcL>GtVWe)Lf-`9%^K%M)77TPZ25d zWwz3)(7JP#RCSz_%gb|ZuBx9>BCbTx<*ospCg;li^ecxfzDRGU1pS!BFg_^8FaL3Y zPk7$9XQHM*ovOk6s+z8jjv9QSX3v^AHQ6y0pP*s(HI zI3&d|tf*ZY2B}ws;qqFq_v2OK)UqMhNGtObeCspHde%B%R+gF7_^rUrleIX=G%n_t#0$C z(tkFaO2IK7qMLOa-Oa^*e%w&Tv@W{4{$-t6WdJe?ZzR0Xb*B}OXrT*lUjy%2W&isX zb6#0%uhsnkGu|$IZrIbBcILS~*2GrV#7^tmLti#8D>Elw=<=9X8rG2=$E_Ju&)PC; z$3Gr1ucF5vT3tn*R@sA}u%5G?vJO~3xb1ms$L)(G=1N)5x++e!x<+sIT)^GBWn zFy0K8!yEyjG1hp}YtA3T>igbywEL}Vz@LA?b7oW71Fb*%%rn=yE-|MX)Y=cbimeuN zLD{?35wo%!F`!7yDkHbw#q~6*@dp20d&D(zn|WjTS1f89*AU&z>li*puxqt8H}dtb z!ryn{;w|Q+@}19_bIYu6-eOjjL(cpnZfMCuN+)wlnIm$`xkNlxz$$y*TvldIFNdw6 zp$Dz9q34hM!K^DaYv5_=5^KnM!@PNnIjIzB+TilGn2#EzZ_@48W?dPww9YW=%gj15 zZ4{bwM;C%D;p-Uwra8ORoLG+U2xpgC?^|ccMYEE#(B(CkP~LF{&d$rgw~%(s`NBF~ zYf43PG=?JiMYEunob=h9%^&sg7J-j$clKs8d#^QabTm;k051n)Dy=`tPZW0zU$)6RIXmU^Kp*oGwEtK#Bx-*6_1_Ry^)=nC%anZ1!dpYaz`eNT_%=txy zi@s3c@)y2Q{}RfO+(aep$T{nFN8Vf}B^|lpTL~@qy?nXhJ2gknb=+VKr+W_iN_a@g9C8 z=*ZtB)rxCWy9__7X1V%e{|48>;@xK0Te#a;_=Q5(TSafwzgYj$=n>;pIdK&dVgZl> z(E%vL(6fLhGxWNKW@FNlv`01c8x38Hxsx`51)@zu>j2%rv@M$Ugr@x*5bY5P#9=^E z^Jxvy3_XdtyKD+^8Y?Hd0#a*0ir_=!Y75Q$zouA!CH2uLzK|F;+tZn)_)DJqD3}46frf6?P!dp>%RU51miV@YUewUfHS}9R3%SI9WobAQ zl4YO3nN~e>PXUrN7Z&W&gG&I(deBx%(&$Uy1zh4YfIiISqDg#-L9c` zG$cwK?j%5A*1SMNw&tcZv_V7kS%@t0vWDK&+%uYX3Eo>txq3h{8aD%~VQZhzTzbDH z<-QGQ5p(~n<<09@YK+?uv07)B#7dSRb0m-^mYUpY} z(%O7LH?c1_YT8X2x($%D8wDio_G#`QAZd3!AZhn=fMnes10?JAEFfv^7aIBlAZhnK zO)JJBkF-_+NLm{YNS2)jNLr)U!qVD&K+@V0K(dzG0m%|S*U;nGU`yJw8hTwrH{pas zYQ7VYlzT+eb_0_A?I)UcM9YoCJ3Lw9DL``6dr8x7!K*7d{v`oPxm_AMtfBFE;Ul@N z8hTJeKhe-gyg-q1ivb0>-7M|tu+c~|1y$L zFDBs6MH9AyEI?IE^J*G~9L1#tspKNKiU#Ku6wus>npUZyD>b)P(`IU@PIIr(w4jD4 z?$R!`l=PW0L_(AylH06lEgHIAbK5j6s-bSp?a{QPhBBHvsA+37v|e-X)wE3-x=(Yl zj1V*M=MkS$Ttf^2n$5JR3tvGqT0_&&NRl=Z zqs2UEHpD_eL5BJOQEnLG0S&#Zq0y?Z8e*!3mTBm24gG_LzNI0j?;GOJns$lmp@yK| zE$!Z-p>9C5@-Xo2ylhKD432fY-=v|_fT*RR1wBtfVL-B$n>6$V4PAXnB;}|# zO6VO8l_GhiTqPjMjcaI7L)!sK?~ehJB?Nk&EK#DN>ov4XL*LcVTY!``^hQ}Wq@k#W z(i(aIkSy_2O{+)$m9!Q>Qp*+%4F(+UR{%-wGaCAx=F(xV)bb7>DOY@%16{14xQ4!? zp?w-Mpn}%khPVKbj1Xn}H2it^SJKi($Z)KIA(r6}1YFA*LP*86`16PW?y>~+Zb_rw zEosynCC%xLhImxdc51nOnzmolPH5UmO)G~JQXl=cUD~)t(}J27)inD5VM)3BG;N!v zJ+EndHSKLpJF02b{9zN^!1du7A-arV2cwcD?pIG2p7bZg4*Xg8%Tm4t%%4Z0ng;SH z?*psMql|(twRseZKwTbXA}GN;3eIwPrp#KWcBS^oQf>vMMN?=zrmGEqQpWkNt6Nh5 zh$KqUWe6M(DOghE9dDAVXU>^ECka&Q!3!s%kjuhURQ!rdl`SEPL6OQTBPqJJD^fAw zbIHgTo$@6o`@^N4MkzWL$x`-$LI)s{;=FRFqmHc1&mlvH8d=J}fWE2N!2^+x`rx*i~Em3|B}WcNyxc}l?)mYdgx!#IQzO;$)H&wL%m25aztZ_ zgocL9V+g_7F5RtU(DaabHK%kIT%%;r^pHVQp5ywglJT)j0GK~diayDui|S36Dd;m> zy0R3)=-*qRwxUa_xM44%&mC z5I8VA6EeFe-gD1A4p^Ovb>equ)j2=iP(ltT`x6djv6J3)iQ%765C7yCUrhh)-#TYM zZ62 Biass, S., Bonadonna, C., Connor, L., Connor, C., 2016. TephraProb. [doi:10.5281/zenodo.3590721](https://doi.org/10.5281/zenodo.3590721) +## Acknowledgments + +- **Zohar Bar-Yehuda** for the [plot_google_map](https://github.com/zoharby/plot_google_map) function. +- **Alex Voronov** for the [plot_openstreetmap](https://github.com/alexvoronov/plot_openstreetmap) function. +- **Wim Degruyter** for the [get_mer](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2012GL052566?casa_token=OBXKwSpV8vcAAAAA%3A4HJhAV0JgGCN2sLJehl8FMWA6PU8oxxLWwEQNBJbA31-M-iH6iz5ayjHwT-GrnWHEjIkOOUTFYYOKiWM) function. +- **Francois Beauducel** for the [ll2utm](https://ch.mathworks.com/matlabcentral/fileexchange/45699-ll2utm-and-utm2ll?s_tid=ta_fx_results) function. +- **M MA** for the [wind_rose](https://ch.mathworks.com/matlabcentral/fileexchange/17748-wind_rose?s_tid=ta_fx_results) function. +- **Scott Lee Davis** for the [Google Earth](https://ch.mathworks.com/matlabcentral/fileexchange/12954-google-earth-toolbox?s_tid=ta_fx_results) toolbox. + +Apologies if I forgot anyone! + ## License TephraProb is released under a GPL3 license, which means that everybody should feel free to contribute, comment, suggest and modify the code for as long as any diff --git a/tephraProb.m b/tephraProb.m index 1ea7c57..b82a5e5 100644 --- a/tephraProb.m +++ b/tephraProb.m @@ -40,7 +40,8 @@ %vers = '1.6.2'; % Oct 2018 Checked windows compatibility %vers = '1.6.3'; % Oct 2018 Fixed problems with ERA-Interim offline %vers = '1.6.4'; % Oct 2018 Similar wind profiles -vers = '1.7.1'; % Feb 2020 ERA-5 +%vers = '1.7.1'; % Feb 2020 ERA-5 +vers = '1.7.2'; % May 2022 OpenStreetMap plotting % Check that you are located in the correct folder! if ~exist([pwd, filesep, 'tephraProb.m'], 'file')