From 2bef32e47c44fa387a9bdf0a0c6dc366f6684836 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Tue, 25 Sep 2018 15:47:17 +0200
Subject: [PATCH] The component extr_map of type snapshot is now defined with
duplicated values at longitude boundaries, in case of periodicity.
Add argument "periodic" to procedures get_snapshot and
set_all_outerm. Duplicate values at longitude boundaries in ssh, u, v
if periodic. set_all_outerm can now call local_extrema with periodic
true. In procedure set_all_outerm, we shift the longitudes of outside
points if necessary.
In program test_get_snapshot, define periodic from step and nlon.
---
.../Graphiques/periodicity.odg | Bin 13595 -> 13587 bytes
Documentation_texfol/documentation.tex | 13 ++--
Tests/test_get_snapshot.f | 27 +++++--
Tests/test_local_extrema.f | 10 ++-
Tests/test_max_speed_contour_ssh.f | 16 ++--
Tests/test_set_all_outerm.f | 8 +-
Tests/tests.json | 12 ++-
derived_types.f | 8 +-
get_1_outerm.f | 12 ++-
get_snapshot.f | 31 +++++---
good_contour.f | 8 ++
local_extrema.f | 3 +
mean_speed.f | 1 +
nearby_extr.f | 3 +
set_all_outerm.f | 73 ++++++++++++++----
set_max_speed.f | 7 ++
16 files changed, 168 insertions(+), 64 deletions(-)
diff --git a/Documentation_texfol/Graphiques/periodicity.odg b/Documentation_texfol/Graphiques/periodicity.odg
index 36b9ecb48d5d7003dde801837ef641061a348dfa..0aff6989bddf09d78b7c73784e04df28c3d41ba9 100644
GIT binary patch
delta 9663
zcmZvC1ymf%w)PO5;1)c%YhVZx+}$BqaA!gw5S$6_I!JH=!yv)kT@qXd5AN<B;7`sy
z=e~R2`>VUl_O4ys)qAb__O~kC9%7HBu7rq02mqh~0J!f3;;~fmeitg;1Y!JZcmTkk
zhXlgSLmnSW2Sx;_D{IT3AS0nY$0fx;C&a@hCM3irrN?|JNJ>FXMaM!x&r6MgLWhe@
zPeJsG1dE>nmxhv-lb({FmW+{sfu8*}KRct~E7sRQK9<+=oMM{1Z)5~S)dAK35?2x?
zA68*06(I!!VHI=nS5y+OSfzyoG(;GsgeBx9-{{NonSqss>~%zBW#nX)^>me_banM~
z^z`5A=@}Xr$Q#(n898ekS-mqceP{Mg&h*1O3tK%qPX|32`*&I&jHR7TG;B=u-AtAJ
ztkk@0bR8V*OzqulUA){Jtb9L#+}+*ztf2a!KttzXNvC4jkKfchE3|wn?L5OB!4NNS
zFn1ugR18?kAj80^$mv7T2ffeL4*toYz%M@jp)Ntm))7tM*c^|ze7}T3@6sTjjGoW7
zn!(OSLC$6o?xxY6mchOr2|jj10Kh0N?j#G#ptww*rv99?^pcU$ij~!xi%Ud+CnVS%
z65^c{=9d}doDt@h8|4L!_YDgT4UK{%#3rSsrh)_0N@LQX!Ku{=8L-g$(4_XH=-hfp
zZh2bXx47J@go5V8!j`n+x-TW|NlnE$u>rYB5d~=xWvL<1jD-5E$kMNg^@S<0{QTUC
z=K9jS)~dpW#>UUtePQ`OAf+Q|C7mhdW7!QqQYyyss=M=Q2J#z+ztwe>HujV>4Z@nI
zs==-OUplf%M@rjf+e*7;8oGYe_RY6XEVqrUwa;u$XFxl?B@gA5^*0nv6<5RS3zut~
zc2iPTYD)Jit54h7=34vudm3juOQ-uA_Gj99d;5BZrn(n;`rxzU!y~i9@R^CBnVGqc
z(Y3*u?TMeeLn}kG+p|B`Pv_PSr?xJpkHGg!TPHtvFXm2u?eu(K9T?x7?%kaEzCJg&
zKij|ib7&R*b8`{Cvh;Ird}?>*=kD_O`M|){^z_>N>h}8b{>J?I_QK}+%G%c6`qBE@
z{{GSK#{R+K%J%8{-r4HW&F;~~-r3&%+4brE#_7?)&HDP?;o<H1-r42Z+2!r;3w(KV
z|Lgqh;r{aR@$rYa(lh{oov9!#q3tn$n2F;-uu0iJ*i-v{F)@sJO<lqt;V!XJC0MS8
zc9=cv<C~H2?$aDwtB>TdBgHpP7LTdk;_NzK+{g!{1$Pn<ykS3C9m<pVQ%I6an*=Rm
zU%VVc!&0FB!O+;!+)U30N@%`ry#^l*8&B8Qc%E;yZage)U*ESr#V8RK&T0c`LDKVK
zH3LbXS~xxW9Rk&6MENf^w+eV;Lphb+s?PtE>#RwszTM73&h|mEzCU}{wsqxb-5ae7
zxq7ri=RPE@aMI}}uDC8aHvyksCTo>lE*&GDhSHi^bVKyGkWVkM)J|4DxPiyaC#Y_S
zaf>L+{9GnU%c)YEHXy>u2j3R)>jh<4;unGC<=>1kV8-K~!!w^<pOr8i4gk&bvQ=IR
zcwOuncIeG~7EbPQxGw1bK2fGm+xF=mwZQMw#TZL}$6_QqBWSN0Kek(o@tsRwLI$6S
z?tD`0vk#XOt$DCWJ2Bg~E_3h%OoT(jAq%-9qol_5(2pDT+=u&a_UEwzd`S$*F!O8f
zk})w^)8lZ&klubYbA68c3d3)N_BfNamP4+9bFcLX{(XeEnSYv7#l&^J(``HQ>&2CM
zKhpp+$EHD@xWvNIN3EPjk+!`FCa?YZ>*<H$P2`HQk=pXdvlvgjiYxGh;Y>Rq=BJxN
z;0t+K1A2}IkT;Zlg5LO3v$2GMm!|Mf!^4&wTj#vRCA6<)icZGu@vbaNh~7+H8Ezkb
zvDYQ?2?sT?8oRY7PBAdbgx(+x7_Ys{p;pl`PNHejE;NJnwZv+u7o^u*9nXi-7esR0
z(lP4Z$5uOkN9DL_Jw*YN{!A<9ln9y+;7?X3`Jp8zZm_Q87+`raTEuIRR`?CmxUOVG
z-Hx>dlhPKs-d^}fnoT;NpP2Pmnj>iAr1Z{>%u8uOy_snH#I@Wm#h0q(N|N2~NeC;;
zfNMuYUVVNXE(3+;i|LuzG{-~B;~`hU!z+iAVT#_!PL<Ck=9j^p%c3S<;g5DKcE{;@
z^Av_UwZjVP1`!ms^ea3uG0X1!J|}91_IxgdUk7W9>ryAeYMvbI@TO7U@*^2Jvi6|(
zr$hs$kpoKfhERI-HZVl_6U~hg>e8LPn^|xM)wA7rG1;t*&8KBHlQhlw<5;W$T#A_J
z+alE-N#bz9ffL~SCVcc$1SPU>>$}hzN@8X~!JdS^C0L$oY!T*JeYAv&wi5m8H5(~K
zLsKWFG1=%(2VL=7wF^4zufoP;W%tDkq!m}>9jX_|37Sah%d7Lv6%1|b7bpjeqOXAi
zXeQ^14rzk}l)p%OCJm77V(p24kR3aOAUy@2Y1T_EuRMW2!V$YrT<4aCj{I|FLe!fN
zdS9n@&dgVli+brt*v9!ED>sl{a>pXTiLGDz8ezB@$%9;%`OLcb32xF}*IP})Qpz<5
zdI<SD#lIbH5JM#y5{<*Ml?LLBhtW(7JzTtXI`?Q#RE*eJq}*#KAZdx#c87<D3kwT9
z%U1{@KG9%7wd9Z^dxf4!3}GT2hdV&moJq-#DaGrHlFA8n&e)qvoRRdCBpL=Mbxh#J
zGPw*k+od12MgCW$CzBKZ-q+ux(J*K<8&9B>U&^=L9kK$#u@LDv)ENV%bNY~|gVdqc
zb6%Ps&$%T@g8Y*;;G#Yocw88@70vSGq#n{6g@@qMHaEZ87ye|3ug7u*)JkuC;;7Dr
z0Xf%j<_5_=!Zt>u2Ffflt(rM?6jDYin&7#a?t3x_`+Ofn@)ytf4Vi$bOxqP&jKOT4
zd45fRL;{s<FOx{{Sh)tn{1-Mn0fI{hwh<5rbQl8NLNDx*QXRwpg_S!!K@AfV*R+6l
zF@g)UfsSkAXb|=(37&4_Qj%I16?^K!)AzONFE7~?KBlLT=`=;dKiJB%4OtZspA-tw
zpq3~^Yuc{!d)hBnLI#Gbia^9?K()8E(HUp?3RDnHrB{2CR%EtHR<JGAVlYan*@(KI
zeT!jh^kZQkc6$q)Ic##Y=rR%pO{Tyn1ViSHI2;x236;sGpgh)4enJEEz#(!g8k8p&
z%{<e20fDa)Dpt_spe7fET0LGmAVzU?@r3ehIRohz{yw;QJL$|X(n9_$qrs|#@=ldj
zIBMMcNsV))ocTh|VV^eH&RUJBDaB)}lt_d4?<OVDe(j&^Hgk0yPv2!EqrCY61^3fg
zBc2TeY;$x(=2OLLe@T44X`{)1u$rH#?Uz1RNQ6%N1Ax<oa7qKX5q~?xD6ztET9iR#
z0&s~z_5W=2J+fSBYPb%27rmpkeVUY?wKkA4y(wl-3d_Y`c&k4SMGdq#o96stq1eWC
zR|3mLBt8-^hdA3)KZP$lx<2Y&6f6p41ht4rY6;gcq$t&|X$3BRSOw6SEX<C%ZD<+$
zV14pUM+5SRaQh9Uwc4P#g~e42m_rLyu=&do22G~?$!bI}EDQ6Q<fcV|Mrr)^gy`PV
zU+a9P4?Mrw$E<1vE-7?X$q||;EsI(G+FTThlFfMy4cC|iVJhH-(~iMK!BumPQ){sx
zUP1du5<L6erKRQ<mj{MobO9>SG#0VUcrHz@w<8gUvv9OVvVih@rf`~|>?WdgvKqb^
zLSV&eBH6+YDT||h)NWz;cOhZ1zSy{|id`S9nliAoGCL>G)|*hkyzSHC%>G_q^l@K4
z3I~=$M5!6bvMz-LsO&8brix}?R*oVo<3qA{pPHF50g}zkbWxpaD3yafF88Ju8}Jeo
zy7;5%<Uq9-1-;=4Sda4p=5SF1w>v3wp>?3LJt22ydFERa5vRwNTb%V6*%b>hxuF?x
zoR#NCC|^5gObv-I339Kp!*_oaDDl?K$ZK@h%6Zto3=C}5TstiUclP@kc}s`KQ3?wQ
z2?+~Hy@mfaAZSl5Qa!-`UVFG<k5T~sa|24F4r#mH_Avf%jZ7qmHI$evqbsmRz!p0X
zm&p|2ntv3y_M)&Kn?Q=8uEi;@n(KSQFNt1JxA?9=*P4RW&diAmt@FN{<?WxJKfMFG
zi11_7D#W^7xw?TvP7Co8L&nNgRK>GNkxl6MNhH439aYL_Y)dd`k{bfO;5!=WY~wQ?
zk4Uj9(L8j(55I&AhodMP+LHO`Y6Z_fxz^1-`paQ?@pZ*t)JEo`d0w+6$%Cqp96663
z$)MHf^hN%tpi3;R5e3t-X!RcG*C-lcDtiz7*OK^&C@SFSQN|%Og)wTgqK9VI^mmL>
zDu+*Gd-5D24K)Lh6QGan$xO8JsB!E3t4Qgfp4f+zwrQSEYu)MmwIlV$S1eiO3<?T?
zv2k5}cacn^H=l4|<;ZVfc_FUWiSKNN<W#b2@n!|O^yA5N8sT5?Osc7Y4Y+m+_8&(2
ze)nbYIe;}SSf^)JGKSOOIqVa^GpO@4v`?F}mR_r}4dqKjKFYeIZaxa$qhO?FV28Di
zJ?@CLj!EFn=tXiqEJh`iVfxh~qQunP_apd@g31&v`Z0V*US*M%L+N*BH>j!As-9BE
z4%Bg+&WDl7vn#7$rMuBV_lW3qVX%PVt~M~EFK=|N1O_n-_7h}yJikseyB4CzwomCw
zDwSDUaxL;nE-LELDj(KhZv(QA@FfS?{m#>Si?I5M2{ue5yjF4p;)b`96me2ZxW51y
zYJD8N&7W<fYH*_b%D`P0;2e<2XUA2taM3O42%mckOGI=D7Gca!7DG7McZ|1s0{0ON
zaNDrmzs|zsOmKkqr6$i)^Q#(WAHL?O?bJ!I&A&~XMM!Lz!`WC=6g<~tyLymMRmDD1
zMW_ic+7tp97={rJDRYo3Ab}<{tzpU0^kn(tFj#WqbQ3jI$D&;vFR)gbp2gN8YSvf;
z3#x0*L8_RxyaA_o^aJI|L+;bS6byEt31u}Ez8_uuYm1feH_Iv}mgs1#F&|(BC><e^
z_D*+=W@5S*K)@2`Ec0SD6KBsEMBTEM+8V*W@)=~q7Hyw^U0j9vt#Zx-?ni+Mx~J$`
zhq$6~>I15GU1su=(FSVvFQ3O;Sn8g-kEGoRyaJZ35ZAa$G;3Hdd1)wtV{iL5Yq5B(
zvp7b1R6h&G(rK|BAnXKS_0Ue_hG>zWW7>)eCq1q`4owQ~z;mUD{a;CjTss<(`f!N!
zp=V_yPnh^!w<&Ropwy^C76f`0cLG8t(GWRW3$60~xmuf@wB!^CYL&^)7qJx6R@+ci
z70J!x3fpu|{kprvqPfSw6R~N}-#n@Opw71$1m(%jv_w$RM1El|&`%4OrR&Wfw-+)e
zCW-Ee(27i2oH;2<3QIMblcT$V@2oK@>Xno$ASQG&xDW>N;nm_(JUA^aTo52jqm)F`
zLpS0cnYCT81!jY)ov1esB0Vx41Kf}Zd}D(2*xdwM@N7)XbJW<uZrDM06i=5|jT;!|
zb`@&HCn@TBPin_@_V#!2Mgyq#pH9m}ow`;rt+3A(lSJJFdb?ZB%v4iumQpdWLv`)K
zWxY8;T*hoW!D?Vm=DFH^HMVr<MYr)>qDg=M9W=Dv%eQx#D>rq%Dx<0|M-)XLUX_bq
zBM}Ivk&yqm1ziMNe<}a*2xA-$^Vc1rt>1GMQmQ#s(U@<Q_eJqL8P&Ix{~Cky&14Yz
z@-?tjJAG@{AudO4k!rjkw-?2fjVM1!8Ik0{Ag@Uz<e~z;Bf^GGZaGaWCDfSKB_B9V
zKO9P<Qum!4lR_#>*_#lX0Q00aN~VrfU~EF*#e{_iyFL`GS1GF)%AHHZpIjpf^oIYK
z8AIjKj98IZ0v6BGO+fp8LrPc$BXCgrAU>#^v}Sd`#b25nIM1`B&6m&bcR_mWP<X$n
zBNT6n<(OxHLCQ94(^5`adw^JZ(-nZWk0{EYU~WwP-Rlob)@3ECE{l6hL*$GVMua_f
zTuE__4OW^w5Ox(?P^)Git~~yrc?DR^CeCKK69mWQS-kvcv?d}%c2E$E_<+94E9r<^
zQ6??<7?pEFijK=?khc39jpNzTO~^u{<;ml96?ValXvg#meEeSkzmtNDs&tv8X~4g@
zAT%_zKU|Q?AIU%IH$4Pd{*8ea*+`@UI%rGaQfh{^np#EmVPEv8IcDMfN_W?_Dmi1(
z#T2~fOGD@BauNMXGhXG?L<GN9*169_R<z+!uTh^CPA%eYV>g>6--z*IK#5&<KJVus
z8heHe5TzQRg5xs66nSJkru=Es@W&CeSuCR&-;ZlI)<iwMx#2RnM*0eNrs(ED4|c>I
zIj3km!Ythl*hifc(pJH%F@4P(Bo$>{qT^A{9BHhy$G|*UTEn~9M(Qx)m5ghS2eN}8
zsCrSz2=hH;jTc)Swvu!3=j84Wk9=%CtywE+c8}_9(^&l8^xlcG*H$&oxphM4WWdzt
z*vb6lX=ZECp`7O}GS?M&P|)g%!fG`}l0g9Eua9O$vC1BH24IrQNT8>2rht*Gdt}4x
z{NHoKSk!qz&nV%6uB~U0t6;^G8MbW|uu@*T-Sgt1M+uWbf;t535O&Mko)LV}=9C=Y
zux`Gnmsun)ISY>Uy%gP-Fh1<%-ymqyQMiqFhx*sL;hN0+%6w}9Pqc*@GL!3HavjV7
zip7=2>t0{=G}A3XjAij1)bKa;(_~^Q&Ma>elZHzsqP)sDM|9(UOm);YDjnrjiKlnx
z7=lqYU#csiqW$)ms1XDJKoJQ5_%99+B1%L9Y3IR!JP@KG{`P~v+AKu)f7Q!=%ia+a
z{=1IFhxz<J_(O;rF+b|x{sY+}(*9eI2la1$5#mbp`fr<D|L9fk5vQVq|F^?E|Lu7e
zPVZf<-n(*ofE*8vlpyPMcx~Hy6+4Op-YfSs8grz2nudmgL(qc~O5w%^YJ6ye5<aGU
z>n|pki>GrgXGD4$O?-wu;GBjo>TuQg;=$vXtt(=&@BJ{V3KSV8qg;hLr|zcTm8upX
z5DBK*c-aOnkJ@yidOArCGCNpT5dY;FQLTBZpxb4F4y`uVj28X=F---Hb<I5CE8z>s
zM0Om7O*9-jEV6oB%q-KWiW2#8tIjoCZ&Tf@Zy{c0a@A>T2(<e-_fPwIPF=6mD{T)Q
zdPRA_C7+yCi&Zacc(oeNSqb}jIE*6&tmI}_6eOFr>Y3Fi7u?ylGPupRYUE$Oid3CK
zl=hX(pgVmHAFDR@v`%GUYnO1cb1Vn{aL8WO8G(I;^ZAfj#Qdww8w2UUJO{kXj)PHL
zN(!5dMsikN?l;;cIUMP^F4Mg27KfxLw=?IHjAgb$7sgTbh)HVeJ{9Ri?uJ>oU|MCX
zXH(|eTWL-29i~1;dV9wCDy-vA5;ue67UjaZaG7>R$3BmaeU|EcPJCvL$3q=Kg`o8<
z1o<Q0Me*bA?QWo0H_JSCe}8I{0GnbWD2`=9fb`Atjs=1m=1A1gXX}glLW#&yTx&71
z&Z@R7jLZ~m=pBa2{^KRekVnEkSupu@7)s|mUGDw~>P(4O8pJ&B{21In6?FUWf{+>2
zF8iegw7j_lm`!)c*70_cB;Bf=$K?-7X;N(@fH*{-J(4WQs%a8GR4T|xdDr2)gba9z
zy|>$O)Z|f38)7h;{}zJFz(E~n)|uclk#Kb6$!CnE)lbZV-UW&_L?7hB_y+#@0%!S2
z%EAI;;ml0y+cKQhP5HL6$f5X+u5AFWq;Fo(_kP9C4^R&u!$Y+IcAe3oBh*`#kl~na
zHeYd4OY~W3GTpan`>3~S%uK!6^*5}!>~7k;ehgit0ED8L1VSsym^YlhvzWwTeQ)k+
zR$=}FAL2k!rwmcvp<=n<5gf3^^ctzDW+kU$#yd*_2ImeFrdBzHS2B+UehZ&GO?vvS
z3D>E{45aF7R>FtYG`&J!FiYWmzX+WoUXNTkryFTY&pR-wyr8&Glm~Y|pd78XrbQc=
zyeitYeTQGSeuptbSB3d><UzO*;-L@nTU28fj4=3o5FzU-$Ct~T-RlRI8o{&TpaiS4
z(&$Dag`v3Q)e~#D$oG6n=HaBb25ej-@vlc&R^Yfi8?4L0Vd0frH$hX{#86yAdv(<B
zsnNIUmG0bmaZ6su7lPHM(NC>OtC&21**}hn?;DozOjhcZ!n(3LXzKEf7}13lMMc`2
zCI}KgXITRf3QG2rs?h{Xi<x0mwLos3Be!oX#`k#`*M5;ZbY{-G8bI*PhVw+JeeP|d
zVV6<(!Spi-&#;gP*?vAq{Kak3h?tAocb))KtluFIH=_W%w?zr`d=XG`TH06Fg(L1R
z+q{Yj3=t9H_;tV~kR4cU@wQ;~z^C)pHN$O+Bpapr9(n7dDZerJo-aer<bx=xjjt>1
z$oh*8LGfPY_Ygu^&p5EWVnmlfw^;h!=<(;EK+&NmIh`+MLU%W;I$zw$Gp<p}_peZo
zw~;NSL=ohrz&!N6PwHeI5AsV@u|f`q?a8(wXG1mDx}dG_-I)d}_fMsQH|#T=MB6J-
zx^c>4=Ius(n^8kxC6pfw3QgrNjzaIwu4V-CxSHX+suYwWLMtDBQiMCXCr}H-G>EMj
zyLoc2l*ni3q3`m)4r1nQs2uM<c(#H>?w%mM&_~|pZex}atOA~Du|R(UO5wk7`{@>I
zi`6(Bw}Py?cakvzK}yjCe&2~h-Po%n)11A$4E(iF{zUW{T$@<qxc!Pflin7C&itd>
zw?xfw7v*MdcWiSpaa+H)1zHJwH^a{{o2md+WZNaqxV0io#HU;<u)|;R>n`6X2|jRh
zb$tCkhj&aE8nlKIbT9S2iA*xYeu+WyiU^NvnL%=73ip-HJtv+^c}SdF%1u;R=xFBF
z?&%b6)iIx%H&{=4NA#Re_)QbO^vJ#<ymYl`dfJQYaNnN=eR}O{TPE0NFw*mroJrE9
z4rn9`%ekjwm*f@*UW5~kd7ljW=I=`AZ=Ylfa&8zn<Sk~uFip@LQjc+WJG%#uBYh%+
zJEg}XmClvBt>U;APo?%vk8>0p8rAWxMI^Cg*}NWg0Q(hkTeIZF)EIrk>Ya;v=;2YP
zRn|0`0KPD7pSov505(CDU<kf`my2&(M>@4JGWiB*oeyZVH-FVvj^Y#e7556;Lx)X_
z@9HnyMQ?>0ridAx>t0A!|Gr;#B~@ZYYf!!+dnOFw!komX*b2!D&SH#e&dj`<1ugtQ
zvtAehbE4uJxmi+H^FmLqd_xRcLk#3Bsv32Ud!P+#MHUm9(UnsX3dU0+k%rBl$B89m
z>o40?O;u9P0qsRIeoGD2W>fs4l;@*u-`n^2z1zV_1`$=ccd4%#7+5PWfb!87@!JM4
z;WNkM=pIjeQP?;BN#?oYRUyB-EdU06WEmPb@FN1|rPX<-KEF;kmAmyEupiZW&V%Tk
z;7XEbfZp7<lL2+PMU^f$91i(&R@KHAf}695EmiaH1oys%D<<tDt?ba2Uktywuu&WR
z)P(o;SmB=X-Xp@We=q5@cAv-@b=I(4JFcb8YqL6VLT6Ywc%m6FfNk%qZg5z9r=vFx
zmLtc1$ggyTkGGn8D!4Oo1~0t*6h@&m$^FPxl|Q~Lx`?JE2Q*&SG5=9>*jv$FEgVB5
z{C2jwPR`=^=SY80m(0s>9Oe`wo!4?D=<~hskP@2cH?3qC0FF$NJDA`nBE%lHm^6)L
zkS{4~KVZ)ZYZT>s@zu&lK$SMCc0_Ey9D)n>GfPiZa%ztYmF`o*II0k}%VS%qD1&W>
zs2;VEnvuq|>2opMH;>qVz|#~N$C}q~Oxz{@4|w_q^kB14f)@>poz?^heV0u%zUtfQ
zeoGc#m_tz^MOI{HRTss#_0&(DMXry@tj$E(nklE;Q4Yl)ptl4XP%*0<e2h6HUtWKn
z-OldgauL<g5%rE(kx6zHO(Rx9*(L>frZScb!=$*~*BS?BR*L9*%b1)t)pPoMd)gRn
zcUIPK3n{^IlVA+QXOS4lqWDC+I&VMkDAf<_PbYnGbhDD+_HAMx%Ex5wH7Ar|8s9}<
zuQva|H*D^6op09>udXiM7jRAIv$zc`?cYJnI*hEce`pf3I7qi$Gy_(=^_b%xiTL(p
zV(|PW#JACP!6Nq7f8q^;*bq_9?D7lQdSBkXq~gW-VJGnXOyni2^g!^Q(wmkbiOX{@
z7hP~nrF2_NeBPy1>1kCCGutSM%yMo)Md@;^uk70$^a_Dz*v~PzUT?lqf6bE=ymwkK
zc|O`(FCZ4YD2wC5=24_AIv6F|r#@2OX_Ck}G_*5=Ho9dvWs;F)eDEB65V%636-JEq
zUh}K%y=E`iSG(qPQ#3#W9zf0-yupf%p+iZO?@BPfK~mS_CB(l{5QT?i`}GbGCa!-1
z8rj(*eEF=EA%qYKdy-w@Bq5Y^+}kun=<_Hy(QW+BH9w45E@EMq?hHCS)k`sB-aBEj
zho?xP#F$2f-V4l$qc)32WB(ZSj<RvDn+Ub~tJefr{%p^Z+79=?0)u{!rz&c@w3Jg`
z=^V6>thgN#|4G29@r-2LOJd%%Xkv*`6nAG*UnP8HP_nk+dVlRzW^&d%EyIltye9QJ
zXyslr;#kTpm0{{pY}f*!xnu;*n~mrG9UmI{+A(9s#f?{an}(+cuhs~uUaacCwslz8
z`~)r-zJiMUK3srM2y=^@h{L&868OU{pp^-g!xC|=hdv{ik+E<TA@eJOAteE}8ZSKw
z?YNXrNA96R+=T6({#+H{R5#s4|LM4FVMvrq*w|fBtK8T$6M<7b_e%4&JaN7~kV$DY
z<UQMUy^;t2MrQG(l*8zY$hQ4r2rp{FDTN@o3@#yj#f8p}y3+83)90K$Dqui+6_o~)
zsNGa}Lts^bLMDRhAY?Ywo^ObMAevTh<l|eX+SDpq=Hx^j_p_b7x(WA8;QB$?w-4i4
zl`Ha^o_Ihv(=b`GeNIdaviq27A(hWfA$HwaZ0;Atm<R-FMbW#5nyXcd!k~vrMotWH
z(TLIKZh}tOE}S{C!~aa9WyHr=3r@r%ryy9E9}&IZOs7hgni6w?5q8^YLkfcjP>IV!
z6L=FBhUM<~rcMgb<`2{h`ZjV?hi?d_!glJv{CFqnz!pgM-uYx7S&|Rb5u7H5FhvQg
zQT%Li0<*4@U0ip4E5)nqCc*-g4D14DK2*PNI<`>;R6Y~F`3O|1r=1b^3I%`Pw$pmR
zT7%=W%6idUS~i@+lS+TWo-+gcvf!pJOoa%?EU|QZ4l$X%J45uW@!MMAQR5(XEib4D
zdfB?M-p38iuoyu;uh-?Yu0t!P0cV*^sN1<F#|`%DA~*aHU2Mb;p8=MMj*Q?2si?ZZ
z&uD@!yXj5vdBiLj*Rx`apM=hvbUqcA@boJ239V0-eG=R+L)x-uVB~Ud#gLE&hwVl+
zlOxRVoQg0@+V^!EYq65qmAOfEx^gWxCrHsRxp_6-d(NA@CEL{8by`JG#2fA5k$&~F
zrdcfc+Un_UG<TI#IObs`srejS1XgWHy&P=2%hy<bAe+tX?^`>nx%Utg#bI(U1*!BW
zPgN=|_fJ%m+^RnX?}2vGJ;rLH@4wa-PR`>;jba5`U}sF{^1f<cN;&+n7n|(1tbHv@
zqYdhu#q`dbua-j{yyxNB^ShpNz+dwP65pNIG=AxZIjH&=dRSds3)O*BAg|DwRPJAT
z%m86z+XG?qddWNRJlm;$Bh|D#G=3ZF?Mm>vv3u{vdr4Sv=ENo6DJs(4hVxp>hn~Yh
z)1M3?{Y}n?26-m@Np0i$Z?gQxKDz^hdG7fdsARO<7Ee`P#3*^pyOyrrya{kelv^dK
zSNnv{&&6JEH~FKGoO&7j0C&Gj)^X7H?LKJ2T4?1(p{uW~=&ze|i}Lp+nj*{7DCH#?
zewh*CDa<V!?Y|29wP%BtHEw&!8)s<1q`VO?7Luiia;u$UzO&d{>DuugkYuA-Hsl`n
zGLG_*WzLgc_d|q}pS_my9gps%$)D}W@N9J03<UrfVfeEh{S}b^y&qKvap@ua6>%?d
z^J4$^(H_W49s%?}U4Kyl5K|r^ayxtg0A%HA_9x_4S3-D32>AEWAV?e!`7;n9WR;Kt
z>;DK4@*rnCO#hmPm;BirAw-#&>VK&%{TAR6K|24BApW;Nfe2E@%R={`)YG%SxhJN-
z%YPLi5L{mVe^0bX1ZgFrV*Ri8`7hpz;=cwH5&xc^0Fv<P#a|Zwe32S4M1_dL>|YlC
znVu#eQOyz^0C*t?0AT!=g<xU;;G?Ujqt(a%T%v!*{U<E=Lty+*qJenvk%|0S?Em~%
nnE#{jy9Q_eZQIJ#)$YB`|Ia#PoR1vA33ABCjMC5jXW;(>5+@;f
delta 9616
zcmZvC1yCH%wl^-pVQ~l&5<IxOJHg!oL3Z)ry0`^j+&wr1cL?sfxVyVU2w(Dl_rALC
z*4HyNJ>AE;tGni$^P7f%(m@z13h)R37#I{7m~9v4I1FXXKb4ZAj9;z*E)0wSZX8ey
z4ELm>s3C=lf{coRhmVR3Kqn!^#3v=dA|fP40x+VHbK#S|$DrmRreYv?&rOAnz(|Nn
zM@Gm-PC`lXo|c-7mzI=|m4g1mduj$ACN^F+IvN2ECJs&>%8#Pdyz<P#n(TaHuokd5
z4mh+vv>zq4gv8~zWc7HIKJ%%XNima2gV{+GxM*aBxTHk|w8TGZi?XUq@@mQQ>qv7L
z%5xhi2$?AITBr$0%E~K<%gW2k$ttTEDk~UhDQjzMi)#TT^&PeJOmy_j6wO=>&4AiK
zPfIOXGktYiW0kLl@|I=>?w^%CEHun5EKPw<pKaZNwl1!=7VeH94^K}{^8$5iu%Es|
zpqPEBoJ)n8*N@N6VU|7+mH-y+AQ9;pY2`#^*D`OPA3pkl^>%?NUP0fS!!xWRTkYdI
zyrT1b623bpb_LoR1Uu?SxS0od_{6%Ij>5o9;Nnga63%@1FvrU~A}KfWNn=q)W<^!C
z?~C1rwe^m-ceI~#l0PUV$SEhn9b6da9UdH*8XX*xQWly~6W$aaT%8tE(3G54mRQuB
z+E!VT99o!}Sd$jqkQdcboLuz%dwzL+(T|pfA0>X-{eijtkc#oN%AXn4!)Z0sMNI>R
zt+RC<gUMZ~nO)^Y!$lnn9pycX&D|rl{mUKWo2}Eks}a$?WvK(znbT$2qu~6?;rfE!
zw))wU;_0I5`KIE-l$4XAqT{-{%Z`r4_Wq&n_O6-Mg`vjH(T3Cc&Vhlx{^8}`<-Y#e
zdFV*@O7H0M#LwmN`K9S0=<ICQ_}0k$e*fC#^vcob=IH$X{OINlbn|3wV{>uqZ0X>9
zZs%g|;AU;_V&(W|@$%PUA9!G6bmm~be}55rJvewbKfkrKwY{-+vI#xgT|PKkz1<nR
zIb7QMQ??Jb&$f3CPY(9Bj}CWFj*mC?FSn1bHcuW7&u)&cj`lC^c5i>3T-=`AJe@qg
zJRI!bU+mmo9zR@cJ>DL_KA&CPUS8ck{JCx)pMPCl{d&6n{rfi<1_tK!^;KitkPZe0
z^Ou~Yn1=h(X{MJqU>1Mio%~MWgSh^VL~KBwrif;)PHeGexin*Mz{1>|%rDC6fEfA0
zhb=SxZ(;RBZ>L=2m4e?04{87e#AIDFn5br@B5NODR-w+TT;b_{0?)j!H{&N3i>u9Y
zXQShznVA;ct*M#dW&tnowb1^GljYecfxp@JJTYM4G*aP<vF!8C6-&pf%ah?2^4MY-
zi!aM#=jM~nPq|hr!+S^alE*UmNuoevLz(#x$}5nxR#)>}Ij?QSgO6u^nNIt2U8P&c
z%Wti{Zc->3UrNz1S79#)_U9;p6vEv19)5FihKL_dvzB4Os@6h}6JX-gCC#=U`vq=m
zybXKl+jnBe#S1QlTuvb;kM5%sjlN#wT^^)pic5&cInN8xgKKwpC5)}#I^{EZDM)^6
zs_AV7xoup;3dM|q?2Uzpb$7bqAl!Y0DgyRPqv~33KmKNS1Xb0Mzl#R6vpg>e(SAzj
zJ`IKG>?rI47i8)qu1W|HtOsqO#?QqOjZrY1Jy;BSUGxuHMJ(adr41H&ea9%gC2XkD
zrC0ZQsA(_jXF7Nr;T<Jn@XPxqBN(GrrVC0+LP&TDWL{g~B@G@>$m;bm%A>0&YuP&U
z)e`aHCZj(-N>zW_shob&EbNYYRTeZ_YP*)v?r@O+pUm(TPK+-i3yf46rTU#soo@+?
zoPSJ~luFYuWUFdx69LL}ypOd7&Eo3>(y_2)Yfpu?P(xgDOY%C^^gk%jkXn<kPbYKg
zk|!C)0bNs|!9}f4%#S!zcb-cZBVL{G1C@hc4B_QRSj5e}qard<rvpg$EPc%$ZQCaI
zB~lT<t5|Mxgo*2IKfMD14Ey@x80^ktRe*BTcWlBHmECZG8i(vW4+BB_d-_|}yoP%`
zi#7aCg}1pQ@sbJ-htDhHBmv@y;#6-*7oW;)W1?=DNbpDYuZ_L!^nrvTuHHt|Z;eDd
z2=ApyL>g7}p2bEAY)M4&o$E%#jRo8Vb!cqCBh{2lEq!W0E!AZ~9Fx8PF<`KHY&xwj
z!|6r}n-uM*uLYcbVF(^NU$&)4FhBRbVIW=fl$mT$;=To{EalYg?BKUiqEMoaVPnHl
zUZ@eWi29vNDwwg=^`{#|ym1e>jg;9S!kTdK-?1T`u$8G*IxNhPDYpD2*i$#_WKiP=
zrVVSDq?Wy!M+hRgVL>@y<PXmADRz7k-aEIJ>M{TZgw$kpv&*h|FDTKpCf0~?A3hUF
z)aWb460ok%=lN8`o@`s=1a>Ze0W$9`v#G;$^b@UE)GZ89^5MCWngfRpQ7b4=R6#ZM
z&eA&{1U!!?G7qc=DTOQ(t%~+?C8OrS8Uv?z^76bgPfIKc^70=f(cVu@JBZSFGOB6i
zEzTA{>d`jyWxZ#dIRp8)flOW51_=S?{07YOf5)&li!J5L(3*B98DHB=IV-*bcesF(
z_4>l{A>kL3^&%R2W9E?OZ_SsVBDeKojsDXG5pC2vcpJ^P;jq*Q`9CS63iiMZ=fnAi
z*AJTHq|`uP#+E#{<kqX4{s+_`R_B0D&QC&H!SmP{vOx|cxXkb+&BG%jGlvep=*953
zI+4|7LdfEpp>x#0*^<YQb4<gka1sfEd18JjLu}pJrVZnVtWG05u4qLTnZEh>z`JSs
zu1CFMsdN-#Hj<XSepN5!xFm3N+`9gEe-fv<=+6`(akBICdzLY8pyFfMn=DMxILk3&
zUyIc7)beox`HLNitEi>Zh`ujXY8Dl|H`|%BRMx*x^ova!0M|u#up*(9F+Pw=Gkawm
zsXI+U#ig0!7e84ziuV34XZga-W$(FSuVT&JZO*O23nZo?e$%0qnzRR=PLU`w3o{sK
z3M=7St};vYwdt+@EU_C>pLoz$7)JupmBDF@crW_p>s|~Fb!<O&Mp<Ejd3=Fb1r|{7
z$M?<RmnaDzYYWIqkJj%=v?y1M2`f>M0cUNrDUo^*IZUfki4@1^;f%0GH4-!gj_V^^
zPQ++Ie3O7=ww}dR*ij0Y{7F+cDdY^>2V`(}MtHkm)|&*a$r01YE(^AiiyA6<ryR_J
z&W{oQq29)-i>A7g+a%py2E_Ng9Ed?wPF4RSC2d_IFMwGR)U3QUHB-7+V~DnUc6REi
z$Mr4x(K^3>H_H@wL3{)Pl(3nEmEEqE<Vsh_B}%oGj?CH|SRa8&D&}@UKH2;QQcC!b
zR2{~j0&S;sRhgmCLmbkppssEIuz{F7OUWz@4<qSkQ6mz~wZ}xEyaZJO2J?pHNIk$i
z=B=23bu3jdintXw&BR<T7{`PD>05zT4!DF}wb^~p0-Ba5`D4;;A!kWPHj2uLajtkt
zM}~Qa>bR|1#Q_ZEW%pu|1fl1?`LR*_`2^>t@TDO?$tkSzlPE%S*lIDGi5MVPoEeHh
zAc?~O<*lZUP|5B0Q!n|mjbJg{QKns?g1+6X*~WeaVq?@^1!7TD7(~<nBChf>zCe=k
zS;aUDA+8V&DPLk53=mmKogjR{Ei{h|*o&RpRlLd12{ySnq?L!JQESQ7W{@}n*n*cr
zZ!s1kUv4w>=R2}mNNQW9;kitKj(FY`UQC?_t_F+ivppQC=Mzd%M={1taesw|IQaV2
z^u5uTTJR{RKMwJPM>LLTpV!M!q6ke7<S*qK`hL5q>q!7L1N7RDZ&EIyS4l-hu3zdX
z6}-SZk%wm18)-F%X0&~{+EJO@UIYf!zn9EDzTPZUZO%GYTaVt1e5wcu9Z;q)u3CYv
zIvAV%bO6~$!_ic@JP0p%Tb#6vYz&5j*S{7a5RM02^PCaEhJFF!nW+<Y4@i`*3%)SA
zd)C0E#gcFbqCbzu78NxB<C$#7B>e(=;A_Ez!y;6a>oTvZG`26%I6uy4xFpeucaApA
z%{fA(Om3OfKye?Y7g$I>*^ccP0kZ*M`q!qx1d#4hDQ7(<4YMOjysPr;QFxR+7iFNL
z+`F7T)VNYse-iOLqAxJqG2aZ!RF;DX3<C;m7ucBQgl>d%*f2lpp;XXnoRSdxg;;`N
zwSEO9q`8cdge>K5Ve&<m6NxyTTsvekRq*?+xmg%6s|yG|P2X4XZBLU9{5<|{z<*uG
zpc5&dJ1BB1rv{qivTkxcEfo><!ru7Sj`RMY2Tb#DlE0C!-DRn2wT^b7st?zGM;7*2
zs8JI54n9PEMBr|Fsb0oznA!P2-p>pC^q3gOj2B*ubF61$gV`*7i^d~h&vd+BG@QLc
z<-ZBNYLFS>0RJKWd5)KUj1nqNr5OW(dXcHnVn}|7MMUTH2blSc(EuC5Vhdi#ELBEC
zh<{X65mSy~;!OX9N~-6D+4_vZ0JYS6G7Zf##hAr;$)uSCK`FeZ<GK8Az?LAQGF-4g
zMhTljnqdiF&perg68Q`Hw)N39add8xD>nxcn`*Abf{QhGLBH<%F(k9P3<)K!l>xqv
zq8izK74<%x{fD&SMKVd|#wpj>s$IQ;z~WtbCq+!uQnchuh(u|3r`eM==$kX3{^a7E
z#&`tvH99ZJHZJ<A{^Z9jcswsD@g?2SpHK+2xVfI;s91OAKM5e|dU$j3i{ztpnqvjF
z5J;D8I%6{ydWH*G1POduLBULyBj*OWeHpe|%8+J;PnVE>-bmvPjSTHIr@<uF4ehb$
z#^5_p8|mujADwsTeA)znWPYSVB-Pr_*(^c0a<?4=s^g;}QD;Kne)UERcFRHDRj<hJ
z)u}VpsvX%P$-CN@?qot!UXQxpM~l`h%H;^dzN0O~>9KJ4JseRoh5QW7%EK!0!SEJ-
zWY#-Uw@54$E{r|+Hc~6^&!4WQ-~s~x3J2zkoUV_{y22>(T}ZqTeo8pwjwS(yeX&(g
z3a4SJj1HRaZYV4QXV|LeffqIoK4qsxGH2UGSs8=m;k!(abRIg2IANHIqJyZ|?aV-5
z|MsQWpV=dWu&;LFo@*P5kRolo%Y}h?6(zZJL#9V`B-*5&(5~y6GmQ~TC(s3!0l6;n
z1If{pLkfRQbk5%L?!4UCSx3kh>^0=h4QiTr0{wiw25=O<&ti-N>=Q}pW?89%+R@S0
z<_U6HuDTh3OB+8e3t|j%Q44%~>q1IX(R9eSvvh;Cs&$l4aMw3$HlzzTz`i=a0t6cQ
zK);YTix~BkVye_|dfOVwFOE&lEL>~Rp4dU9F%DH$3Dwb$PkM++spgAL9Y=^cEF;s@
zU5Dju!BBv>Vl*}!P<W<f_AjZrvbThB^RTvFHPE&l*`GbAn)`w)*7H=gxoi$VXv42h
z%H>{hPpJ}8w{k{@4{~s`_f9Cx{8peeCtFc@ZzcW#dk(1ssdDXMLBY008gvAY#1zhn
z7HuU*@1>|db7Y+t7MBFkQ@yM&kkUFDrvf`-)IISgCrOTUk<~CXw$I5jSIYZ_7w9@l
zGRD@)tU_T(k81I7`8$vt(ud!2sha(u<*XU-j`$OfynD9RWblGYDkY$r(|=Y6MfWQ+
zBLQjZyT)R%xFiqi%@vY%<;{$}9frONDK+G_pXlB01?!1@Z7uZEKLC~3#(`FUih?Wj
zaaWA$>S}9~l7p5kI7v~OJX{~A-zy{uV#^Z5if9o&%r}@{b*~R2uBg>M<A+yaxVEga
zRjE;8m%Ge+Oc4Yv!7OnzP_nENCZIuEq7sU-+Qyhm=!IqkJe-o=Hi5$hg<@5A5|fP#
zN|}C*DU=+Zu2j@<zFx1MoVt5i(tw=`fmd8o{WTfAYH<yM>z*n}NWzhv?pBn*xeWeM
zd)&{YYL3QsRvwc@(``#6HbGV;5xYJw3TT9WV#6+N#>^PKzYF{1T+2|W;q>&lbl2Jy
z*nBjQAJTq(Ww98O7~!JH(o)5%dn7q%a-}nZ=_tD9?_qr(7)h)uV=F`@@B~g9uLb%P
z4P@@gZ$}m9!xBgcOtGh=MWIhyR+C@v$zP|PTL;Hmr6!`P7A1YcmFSr_4S-=5he0rk
z53)eXS+o#4+Csc|V<gJ>{sx@^QTU77b@_t3-;3Y<k5<IEUANlj;Fl9*xRJNmx(|N0
zFK!~F_MfE|#-8}^Yd0r!>u%(1>!$X0+)yU_8vjJ@2u5Lufp@6yER_G|pcwL`3EXtH
zDru0SSJ*8ItM^~fjn4k$%P?|M%954hCISBgFQK5I{DqewuDCSdR&Cahay%xJ!ndAF
zrHN_s?NYs2=BoipwzUOdV^;roySm}H)O(AyQ9+1MR@#ymqHR!ob4<_Ry(oKMaQWBQ
z=ilnp*A(BdQaybM(CM_xtwon|VD&xWHA(mjzjn^VmrkQdm@oR}CaKN_WY9ZD6;)0F
z_-4g_55Evv>qn@7*Y7`P2gz%rd<RAtxQ2Xf*jQ4Y$Lu2xAnC_M+QiN$7GP-8goYU^
zfVfhq^z|AOmO<rw$9nQ5Qa7^c1RDhds&!&B_mSLD#rCJXrtn2+t_L3FEUtR;Q!9EF
zy&bBC!0}^i8DdUas|8Qs$cCmUc;9$o;{`ZlzEq2$)&C9rV(Pl+9!>2e^O+>0!Ndwh
zifq#myAW3SE@hG?*O3f2((#ul%FrZ!zXoDcbzm<18Zeo$<IM5OBd{}*&;BZz^DcDO
zB;pi&;N;$|MfsVM(NHivKZL~&_Xc|_Ht6yN61#eMVJ3t5WxE0<RI4nzi>X@`ZxwtZ
zITrY^vS)1C7s&w%uE;HRUhD{4Cr?3V{>kQJE|kt3B`X{6tl$8!C;waz*`-gPs_lQX
z*WAP(9e`x)GQH<Z_E<bP+4p&}K3nXPP(=Y5g(NbXc?=c?MjinM=Kr9GkZJ%GBprYY
z|98o&0UZCTsUSWCnE3xf8kPQTl7AO6NG?7eq@4p3?LV+c$Xk3W2pR!4^1n9^1YhGH
zO@I55@d5w5YvFSv|4ZW^ezt$}*2WX0qW;71k&+y&uK?L@!0Fi6u0D_-_S$%+Qd=X|
zR@2uP9&O4f&kSYx7#2i#5|4nP#_}3YcKv##FbzXNEv1?#oi-GSOK$PU2TI?EO1r^B
z+g@EEz0+@MGbV48xit)v97rZ(3#KxE&BRdPfBJp*+wJ_&TNzFe16Zx-phpL*g!aJf
z@0<qyG%_XNF*Hpt#W(bwUfw~NIe>o5K~F(xc}1N~TW;u^t)QL(Yhxi=8-4#9r?qIN
zw5V*|KW&MSIkR9s7QZEe?3o|k>4I8&0F_Vd8nHr*+-!24D?{J%28+@fjH+hq3)KjK
zA?MF&M7HRdBpo(Gra$qKNTS)YndZ@FCqdvf1VXEFpNvLZwH8H-l&BCS{Nlz^;vT+g
zzg9uk0eaAg^OiM!RCJ{#CyZ0IsEu<5>YiP(&gO&4kM0JWfQY`XC(osK=n)h#SglQ_
zLkrEk_Q<7vaW-}7@7?sX)Rt0)5r@Y<j&ixI&JFJq+&|Ghr^Uh55_mZz2#GpQQYwQX
z5;?Hhxv-%D+JZajO1=Jr$OwuRWOgan4l1|n6coaUXc;-=9DHNtv79U_EfOC-$qIVn
z`?*nCsme=-!tI~hlJ^IG#&(o3<IOL(BCrzfK<Q$#2%O53-P#27%UI8a1#gOJu?|dA
zW{R<^Y9Um0#tY0Ix|=+C`O~Yt(BpuALMstv5YY{{1(o>yRaWso5^&YKC?^ATDkV@U
z*S@6S2nU%`ai3aQRcR?^=hT!AFk@s=OF0-=R`wSLg6K5k$>4~&nZ}Bk5ADXxhAcAp
zJLG`8>+?;?)d7Ny5fYad+s2&2MHlf{QWa?>vqvdcKCANy9!&Ej-DFRs=Fwod@*v)@
z>*srOAJCxM<LGGzgaQR109$o~S{VC2Ck$F%&DH33MWG^nZGZP5pU32);pF?in;r(X
zBzg*=kUTmP01F7jH<gsk=3e7e;iE7owe&{a7T%aRs2bT$sdtx7NwB|$qa|GIr9<s=
zo|XX&Z~(lC8!9HAskwc>5P`A$=14r1(cT^e0H%@K#OWhMEqKfCh>Q`kJ2!786cT<1
z%Er)`XgpzJk-GccgSbHwLj`&H=3*Xw7=YaGL+o3Y1Np2RX+sRxW!FLRxy6XG#+;@?
zQ(?I_VgYFxRQ)l|TVnLV?ku>}1jF@f@{&W8IA&Vd844>fgXe63ZE!-73NoU<rT%1S
z0;2g2&#-!PfGMt6|4Z8_8vFVgM`cro$v_e|-9<e=W*jbCixxqOB%17v6JVc}|NBf(
z0miJB;W}1)NXG}N^Hk<pZoX6&X|x<OtJH}wt$+o>QBkc7n%{$6EX^Wnq*Y}E_r*e(
zwooy?9Ey1U;7j6a9q?=DE;iF0M|V77r0>t$v#rOn0?{%GbzTE-;kOvUYCOK&&zavQ
z4uwdzJ=mIV-!z<@kGxs`YGbc`^g>>e-AI~A*u^9HC~gM<C^tjE5ojTMG<}^ZPb?IE
zK|NIEuL`=Kw8>Xn<>#ETk@o||j$3-fo?HCNhb&?EJlH4$H^2*L;A>IOD%<_eL>+F8
zfx2e}7084R+$#pt&$1Q7d3`tu5irKV2?u5b7d?Tqu7m8yp-BnyIaX*gk7hrfGE*}D
zYIr^bqr`7}G43HZ6rMLaNMys%L!D>KWR47z^pKk;h>oYJfvxLT!zntc*M4C`)C8~n
zx#$&*O^E((7Wn<v>}SO*hO?G?s*n|nu*)|Cn(pICj*)03WIN7-@LwNwFzJ68{`rob
zgI{W4%m08Ut+AVjG}0VIK>0HLq69S!lcnHhr{hah@#`7x)TyR^`&6^Z#E{ro1j}HF
znMpKLWnXw6+x>gEWV5U+|3D8|OppF^bqjZ8;!3e5A#hN68P%8KJtgWd#yC{xA8@`_
zhQyXE@r1{EaRrs9_ga1x^0wpr<fSb!*fAw|YK)V+JUa8M*zJ8-Mq^^Jevf8V1x5!|
z8IlF(3OCklS@rXMUg|ze?{GuM$=GA<lmXj`a^zN~o=fq_<_XnW0sSp;LN=MN;g&&V
z-P}^3IVSi6{Ue-tQ<Zv~jWwo!Cj1fBt!!k&uQJCAH09FW%>K3Zik#au8^^5>!p}2}
z$NLpc8W--v#nj9*3=)nl8|$=h5>aLl2yGY<9th2TONvE(Jd6i+@pb;v)o&l`(m<jG
zCP<If%j7fQ8B2O4G+ju%UsHW|>U%srF@0}(>D2_*B<QCRLdoU(bE9weFNF5fRgp&_
z2|(;4!+~66d(rvQRZB{%a-_6PeYBB{S)jIDWx+;U@Am4bOG8Y`rqV;=)xz^dK1<eT
zDv47ui3*7Sn^07*8j9i;5FbP<N1W43i&4cz^Pc(hsofg!DztwcOx`2Sl0@aJ{XVGN
z6>K?Bdg|?I^h(ylznfi2k~DKQp{8`)1m=t~3J)}CbsH8rxU;RQ^{qT4znpN#GZF*V
zjd8TeuLo%tmAY969KS{Do3GX=vd+F5RuM=XSB+uTxO(F~O-t(&@O#a*<x(3zr^n7m
zc%VXkEZpgIufSuKXCWgYwyk|{XI2cg5Bys?*O|qyCU6NmA?fH`?|^5OJv-p*)wGJZ
zgipSTC!y$K+pDd_BGnr@@}7t6!NlJR5{II0pWEc<ni~S=-#j4_T`>m13wUikwpK~9
zC+zLMMzQJTz$qBR590Z3wI&@ZdoP1&QL}U+Q%kFkA$#z{@3m)@ten?|d=zcx@4<!P
zRml=mIfd5)U7SMNY@%&D#YDHiJ(BXjzR6MeXoLh)GYt0&;J3kgAYgQTi6}a*>c_MW
zghz<+2l95nO*e+cATt;fWM7?I$wcb$L}k#NVb_mi<a1xdOx(Z1{tbm@MsoZ<{rQE;
zC-m=7=wCt0-!KQVOGgae)ibo;;syAu8>{B(0;L({bIcoZl2jj1a^QMo{JKYmX^>X<
zUFTfij8+*7h&C3Yn&Z}x06qNT7xu}Pc{&CfSA6va9=3?tqv@Iz%vtK8h%^NCIJzr^
zAw_W=az2bNt_u53Uw-B(F;DnQ7VrVYxi<-G1UCClPg@fPWaGea(e80D$X9vQ5zyc-
z>lF_zE>pD2M!IH_S!ZerR4b+*YV|5UYHKGaFtaa;=f+;1pQX2;e?>(MVtL@~p7$}D
z(MuSe+R{j{J)VD|FJ#Hn3_b1g-sDqi9k=N`*LOnRRBzT;u9G6y24!vXI3c2ZRSs6Q
zU61lCJ8HV?{&ojmHEn+rXHZ%f6Y4{hVCM#@eUW_!OOF4JuysU!+fZJOAREV^UGft{
z@%pmxyWUxZ19Uh@OUaAU564|4Ny*ZL@jbgrJC1XI-#2ZD`_$M%Zuz@r`7v?KNn#g^
zC)&(}q#dwRro?oeu@jb`f<Fm#DZMl8;i|(b*<f)0@k^K;ytc3h_~|fPho$lL=NFm+
zd4`Y5T2BOx<FTtdt%tmR*{Ts#lqJ(-CYakBq>=GG^uoxV5G)i2rBa>l4a95HJhz{4
z4Y2x@-H>;x=gg8O0f2SOm`PkD2qB)qH)e(8VHrASrXN-<3)>i?Ht5aK+9nW_o6^BV
zum-UX`VJ%n+@np4B!}@18Hds&c`Hi9p4mg}AcDnVdBj=|y4kibZz|b}2=tFq^bURS
z>#3SIrw_1ID;yoFN^iB&Yx(*Ou%WFwlQ)M~Cz2HZNK`7}LUsSDi;K39iyv|56@Gc2
z*AYtf%(^K#pzU?rKH3wS`c*E~+nHLJ2R>{!3CR`_?5nl$;a(zF=@g+TZ>S6kff)Q8
z^km=EgW3Sng0?3g5E9>0qAh%14(xseA`vadx=|g^XkL&Q$gKG8mLo9vb<Ro~%)R=b
zN3;AcK99CKPQ)r#>EcBv%zozBD*cW2_M_0-7OIzZw$01fn9^v>3WIExr`&g+CeTra
zad(P`!95M<CgM8dU2Kn7{cLg;WJ+hP(o*X0hHrwhB1}`@gR}hjGS!Z}X<_N_0AsL(
z#@UrUcITKUKQVX#-;*sG2nF{nJF`os4l2L#NMP#Xqd5Go)0Gy-glha?J!up_Le3T*
z$**{b-&E)_=||KV4(&8f6!Ez$sAdt^LMOa}f)P`Eh7Ot1NPxIRvC0ve#ra1Or3@DG
z^p;$j9WE~py`SIJ=EZ5j8fW<=ordXckh$x~^1<WSh~#}QN#h7?5SwtZLuW}RP8f@|
zmAJy?xcp|t3{-80!Rpkfl+8oR`2HZF5*iq$h02Psra9(XOq+o<;@}@vN+%^<jPq>u
z3T{e28z9D8kPXOuY1Wg}k3+(qd!)EBhiZVah^Lw(<0@sCHpzwD`WTx^aOG2{Y57um
zQHQ%c<+oy=(zdPm__@OscZQ(UKF;vrVQ&v51425F8KNcSr-Vbr&uKJ6?uZqGKG$4-
z$&<V(BiFWTvs|~x8_w-IlQ2;az`|SM14}Q5_<5mY4z&(L$Kj>RIUZp5<f^(0VjM&_
zN_b4$iy2P`17O|wF34z0l;&fqvyfSFLucN8XWi`)pGzRtdFj;On<H3IVd?#nz-m}n
zsPKJ3y(dzN_@eGy7u0NX$6^O7PWi2m&97UauGJpx*C)^Im`LuGIwdo^UH*|&2zX#l
z{MFpJ{iK-Z*Cxs7z|zxZ=`}C$e*2*WNAaP~{4)P}yLZ*$)7Knf+edDPH^s(0bwi@Y
z+Oi(-{T@nNfF}#5)wlS*ieE!`oTSqu)MM|ZvzV6y9Dg214^I#r+J7)B*UT8N)f`xO
z%B`4PuYK|a&W8q<m^m==v}H#{je$8s+kvYO%f1Iz9Kq_dS{0?^%Zam~<v)1loD5@Y
z2VEV7*S%JL&vz_!LC9%A<@b`Afrdbe?|WH2_Sx{;m}hBEFOzihp3MPJo|e0T87D8u
zkEl?-UoN1Q-Dr}FJOi$Jn<}28iRb=e%Q<NXGV}83e28KD90i88wV4fI$-F4Z{wDNA
zy3>vDGMf+n0$pY~+tr&EU(o~i(ZpnEOehr;hvn0@fREvRurv)Xm%N)Y>3E4F0yv`)
zXJk)x6UiP5dg&^?!TQJBfIcS+BEi7qz5g$7Q|ry54f{{pJIBiLUk7<0glxPNf6w;(
zXWIV9W&KC#hS;$Y{v$rlMhe#lfP@3cFks;TF#mf#2=bTk41je0e+g{(5Nm?}mWbvc
zg|ozmun|x~MmXO5quqd=6fPGZ@`(TTU)ukj;a{?o!Y)EEIJp0pjg<hx#!g51cijIm
zP#B8;Y`|C%Kz_2*|1S>kSAF3hIRZS0C_dVMnRv@d3SK3EfpK#7urYV~=x%GHqVPXc
z{I|#cPhXDqf9!Zr!oZl>+c}%tIsbnJ|6^Iu|BoOn{U3wo&dxwP%ReR|O`MqjOmK*k
R6xI#0%Snec!TOi~{{cHKHu(Sm
diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index 0a29d868..9c912b3c 100644
--- a/Documentation_texfol/documentation.tex
+++ b/Documentation_texfol/documentation.tex
@@ -1100,14 +1100,15 @@ Cf. figure (\ref{fig:periodicity}).
\centering
\includegraphics{periodicity}
\caption{Cas de périodicité dans local\_extrema. Le tableau
- représenté est extr\_map. Pour définir la case (1, j), il faut
- avoir copié la case (m, j - 1) (qui peut être non nulle) dans
- local\_map. Pour définir la case (m, j), il faut avoir copié les
- cases (1, j - 1) et (1, j) (dont l'une peut être non nulle) dans
- local\_map.}
+ représenté est extr\_map. Pour définir la case (2, j), il faut
+ avoir copié la case (m - 1, j - 1) (qui peut être non nulle). Pour
+ définir la case (m - 1, j), il faut avoir copié la case (2, j)
+ (qui peut être non nulle). (La case (2, j - 1) doit aussi avoir
+ été copiée mais cela a été fait au passage sur la colonne
+ précédente.)}
\label{fig:periodicity}
\end{figure}
-Cas de périodicité : field(i + m, :) = field(i, :).
+Cas de périodicité : field(i + m - 2, :) = field(i, :).
\subsection{get\_1\_outerm}
diff --git a/Tests/test_get_snapshot.f b/Tests/test_get_snapshot.f
index 15b54757..cdf6a7fd 100644
--- a/Tests/test_get_snapshot.f
+++ b/Tests/test_get_snapshot.f
@@ -3,16 +3,18 @@ program test_get_snapshot
use, intrinsic:: ieee_arithmetic, only: IEEE_SUPPORT_DATATYPE, &
ieee_support_nan
- use derived_types, only: snapshot
- use dispatch_snapshot_m, only: dispatch_snapshot
- use get_snapshot_m, only: get_snapshot
- use init_shapefiles_m, only: init_shapefiles
+ ! Libraries:
use jumble, only: new_unit
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, find_coord, nf95_inquire_dimension, &
nf95_get_var, nf95_close
- use nr_util, only: assert, deg_to_rad
+ use nr_util, only: assert, deg_to_rad, twopi
use shapelib, only: shpfileobject, shpclose
+
+ use derived_types, only: snapshot
+ use dispatch_snapshot_m, only: dispatch_snapshot
+ use get_snapshot_m, only: get_snapshot
+ use init_shapefiles_m, only: init_shapefiles
implicit none
@@ -32,6 +34,9 @@ program test_get_snapshot
! Just the first two values to get the corner and step:
real longitude(2), latitude(2) ! in degrees
+ real step(2) ! longitude and latitude steps, in rad
+ logical periodic ! grid is periodic in longitude
+
!--------------------------------------------------------------
call assert(IEEE_SUPPORT_DATATYPE(0.), ieee_support_nan(0.), &
@@ -52,12 +57,18 @@ program test_get_snapshot
call nf95_get_var(ncid, varid, latitude)
call nf95_close(ncid)
+
+ step = [longitude(2) - longitude(1), latitude(2) - latitude(1)] * deg_to_rad
+
+ ! As we are requiring the grid spacing to be uniform, the value of
+ ! "periodic" must be consistent with the values of step(1) and nlon:
+ periodic = nint(twopi / step(1)) == nlon
+ print *, "periodic = ", periodic
call get_snapshot(s, min_amp, m = 1, n_proc = 1, k_end = 1, max_delta = 4, &
nlon = nlon, nlat = nlat, k = 1, max_radius = [20, 20], &
- corner = [longitude(1), latitude(1)] * deg_to_rad, &
- step = &
- [longitude(2) - longitude(1), latitude(2) - latitude(1)] * deg_to_rad)
+ corner = [longitude(1), latitude(1)] * deg_to_rad, step = step, &
+ periodic = periodic)
call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
diff --git a/Tests/test_local_extrema.f b/Tests/test_local_extrema.f
index b5792c0a..fe893ffe 100644
--- a/Tests/test_local_extrema.f
+++ b/Tests/test_local_extrema.f
@@ -3,6 +3,8 @@ program test_local_extrema
! Reads ADT from a NetCDF file and writes the list of extrema to a
! CSV file. Also creates a NetCDF file containing the map of extrema.
+ ! Note: the grid spacing does not need to be uniform.
+
! Libraries:
use jumble, only: new_unit, get_command_arg_dyn
use netcdf, only: nf90_nowrite, nf90_clobber, NF90_FLOAT, NF90_INT
@@ -21,8 +23,12 @@ program test_local_extrema
character(len = :), allocatable:: filename
real, allocatable:: longitude(:) ! (nlon) in degrees
real, allocatable:: latitude(:) ! (nlat) in degrees
- real, allocatable:: ssh(:, :) ! (nlon, nlat) sea-surface height, in m
- integer, allocatable:: extr_map(:, :) ! (nlon, nlat) map of extrema
+
+ real, allocatable:: ssh(:, :) ! (1 - copy:nlon + copy, nlat)
+ ! sea-surface height, in m
+
+ integer, allocatable:: extr_map(:, :) ! (1 - copy:nlon + copy, nlat)
+ ! map of extrema
integer, allocatable:: ind_extr(:, :) ! (2, n_extr)
! indices in the two dimensions of extrema
diff --git a/Tests/test_max_speed_contour_ssh.f b/Tests/test_max_speed_contour_ssh.f
index c4a40f51..10a18029 100644
--- a/Tests/test_max_speed_contour_ssh.f
+++ b/Tests/test_max_speed_contour_ssh.f
@@ -2,18 +2,20 @@ program test_max_speed_contour_ssh
use, intrinsic:: ISO_FORTRAN_ENV
- use max_speed_contour_ssh_m, only: max_speed_contour_ssh
+ ! Libraries:
+ use jumble, only: get_command_arg_dyn
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var
use nr_util, only: assert
+ use max_speed_contour_ssh_m, only: max_speed_contour_ssh
+
implicit none
integer:: ilon_llc = 21, ilat_llc = 215 ! lower left corner
integer:: ilon_urc = 49, ilat_urc = 231 ! upper right corner
integer n_lon, n_lat
- integer length, status
character(len = :), allocatable:: adt_file, velocity_file
integer ncid, varid
real, allocatable:: ssh(:, :) ! (n_lon, n_lat) sea-surface height, in m
@@ -28,14 +30,8 @@ program test_max_speed_contour_ssh
call assert(COMMAND_ARGUMENT_COUNT() == 2, &
"Required arguments: ADT-file velocity-file")
-
- call get_command_argument(number = 1, length = length, status = status)
- allocate(character(len = length):: adt_file)
- call get_command_argument(1, adt_file)
-
- call get_command_argument(number = 2, length = length, status = status)
- allocate(character(len = length):: velocity_file)
- call get_command_argument(2, velocity_file)
+ call get_command_arg_dyn(1, adt_file)
+ call get_command_arg_dyn(2, velocity_file)
write(unit = error_unit, nml = main_nml)
write(unit = error_unit, fmt = *) "Enter namelist main_nml."
diff --git a/Tests/test_set_all_outerm.f b/Tests/test_set_all_outerm.f
index dfc52246..53f22a90 100644
--- a/Tests/test_set_all_outerm.f
+++ b/Tests/test_set_all_outerm.f
@@ -1,17 +1,19 @@
program test_set_all_outerm
+ ! Libraries:
use contour_531, only: null_polyline
- use derived_types, only: snapshot
use jumble, only: new_unit, get_command_arg_dyn
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
find_coord, nf95_inquire_dimension, nf95_get_att
use nr_util, only: pi
- use set_all_outerm_m, only: set_all_outerm
use shapelib, only: shpt_point, shpt_polygon, shpfileobject, ftdouble, &
shpclose, ftinteger
use shapelib_03, only: shp_create_03, dbf_add_field_03, shp_append_point_03, &
dbf_write_attribute_03, shp_append_null_03, shp_append_object_03
+
+ use derived_types, only: snapshot
+ use set_all_outerm_m, only: set_all_outerm
implicit none
@@ -82,7 +84,7 @@ program test_set_all_outerm
call set_all_outerm(s, min_amp, max_radius, &
corner = [longitude(1), latitude(1)] * deg_over_rad, step &
= [longitude(2) - longitude(1), latitude(2) - latitude(1)] &
- * deg_over_rad, ssh = ssh)
+ * deg_over_rad, periodic = .false., ssh = ssh)
call shp_create_03("extremum_1", shpt_point, hshp_extremum)
call dbf_add_field_03(ifield, hshp_extremum, 'ssh', ftdouble, nwidth = 13, &
diff --git a/Tests/tests.json b/Tests/tests.json
index 8d6c9f5e..670164bc 100644
--- a/Tests/tests.json
+++ b/Tests/tests.json
@@ -5,7 +5,8 @@
"$input_dir/example.nc"],
"title" : "Good_contour",
"required": [["$input_dir/outside_points_1.csv", "outside_points.csv"]],
- "description": "3 contours at that level. 2 contain the inside point, one of which contains the 2 outside points."
+ "description":
+ "3 contours at that level. 2 contain the inside point, one of which contains the 2 outside points."
},
{
"stdin_filename" : "$input_dir/good_contour_2.txt",
@@ -13,7 +14,8 @@
"$input_dir/example.nc"],
"title" : "Good_contour_2",
"required": [["$input_dir/outside_points_2.csv", "outside_points.csv"]],
- "description": "Select another good contour. Case where one of the contours tested does not contain inside point."
+ "description":
+ "Select another good contour. Case where one of the contours tested does not contain inside point."
},
{
"stdin_filename" : "$input_dir/no_good_contour.txt",
@@ -70,13 +72,15 @@
},
{
"args" : ["$compil_prod_dir/test_max_speed_contour_ssh",
- "$large_input_dir/h_2006_01_01.nc", "$large_input_dir/uv_2006_01_01.nc"],
+ "$large_input_dir/h_2006_01_01.nc",
+ "$large_input_dir/uv_2006_01_01.nc"],
"title" : "Max_speed_contour_ssh",
"input" : "&main_nml /\n"
},
{
"args" : ["$compil_prod_dir/test_max_speed_contour_ssh",
- "$large_input_dir/h_2006_01_01.nc", "$large_input_dir/uv_2006_01_01.nc"],
+ "$large_input_dir/h_2006_01_01.nc",
+ "$large_input_dir/uv_2006_01_01.nc"],
"title" : "Max_speed_contour_ssh_north",
"stdin_filename" : "$input_dir/max_speed_contour_ssh_nml.txt"
},
diff --git a/derived_types.f b/derived_types.f
index 7790e056..8e1cc019 100644
--- a/derived_types.f
+++ b/derived_types.f
@@ -42,9 +42,11 @@ module derived_types
integer number_vis_eddies ! number of visible eddies
- integer, allocatable:: extr_map(:, :) ! (nlon, nlat)
- ! At a point of extremum SSH: identification number or this
- ! extremum. 0 at other points.
+ integer, allocatable:: extr_map(:, :)
+ ! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is
+ ! periodic in longitude, else (nlon, nlat). At a point of
+ ! extremum SSH: identification number or this extremum. 0 at
+ ! other points.
integer, allocatable:: ind_extr(:, :) ! (2, number_vis_eddies)
! List of coordinates of ssh extrema in the global grid, in index
diff --git a/get_1_outerm.f b/get_1_outerm.f
index 0689f934..a4d51c50 100644
--- a/get_1_outerm.f
+++ b/get_1_outerm.f
@@ -14,10 +14,20 @@ contains
! at innermost level. If a contour is found then its level is
! farther from the extremum than innermost level.
+ ! The length of the interval of longitudes of the domain, step(1)
+ ! * (size(ssh, 1) - 1), should be < 180 degrees, so that the
+ ! geometrical processing done here is non-ambiguous. The
+ ! longitudes of outside points and the target extremum should be
+ ! in the interval [corner(1), corner(1) + step(1) * (size(ssh, 1)
+ ! - 1)] as there will be no automatic shifting by a mulitple of
+ ! 360 degrees.
+
+ ! Libraries:
use contour_531, only: polyline
+ use nr_util, only: assert
+
use derived_types, only: ssh_contour, missing_ssh
use good_contour_m, only: good_contour
- use nr_util, only: assert
use spherical_polyline_area_m, only: spherical_polyline_area
logical, intent(in):: cyclone
diff --git a/get_snapshot.f b/get_snapshot.f
index fbf9c299..8f693f3b 100644
--- a/get_snapshot.f
+++ b/get_snapshot.f
@@ -5,7 +5,7 @@ module get_snapshot_m
contains
subroutine get_snapshot(s, min_amp, m, n_proc, k_end, max_delta, nlon, &
- nlat, k, max_radius, corner, step)
+ nlat, k, max_radius, corner, step, periodic)
use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
@@ -44,12 +44,15 @@ contains
! corner of the whole grid, in rad
real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
+ logical, intent(in):: periodic ! grid is periodic in longitude
! Local:
integer ncid
- real ssh(nlon, nlat) ! sea-surface height, in m
- real u(nlon, nlat), v(nlon, nlat) ! wind, in m s-1
+ real ssh(1 - merge(max_radius(1), 0, periodic):nlon + merge(max_radius(1), &
+ 0, periodic), nlat) ! sea-surface height, in m
+ real, dimension(1 - merge(max_radius(1), 0, periodic):nlon &
+ + merge(max_radius(1), 0, periodic), nlat):: u, v ! wind, in m s-1
integer i
! Window around each extremum:
@@ -76,7 +79,7 @@ contains
! and v when we compute the max-speed contours.)
call nf95_close(ncid)
- call set_all_outerm(s, min_amp, max_radius, corner, step, ssh)
+ call set_all_outerm(s, min_amp, max_radius, corner, step, periodic, ssh)
! Done with outermost contours, now let us take care of the
! max-speed contours.
@@ -88,11 +91,12 @@ contains
llc = floor(convert_to_ind(minval(s%list_vis(i)%out_cont%points, &
dim = 2), corner, step))
- urc = min(ceiling(convert_to_ind(maxval( &
- s%list_vis(i)%out_cont%points, dim = 2), corner, step)), &
- [nlon, nlat])
- ! (min should have no effect except because of roundup error)
-
+ urc = ceiling(convert_to_ind(maxval( &
+ s%list_vis(i)%out_cont%points, dim = 2), corner, step))
+
+ if (.not. periodic) urc = min(urc, [nlon, nlat])
+ ! (should have no effect except because of roundup error)
+
call set_max_speed(s%list_vis(i), s%ind_extr(:, i) - llc + 1, &
nearby_extr(s%extr_map(llc(1):urc(1), llc(2):urc(2)), &
s%list_vis, i), ssh(llc(1):urc(1), llc(2):urc(2)), &
@@ -119,7 +123,7 @@ contains
use netcdf95, only: nf95_inq_varid, nf95_get_var, nf95_get_att
real, intent(out):: values(:, :)
- ! (1 - merge(max_radius, 0, periodic):nlon + merge(max_radius, 0,
+ ! (1 - merge(max_radius(1), 0, periodic):nlon + merge(max_radius(1), 0,
! periodic), nlat) ssh, u or v
character(len = *), intent(in):: name ! of NetCDF variable
@@ -138,6 +142,13 @@ contains
! Change the missing value:
where (values(1:nlon, :) == Fill_Value) values(1:nlon, :) = new_fill_value
+ if (periodic) then
+ ! Extend in longitude:
+ values(1 - max_radius(1):0, :) &
+ = values(nlon + 1 - max_radius(1):nlon, :)
+ values(nlon + 1:nlon + max_radius(1), :) = values(1:max_radius(1), :)
+ end if
+
end subroutine get_var
end subroutine get_snapshot
diff --git a/good_contour.f b/good_contour.f
index 69f41f8f..491220f7 100644
--- a/good_contour.f
+++ b/good_contour.f
@@ -13,6 +13,14 @@ contains
! first satisfying contour it finds. If such a contour does not
! exist, the procedure returns a null polyline.
+ ! The length of the interval of longitudes of the domain, step(1)
+ ! * (size(z, 1) - 1), should be < 180 degrees, so that the
+ ! geometrical processing done here is non-ambiguous. The
+ ! longitudes of outside points and inside point should be in the
+ ! interval [corner(1), corner(1) + step(1) * (size(z, 1) - 1)] as
+ ! there will be no automatic shifting by a mulitple of 360
+ ! degrees.
+
use contour_531, only: polyline, find_contours_reg_grid, null_polyline
use geometry, only: polygon_contains_point
diff --git a/local_extrema.f b/local_extrema.f
index 99063172..aba4ecd6 100644
--- a/local_extrema.f
+++ b/local_extrema.f
@@ -17,6 +17,9 @@ contains
! in extr_map but there is no duplication in ind_extr,
! innermost_level and local_min.
+ ! Note that no coordinate grid is used here so there is no
+ ! assumption on the grid underlying "field".
+
use nr_util, only: assert_eq
integer, intent(out):: extr_map(:, :) ! (m, n) Map of
diff --git a/mean_speed.f b/mean_speed.f
index d3737069..5823a291 100644
--- a/mean_speed.f
+++ b/mean_speed.f
@@ -9,6 +9,7 @@ contains
! Interpolates the wind at each point of the input polygon and
! computes the mean azimuthal speed at interpolation points.
+ ! Libraries:
use contour_531, only: polyline
use numer_rec_95, only: bilinear_interp2_reg
diff --git a/nearby_extr.f b/nearby_extr.f
index 45d66bac..6f69677f 100644
--- a/nearby_extr.f
+++ b/nearby_extr.f
@@ -6,6 +6,9 @@ contains
pure function nearby_extr(extr_map, list_vis, i)
+ ! Returns a list of extrema that cannot be engulfed in a good
+ ! contour around the target extremum.
+
use derived_types, only: eddy
real, allocatable:: nearby_extr(:, :) ! (2, :) longitude and
diff --git a/set_all_outerm.f b/set_all_outerm.f
index 5d6f8573..893d14d2 100644
--- a/set_all_outerm.f
+++ b/set_all_outerm.f
@@ -4,7 +4,7 @@ module set_all_outerm_m
contains
- subroutine set_all_outerm(s, min_amp, max_radius, corner, step, ssh)
+ subroutine set_all_outerm(s, min_amp, max_radius, corner, step, periodic, ssh)
! Extraction of eddies: find extrema and set all outermost
! contours in the snapshot "s". Not a function because "s" is not
@@ -13,6 +13,7 @@ contains
! Libraries:
use jumble, only: argwhere
use numer_rec_95, only: indexx
+ use nr_util, only: twopi
use derived_types, only: snapshot
use get_1_outerm_m, only: get_1_outerm
@@ -33,43 +34,67 @@ contains
! corner of the global grid, in rad
real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
+ logical, intent(in):: periodic ! grid is periodic in longitude
- real, intent(in):: ssh(:, :) ! (nlon, nlat) sea-surface height, in m
+ real, intent(in):: ssh(:, :)
+ ! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is periodic
+ ! in longitude, else (nlon, nlat). Sea-surface height, in m.
! Local:
integer nlon, nlat
+
real, allocatable:: innermost_level(:) ! (s%number_vis_eddies)
- ! ssh level of the innermost contour around each extremum, in
+ ! SSH level of the innermost contour around each extremum, in
! m. By construction, innermost_level < extremum for a maximum, >
! extremum for a minimum.
logical, allocatable:: cyclone(:) ! (s%number_vis_eddies)
- integer i, l
+ integer i, l, copy
integer n_cycl ! number of cyclones
real, parameter:: min_area = 8e8
! minimum area of an outermost contour, in m2
integer, allocatable:: sorted_extr(:) ! (s%number_vis_eddies)
- ! identifying number of extrema, first minima sorted in order of
- ! decreasing SSH, followed by maxima sorted in order of increasing
- ! SSH
+ ! Sorted identifying number of extrema: first, minima sorted in
+ ! order of decreasing SSH, and second, maxima sorted in order of
+ ! increasing SSH.
integer, allocatable:: selection(:)
! identifying numbers of a selection of extrema
! Window around each extremum:
+
integer llc(2) ! indices in global grid of lower left corner
integer urc(2) ! indices in global grid of upper right corner
+ real corner_window(2) ! longitude and latitude of the window around each
+ ! extremum, in rad
+
+ real, allocatable:: outside_points(:, :) ! (2, :) longitude and
+ ! latitude, in rad, of all the significant extrema, except the
+ ! target extremum
+
!--------------------------------------------------------------
- nlon = size(ssh, 1)
+ copy = merge(max_radius(1), 0, periodic)
+ nlon = size(ssh, 1) - 2 * copy
nlat = size(ssh, 2)
- allocate(s%extr_map(nlon, nlat))
- call local_extrema(s%extr_map, s%ind_extr, innermost_level, cyclone, ssh, &
- periodic = .false., my_lbound = [1, 1])
+ allocate(s%extr_map(1 - copy:nlon + copy, nlat))
+ copy = merge(1, 0, periodic)
+ call local_extrema(s%extr_map(1 - copy:nlon + copy, :), s%ind_extr, &
+ innermost_level, cyclone, ssh(1 - copy:nlon + copy, :), periodic, &
+ my_lbound = [1 - copy, 1])
+
+ if (periodic) then
+ ! Extend in longitude:
+ s%extr_map(1 - max_radius(1):- 1, :) &
+ = s%extr_map(nlon + 1 - max_radius(1):nlon - 1, :)
+ s%extr_map(nlon + 2:nlon + max_radius(1), :) &
+ = s%extr_map(2:max_radius(1), :)
+ end if
+
s%number_vis_eddies = size(s%ind_extr, 2)
allocate(s%list_vis(s%number_vis_eddies), sorted_extr(s%number_vis_eddies))
@@ -96,19 +121,33 @@ contains
i = sorted_extr(l)
! Define the geographical window around each eddy extremum:
- llc = max(s%ind_extr(:, i) - max_radius, 1)
- urc = min(s%ind_extr(:, i) + max_radius, [nlon, nlat])
+
+ llc = s%ind_extr(:, i) - max_radius
+ urc = s%ind_extr(:, i) + max_radius
+
+ if (.not. periodic) then
+ llc = max(llc, 1)
+ urc = min(urc, [nlon, nlat])
+ end if
+
+ corner_window = corner + (llc - 1) * step
! No need to consider contours with amplitudes < min_amp:
if (abs(s%list_vis(i)%ssh_extr - innermost_level(i)) < min_amp) &
innermost_level(i) = s%list_vis(i)%ssh_extr &
+ merge(min_amp, - min_amp, s%list_vis(i)%cyclone)
+ outside_points = nearby_extr(s%extr_map(llc(1):urc(1), llc(2):urc(2)), &
+ s%list_vis, i)
+
+ if (periodic) outside_points(1, :) = outside_points(1, :) &
+ + ceiling((corner_window(1) - outside_points(1, :)) / twopi) * twopi
+ ! (Shift the longitude of each point to a value close to the
+ ! longitude of the target extremum.)
+
s%list_vis(i)%out_cont = get_1_outerm(s%list_vis(i)%cyclone, &
- s%list_vis(i)%coord_extr, innermost_level(i), &
- nearby_extr(s%extr_map(llc(1):urc(1), llc(2):urc(2)), s%list_vis, &
- i), ssh(llc(1):urc(1), llc(2):urc(2)), &
- corner = corner + (llc - 1) * step, step = step)
+ s%list_vis(i)%coord_extr, innermost_level(i), outside_points, &
+ ssh(llc(1):urc(1), llc(2):urc(2)), corner_window, step)
s%list_vis(i)%valid = s%list_vis(i)%out_cont%area >= min_area
end do
diff --git a/set_max_speed.f b/set_max_speed.f
index 5a458aff..42ede51d 100644
--- a/set_max_speed.f
+++ b/set_max_speed.f
@@ -9,6 +9,13 @@ contains
! Defines the components speed_cont, max_speed and radius4 of argument e.
+ ! The length of the interval of longitudes of the domain, step(1)
+ ! * (size(ssh, 1) - 1), should be < 180 degrees, so that the
+ ! geometrical processing done here is non-ambiguous. The
+ ! longitudes of outside points should be in the interval
+ ! [corner(1), corner(1) + step(1) * (size(ssh, 1) - 1)] as there
+ ! will be no automatic shifting by a mulitple of 360 degrees.
+
use, intrinsic:: IEEE_ARITHMETIC, only: IEEE_IS_NAN
! Lbraries:
--
GitLab