From af77c95f29001f9fadf6486ca35cb0b29c3cff63 Mon Sep 17 00:00:00 2001 From: Rowan Orlijan-Rhyne Date: Tue, 6 Aug 2024 16:44:58 -0700 Subject: [PATCH] P3 aspect ratio as option, docs on regimes and density, The Duck --- docs/src/P3Scheme.md | 30 +++- docs/src/assets/p3_mascot.png | Bin 0 -> 39214 bytes docs/src/plots/P3AspectRatioPlot.jl | 105 +++++++++++ docs/src/plots/P3SchemePlots.jl | 46 +++-- docs/src/plots/P3TerminalVelocityPlots.jl | 210 +++++++++++++++++----- src/P3_particle_properties.jl | 38 +++- src/P3_size_distribution.jl | 4 +- src/P3_terminal_velocity.jl | 157 ++++++++++++---- test/p3_tests.jl | 99 ++++++++-- test/performance_tests.jl | 15 +- 10 files changed, 576 insertions(+), 128 deletions(-) create mode 100644 docs/src/assets/p3_mascot.png create mode 100644 docs/src/plots/P3AspectRatioPlot.jl diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md index 8752d9b49..0007310c5 100644 --- a/docs/src/P3Scheme.md +++ b/docs/src/P3Scheme.md @@ -68,6 +68,12 @@ where TODO - Check units, see in [issue #151](https://github.com/CliMA/CloudMicrophysics.jl/issues/151) Below we show the m(D) and a(D) regimes replicating Figures 1 (a) and (b) from [MorrisonMilbrandt2015](@cite). +We also show the density as a function of D. +Note that because graupel is completely filled with rime, + the density (``\rho_{g}``) is independent of ``D`` between ``D_{gr}`` and ``D_{cr}``. +Following [MorrisonMilbrandt2015](@cite), for nonspherical particles + ``\rho_{ice}``is assumed to be equal to the mass of the particle + divided by the volume of a sphere with the same particle size ```@example include("plots/P3SchemePlots.jl") ``` @@ -115,7 +121,7 @@ N_{ice} = \int_{0}^{\infty} \! N'(D) \mathrm{d}D = \int_{0}^{\infty} \! N_{0} D^ ``L_{ice}`` depends on the variable mass-size relation ``m(D)`` defined above. We solve for ``L_{ice}`` in a piece-wise fashion defined by the same thresholds as ``m(D)``. -As a result ``L_{ice}`` can be expressed as a sum of inclomplete gamma functions, +As a result ``L_{ice}`` can be expressed as a sum of incomplete gamma functions, and the shape parameters are found using iterative solver. | condition(s) | ``L_{ice} = \int \! m(D) N'(D) \mathrm{d}D`` | gamma representation | @@ -175,6 +181,14 @@ V(D) = \phi^{\kappa} \sum_{i=1}^{j} \; a_i D^{b_i} e^{-c_i \; D} where ``\phi = (16 \rho_{ice}^2 A(D)^3) / (9 \pi m(D)^2)`` is the aspect ratio, and ``\kappa``, ``a_i``, ``b_i`` and ``c_i`` are the free parameters. +Note that ``\phi = 1`` corresponds to spherical particles + (small spherical ice (``D < D_{th}``) and graupel (``D_{gr} < D < D_{cr}``)). +The plot provided below helps to visualize the transitions between spherical and nonspherical regimes. +```@example +include("plots/P3AspectRatioPlot.jl") +``` +![](P3Scheme_aspect_ratio.svg) + The mass-weighted fall speed (``V_m``) and the number-weighted fall speed (``V_n``) are calculated as ```math V_m = \frac{\int_{0}^{\infty} \! V(D) m(D) N'(D) \mathrm{d}D}{\int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D} @@ -188,9 +202,21 @@ We also plot the mass-weighted mean particle size ``D_m`` which is given by: D_m = \frac{\int_{0}^{\infty} \! D m(D) N'(D) \mathrm{d}D}{\int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D} ``` -Below w show these relationships for small, medium, and large ``D_m`` +Below we provide plots of these relationships for small, medium, and large ``D_m``: + the first row highlights the particle size regime, + the second displays ``D_m`` of the particles, + the third shows the aspect ratio ``\phi (D_m)``, + and the final row exhibits ``V_m``. They can be compared with Figure 2 from [MorrisonMilbrandt2015](@cite). ```@example include("plots/P3TerminalVelocityPlots.jl") ``` ![](MorrisonandMilbrandtFig2.svg) + +## Acknowledgments + +Click on the P3 mascot duck to be taken to the repository + in which the authors of [MorrisonMilbrandt2015](@cite) and others + have implemented the P3 scheme in Fortran! + +[![P3 mascot](assets/p3_mascot.png)](https://github.com/P3-microphysics/P3-microphysics) diff --git a/docs/src/assets/p3_mascot.png b/docs/src/assets/p3_mascot.png new file mode 100644 index 0000000000000000000000000000000000000000..dc785fb208b4d635955e464262d4782fbc4f30a0 GIT binary patch literal 39214 zcmZU(1ymbAvoMJ)s8fpqy=w#>s008TqqO8_yJ^flP(NJE$k8ILd0RW)9y^M^; zI~f^>hP$h+y^{?9z?5cfVS)FKm1)G%(!ye7@+~vEyN_0EY@(J$$WZST1bILnb0ThZ z5IGqL6qlsO>BDtE8#c}pU4f-7d~%euIx=NYzxD<4kY1PECqbdSlD_Gc|why z|A-NFF_kAZ=oP0prHGW< zcZV>t^iCR2UBOD`6d6cKtWQN1r9tZLfjkaF&i-p1~vFsUsfX%!nv zM8psb#QruJLHM1HIe|V$c2`P!e`SZ(*VMJJ15DjAgM%>HGi`nDL%f_JO zl`9V?xC#mI;iI|a4q$Qi)v$Ut2JdWDRRJuoFd6_DV-G-jg@CVx?6m*@h##W?sIN8N zYmqHN_#fBQD5O#(Amb*65{LZTLm&`Q_xHBKTC(#0&HnmJjOK%Ntt-;AOh{~7uJ zvc&%}^MAFz+F2Z3l;eMVCXTLP;2Qt{NCMu;O6m9lPrCi;Y=(pO^?iC=50~HmqG48Y z()vQD45I(_g}s_u)$%vm zgB^8s8On4jz95egb?gyUnO#xrV%cD1uz}b#!Yk<_zaX30Vr0Pmas*)?iA%e-`3f_V-@OG>wJ`$f6s0E@x4$h2vR^(E19Emf9ETx(g8vH;5A^$2-Y%(V z&KY&_KDFkPyrQX+VSelU;lJ%+puAYC5j5plP*LH=%z(I#3XTDyBRQXcPnnlredS-1?b{dz?M1GQ-e7xb(yi8w9e>dd-}kcL zJH1F%891VUno*|ozuV-E4w&GWD7QuIM()*NaYhTqL<69f{32nW2RMznt^v^y3|2Gr zM98ZtS8ZVaLF2X1?m~7}SEOaX^MLido3(2cmoAe^|EJ&g_cPzjdT!YYugTgZ`7tr2 z3;j+m6c9O2JN;t23s3;1V{Nk^Dss^!!15TseL3HW0p=;;N#d7lWw$QtXHkQpowbrhFni|nSXPV!Ls~&nfTgtoYckRzYFqGHB!Y-mI7EK zBGM9b0~*eHb0z`oPk)J@-HpUP<<-)o2V!ZHTK*y9KdU3Iw^4DiCc{KI{AE{&Mbkj- zqB>OCn7dem)yw^k+>wGe*RrNp)xZC5i{J(V+7{C1AFczfxJD_}@zi~5xF(Fo>v-S~ z8AL?Hb07_tt#QN_~hATM6rw z{t%p+xHNPkqhI}*SC_~qA|}=jX(6&8#rC%k{nZ(r253-)7LW}`^hI5-^jt|u-#zj6 z{WX939_=mQxN;h=rh8WEbn+bK(6FNJ23MROrnsQU$!M%&{e; zd#-4`?e^2B?f!G!$GcHemi9b3hPKRXm6wLn??6@}ve-nT!;%qp&R8#A=$$-WE_Rx+ z=K!s=WMn(bM)iEUZf#g6*3Bz(u@#3Fk2QuLMapr=3b)ViI5O=vO&aOBZT@a6|LN^~ z+FI}72p%ah;&X<_2zj3%3nD&Ji$2OC&ItFSFaU-n_tQudn^IL{i406ubBOkc=Ope< z?B+yE7ORP;DsZ}Ax(}Zu{_;TSMJ36wm?Xuj>(jzQ2w7eEgbXp(92zrG-d4`<{E#PG zAe#R%Wo*Hlv|v_OS?UAatV1sR3l?W(kks#d;cO@_KWqi{Wsv%L>jxA6gi(kP z!DVAj)^VnQJb$-Z+ug5*+)RfJ!L(?hVGRJ2Q|vFLRb)Qa zlKY_Vz!Am+LsGrvXG7&Fk5lJ1wsR0KUYY;EY^j1(tN8t~U zaV{ibvcFoD>U4c*jP^n5zkMFMnue0JiAxG};3lC>W`+B|MR8F{@KGoLNVEO%oWyVp zZk@{LXW6)c!Bb<0veHWRQi&j;MelDyr`a2Y{IRblOlzp1-3T%N`gEp_Xe%XEuIhrI za&k#y^`=8X0qRzQ%VTED(&zE&XlD;~PpGtSkhD!^mOjMOY96JOS3rPxpKza>n0p;F zfdo(^u&aosPXc1VmlI;M*Sbtz*gh{qo*5!P5HdJI0kw+5$aXDAW}1vj6S<<%-V3n0 zU?Iq4&7liSE&G%SoH|y1x=QJfD~L|{#>J<{UKZprJcnq0|8!l^B{9|Mhe?2x9KPuw zHz}`0tQKxm+D~8*6_w1hoXc9RVuSH*hiHbiUS4zf!5o!9j;?tHE{x)RkE2;zGngnj zMfPu^7P&s$##!sMi5`ra`#HrM(#8g7b;eiKDBS^D=EHfTFAg5Sz$)sY$gf=Pr$qAc z9@#(i=OUscZCK}N$DW*cPM!fEsE7vv!F~J|snsW#9VXP%B6q>>(ca1h*I~YQ z$AsP+s|#E<-ouit<4H@B>-`}XJY$)b2dPOKmIrs_gPH&eWTR5Kk-H~~_(5Vp9lG5q zg)zlA9s%=bPY-jC>%qU<8!}zY$to**E)))TZl9U>01k+xs_9`y$ujgU2b$Le;!7<^ zDfJo^%N@v51!6c5jW5*V){CnGx<6T_4v#faYpjeuOPv_4S@6s`fH*~?f_ixv3_Lq} zgWxEj#5d@lcr2+>6%_X}mQ}l3TS*}LeQ?f04B*6ktcUY;`O$lsLIs1-KfRumxNBf?>IaScpW9h2lj8i)A)QhQPus1 z-j$tEdS8$Cr2IEULb9LNUY=G&1Nr`N$NRsYXM&h-{rsr7uYUelr!dwOOTEd-$}bh|C`=NRV4Xurhz9%ae`?0ij#>X zij!<>H*7@*?Qm9u+DA*)QAL=bsSk|WzwpbTBGKCj=2*tFpHjZ@^b}#gY@$Tum6MFW~UI!nx}3_m@`B@wLL;SwS5l#yY5nL>f` z;Rr?$tVv}4!p5^cNcnKshoiFLX|zZ2SpAuA+6Xx`KDkl2VDY?j5)`YChMW8~wjz+m zH6h5{LzyMN%$Tfnwc8(K#Pa#<{Fw;e@^6wSFG<|rlRUo^q<8lo$*|823HLS|{fy9W zR6BIMek(6oN9Yqw{-0{E_rq~t*g~E|lJ}9NpW}`G{C%^z(te6~wu|Y`_r?9Wj$i5? zS6$jcjp;eQQ&wlw_1$tYfQ!0YhvO7K9pFyroP3=t9~n@I;M`OILeJkddaEtcT~)Xs1zp2?lv}^{MB` z;PzHeECbEbsCY6)Q(w+X{fsuI%hSD7K$~D@2XVVuQY&8QXM`Q;p)(6PnT={oDH2Zm zSRVCB>RhplQwddpb2G&NBICo8E1Zh*8~2*YQ`8n{O2q}#G5dR57-DEr6#W@RMd#oB z`OABfx0{0b;_AM$tHyI8!Ow-|h*WokyC6spq+#Tf>txoEYx6>T@DO`YfwqQ35PU?DA0z3 zxx!);ECzH##yNk504^-G0;S=e7rO~U zcbe1gZOPvOSM`V2#iiAq7p(ANaT02%h-Ih)A3LmfqPt*8gjmat zaQ;?4&)ui;=2Rko{nqit)Pemw?V}XQJXhBC$Zlb;p6nVXd48xF4u>+mFS}#vC1-oD z+@CVfm#5q9zr>2y0|W7kjZk}D5G39?N|6B|tx{Tvja~kMA^Xlw$LsqLM1bavd^UUt zin8zMIJv0ZQyPV3$5bpF?>=`{#ig>3Bk}S15$AdOx~cv2U4VD%{Hn7+`QQ`k95jgC zD>*3PAnNLH$=KLn19XALuBwWlyuTv*cH+v3)qk6CYe%iz=vK5?EP_P2MT)mSDwwBU z*G2XGvA3ZDk?Nsg6#X8$fbC1&nKXo4Qy@kA`VXM=)697dcT127%y%%iyhfz#fGZz{h(x4OPEf6+=|zbcf)ug6hFbgx<3je+bG zW;eB|VJi7}M4QOf(S^@mHEobAPEHg z5HirYMR&Ry?5r~6p91Q<>(g9H`98Kp`W(h|MR{IwC9>82o;Sa?M{p_%w9XC3!fK() zx~(>uVp%(eXDinRyYrn5<6nBKuCWjLD!V~PZSAJ22?x@40@xMjK^MO~PiD3cgL?j? z%&)4t(sQK!n&QfWZXI3kZtOBIl#Huq4%=X$fd^TH)C)d{6>$^Ja&migqf0J$wX}8C zSKsz}A!!vBBgnr|;nc^;Hu4#GCh5flsKxw4WiVI<{|3R5QVY@d_J!E|4&@4Oc2xvjsR!6yb#^pgV){xvJi$8A^6U!s@l=G4^m)d?fAdWp%l zj7HkQ4fq19enEV_5 zWNFnrsZ=kZzvb&!bhlgWJZek~DB7S*T993sp(LBn1StAV7|fZ0p!0aDy`tQ- z-gEf0z5Mrzi~5G&!RYnk!*_)?hx^gnr1?tm>txRxPl5UZ(M<0%r*eUB7kUo0nigzW zD-kXJo)q~X!X6`9&at3;g?uh$LXz9>X;lnz@Bh-y43}C)l+^V^@3WXQ09_F@77lhfC6RI4NvY_> z-tEbrA=cT-2u_9iEk{F3YG#p<{A%C_=7(Bs)yTwEJUCv+KlI2}TSronWjh&U>KRr3 zT>RL6VF~(|pfDE@a!-jMqkVtlL-kIbLbzZGMl%*llSNW~KN>)GxG(;VYPk9Wy5_m! z{Ls{hOJXEa`1OE`95E3{g6{Ggb$`?E(~&4ZPNZM-``BHvR8PTF!BJ&;OwN)I>6Iwr zN${D#%hmO>`U`E{zHH#QEjyy-Y5TkHlCP@(lTp!(PeJf~H6WQLw|5k-Nh}?ZEw8Su zsq2Y2iNuN9o6)JmMe=!k6_K2k7C|{&GRd4)9IND%wgxNF8j({zatxU!I;Nj8{#kgn za+u|fLgBs35!dqHfb~`ZvHO9GmatDJYd$Z&_i?w~_@fD2Vk!!)2nqp5P!2JVJCj_- zukgL@gn=%$QvLj4y5t7bo|*N$+p_XhpZu}63dUtYDL13sRj?soirUt7azW>9)0GGL=UoFj&ZsX8o@jl65Q70K) z`ML$Y*Hm|^|1jkbX>Y;P^-dm&K!3LU@@!n#){?-vxBs(o9ctJ0_a#7P$*1xVv z@;Fa9v?b@z1F<$ooZ{RIx+;1n78_Y(}L^mP?|%V_jg1u(D-K%rm%FzN+^%u zwp!FCCRViAwE^5=`3LIUSTn|94$pKY1=FqT*YW~A^?Wt;d`kSxdHRq3)TX_Sf2O)5 zvX|Il%F|$f)~8KkwX&#ATmNJ7{rNSYQ>B3Y-v=fXi5Lbu+3%aerWR0N+8=**Y5zDh z>Iw>{KYzQB04*N+d|e&#Qdkz``Ze$NXxIcj=>W{`{GS!#H|6eZ3vwb0)D&EA=Z*-Z z`>;v`(R!TIp`an#SHRTL4kZVY?K{ppnab_VazZV4MZ@w{S7l570&Q0x4T^)RvkxpUWlWG?0TT@{S=>l#(cn`i4&B=k*`q3p| zKK~H$HjtHV@?7y#=I%+XUDsL4*BHmUD!Y2HQ16$E$&Su*M{^T0;le6_`JM&x%0B86x<~RUoenU>1=-oXMwFP~+;R zSD6G^eWFHU`B&g081A2}HYeWSj+~r4BaHN`A37t$@li%-GzNyD_@kcTa7kc2P5goC z5k&zC=>-Ls;wG!t6paQ(o-NYJ#zrh4N@5~0s!NF7)HhUw{C#*)@W-xnu?wOtX!W@@ z`Q0kJ033=l_C8hePz_D&Xm=0f%NPGTrl>oNYAKC6fziI6@jUqQ+;wOw@uK{#s#xsE z^gs^RIg~1)G=Xg(Y?v1|!#>aN5{4i;cr#NhbGS2_-sha?%P-Weqh?SYD>nbD+OQLl z@1MNHBefl_C#*Q^P9D1b=eU5u5zqqMCk0vS^%>c*!r}zAKkLpte(pY$cvR>a9EgGf zs)WNyFc*gg>bNEvsbVaNrwO3cD)0K-rIG`V_vf~8g6V^vT%YO0ZX#mv(?mp@Be+V^ z)}Qy*4jC^9Dg5uC!gu-o)v7zaoz3IbGGeCoE*+)$!50$G>T?Y*?L84r3Zj(sN0eUZ zLzwx>o7&+qw!hA*0jkCGBVMBWZYZ%Pj5uPc+at_7;@FIw_HhYeuBo{^JJ*}yrjpim zM)a~mw*THUwccXRM$rw5Kd7QCVO^lMW-IROJ+)V=2kz-NCgRHDTa){W!;xx2aq`j5 z=KhGLN2pX0_=ao_pkGxTE~B+0?kGu|z1c*U+8j3+dJyE{rFd-UxctDQX~p ztT2fsZ|JK4FtZW?-|+FBwy3#OLw3rEhKZzPFQBe^ky z1*=-Xnv`2!j#!%Zn}QPhwhMlIFC|7O=a>zhf7pVp*@wwFqxoo7b@c@lz_ zMrY3@PAZW_b~S}?Ns5^zdIAM)@9V!Agu9lkjWF@NeSuRguvnn$= zAI%3F+|nX!o`#J~Gb;?YgnK);{KGoNUEcP27)Chl_op+aDAeT^WfvAIIC-4eJL|GM z2#a|oa)+V9=m+>Dn*{}^jvShjIFCnvcIDajJfMR$ZI5(y`GXmuftWCLxK^m8k+(AIU@?e+HOYC^*>kT<1L$U69! z{F}7GH(v@LC#GfJ63entbwu(u$d@<3yG*oJ=}{4A5#=e^X2Oc$R$2lq@a?3*bE{-L zQQ{RV<8#uF_@(j}C*#gYNwA?f5wQH}z?qs#Hu0Aq*`ekGXyKneErPZanGSq1Ur%m+ zp0N(1B;0Q+n@`a`3#>m^T<;-D4CyMm)hx?C8h~sWOhy+yN%Qy9IgZaWx(IHG0Gq(P z>6y+Nff0N?NZVshykp+KcOEQQvNSDM^lg0o?-z+s|5zyyEoJbd=YP{f0Fg%1@GY(7W*y2+CEwn zVI+8O_?HXbjK9J0`?l-+$Y3yvg`@_03i_z?z`0ia~%`wyo~!gT0ls6wz@KkmH*Y;b$5FEigWG zN*#`%ONP>j*_#8yiUGQqmQ-*N2%c38L)yJM<{AReBubF>1+ppNcV5a<)^DAf=`Qk+ zCW>WZ5HRuD(vXdQ`x}hxkO_+?5ms&b`Ba7BcZ=93n+zJ3!5qdcoxVD*67eKBY;dnn zY8IZ`eXYv7+4#W&+A#~qqa^H-dd`6Dqx{oo(Bnf^>$5=$i&#fe?IVnpm&uydsi#ud zXUOF1Ofo@r!Mk-dL)2sEDw)?g(?OQM`t!}Gu!(hF zVHN?;4G8g~@FMF>f>q#I=MxU2P3I(|VBv7M!ikY#^-OQo}EQg}Botzh^M`;v!GS1>9b zU*N87!O!Inj&J>9R)01uYx1sD)s3wpKKDrahzrnfBT1Hct)p}-^1@pxulpzk0Hd*& z$uP1)cXpT1^@JFqKCr1Dc;kbJOs&NUqbkNJ9A=%kw@QCjW6e(q7YwIL=dVtFH3Lw4 zR4|EW(BPrltp1~NNN5SaAq9jbMns*`SwP8XIEiT}2(?O3YgiM8s1qc4UC`bwHkQEM z@TJ}F2xXnQ&jT^!Tb8z$(DZhSM10k!$sz1=N?fHf^btv?{!6S?XQS(M+PlW#Q=bqOG zX>d_B|IZ`4D*0GT4WkwIcBd1o3QejY%p6d}tZQ~Tml&=4KM(N;gR1z6ZlL3zKdBvs z0jtCxL~hq7`q`Wum8X(N_KYW&R-B&HcEuxunW05{1}YA2s8i~p+i2?hIhk@wYW;ite5t`3EV zDVUT2=c-WEuz`^FR7^)(q5`Z%EjC&7+RTJPP%JFw7~-8d=k~WjNMmkH2o`mkYBIEl zp`LypVPDkkm+sU@{VEE5GWm^0mzI|1nx`qlL-UEUNL_RNYK$ur@xbzGZ`%kZ<2(uk z8~|$cso4pJSESm9c)<~Z{O^y}!0$}$46!qV70wSoH*mRz)%J2#0*_OzGGx?C)Ng`q zy4eorO4*2AX4ddT1rB1mj9QW7A)NW*LBe~cZy&ki{7}Aqjx*XKLZTeEho#@sWTl>w zN-2rUHBqA=VlTeX_AEtN`otbA_X^HLmA(O&G?HR6lJs+<`i6meei1p@KAheKt}qk! zx|JjT#TR=@UgScCW{b0szd;A0BLXdj7-ZZdV}JICpdXqqgvFX+jOWpS`%r-FjDv3L zh*&u5cG)S@f4+b73`A+NnRa-+=0yxKjB8aK{0^2C8{S>bkC9AM)ew4(r{=9BJn`h- z%0S&o=fuU(}AMoqoSls;5Sh(XH z7<e1VGRr)bX8cI)0^Z z!i|x7-6gB8G>2QXpcxD{B!2aWD-Vs6yKW*R*w?2jBl9*FW`(-jjJ+jJmoDP(qkE)A z{aoNF#!#vtVY}&q$Q9l%W)bh|W<5XfDoLE{FY@k{fk9sROraCap$Mf4wb@H7tm`|e z4q-d8R28g54xFLGATGFX0WQN4*RGS3pa{11wfl+m_ah~F$K@L5N6`uWHPaOCORBL1 zH?ihr_33qT=YLMk-;8685Dp&yU45%@`(X49Hc}eZX(*w(#}%)hllR_2a+2E6`LamW z&UPsXDaQ{_-?vZ;__dRKP;m7W2Q^=EpY~F5{lGZIW%o zoPZ>==P2R=oSB1VE&_#c@5m1`4PZMKnNjLQO3tcz(WAE?yCH5|YX&3n9Rlg#D)n`0 zYEcAr{EmJSSmAH?)&bRQB_)tuqdjlxIMlg@QBI2I@pVJMdUQknfv+w~`ZWpP@BLaC zr@1390ED)8k%qfV7a*mbow`eSZ-<^>9EuG&q+c8}X|TJYPBprVPW0j-P=J2OPk)__ zA-b{>?a7CRa*36?qSlkG!?td>y+V}>vAfrh5BrspLju%{Xc!8!lQI1bI`fX!( z`{Fv1`S-3&OQTIbDvuzMaDJtl3rQ_CAVa+84H4flYA_AN%jjCgRjKI|<=;2FiqUm5 z0EPd3&O$G*vy&kv5S}Z8Mt6h`C!uzavQ@uDb45Gj3R}tM1xOyk7&%~-I1tv2!K#JS z9d?<_U>|_VglXt(FPDVTYG50<_76p;ZlcH9tfqdlj0Mvbxl@dedja{$TU01E#rlF+ zqz=&`Cq6!*1s<)U9n2>*m&f0By5}t^fxkwq2Y^AAm8PT*GD?r%h5Heom2PYdY$Iz} zdMJB37%C)Xt!2Vv^I6c=!23VYq5sr-$>|k`4osw^A_2h|l(c1XKEUVIMj00SkNe;1=PB{n}C!r{_ z&Lbv^)(DCZvII<@kt;7dp_Dqz(JTkhT|H-@SOG)vbT zu03L^r>{>!;ZIWXJ{Tx%d)=i7!wXLPwv6o4Eif?(MHUn1%jD;6LaiKG;!ElO03b;u z;PeY2Y?>PM9md4)dV`BB3d-7EYpgkBO?oFe;`nUkE~_kX{$oozqnzF)^EE6#{&Q?w ztIJ?ZPs3<=&Kc7U8&!v2@m_G^bv9D2LPQv7`!kA~)k;R_4_Zp>*)YV=V2TS#Ec2>S$dm4qh3{6r z{7^}%D#x`%YS(N?W8&4sprmM`h^8I^?#d&`CNsupW^p_8HA$2hO8SX$h!p;X6bjZ> zx3rnYn73J=`CqR|Tt0f(CoNVGlW3ii2_;3B81Z9t@!?5MS#bMHqJT_jrbvMn21zuT zCKKFZS(`b9<$qPS!-B~%F`?Y5^Z0;6i16@rqA!bD)b2vbxRb&%h5kp?Egt(?pvTo!(MXu;_z}SrDO)iP@KY$4Da(zj%!eIr zXgO%G0xUOf)NsB!kGY7#8*iCgt?)SqStqD>;b*WcND3cyb6Sh7@qULZT1AZ> zZPqU)Y>up4Cy@I~A4>hUwWBctJdfHeQ2Uxd^SZO4_m4XtWu!zD&+$+JBxxeAt^4Qp z&qYT@-Iv7feVDEO$q)c}L3?Z1T~DbUicEo~5li6p3B828nZ!MH_*}5M82@wBdRaAD zk*7cTPTDq!&}5Oo4>wpl+0R}ZUS&H|Wtj+YO6rWxz6o*d+Eh6`RYoBnawcD1>+7pX z;eRDjQy*)2RHiA!6D9ToDG>nFEG45hN0@@-%HL!;R>mDO_4c`HdM`w9{dqQ<-PQN7 z5ojX@I+tEA-9IoTU_dEUXEU%QF-tPj4>-sUrqhzb<}_1+F2Ck_g?zm?Uiq~3-k7dX z+hYa>0WWjvc8RAf=Gn%U{(5+QpIoy+K=V^22T(G+F1kOpxbv_6*3_>XpeUtdoN~I2 z6*a4L^RqOiBI&PgpwQ5lUAd}oX8?9-!7k)|`Eb^@WJA1BsHDDHl!R9t5w^ylW=*vb z$bcM(SIoCsR*FwVsJZmiS0h&;E+0XXQ2H0^Gd~G{^g1tFH%9;J=U;i`mC&!eobP-# zwVdnxbQyHP3PB{_dnd#YR%HK{oCsz&NznU#zL&a+na+Y}c()aF9^=@d?P&;Bt|G(-!9Bk^U1 z+!ioVSqx+%T$3Z~&wZO*5kl$!GnS5(t@_y^XJk>ebn^yG1HO=F=0v3Mc zjXZTgoajWkPzR(V9Y;omSB1lnb^vdd;q=&&)=%CHM<0L;ULZ8GH92gay=UGXaY-Ue z0WwrpoGra;^d9Ehx@MzQq=dWyS5qW6!bWW9*wL+_L3TnX{7V{em}%LA!bDE$xL|7G zQDx&@G`@I6J%>o?!|U_roqjWU36P2%Kyn>N_}eI#3cS}nPNN|yI^vY-lBMz|0330T7C%F>sNCGsaB-I+NYE5c$U;dGw}h8jMUZ%s zk+DG~qQKY>omg6L3W)!W=wGH0j|P5HHk#Z)InL*!N#fwVpD&uk{hsF)tCElGhKm2? zITQI&Ku8aZMxk!*J1Z(~AUh{v=mn{gER;edNWp7}?lcq~6FCbQ@atM-mxp-dYZ4NL zXeoNj$Xism5l?ds+4q3EIo43LtU>l0V2KGkxD3k8_*7-7%UNO}-I21?4u@sHXzQDk zOAp#9{;OAG@z;13MBQH$t-Mn(f_5)g2o3pTQJ}w9b}{xfAD2ss+bfgZctKSVd*0u@ z(%2S?3F+P=nq*wJx!ql2J=*UGM}yaR88THhgij1x=_;bbd*4EVx>!Q`U1F#j0yNhZ z>L0M`LNDWj_u^dYd(?Ed*La0n4DZJ5Yy8z}MUYc5EU-Z5_ea_q@9bI)M%r_Pjqa*K z&{Cj^<5&4xrf{dv1&yyArHmb17kgoU+37?jlxMJ+2=-A^G_UX; z#IacI&s~1k2eSrWWu|7FCQYS9$6oMe!I1}gsU!K}5f*>vhD-QVx_iqGcqOyYHj;Y7 zn1ULSY6W=(VyL3NMQ@;SZ~yZ2H)@-1uq0cu;~7E#{X`+D06X;43t{58F{tnJ&s5t7 zSrxMgs-fmr0wv!juFbwVc>Rcif>ao@F#o>DXMSooUp4cP)icr_F`Me6 z-UaR=$Rk{0zd<AON_l!EY)FWH2= z3OJ`=0Jh=ht&^@sca|SkY6YTcWJrQvb%s|()g=6nF<~mWo*h7fIQgwiNqTTPocj%f zRBHN|fNYGN8H>KaMebJ4fZal9bzcXpM8Rh(2}QYOCGQru9SkHRVjr8;+V@zPoaFE& zD1WR*jfzTAcyJL{##QAMz}@!`)iM7F!|20MsE)=M%Edchk~v9CrSicSKs{d5^+Q0Woq$@X-nJSeRB0No9sS*sJb#15c-71OEvIz0TPvrw84v@Vpncdm!`vTrLW>4Y(PsF!~f~1)7 zig-=kov#A}r!VoRKW83(f&>pQ0zT=;_>3xdipUy?vV2Z^ef0R}g-Pur%__sqHxm%! zgY9xG9OF{d?ps_OdY|KlUl}2;OkISGrgwI`AI1HFr5e3pLL7dAPiq4A=iuuXB5ia;5UW{Z86c&2X5zG=XolI z57NK8JrsVwL_clrhoV9?hpCnL;WH+y(w!X>e_#gvOr<>^4p8)Vt~GYz>=k zKPHbBqx0{lRfR+S&~|_Nm8NhGw>ZY20ImIMbK;H}I`Q9=Kak#uC+)izN@@;Y2ujkW z-X$FpER-*JICubn$;A;1x0;7p<5)AM1jLJ<{GqvAkXUZ%6roOrqe z@RaIoUl8~*f4`Z1fWKX{d@J*i-{|m*hrr%u8-n^-KK&P!yO4m5?!&taT97l}=|($w z&k&lUA>(8#u{7@TTGa?G8n;FZr<5Wzj@5hAkPbzq=laJF7WL1c&8$;gzy=J^(FhSB zHe}3G0-1vjCm#7rC?yWHFdAPp8XD&FXg3S(I4peeF;)c*e}Z zA!!q6@{R8=Fb(A=R8*wL}h5Oc9S81q_AKmf#lh3+2?RNOQe=lu-UK${Y z6Di_P;WjhHupZ&F&A{Ry$B-J;D=GpN z7i9@|du=*u+A$5Vy)g?u2t-74*(_=i$U-a4LhVgLSN4U7`y2HnFC!quy8Y;ns^SHp zk;v76N&(4!(MVou9jLKIU6PNYUE__2npmA-m#t5J@W&Ui# zDHszvgyMxL6fyen!^s(w9OribXGnW|_{oBXJpm@^B*OXi!dgLvab7@9!b@GCM?`?YA5-S< z;U!ghCfPr}VAI>t^f+5RA!RmsX8UhdB9`arlOagNm zB!>U$v#{PEZQA>{R_ewWScxjJE&_o{o@Ofqt)i%LH&C#&V`!uVk>tle-5ju7kY?Q{ z9zIq@PlO6-P$gI@RLN=G8jmUOAyxySR`A9jN>y<3GSe*0yB!LHjP{)P<9t=zncl65 z`J*Me0+^+o{e5-88fG)hC2h;2^jmn;%U_-@b*GMk)l=SEBrEJM&Rz*NQ%tF^P%g8q zE*Y(~C;WXlv(czZlPi0exbdq~UN4)Fw#FuJgSU?lCw}2W$?FZ{EkFXWqGGwE9i7&z z?DB2}uE{br_#A~QEas*W`wlPiUg(=`?p*MRa@oBU9m#ka4_Z8epip`iJsmP-tiO^r zf?Q~1Fg}tUDJh^JowL;BtPbVQ#3}wBnj)+2TyC{atAV>aT!Z!_zrwtRTDY5ds~Xx+Rb;IKPGsp&oyF@R-Js9aRVU zhqMh0hHK|eb+V@gC_uT-Y+@dvz>V=+xcs>HDP@=?QtLF=qkKg4zSXSR8AUNM->g0i zF(98_j#c2O-bIG*2q0e??{#fZL4LZ!!#eB{s+zFwt;mB#8AP-3a29k}%}AwCku|R~ zmrN1+e&kjRQ2vEgN&w7(Xkfz7ZJM(!cuc8*vLoXfZ%{TxGpxyhD~K3PN8H^Isa46U zpsxgXC%94`(8{pV(;%|jQCS+_F(0H=2<7&NjZf|8L z{DLn5ooO~lh*|<3k7ZyrWbjVurUg-+FZi(G znUAQ|Yb3deneUO8X~guAY?9)%#nIkJIOs-DCfcwE-(j05R`GW@FKCjEVXD_zqvFs_ zH+wOts1=}=e_67O?1{h;C@^&{O7H$5cAc_!wyX9*5iYM<8P}g%UhkBTvB}nft}%;l zYx+1$&X&$_9FV%JV_NwgjW+y5+e@=d{u}$qjUUS!oud46)AmE1 zhyw6d4rJd-DEY>b406ED`Tqc4K%l>mVWC+(IhZaH5=Clgn2KtlV^N0X zM45mw%+_#6fh5>9uWPY1W#Gi_O#YdnEwt9})XP#sT@qBzVm0^j>o=^hl zYP#jk1x^FA?CzaA(;QLK0>|M4j?5H<5};d2tA;Kha2y(~qKqaTAhct(OhaH4fkFrI zWPlmbXnH=lhdr|oWNL(%5WwMmvgta5H%5NK@md&8RQeq`YAOidvPam!eP!!vH55Aw z(b5R;ieQdSta#$002M$NklKnp*042LEFJzzwI;jo6}2wdl0bF`@Q9`uRS*7p-X#7KcaP2;(89yzPwIf`^T zMCEAFu@I8f*KNi7pvk_DTwBE+qmh9b z#GdC`E6muwmxP*=%6qnec^tYH`io^t>=LKq9<`@d8I;g3#t*`9;Fc7`d9Vg>Oe26a zqGJ#=0S-VOA|8P8U7I+8En1)}b`BjOs@TUl!*=HAJg5#_l@7%zb}3xI8%+?pme1hb zOme9FE~9ca>}Zke{6$sbE`$JjGghGu?%2fp{Mw6sKc2?R68Y?(bmPx^>w5?i-t z51BpufN!cHSFsJUcMRbv% zB}xfPkjUZOT2>bizt>C>V)1Sgcm|33(zH(_-$Y6#2|Du%V8{PxZ9^uDNj!6G^9VS} z22sJ7UMDdFjbKqot+o*c-Y3p52t(kSu6sQRUop+NFz0Eer$>yDgHbR< z!Tn(gQ+z_-2qLl+4GcLFSy4_7SH#BIyGwjwy5`vupU0lJ-Vr;;7XD{J5UqI>v*0)o zqULB$1Tm{~gsGYLY4Dsb!JKj<>t+XGxlYKQIYR~=*gT)nx|XBJDt8H4FX~W$LhE7hI=Qm%Z?8eUdUjgYC4QWr(R4$Ar}$XKw&=~i`>p%Z z^H|j;gplCO@bf`qaH92Y^y*Sj1r54ptQOy1Ux&TAejJm-3=~7)b*p>C-XL(L`q}m4~W?>f=T>jLz1c}5b7z!~*?vrPnx!AL3^5ncR za+K-2WCF~F2Z1D!OOxi&pdd#Gu<#B8>M9vy=!=1CJU#;qnwHbbeXL<$P-qkpn3N$m zvd}PV(@DgdJwBS0{YpHCPQ{KRS~#ZmgHwEDQpaCow7Jg~>a5XSoN2X!(?`o?)FV z$4NwdhC{4Uc=OCI{H^oT6Zd{q5fjY6w$vJsv4)wvP zuSuZM82}Cob;yL|=>&wv-3jS7%ZzAh=Mrs!2joa&`7@ z(0)2<{=@!Um*9b1rVTy?5g$t6oW=zGjc4j?coE05ce8=2X6%3zZVs=KI%Ikl51BABhO&*$EVc;JPG z5n97LN|3x$L=#Y6=o;MTTnj_yU5j^+tJ!AMCx$u5390vxs0+E2gx!N<2{cnUd;zuf z&gChdsq@j6?0LPg0?J`*PoJYxwMR#h=fl=lmG7x|UzF$TgJj-j7Sr_!9*F{CjV=^o zi+v8i{*AAsU;50a(igw60U&DJ-#HwQz%Qni6WQ8lfCshn%L?__`YIR=n#rA z)<)v-RTRYp1N)qp@Iu51u%T0MLBN_q>~V)@(eVS6d|n5oMgzGbw^N>hqHtRVu8rU{ zg52ooj6*p_HXpOeF>1bBh+UemGg##u)(cse2hv`Mjv7g-gzRWRn;`00=4K={^@Q(* zsOdPO7zjtXPDm2EByD)pP5UWJX~0%7dtnJrB^`2x_z;>PnAnFaJK!BgE36}9FG+&+ed zcrV_~Ix3Et6a{;BEP^h$DDOpEs!Q65Xr|YQ62Jo&UZX|mEZ~H4y80YPm2}6QcUl_d zY<(%;ck8iFys^`D10EaS4|!>#!6jZk!nt#&)1wbPM92*RnSe&QT%wZr^SsIM8}El; zl%OxaMNhynYx#yM>NK6t$gtI^*7cF1h zZ-HKmEEx>$`lH+C6#(=!m+xZU_(8}5MBW+C*nEG1Q4k-Tj(lQP2PX$`5q&wid^w;4 z-kilWI3r~kj4H4s@`y6vDH=^|W5E-aI@-bGJR(+**pkuIxr5-d63Y1-gL*EmK1LpZNSkirt1bgmjK6)%s}>)XoUd>#l90gdg&o|L-fY`7O0@Zy0lFQ(&1Tlq%^3VNpM9bR> z9?i!dYdH&o56%g&F>(NBisl{wVwnRPLS^6QWDY0jION&!91gO7DeG3|m1B6X1_Ikq zpv5}@6bsTo!P_9dvW2TE#@T7u0T@m0?IWVDE#^Ac42iH$ng2$Fn+*@D250vM*s*QHJ=)B zX&Mdmpd+nCIEoxZt|M(qKpr=0I8Y&?Iw!;TVyQvY5DbJB#0m#7f2NAuYiIB%!e!oL ziN=>@=W4#!_Y-&$nVhkLK||ok5xie@!#c>d%xfwLnJSSbi-3KM56<$#MMSaew9URF zP1E{oD2MXpQnU!fEw}O+mZ>i&FB>9o1dN0J1XtmRqR_6;N~zT+q9V)i+j^2KKTpep?&!bKG`bg&cp3&r5`zXGGS(-v1vV4XE=NIZ2Ib# z{yGTCsbND5;d8_^NqiWFjf9e2e2(BK6jOuCj18ed;WfIzj~n4A>IBI1;3l47PwA3_ z#%-zl=*y7X`1CxtK?b7{0{P=hj8kyQ5ft6@`vh=u2M|^(h)e<%2m(Z|CLzk@A~sGi zaF2jdOq5s+?Fo^OC*||`t)jQ59NWa<;8pUujA3Cp)EsTn9G~YZupQM=s3)SG?;;Sx zS{S4ZSp%97!_#>|+P#`?yZr#q zhnjAroX-SYXC-vrYdRfxebHy!)sx)M+@;H zzD*`DKYn<}?fLMz(ARUE{WSCFO5F`W;`{~qtMa{#6U))YeC;3zxr3rw4v7)80)RZD z9A+QZDtmJXL(7-HfPuX;_XfZKq?|^lXXaeggF>wvSP6g#nBYd%1qNC~Qrl4oDp^Tz zELf9Ddgv-JlDMCq!?NWBwp*;RVVxQsMb&N9(9{eO+CqzZ^T>P$nQ}4_7 zdNeEPM+BbSkU@id9D65MBv4;?^2sOC!TaxpD~y@7pXE9HaC}@;wS@=kLNcDzn1{i zw_Fu^0=H5FP(+a+kdPuNGRZEShs`&z3@1L!J-3CzT`034a0x9;;#>)Ca09So7_cqioPx^<=gGhsW8nD9eKAW8J}-403E%Za z!x!Tu^K>kE=mqbukA{1EkHV0W{ZJ4;FGUl^BERtt0M1xG9k6y>ipmq};_RZVThEN7 zi|I`V-X3+fJz&7ow7FB9n3?VC3?7NfAD`b1DS9Xl6j#3a&99~B4n3RN$W{`(1H@v8 zvmlXJAq7WOHXv}+S&ku2&et)K>h*3DUEtCm-edXM3CTwXK=#hvE=zY1;p{W*$tCWn z246alqLeBG4O7ypSD+_&Y;?BdX%Abh-RZhK5{R57#4gH|6I!OpO}v{mg##3;D&*(` zz%2`z=Ce&Sk9%*3Uw!I9l0$|LN8Sy?j%X7(f-6$jAYRSWwu{`+B-e9d50^t!RJ7c% zNK&`rdDhqV7n~vf$$jNMG(>I}Tg3JuAG>m4>lKC>5M?<^02CZzG5GN*>A-=va9=K? zXoEO0!+&09@Dz;;4B{h5Bo$UrJ4eL8m+|GCqvxl0bM`*Y7ri^sLf7E~>rdcAdChVB zD*jI7SEu)LrL|^r4#H7P^W&2R!YGKwh-*lQXlo|%W6z(n`TLBRW59a=W-L}E+v2#+ zJrUwhU4azcFa)dE$aRMy6mt#8J)+FepE7tnwl_EQCPYATSQj~(Ti5+GrY$U9ypnp5 zVw?mlcaV64e&uG0T9g;j-5mRTeKK-LEo^y3H#i^GINvPvEN<46taxYgJ@pNlLiD8$ zKF>}(g-5aztq(i%9~38pb4tYXS{om1u>;kUFS4=zx{7 zjna;KQIi_RP{SHWJ_oaTtjNhuaGM-i zAZRSZp9KPDouDKgnF1gPpE`#XEUx4~@f=0HKw*R!w)m`zN`PceatXQF@)hlNHHq2JW+n(wX#EfAw!8ir4R?@~RFI zll4ic5RdbK2@oBeeqLq$Iumg?8GFoFLw?Wmt7$j$yj@i4RJY!O}w}w~+iQ=G`3*>Rl(nICRR}8b zQ;FQ8!#U~UfDCn`EC&2BFWj4OGVa1ij|fcNpNAm%6pY7#bofV478V_P zJMM{Z%16os>p618hi62l!Zvj6(18|fcS?vYd_X8FJV#%L{DU4CZ&+A}SOi{#HyRZs@ zIfDuUSUeCEdEhPdjyhJHAO+WpP=xx zpH+^gs+K=-bVaf6`&zHKciLH*$8*8drr#`c#0<>*F8Zqwsg(o>o_*^(-kUyn&qpZJ z&LFo(6uM3R9D&!pfQSH(&Tmo;XRn7Wc5>gK^$(s7MBRKrkFcqE}`g z=D`Qcz70X6>}`CP>&YnU#BX{Z*4Q8rCq1;+&2ter?QJ7gmzIlYG;sRegn8Qu-zAc8 z9KsoL$UT-cK?mol!6hVYhYkoxpQRbs*+(-0Y8cS*nYpaupXF(c(AW?|t0^NpN(}|V zk~bYXppQfx`u>*Zvh6|eQ^J-56`hX68M=;al>AbwelY70z9!BoJ_d|wGwYK@Elu^b z9d(RU_;{ur02SjI1(>pQAdvhO1LC|sNYjp8Hv^oXoyv`L71`3c+R z?`sMkKY$;ZW2QdxL&D9s2E)NsqfH%hWzhZoJ?8<(Q2fa z=h@r?LB$YlNYP;maxLd7bDKK>$|j#)cE!7JSXzI`rf_y@at_fEfJXPv-vEYIKVDjg zW9TkkCFH*wT!RCWAQgC~q75x0#67l(dbWqx)DKM=$gvx3ZuGdo1!-~8G9_~B)xc6| z5~H?>ptL;S+kb*(fQ}UwMc*B1U_`w&$~Y5|62V74`Y|deU}-n;OgY;vy?SC>d_(?d zI-y*w{`LREOr6L0h?ton;IN4CQB8i>{Ns26F;l>Zuh_yQTn~_B zamo*-ZqUSF5MmS*kfRFV25=XB7ZO7d6TpZi1QC-Yjt0c;c`ZwF00S*E(OCz>c-lO# zi9D~dVJZr;pFt@>=P@&y2y`^Y^tGDLDX-152fi?Aj72S$W2cukQC*%D1=%3MAq1WK z>NG*XDaNOA%dkky)~a{2&aiT0{Dt@7zdik zh3XAwwV4YMNk^cHut~7V5vnMzMimz;9O%e2j!hie+d4d0^bQy@Em!C;^jLp}z8;Cv zSS-b=SZV}d1=1lJKfVBPqC@uumN9U7_XYBMN$zdRq_{cDb5ObPLbO=?1o}#B;5i3pFsjX z7~ewvp7S11UdxFu4ABp5TcbdJj$+uxIz2K*j2N~|xrXXkXNlY4%*LZ-fH^ZKT!Eiv zb9s!$14T6XuXUG}1-adeMv5T)>AS1+OWYJk?r4(Ok#TT7+fu6bq-L z+>Ys029FJC$_*=#Kso6I)8Ovevm5aNUk60jVeOX@=HfOc<}qVtw3tQ2afmtaBtDn~ z9U~qUZ%mXmai;!>D+Dk!A$ntaOlQR8L#wblj$e_?DVY4I#;Mt}j+ucX(F43& z_Qhn=a@3#ud`EQ@eJ>L0R$|ZGKw?z)(5jVnn!?;2APG>ZqXT*8jX};3=9hx1S_t!aUsXq%Zu<-w@G^*o1gBLJx`g2Dky$KjIKwEY{vjvBj_x z$GpfTytk!TMdv7S$S50aJ~}=0u!^psasnah=kaO}Fc{vWbu6!dpm70W1(czM$p!f> zEG(tU;u3dc2cu<&NMx5!PKmeAETfLx)WG1em$cMnF2H^hLSu}LdXLF%9qu>mLiqF= z_TYx@XE9NN&mx=z8|!#N)R`mYmh*4!KrmorO3z!)YgU)Q6{Dt2RsUN^^AxxF7SD(Zds4wqx(U)LJ0AfDft2r!`Ct)(qW)8z2^Y zZF|fSQ#i?qEK8$;y>N`dpT4qONyCBGD~ZTap7zoXQ>fviaS{Tq0GnWs3ZQPd!F`LMN`=ryyi=W-xe;I$&#$*>T8Yoh^5>oh{We zNX2LK*q@**9HYhcFrULM3)8>%1X#R3-Y0)kXBHKh*3U2>*Xcocwzs`Ez5fFr zrII?LwLgGkr<`aveT`^O{9uGqS%RXaNp^W_XL>bg2(RRU)p9d*m>R7%@5% zsR_*R4D#XuuTHwb&Eh#^lHeUhHv*N8NwG(V8AcvG`c8Cibn_h7Oh{0aM{f_fEJL8= zo}6>Mhx!dwkNU5zt>bWKShr>JFZ3=*FqRk*e-J07Z2~7c2ii%Z&BSLrzzuk&!v_jL z7T~ZxrUFRC*6XGspXbAg8oU?)*eMlxxF%4;yzOJjW2{~QoM;!?w}%+&ybE_6cx!sk zPkfL9j6^rTsQ}Hi3*Vgwqry9-}s#4WbMJR*@Yr zc^@3Z@<*)pHYF56YBYibqUh#lme}==c+NDO0!M=K9rUp%a;e`~{sKObOqiaBi*R3F zO@QNw_RgZ3eY2vToKLQ@urP<1Hl6;-=P)f+R`S(@04SbbVdXm)$==_&etI3NKOIqkVm-TBTWzQP4Zmd@ zoh9FpewYt4PJjBRe;f@jXSwshUFn%;p20ZK34)0UI*_@M?m{%GLpZKnki{wi;>VP@ zNjEb!r1SyH4Ky>)&Kbsw#1A%;;a2~UtZIiSQ* zhkL2(gvdqRd8sIqtPpmvaYdJ64cwkLLub#c4=z z$Kx_IqVB&AR8!X7do#-eV(ZF z-RaYx{Z;Zv7WvuiTy@X$Thn#}j_-a2c>Fk9IgS>j??3b{&ehXiEL|I&URqzfkRE>& z@rAG(iE#(={`vD4^17ZBs{Wwgz_KJ1b(xUsxuzq~59#Pn^#FY%5=vO2%#W^POj@7b zJcHGboe3$p33)8nsDlRFM%~8fV&kwJwP{akD8G3HQbCvjQ0@y-U@L(Y#0k9O8QcKz zQ%B0~4_pNsk&TYQl?!SttlP%R+wWT69NnOJCfo!c6vuf==b2PDo>Oc?r^7kqN^4`v zj=IAQMVPBPa6;|k zoiiB0JMCc?K|Ar?23R9??pA9{Y3I(H(kDOlY09+?eV?s~o~BG!w7?bM@dNs?qPM2u zckbMk?tahR>Ewyi>CmCWAs+RxoaHe%d=x~0Rm3#c|5QHfaIL|d0;U$`UOzhBca$lB zV#Fio6zMY2NO5oiHY)Jh3Z=rYuHi?A3P^GnfgxCtLNFskr1*JaiNqk@$Ohqr90*+k z%c%&JuTVb39>gR;aV?^y<6c1EdsDEAx7P%;MCx=^?kj3msZ)aq7YA!x58Z=urfFIv z_ww9noo&Ma1a-3a-*BY@%zKKdihevM-pL}2sy4}CVu=#+b9>w`Fbl3K@M3s?;u@to zMKJs}z85`^H;`Ror*2VJZo`0J?%gYz&g$DsB zM~yn9*_1;`6xA}3;I4>XKJJTq>2=k?)KMTIlfbYymf*%%)6@rI7z4VE`3&HUF*i7O z<&F_&2LS9~oVSe|DT#kv%S5RHp35#eq2L_mqz!5p3d*u& z)sfVMwD2h^0zf#9ft!lPSse2l(1x-X2kg{ZwjNFIS z4J01Fd5)aUfD6qIn0+*wF#Q|qBzfXc@ZfY3$MuXqj}i1Qy5oq?usqYKml)qfO^1dhpf^f}zl zdD~6^P<>rbOcu~$%cA6^T^SfElya8sTsO?o~Eg@YjG=p?2_0AAW2 zeZri~**9Z3b`+V!3P1#@1=p0|a-K@~V50h=^X2X%zl$A`L@#k?ACYy5$^aP4&UG$8 zWhqfl?ki&M@{j^c%{oo)L& zSU4&~B#1#7e&m`o^t2LH&Iqxxc8A0wQR|n1j<8CJxg3+0jSDVNtGPSt4CWU&_WBW_ zpja1#MwK#4iK@noE|yEr>iJG*JvdqDFVv$Te0yG>A9dhPjRIhvO1Sz098NAM(8%Ly zlO}5|T-GT^8)QPa%&@TWyxiO8{F1VD$nT%)spw_TPh;e3l<6}pBfkev+XSNLE^;L@ z!FdE(2BAX+tiP|%CCa57xf*R0Jr+$K<8)Z)w*eQY`2bBj^qV)|^-k{11|7et&-PzV z#ddmFZhmMByd>~+JIm?#i5JqD6UWkXN1oyPO--VFT!S)eOq5>18R;~yeoRfxI!h&9 zW7QnWd{H)+*n|^=$cOX0O9(N%nFBwhk;D>oSj+%nd|b*MqQLu+Bh&?Ak+>v0Y0@Nc z!#9L@rzOlu;q}8r)HUT6)+(M7Wjc?5_voAz7;Bj21)9!VD;XR|f+cBw*I-~>Bv#*Z zGX?TI%6Wai$Rkd8mR!mcF86K-j36pv(acw*0IOwB%L}I~WfqU}eGBeMcM@QY&IS-U z!gK{2E@)P3oFZE z9`-Rt4VZfL#EucNj0`-*n?8P=HE2g5QlF4FajB04Jv0rQV7{MV%~DKaZl6 z+&1!aPJKtk-Z~2q>%>bTB27ifxQ}uL3jZT6P~xXv0@;)~A9eCRfF9FO#dU6=o`Q>3 z#mI@#0bn_}PoGxu?Ccjc$*(%n~3cPo7S9-0|knOB^!o{A%k=mGr%De=D704N-Lpm1jSW z_enq!l>}&Sxk-pzBM0R&i_BC+bhQI$lIJ3_K{Q&}qL<|BYfj{m#XKGn z#0TLL^$2v$1z!fWPrz8EF;RHjNlRkVMb^T z?VEuqV)7%@39~g#aSZTS4?35;iUBpI3r>=ldL;&t z$s8p5AR5lvlX{0-M9vf-0~p4wsr#rq=`>wjT!9z2Y8%#FpuVR zKwcg`VEOZy!X0LyP9@Pur0NIi67~Udh~RJ#XZDq3;^kF4c=x?RyP2a%#eF=MYv#Ym(QV5qH1uoy|LvU#keyd`=g-@h-qcdJq}JM! ztj)45Tb69&!3Jz6gJ&p0Dm5?+B!Mi$c%~-Iq-IH_Y8VKlrl?G%rcwo&OvqX@YyoG; z;FK{qB!G=Uc)>BY5td}hTDqm)UccA->-XmOKlgpFTMMzUTHEt?rT6+@a~vA6LKuT2v$s$zXtl#KdigSx_JowyqE`#nlwS zM}#E9g`VLGHUJj);R08-MeP0z-b3IKH&W5{5+8L=1AtJVV|owy&uKd<-UqswqT||+ z`o*u)smAr%_iKb=zsD;%O;i+TU}mahpZ)Bo`61;@%2@~@ z06=jJG%1PTSQVRpqIqLHOe3cE@mP~QZEa)-g% z#Rx5hde(jJB-Q{DK%{4B(p8+%mk|IenaEvdVuaUm0n9lNfC_a(X%KbLbk2_BwR)6xcx>G>);U8~6rO2!+Ofgq;Mu;k&! zn{Kh)`|hA0iN}Bkf#z?x)QjM*CobSbTYPkQ;hradAq zjqVy|+4uE2i~2&ki=0)na6Vb4=9mQ4WTkx9rY9%Z7LfuPz~g!j+x|9<*Q080SCfz( zAbM2bJi1ozAWNA_0Gw}Xk2(beX93|JuP}f{T*WyPMiY1b}mOE|T`mL7hkeH@tA?9+A{?QU0poI1Wq+iv2mbjfh%FpJ^t#n~=9V>3qqh0m zDYF}WA-5y}X$lX5-8**Jmh1OWA;)wPm8*>nc&c_if}$gS0e^W`TrRo?LLXmrHL5^A zNc1^>g<=eHUD1gC#=GvpIkX*9aSmc>Gy@ZDD#J=pQgG;+K!mc2=|%{LA77%DCMzQX zMh1f>I4`1wP;?L4I#Q1a5z^r=Uwj12APUMAUqJR%u0GmJ#0qmko*PHN;rW`sQFYN< zY1`uia9pCIKGG`uk^I%Rj2i-u8x+(dY6L~?buL6Z>+)WZPN(fAhZ13j>o`e!B&N!z z4YwrPElwmr$$MSd!3PFL7=VWC`dz!Zi8tpl9@47O2kHyGb13RxwL?N&XiqpEe_u=w zTp)P#Bk02q@{<%@zUZT$O-qqw+Lm+ZKfeF%EL*^QimoHD$Sc7JDJ7zjIa?#ON<@&J z|3;Jw0r($@F(La3jY(p$O%c>YUd__KqQa8_f6!ONjRcw?n$apQ3m*XnuyCDpvq-nm zd+0f`^l2%CALp6K(j1%jTtXA&l8rOzthbc+_`yR5?d5}qne44vdvRG1)ifUlNF~|% zh|y@A?rB=vq6Rnr>K&rrwRte!_z8Mo33Z+ed8&ReeFP(YBtDhwiO2rVhEAUhj+`2m zlVyk;KtazcPH+&Eu#q7noD2&Y5e&pAanIN4`!Ik$qZ{Bn0?jOsNTDa;$Q_;ruq2MT zNTj8tR5RXEGTuXej97*~h;B zls)sE#{p)|cE9yq_J;f4ZEMzTCj0ajy)e($AQD!THc3Wfl#2IF*o?oIwg)adc+S^+ zC3h01(B$MpUwY7{#!qk=%C|@iaV+C;s70t;L?Ax-YRxBRG0)0lSc~t%x*id+n|}G5 zBXX#7FMy;;MFC2brsWiRj#i24TBR8j2M{5<*aaZx5YcpfhDIf#$ynq>U*;|Nh-txk z#Hpg9XNMG!fm&6cwq$^i!9t;f1=uRUDQ0PZO5z!ejOS?S!~|T4g%y^?0;x;{IYlIe z`#k;lH*DbG0TOv2N;3KRuYUXYZO>bNhI)Ay>ZKLVd7^eS|Hj~HcEzRIxUAN;^SuN# z2hmCw2{8w*+kS&R_tXd+22KG1<}gdk9YrZUB-U6Wqzy@7c~gLO&DSIpN4fz3*Fhl2 z`I>U>NWXyS*p~|qtiX@SYg?j@=5~~JM7or>Ln~rcA7mVH3k-mOOfEw_g6n5B3Cb0| zECKN=;-Y%DH2NARz1^x$^eY`{M!F7GMwKFAeru#{Ou+G^%c-7mQrsund+bej-)rA^ z@Xzhkq5W3i!H0hIl9eY$2*SvFB-^dMyVv&I_9h-c8~7x(eyUO2dtHv8QE{bEJByA* z=i4P*r2y^XR2-&+^8$|mD2Gu!(?vO5$3X!tfoLeu%wx^+ z30O@uvh>_xAnA)X3Ig)!19{=NhJO{+O*ml49HN#;Dfb|ba}R>1b1iYn3fz*2TO4L? z32*e?6>IGJ9XHypefN4bx9ob_9(w#i8yLmQn|Y@k;F_jwX>OV(RiiOUEUuGM8kI#w z&y36E?Sab-9s%b}@Icm%IQ+}Xl%*48hxH5;~ZN!3mrd(obG<_U8A9EPyb9p*Jp zlV@o}0k4HwmK9JhVpa{&Nb864tVWgjHi?)fj{(gSDUq*&yo-pH8bmWW^MFr!jH1ZM z3NT5cVVQ!{a4YvfhGThT9&|8q3UNY-cP^FzcKO0dz*K-kt(%b!rML{?eE=i?R??~r z_g7hZkzqw=4K9>vZ@1fTzuUIna2r5mvKD}4sH=0+7Sb1?%he~Ccftn`fvsI)9^hR@ zrU@ZZ8QKE>TL8?$R@9^M-%H&Cml-_K>jV)d>@`O2=#>A3vk+thZH7@ch)q4Kw%Imt z)ZV_5Fjs>%fjBcO5esWn4P|-l`gL~8p4+W7R`G~G`C_|$>#;}4fjSEDid)sKo%q<( z=pK_36iGbF=;+AWXkjzVMwkrJZ;L@RFPgbGlw#{p2F z2Ut|<5WP-HDwJ_!4$vUAkzS?OHgDc;yLR7f**2Cwlj{a}p-)L9mH6dsi(@HUzIF#W z^a@rbElPQ|gnq?yh~JRUwULdjm#0uU1_ebv{fEn|o&SgLj?cb&4_tPGK_Ih8aEm~6 zHh>@y5eWn)@kT|7m%;{F6K2RpFf$An8UQHW*3sFGF0l;i&mbCcG$b+>)!^G|Xe?|FU(?p>s>GyA!DINn&o9=JlBCwhak zLDM)s2NB8qCM;-B7A5dk@X@A3pKPw*wsmcCVhbZbgv~$A^VgL6fFeczZSQ!eJ^ZE5 z!;u=t5b2Z=VScv4916?C*)F!VyMyT1=$Lg8e_93unVR*i8CZv62ZXgpZvO)`s-0)S ztwii9@(7rPldRp)Z#!QUtdj)6JYK4SOH|!-2HbBtiDT+E4W1qX3CN-p1xgzB807##@tD3CE#VB9w#9!}v z+j`Soge4{_D9;2)m)W;tO?&K zg0En|V@U?>^2LAbdR@B9Yu3iXA>WI?E>#a)cJN%P7dcYY2GLv#+fje8l06s42EktQal2!y`O z#>!K+ZcV@4eeX|M+wwISC5qsK5g(4H6|7qM8mxy+%06H7e(bSSZ!d|05%WN(+;--W z-NaXv=u6NLa9&%@KlYLFiPz8rEe#%dY5OFikIE_A-oM2j`}`3rU?@Ln}pOSexYwMSHTIa@_?F8`|W{Cb%EVLuab-`&cIb$(v>*&A$&||CCZKS&a6*`$4 z4cHT3kx?NMBVWVs5pJpDnuJK0l!*=QCvS!fZ96@1iu0OFLO+UF)@y!9eB;ITKud!s zbdC%MVJq2&xHnU_sZ%Er$tqTx8iCxa*4M@bEGVm>S9EoE+x2(cZMzI^;@X46A=M#+ z0o!i=laen}wzv)V^51D$n>;KM9~W%3 zK_?5IHdrpz77V;1q{B#_lADijN6%}ZH3AWbf3lKbu^${vMRQgHVFsm0V%z`y z;^!LgLL1_8MSGyN!6Q9H-Rp0ut(0xM@pfkCvDl%KD^-ePLqXmYvIzq~J>SVaDC+`- z*onXb@d71w5u>z^F(9c>HB{*ftmC&>5!u71LuCePZK1nd5ji5V7?Va+0W*F>2)ALe_iDB|ljZQ**{XB;>v z=YF@ksbT(BJ&hu01-2FmZH*Hb{YPPp#d*Xs0VbR;o09`8j6)$zwm$c7%5E7}7t z&|J|LT*8OrBjDnca1il_-pw1V6Qe*4*-g1k$N5al`etpmytm&vR-=p~F9ktOnvggS zYP-a;=;K121Ur4W(q+FTBxV6k7qaNsicAFprZ-9d)C*{9}L4XS6 zwcs8L5esvSsN13s2cmPT{?&@tsf*Z<Kud$?Y?{`-M(>;Nz6ZzCoRyiWABRYX zhl$q1L~3F4HYY}1x44|Wzxdowyxr(==Ze7uNEPB?*J6 zQ%^j>#TXL!7U%Q(eVeal31^$T@y3$(K+A)tIW(*5peHSN=vmiqMLYxMQf(Lxio^)$ z#x%^VeQ~)2v~Zy~gEp7A$TKY}ih~gJnkJ$ek0gPxY#n;Om=CA@Iz=8Lkjd#|$C$6e zm8b!n@K^ON)Eyds30@#-eEipCkCyF*vq{GXcgb8-uvi1dwpXv;}XyNRVH9oHBj2>ud&f`?%6dI8fB&JPF!QFiutO*BL9RL6o zM@d9MQ~<}BV=1*9J{9CcNpgb|sxE>UMbNt3znAJObh*NNu-?=w_4@H15a)0mZI1aT z1!BbeVU9(1I;H?2x=tS9if5q2VL+GqoCZS%2jX&xd!V(!1GPn?aOg)F4Z@LC zYx`}qR7L@>+ZY+)=6&R6@{xce3&qKLVRs3(v|?0Fr1FW=Dg*hdC7ATo>e6WK(vyGs z$Ch9L9F_`1H@}F99(+?-m{;)7Tfy5wjt=tfc0T5psl9;%EkRS`Dqf`rS{pp)ygXu% z6bTl)ty;U$CRmm&&jNxHcy62d*~vw0M)NvL51yz_EzOyS1LAnYN?4Y50}3fQW6%Gs z{mCDE$})&(3RR`%il+$Eb0kLA^(pRCWK(fTePoo2>eW^#4bVlq@!uuyftCl)g5VX1 zNji-%k?vS#D_5;1W&v|2EOGYK@zY~aDd)u*n(Ck*Sh9}Ks5;zylEjWv%r<`#+$GKW zq0$c@vyVUU3!E;q=4jkgq~zo}1T!hH3buk}u);EN`QQy6If&6<9DpJIsS@5vb5rM= z%w4j^#}%E?11%1oFuxPbxq*{JI!y{U9K{ToKwO#vKslCAAN$@@U?$QXqJ9#!2M9(- zF6J+BeLRe>5Rp=-D%X1|_O;*pFZR)&|2Zp<4FLr78hqz8FQqIjMYxg&oX6}|D$m*Y zc*VZ)$itkJ-bCy=YViQ`&4vj?`o!f@_Q2Vcaw(s61?pZzKoDWH?nWB8Vh*%+W-3CFPV)2Ks6Oq`l3Cz(uwD& z@v_>50YETOk?CP!Y|*=_P0cqB@7#QE^+R0F?t#>YDQ6#Soqdst$Y%inp7V#1c&apP z2lhW=)2!(k9xGW}9$^Z%{}rp(TF=&7xJ_{apdL@_t{z#MamH7lQm^zG(RX3la=Yc$ zy;h!_vHneKZDR06TbXCsHIt)?g)lgRUX;L~F;@mavKH2uD8=^;@o(eCetYHMpq)N+ z$e#Mj7i|34^R}^njitMmg?A>Hq5?NM280xGAh;fpyPH6f0I6Pz;K5XErZFf74tHl%x-*9DRJ<=F+p#|d16^l6!rdVT^ zw@%>2e09jG6un*uo%yF1+LV6;@RZv5MhymY{l5!+KY@dvwyN z3G%8k;W~2csD1fU|HG1h_DM_gp6mM7*lGfn?6~=MTXn}>mhV~(o)l9~Ibbm~=;L69 zBk@3@{5sQ-6RmQfYBE;%i_4ik(Bj}}{OBT#MqzeE)akLTi(Q>4>m$~+j11cVCr@DY zl`P-++86%N?)$|LavkYRxGy2<2{!Jinv3|mK!q|JH0Ts5k;&4rSBOiEn-FI6NkaY( zkJaHi2|@#lS58`IN6{uL7zbp=Mgf|cDOqW7&hETvy&XJN#mr4S1_GPN{5{7kY+D}D zbGl-O84zB2=_UL2w;#8}AN_&#ZQNo#tJm9__5HSSL%;Ri`F6U5u*HByLFpaUpd}>c zqB;%!PE-`5w6NXl=WhdJePc^jJ=wbea& z>*>nEbqL9eg(C}?#44R@p`#T^C0{F#+ScAatE}!uH)0}{;Ei1bOq^wT`D_)D3>SnN z`G^vDWY6P`efaN3Y~aA(+wmAm!ogAkKbee3{QX>hu4k8UY zaXF(0E`%b^xZ#@G!StlL4^gQbpJ=>dR4WTyoe;LAv2mNJOxl|6c8D--M}Byid^j_9 za^OY#;)j0Ce(pbf7<1ez4-BGNl&JjpkW#M>9@Qc7(c{TAfEx)Z&74n$w=$^E4 ziCmyi?cgZksrq-?`~U0j*p7W~A?&UYC?!7h)%hdjNCJ|SUo1s5Jz2bFXb(O;-A)wU7~ z!&+j{) z1ONHMwr1N6^r4od$7OL3v^3`t*Gc&#C$9_v)!{lylAA+3Q<{PwpbP9>-;R&Yq*cla z`aW-)wqojK&CETWo9vt4{l0zqL+`cy{ny!w?YCS1&h55#(@oYZQ#NaVn1$gZ2*^Q( zf|%#nf*7^ozyK&!hmA{3oCIB3T?{41lXlcX9h^qWiEY z4hbrk#$1iam8o_R@`@V^a>`-uCl630mhn!r#5RyiPXbPqPM@~6^eVR~R&jQ1tJ7AV zrHn+ibe1$Kh=8dACT}ckId~l@OR8igSE8^PCTLgoEVE=~#*V*yz`pU<589oyY16Jd z{T1V{g&t^W@GQLM!j=MR$cmQ4q_lcXAM0aE)>&+Wn8uiZWTDx#$TDj~#IwF*&n_FT zOxx&q$tpj3+D@K+(N28*OU$(Z3??fRZF$QqTW5t;+pW;O-d1hdVXOK#aKDfPk?T77 zp7A%bunELj%*J)h#R7LyS-DAnk~n?T4F#1M2z_GA+6YMBx3Y^UI*6TE#S{-)W&=-b z1qVqg>F~51Oi4cicgbK(@K_5<4vd&LH#5M!on(}2!^BOUKWcr>5=X*p&UoZ++DEzV)Bm>9T~$xSZVsEe#&c)3}=J#Ot)| z6OT1GL>f6riK`N%Hk>8Z;m8=G5hMR}txV3Jtfku#?O3tUg#&3*gt?)QO3cKxT>wv4JtaoLvhmOjT`Y_GvPz6_;4*-Ov9WGmLT+u;F3rhx$~6qzV3GT(&t4TM_`z(MkeTsVzlX{*o6 z@E|cVVN(DwP2KXxQ`*V-nHk)CYB1C(t9HPFFjvn`F*gLamAEzw7po$^mZlNs(4p#Z z(KLCEdUsrBJGN}G4gFhf+f6sYRK*@M*25qXmveidrNQH`E8w^+6<$7yF|E;gbfjXd zmgfkdK1?u>3EQ}Sjf_lKs46xxUbczk6d;(l4Qtm~s;doJnzthZqxhW_$THVy%uS*r z&00I0qcl|FzFAytX6*!0VGeF4({h4I%CX@Qq`a(6z$s2nr0wACN9;3S_?-RvZ+?V; z@jLCIzy4c${QF1k?p3?(z~Q5|YhAakUq@JA97B~|PVHB2(keusMsZJ)tauV2WhLg} z0g-hXFyX;LFm3~yNt+@XP!+4dkuh1@#^IJJ8yT$Gz4yPv?*I8;MbboErG8hGfL&md z*y+_*a7SMU0aps6j~@&j7=Jf=p!LDieATew;c*$UD^Pw4;>i=wAdS4AMd8kOpcX^q z^?Cf7KJPndR!}79H(*+q!Knqc_|N%T*6?50u|$7R5Y|Zk0tuOQ4WTk5PV7 z-P=ZedAP?o;-8Y%r7?-$b?@74?|uK&MoL;row8r~&EK|v@$Y{F9pibs=bfi)^-A=v z0t@7LMs0f123|SDeAJ+wI)2nnpFC+NPaL=DG2$4^knxUo$q1svS5=l*Xp!C2EICEH->e|rQS zCQYZvyE-vG4i~K0(D1ko4xF%;4*$rWczVAbe&uEKpJTT7?tPX-w8UoTqFznq=Qlv5V{tjN2VMu@fjS$%ybQgyd2|mmJu29Zmt;qzJ^O5N@bB#Q(bNpnFm zZjdsImMdbY2whs7EPQxnXR*qW?liXowTtTz^lXz@hYEj6!zb)_KlTw1S9R0gJ81*{ z0|E$b$;b*RnT+e`A(zww@6pYn%iNubZJm$HxjoR*;92~-5>T8ZwXB&k$*g<|ePjj; zR1y881gK`FnS?~)SD7M)Ko+y`#H?&$fGms70U)yELp5}t3L?{336TjANV$D{Hk0|7 za*1^{aIykcG7;mf7PFx{sNOtvD7%6O^8en|^ zwFN$;c`}v8^}Zc{wS13=WX>+49WS*Rd*#RxpD0wwsV@B7`bgD2O=^0DqskDO)00YR zW?fMSMbfExKt+-%Do$nE@Q>>R5Q=|*lO}pc9gdN~y$0VJI8YAVW6q{|U>T-ThoiZ_ z+3M?{Rsw%$(ntPm-Wb&%L_aYEVF`<~a|q~Mw9YY+GQ87)qno%079vIy-ct|3%R?_4BLaa0mv`hm19lD2v=byG$2`eAF%B}C znV0K61ee7Gok{006ij0{z{Cv@rZvF`YY;cd6RhLp0ClbCao3&_TxtF z-ub}41~^11O`{XhEErOwebHY(pq;o}#Xayk0Z$m;1srfzevWdP4nm^a1!*O?>zFL{oe0<)aENgWFwfjUgmYGxY4vT!Pwo? z4Ty3!G)OXFI1j0PeDZXtGH1t!m~p|HwsG@TmPYruSMvlGwpUIL*da{HvZft9F#zzg z=s{$BWAgKjZ@$Og^UDv|)~&l;l#@YWp$I)Of)@JY>fL%BfJdD#^ymY40E=fFU@ry7 zfM+N=gsLgX4^BpfLfWp|yw=L&Q)oX$`)BWcpY?3s zWh)7uXhr5ul!Y!lf|?f>bRqPiffFZf4(nStnY24r;FHG1JQ5uQplCxyn|#I43EgqE zf4|Paqn{(#utcu<;GZ0b7ZFY##z`%DV|)DcVf)fwe8!%6^1C)&9uXp zd*F2f9{mO)<)|zu zK`P7uIpz9633Cy#Yg+_fgJYy+Kc3@id1iBYFyIn{fG23;5g$|P3PobzIRaSY5rHQG z6z)-J;wtR&9o5o25N}+CJrA3ymT?|~RR*Z2<@H~-4@Kn$KM*Z1+0ukjv;!E=rO zfN}G$To1(HxpI9UKlvK(ffzj3_zxI2|H}1144y03_wkdj@gDg9j8&K@LI`D!00000 LNkvXXu0mjf{H{Xb literal 0 HcmV?d00001 diff --git a/docs/src/plots/P3AspectRatioPlot.jl b/docs/src/plots/P3AspectRatioPlot.jl new file mode 100644 index 000000000..9ae6090e9 --- /dev/null +++ b/docs/src/plots/P3AspectRatioPlot.jl @@ -0,0 +1,105 @@ +import CairoMakie as CMK +import CloudMicrophysics as CM +import ClimaParams as CP +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.P3Scheme as P3 + +FT = Float64 + +const PSP3 = CMP.ParametersP3 +p3 = CMP.ParametersP3(FT) + +function define_axis( + fig, + row_range, + col_range, + title, + ylabel, + yticks, + aspect; + logscale = true, +) + return CMK.Axis( + fig[row_range, col_range], + title = title, + xlabel = CMK.L"$D$ (mm)", + ylabel = ylabel, + xscale = CMK.log10, #ifelse(logscale, CMK.log10, CMK.identity), + yscale = ifelse(logscale, CMK.log10, CMK.identity), + xminorticksvisible = true, + xminorticks = CMK.IntervalsBetween(5), + yminorticksvisible = true, + yminorticks = CMK.IntervalsBetween(3), + xticks = [0.01, 0.1, 1, 10], + yticks = yticks, + aspect = aspect, + limits = ((0.02, 10.0), nothing), + ) +end + +#! format: off +function p3_aspect() + + D_range = range(3e-5, stop = 1e-2, length = Int(1e4)) + logocolors = CMK.Colors.JULIA_LOGO_COLORS + cl = [logocolors.blue, logocolors.green, logocolors.red, logocolors.purple] + lw = 3 + + fig = CMK.Figure() + + # define plot axis + #[row, column] + ax1 = define_axis(fig, 1:7, 1:9, CMK.L"ϕᵢ(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$ϕᵢ$ (-)", [0.0, 0.5, 1.0], 1.9, logscale = false) + ax3 = define_axis(fig, 1:7, 10:18, CMK.L"ϕᵢ(D) regime for $F_r = 0.4$", CMK.L"$ϕᵢ$ (-)", [0.0, 0.5, 1.0], 1.8, logscale = false) + + # Get thresholds + sol4_0 = P3.thresholds(p3, 400.0, 0.0) + sol4_5 = P3.thresholds(p3, 400.0, 0.5) + sol4_8 = P3.thresholds(p3, 400.0, 0.8) + # ϕᵢ(D) + fig1_0 = CMK.lines!(ax1, D_range * 1e3, [P3.ϕᵢ(p3, D, 0.0, sol4_0) for D in D_range], color = cl[1], linewidth = lw) + fig1_5 = CMK.lines!(ax1, D_range * 1e3, [P3.ϕᵢ(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw) + fig1_8 = CMK.lines!(ax1, D_range * 1e3, [P3.ϕᵢ(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw) + for ax in [ax1] + global d_tha = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw) + global d_cr_5 = CMK.vlines!(ax, sol4_5[1] * 1e3, linestyle = :dot, color = cl[2], linewidth = lw) + global d_cr_8 = CMK.vlines!(ax, sol4_8[1] * 1e3, linestyle = :dot, color = cl[3], linewidth = lw) + global d_gr_5 = CMK.vlines!(ax, sol4_5[2] * 1e3, linestyle = :dash, color = cl[2], linewidth = lw) + global d_gr_8 = CMK.vlines!(ax, sol4_8[2] * 1e3, linestyle = :dash, color = cl[3], linewidth = lw) + end + + # get thresholds + sol_2 = P3.thresholds(p3, 200.0, 0.4) + sol_4 = P3.thresholds(p3, 400.0, 0.4) + sol_8 = P3.thresholds(p3, 800.0, 0.4) + # ϕᵢ(D) + fig2_200 = CMK.lines!(ax3, D_range * 1e3, [P3.ϕᵢ(p3, D, 0.4, sol_2) for D in D_range], color = cl[1], linewidth = lw) + fig2_400 = CMK.lines!(ax3, D_range * 1e3, [P3.ϕᵢ(p3, D, 0.4, sol_4) for D in D_range], color = cl[2], linewidth = lw) + fig2_800 = CMK.lines!(ax3, D_range * 1e3, [P3.ϕᵢ(p3, D, 0.4, sol_8) for D in D_range], color = cl[3], linewidth = lw) + # plot verical lines + for ax in [ax3] + global d_thb = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw) + global d_cr_200 = CMK.vlines!(ax, sol_2[1] * 1e3, linestyle = :dot, color = cl[1], linewidth = lw) + global d_cr_400 = CMK.vlines!(ax, sol_4[1] * 1e3, linestyle = :dot, color = cl[2], linewidth = lw) + global d_cr_800 = CMK.vlines!(ax, sol_8[1] * 1e3, linestyle = :dot, color = cl[3], linewidth = lw) + global d_gr_200 = CMK.vlines!(ax, sol_2[2] * 1e3, linestyle = :dash, color = cl[1], linewidth = lw) + global d_gr_400 = CMK.vlines!(ax, sol_4[2] * 1e3, linestyle = :dash, color = cl[2], linewidth = lw) + global d_gr_800 = CMK.vlines!(ax, sol_8[2] * 1e3, linestyle = :dash, color = cl[3], linewidth = lw) + end + # add legend + CMK.Legend(fig[8:9, 1], [fig1_0, fig1_5, fig1_8], [CMK.L"$F_{r} = 0.0$", CMK.L"$F_{r} = 0.5$", CMK.L"$F_{r} = 0.8$"], framevisible = false) + CMK.Legend(fig[8:9, 3], [d_tha], [CMK.L"$D_{th}$"], framevisible = false) + CMK.Legend(fig[8:9, 7], [d_cr_5, d_cr_8], [CMK.L"$D_{cr}$ for $F_{r} = 0.5$", CMK.L"$D_{cr}$ for $F_{r} = 0.8$"], framevisible = false) + CMK.Legend(fig[8:9, 5], [d_gr_5, d_gr_8], [CMK.L"$D_{gr}$ for $F_{r} = 0.5$", CMK.L"$D_{gr}$ for $F_{r} = 0.8$"], framevisible = false) + + CMK.Legend(fig[8:9, 13], [fig2_200, fig2_400, fig2_800], [CMK.L"$\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) + CMK.Legend(fig[8:9, 14], [d_thb], [CMK.L"$D_{th}$"], framevisible = false) + CMK.Legend(fig[8:9, 17], [d_cr_200, d_cr_400, d_cr_800], [CMK.L"$D_{cr}$ for $\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$D_{cr}$ for $\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$D_{cr}$ for $\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) + CMK.Legend(fig[8:9, 16], [d_gr_200, d_gr_400, d_gr_800], [CMK.L"$D_{gr}$ for $\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$D_{gr}$ for $\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$D_{gr}$ for $\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) + + CMK.resize_to_layout!(fig) + CMK.save("P3Scheme_aspect_ratio.svg", fig) +end +#! format: on + +p3_aspect() diff --git a/docs/src/plots/P3SchemePlots.jl b/docs/src/plots/P3SchemePlots.jl index 8730ebd96..77eb66a05 100644 --- a/docs/src/plots/P3SchemePlots.jl +++ b/docs/src/plots/P3SchemePlots.jl @@ -9,14 +9,23 @@ FT = Float64 const PSP3 = CMP.ParametersP3 p3 = CMP.ParametersP3(FT) -function define_axis(fig, row_range, col_range, title, ylabel, yticks, aspect) +function define_axis( + fig, + row_range, + col_range, + title, + ylabel, + yticks, + aspect; + logscale = true, +) return CMK.Axis( fig[row_range, col_range], title = title, xlabel = CMK.L"$D$ (mm)", ylabel = ylabel, - xscale = CMK.log10, - yscale = CMK.log10, + xscale = CMK.log10, #ifelse(logscale, CMK.log10, CMK.identity), + yscale = ifelse(logscale, CMK.log10, CMK.identity), xminorticksvisible = true, xminorticks = CMK.IntervalsBetween(5), yminorticksvisible = true, @@ -44,8 +53,11 @@ function p3_relations_plot() ax3 = define_axis(fig, 1:7, 10:18, CMK.L"m(D) regime for $F_r = 0.95$", CMK.L"$m$ (kg)", [1e-10, 1e-8, 1e-6], 1.8) ax2 = define_axis(fig, 8:15, 1:9, CMK.L"A(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$A$ ($m^2$)", [1e-8, 1e-6, 1e-4], 1.7) ax4 = define_axis(fig, 8:15, 10:18, CMK.L"A(D) regime for $F_r = 0.95$", CMK.L"$A$ ($m^2$)", [1e-8, 1e-6, 1e-4], 1.6) + ax5 = define_axis(fig, 16:22, 1:9, CMK.L"ρ(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$ρ$ ($kg m^{-3}$)", [100, 500, 900], 2, logscale = false) + ax6 = define_axis(fig, 16:22, 10:18, CMK.L"ρ(D) regime for $F_r = 0.95$", CMK.L"$ρ$ ($kg m^{-3}$)", [100, 500, 900], 1.9, logscale = false) # Get thresholds + sol4_0 = P3.thresholds(p3, 400.0, 0.0) sol4_5 = P3.thresholds(p3, 400.0, 0.5) sol4_8 = P3.thresholds(p3, 400.0, 0.8) # m(D) @@ -56,8 +68,12 @@ function p3_relations_plot() fig2_0 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw) fig2_5 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw) fig2_8 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw) + # ρ(D) + density_0 = CMK.lines!(ax5, D_range * 1e3, [P3.p3_density(p3, D, 0.0, sol4_0) for D in D_range], color = cl[1], linewidth = lw) + density_5 = CMK.lines!(ax5, D_range * 1e3, [P3.p3_density(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw) + density_8 = CMK.lines!(ax5, D_range * 1e3, [P3.p3_density(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw) # plot verical lines - for ax in [ax1, ax2] + for ax in [ax1, ax2, ax5] global d_tha = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw) global d_cr_5 = CMK.vlines!(ax, sol4_5[1] * 1e3, linestyle = :dot, color = cl[2], linewidth = lw) global d_cr_8 = CMK.vlines!(ax, sol4_8[1] * 1e3, linestyle = :dot, color = cl[3], linewidth = lw) @@ -77,8 +93,12 @@ function p3_relations_plot() fig3_200 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_2) for D in D_range], color = cl[1], linewidth = lw) fig3_400 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw) fig3_800 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_8) for D in D_range], color = cl[3], linewidth = lw) + # ρ(D) + density_200 = CMK.lines!(ax6, D_range * 1e3, [P3.p3_density(p3, D, 0.5, sol_2) for D in D_range], color = cl[1], linewidth = lw) + density_400 = CMK.lines!(ax6, D_range * 1e3, [P3.p3_density(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw) + density_800 = CMK.lines!(ax6, D_range * 1e3, [P3.p3_density(p3, D, 0.5, sol_8) for D in D_range], color = cl[3], linewidth = lw) # plot verical lines - for ax in [ax3, ax4] + for ax in [ax3, ax4, ax6] global d_thb = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw) global d_cr_200 = CMK.vlines!(ax, sol_2[1] * 1e3, linestyle = :dot, color = cl[1], linewidth = lw) global d_cr_400 = CMK.vlines!(ax, sol_4[1] * 1e3, linestyle = :dot, color = cl[2], linewidth = lw) @@ -88,15 +108,15 @@ function p3_relations_plot() global d_gr_800 = CMK.vlines!(ax, sol_8[2] * 1e3, linestyle = :dash, color = cl[3], linewidth = lw) end # add legend - CMK.Legend(fig[16:17, 1], [fig1_0, fig1_5, fig1_8], [CMK.L"$F_{r} = 0.0$", CMK.L"$F_{r} = 0.5$", CMK.L"$F_{r} = 0.8$"], framevisible = false) - CMK.Legend(fig[16:17, 3], [d_tha], [CMK.L"$D_{th}$"], framevisible = false) - CMK.Legend(fig[16:17, 7], [d_cr_5, d_cr_8], [CMK.L"$D_{cr}$ for $F_{r} = 0.5$", CMK.L"$D_{cr}$ for $F_{r} = 0.8$"], framevisible = false) - CMK.Legend(fig[16:17, 5], [d_gr_5, d_gr_8], [CMK.L"$D_{gr}$ for $F_{r} = 0.5$", CMK.L"$D_{gr}$ for $F_{r} = 0.8$"], framevisible = false) + CMK.Legend(fig[24:25, 1], [fig1_0, fig1_5, fig1_8], [CMK.L"$F_{r} = 0.0$", CMK.L"$F_{r} = 0.5$", CMK.L"$F_{r} = 0.8$"], framevisible = false) + CMK.Legend(fig[24:25, 3], [d_tha], [CMK.L"$D_{th}$"], framevisible = false) + CMK.Legend(fig[24:25, 7], [d_cr_5, d_cr_8], [CMK.L"$D_{cr}$ for $F_{r} = 0.5$", CMK.L"$D_{cr}$ for $F_{r} = 0.8$"], framevisible = false) + CMK.Legend(fig[24:25, 5], [d_gr_5, d_gr_8], [CMK.L"$D_{gr}$ for $F_{r} = 0.5$", CMK.L"$D_{gr}$ for $F_{r} = 0.8$"], framevisible = false) - CMK.Legend(fig[16:17, 13], [fig3_200, fig3_400, fig3_800], [CMK.L"$\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) - CMK.Legend(fig[16:17, 14], [d_thb], [CMK.L"$D_{th}$"], framevisible = false) - CMK.Legend(fig[16:17, 17], [d_cr_200, d_cr_400, d_cr_800], [CMK.L"$D_{cr}$ for $\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$D_{cr}$ for $\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$D_{cr}$ for $\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) - CMK.Legend(fig[16:17, 16], [d_gr_200, d_gr_400, d_gr_800], [CMK.L"$D_{gr}$ for $\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$D_{gr}$ for $\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$D_{gr}$ for $\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) + CMK.Legend(fig[24:25, 13], [fig3_200, fig3_400, fig3_800], [CMK.L"$\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) + CMK.Legend(fig[24:25, 14], [d_thb], [CMK.L"$D_{th}$"], framevisible = false) + CMK.Legend(fig[24:25, 17], [d_cr_200, d_cr_400, d_cr_800], [CMK.L"$D_{cr}$ for $\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$D_{cr}$ for $\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$D_{cr}$ for $\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) + CMK.Legend(fig[24:25, 16], [d_gr_200, d_gr_400, d_gr_800], [CMK.L"$D_{gr}$ for $\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$D_{gr}$ for $\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$D_{gr}$ for $\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) CMK.resize_to_layout!(fig) CMK.save("P3Scheme_relations.svg", fig) diff --git a/docs/src/plots/P3TerminalVelocityPlots.jl b/docs/src/plots/P3TerminalVelocityPlots.jl index 0a2bcf3cc..db7d52cfe 100644 --- a/docs/src/plots/P3TerminalVelocityPlots.jl +++ b/docs/src/plots/P3TerminalVelocityPlots.jl @@ -3,6 +3,7 @@ import CloudMicrophysics as CM import CloudMicrophysics.P3Scheme as P3 import CloudMicrophysics.Parameters as CMP import CairoMakie as Plt +# import ColorSchemes as CLR const PSP3 = CMP.ParametersP3 @@ -24,24 +25,58 @@ function get_values( V_m = zeros(x_resolution, y_resolution) D_m = zeros(x_resolution, y_resolution) + D_m_regimes = zeros(x_resolution, y_resolution) + ϕᵢ = zeros(x_resolution, y_resolution) for i in 1:x_resolution for j in 1:y_resolution F_r = F_rs[i] ρ_r = ρ_rs[j] + aspect_ratio = true + th = P3.thresholds(p3, ρ_r, F_r) - V_m[i, j] = - P3.ice_terminal_velocity(p3, Chen2022, L, N, ρ_r, F_r, ρ_a)[2] - # get D_m in mm for plots - D_m[i, j] = 1e3 * P3.D_m(p3, L, N, ρ_r, F_r) + V_m[i, j] = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_r, + ρ_a, + aspect_ratio, + )[2] + + D_m[i, j] = P3.D_m(p3, L, N, ρ_r, F_r) + D_m_regimes[i, j] = D_m[i, j] + ϕᵢ[i, j] = P3.ϕᵢ(p3, D_m[i, j], F_r, th) + + # plot the regimes + D_th = P3.D_th_helper(p3) + if D_th > D_m[i, j] + # small spherical ice + D_m_regimes[i, j] = 0 + elseif F_r == 0 + # large nonspherical unrimed ice + D_m_regimes[i, j] = 0 + elseif th.D_gr > D_m[i, j] >= D_th + # dense nonspherical ice + D_m_regimes[i, j] = 1 + elseif th.D_cr > D_m[i, j] >= th.D_gr + # graupel + D_m_regimes[i, j] = 2 + else #elseif D >= th.D_cr + # partially rimed ice + D_m_regimes[i, j] = 3 + end end end - return (; F_rs, ρ_rs, V_m, D_m) + D_m *= 1e3 + return (; F_rs, ρ_rs, D_m_regimes, D_m, ϕᵢ, V_m) end -function make_axis_top(fig, col, title) +function make_axis(fig, row, col, title) return Plt.Axis( - fig[1, col], + fig[row, col], xlabel = "F_r", ylabel = "ρ_r", title = title, @@ -52,18 +87,7 @@ function make_axis_top(fig, col, title) yticks = [200, 400, 600, 800], ) end -function make_axis_bottom(fig, col, title) - return Plt.Axis( - fig[3, col], - height = 350, - width = 350, - xlabel = "F_r", - ylabel = "ρ_r", - title = title, - ) -end -#! format: off function figure_2() Chen2022 = CMP.Chen2022VelType(FT) # density of air in kg/m^3 @@ -77,59 +101,151 @@ function figure_2() # large D_m L_l = FT(0.7) N_l = FT(1e6) + + # D_th + D_th = P3.D_th_helper(p3) + # get V_m and D_m xres = 100 yres = 100 - (F_rs, ρ_rs, V_ms, D_ms) = get_values(p3, Chen2022.snow_ice, L_s, N_s, ρ_a, xres, yres) - (F_rm, ρ_rm, V_mm, D_mm) = get_values(p3, Chen2022.snow_ice, L_m, N_m, ρ_a, xres, yres) - (F_rl, ρ_rl, V_ml, D_ml) = get_values(p3, Chen2022.snow_ice, L_l, N_l, ρ_a, xres, yres) + + (F_rs, ρ_rs, D_m_regimes_s, D_m_s, ϕᵢ_s, V_m_s) = + get_values(p3, Chen2022.snow_ice, L_s, N_s, ρ_a, xres, yres) + (F_rm, ρ_rm, D_m_regimes_m, D_m_m, ϕᵢ_m, V_m_m) = + get_values(p3, Chen2022.snow_ice, L_m, N_m, ρ_a, xres, yres) + (F_rl, ρ_rl, D_m_regimes_l, D_m_l, ϕᵢ_l, V_m_l) = + get_values(p3, Chen2022.snow_ice, L_l, N_l, ρ_a, xres, yres) fig = Plt.Figure() # Plot velocities as in Fig 2 in Morrison and Milbrandt 2015 - args = (color = :black, labels = true, levels = 3, linewidth = 1.5, labelsize = 18) - - ax1 = make_axis_top(fig, 1, "Small Dₘ") - hm = Plt.contourf!(ax1, F_rs, ρ_rs, V_ms) - Plt.contour!(ax1, F_rs, ρ_rs, D_ms; args...) - Plt.Colorbar(fig[2, 1], hm, vertical = false) - - ax2 = make_axis_top(fig, 2, "Medium Dₘ") - hm = Plt.contourf!(ax2, F_rm, ρ_rm, V_mm) - Plt.contour!(ax2, F_rm, ρ_rm, D_mm; args...) - Plt.Colorbar(fig[2, 2], hm, vertical = false) + contour_args = ( + color = :black, + labels = true, + levels = 3, + linewidth = 1.5, + labelsize = 18, + ) - ax3 = make_axis_top(fig, 3, "Large Dₘ") - hm = Plt.contourf!(ax3, F_rl, ρ_rl, V_ml) - Plt.contour!(ax3, F_rl, ρ_rl, D_ml; args...) - Plt.Colorbar(fig[2, 3], hm, vertical = false) + ax1 = make_axis(fig, 1, 1, "Particle regimes with small Dₘ") + hm = Plt.contourf!( + ax1, + F_rs, + ρ_rs, + D_m_regimes_s, + levels = 3, + colormap = Plt.cgrad(:PuBuGn_3, 3, categorical = true), + ) + Plt.Colorbar( + fig[2, 1], + hm, + vertical = false, + label = "dense nonspherical ice (1), graupel (2), partially rimed ice (3)", + ticks = [1, 2, 3], + labelsize = 14, + ) - Plt.linkxaxes!(ax1, ax2) - Plt.linkxaxes!(ax2, ax3) - Plt.linkyaxes!(ax1, ax2) - Plt.linkyaxes!(ax2, ax3) + ax2 = make_axis(fig, 1, 2, "Particle regimes with medium Dₘ") + hm = Plt.contourf!( + ax2, + F_rm, + ρ_rm, + D_m_regimes_m, + levels = 3, + colormap = Plt.cgrad(:PuBuGn_3, 3, categorical = true), + ) + Plt.Colorbar( + fig[2, 2], + hm, + vertical = false, + label = "graupel (2), partially rimed ice (3)", + ticks = [1, 2, 3], + labelsize = 14, + ) - ## Plot D_m as second row of comparisons + ax3 = make_axis(fig, 1, 3, "Particle regimes with large Dₘ") + hm = Plt.contourf!( + ax3, + F_rl, + ρ_rl, + D_m_regimes_l, + levels = 3, + colormap = Plt.cgrad(:PuBuGn_3, 3, categorical = true), + ) + Plt.Colorbar( + fig[2, 3], + hm, + vertical = false, + label = "graupel (2), partially rimed ice (3)", + ticks = [1, 2, 3], + labelsize = 14, + ) - ax4 = make_axis_bottom(fig, 1, "Small Dₘ vs F_r and ρ_r") - hm = Plt.contourf!(ax4, F_rs, ρ_rs, D_ms) + ax4 = make_axis(fig, 3, 1, "Dₘ (mm) (L = 8e-4 kgm⁻³, N = 1e6 m⁻³)") + hm = Plt.contourf!(ax4, F_rs, ρ_rs, D_m_s) Plt.Colorbar(fig[4, 1], hm, vertical = false) - ax5 = make_axis_bottom(fig, 2, "Medium Dₘ vs F_r and ρ_r") - hm = Plt.contourf!(ax5, F_rm, ρ_rm, D_mm) + ax5 = make_axis(fig, 3, 2, "Dₘ (mm) (L = 0.22 kgm⁻³, N = 1e6 m⁻³)") + hm = Plt.contourf!(ax5, F_rm, ρ_rm, D_m_m) Plt.Colorbar(fig[4, 2], hm, vertical = false) - ax6 = make_axis_bottom(fig, 3, "Large Dₘ vs F_r and ρ_r") - hm = Plt.contourf!(ax6, F_rl, ρ_rl, D_ml) + ax6 = make_axis(fig, 3, 3, "Dₘ (mm) (L = 0.7 kgm⁻³, N = 1e6 m⁻³)") + hm = Plt.contourf!(ax6, F_rl, ρ_rl, D_m_l) Plt.Colorbar(fig[4, 3], hm, vertical = false) + ax7 = make_axis(fig, 5, 1, "ϕᵢ with small Dₘ") + hm = Plt.contourf!(ax7, F_rs, ρ_rs, ϕᵢ_s, levels = 20) + Plt.Colorbar(fig[6, 1], hm, vertical = false) + Plt.contour!(ax7, F_rs, ρ_rs, D_m_s, label = "D_m"; contour_args...) + + ax8 = make_axis(fig, 5, 2, "ϕᵢ with medium Dₘ") + hm = Plt.contourf!(ax8, F_rm, ρ_rm, ϕᵢ_m) + Plt.Colorbar(fig[6, 2], hm, vertical = false) + Plt.contour!(ax8, F_rm, ρ_rm, D_m_m, label = "D_m"; contour_args...) + + ax9 = make_axis(fig, 5, 3, "ϕᵢ with large Dₘ") + hm = Plt.contourf!(ax9, F_rl, ρ_rl, ϕᵢ_l) + Plt.Colorbar(fig[6, 3], hm, vertical = false) + Plt.contour!(ax9, F_rl, ρ_rl, D_m_l, label = "D_m"; contour_args...) + + ax10 = make_axis(fig, 7, 1, "Vₘ with small Dₘ") + hm = Plt.contourf!(ax10, F_rs, ρ_rs, V_m_s) + Plt.contour!(ax10, F_rs, ρ_rs, D_m_s; contour_args...) + Plt.Colorbar(fig[8, 1], hm, vertical = false) + + ax11 = make_axis(fig, 7, 2, "Vₘ with medium Dₘ") + hm = Plt.contourf!(ax11, F_rm, ρ_rm, V_m_m) + Plt.contour!(ax11, F_rm, ρ_rm, D_m_m; contour_args...) + Plt.Colorbar(fig[8, 2], hm, vertical = false) + + ax12 = make_axis(fig, 7, 3, "Vₘ with large Dₘ") + hm = Plt.contourf!(ax12, F_rl, ρ_rl, V_m_l) + Plt.contour!(ax12, F_rl, ρ_rl, D_m_l; contour_args...) + Plt.Colorbar(fig[8, 3], hm, vertical = false) + + Plt.linkxaxes!(ax1, ax2) + Plt.linkxaxes!(ax2, ax3) + Plt.linkyaxes!(ax1, ax2) + Plt.linkyaxes!(ax2, ax3) Plt.linkxaxes!(ax1, ax4) Plt.linkxaxes!(ax4, ax5) Plt.linkxaxes!(ax5, ax6) Plt.linkyaxes!(ax1, ax4) Plt.linkyaxes!(ax4, ax5) Plt.linkyaxes!(ax5, ax6) + Plt.linkxaxes!(ax1, ax7) + Plt.linkxaxes!(ax7, ax8) + Plt.linkxaxes!(ax5, ax9) + Plt.linkyaxes!(ax1, ax7) + Plt.linkyaxes!(ax7, ax8) + Plt.linkyaxes!(ax8, ax9) + Plt.linkxaxes!(ax1, ax10) + Plt.linkxaxes!(ax10, ax11) + Plt.linkxaxes!(ax11, ax12) + Plt.linkyaxes!(ax1, ax10) + Plt.linkyaxes!(ax10, ax11) + Plt.linkyaxes!(ax11, ax12) Plt.resize_to_layout!(fig) Plt.save("MorrisonandMilbrandtFig2.svg", fig) diff --git a/src/P3_particle_properties.jl b/src/P3_particle_properties.jl index f88d15e06..9401301a6 100644 --- a/src/P3_particle_properties.jl +++ b/src/P3_particle_properties.jl @@ -131,6 +131,38 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT} end end +""" +p3_density(p3, D, F_r, th) + +- p3 - a struct with P3 parameters +- D - maximum particle dimension [m] +- F_r - rime mass fraction (L_rim / L_ice) [-] +- th - P3 scheme thresholds() output tuple (D_cr, D_gr, ρ_g, ρ_d) + +Returns the density of a particle based on where it falls in the particle-size-based properties +regime. Following Morrison and Milbrandt (2015), the density of nonspherical particles is assumed to +be the particle mass divided by the volume of a sphere with the same D. +""" +function p3_density(p3::PSP3, D::FT, F_r::FT, th) where {FT} + D_th = D_th_helper(p3) + if D_th > D + # small spherical ice + return p3.ρ_i + elseif F_r == 0 + # large nonspherical unrimed ice + return (6 * α_va_si(p3)) / FT(π) * D^(p3.β_va - 3) + elseif th.D_gr > D >= D_th + # dense nonspherical ice + return (6 * α_va_si(p3)) / FT(π) * D^(p3.β_va - 3) + elseif th.D_cr > D >= th.D_gr + # graupel + return th.ρ_g + else #elseif D >= th.D_cr + # partially rimed ice + return (6 * α_va_si(p3)) / (FT(π) * (1 - F_r)) * D^(p3.β_va - 3) + end +end + """ mass_(p3, D, ρ, F_r) @@ -155,7 +187,7 @@ mass_r(p3::PSP3, D::FT, F_r::FT) where {FT <: Real} = - p3 - a struct with P3 scheme parameters - D - maximum particle dimension - F_r - rime mass fraction (L_rim/L_ice) - - th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d) + - th - P3 scheme thresholds() output tuple (D_cr, D_gr, ρ_g, ρ_d) Returns mass(D) regime, used to create figures for the docs page. """ @@ -175,7 +207,7 @@ function p3_mass( elseif th.D_cr > D >= th.D_gr return mass_s(D, th.ρ_g) # graupel else #elseif D >= th.D_cr - mass_r(p3, D, F_r) # partially rimed ice + return mass_r(p3, D, F_r) # partially rimed ice end end @@ -201,7 +233,7 @@ A_r(p3::PSP3, F_r::FT, D::FT) where {FT <: Real} = - p3 - a struct with P3 scheme parameters - D - maximum particle dimension - F_r - rime mass fraction (L_rim/L_ice) - - th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d) + - th - P3 scheme thresholds() output tuple (D_cr, D_gr, ρ_g, ρ_d) Returns area(D), used to create figures for the documentation. """ diff --git a/src/P3_size_distribution.jl b/src/P3_size_distribution.jl index cabc64d7b..a20417ab0 100644 --- a/src/P3_size_distribution.jl +++ b/src/P3_size_distribution.jl @@ -124,7 +124,7 @@ end - F_r - rime mass fraction [L_rim/L_ice] - log_λ - logarithm of the slope parameter of N′ gamma distribution - μ - shape parameter of N′ gamma distribution - - th - thresholds() nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d) + - th - P3 scheme thresholds() output tuple (D_cr, D_gr, ρ_g, ρ_d) Returns L/N for all values of D (sum over all regimes). Eq. 5 in Morrison and Milbrandt (2015). @@ -191,7 +191,7 @@ end - μ - shape parameter of N′ gamma distribution - F_r -rime mass fraction [L_rim/L_ice] - p3 - a struct with P3 scheme parameters - - th - thresholds() nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d) + - th - P3 scheme thresholds() output tuple (D_cr, D_gr, ρ_g, ρ_d) Returns estimated guess for λ from L to be used in distribution_parameter_solver() """ diff --git a/src/P3_terminal_velocity.jl b/src/P3_terminal_velocity.jl index a6c0729fb..0621e1121 100644 --- a/src/P3_terminal_velocity.jl +++ b/src/P3_terminal_velocity.jl @@ -1,29 +1,61 @@ """ - ice_particle_terminal_velocity(D, Chen2022, ρₐ) + ϕᵢ(p3, D, F_r, th) + + - p3 - a struct containing P3 parameters + - D - maximum dimension of ice particle [m] + - F_r - rime mass fraction (L_rim/ L_ice) [-] + - th - P3 particle properties thresholds + +Returns the aspect ratio (ϕ) for an ice particle with mass, cross-section area, +and ice density determined using the size-dependent particle property regimes +following Morrison and Milbrandt (2015). The density of nonspherical +particles is assumed to be equal to the particle mass divided by the volume of a +spherical particle with the same D_max. +""" +function ϕᵢ(p3::PSP3, D::FT, F_r::FT, th) where {FT} + + mᵢ = p3_mass(p3, D, F_r, th) + aᵢ = p3_area(p3, D, F_r, th) + ρᵢ = p3_density(p3, D, F_r, th) + + return ifelse(D == 0, FT(0), 16 * ρᵢ^2 * aᵢ^3 / (9 * FT(π) * mᵢ^2)) +end + +""" + ice_particle_terminal_velocity(p3, D, Chen2022, ρₐ, F_r, th, use_aspect_ratio) - D - maximum particle dimension - Chen2022 - a struct with terminal velocity parameters from Chen 2022 - ρₐ - air density + - F_r - rime mass fraction (L_rim/ L_ice) + - th - P3 particle properties thresholds + - use_aspect_ratio - Bool flag set to true if we want to consider the effects + of particle aspect ratio on its terminal velocity (default: true) Returns the terminal velocity of a single ice particle as a function -of its size (maximum dimension) using Chen 2022 parametrization. -Needed for numerical integrals in the P3 scheme. +of its size (maximum dimension, D) using the Chen 2022 parametrization. """ function ice_particle_terminal_velocity( + p3::PSP3, D::FT, Chen2022::CMP.Chen2022VelTypeSnowIce, ρₐ::FT, + F_r::FT, + th, + use_aspect_ratio = true, ) where {FT} if D <= Chen2022.cutoff (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρₐ) else (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρₐ) end - return sum(@. sum(ai * D^bi * exp(-ci * D))) + v = sum(@. sum(ai * D^bi * exp(-ci * D))) + + return ifelse(use_aspect_ratio, ϕᵢ(p3, D, F_r, th)^FT(1 / 3) * v, v) end """ - velocity_difference(type, Dₗ, Dᵢ, p3, Chen2022, ρ_a, F_r, th) + velocity_difference(type, Dₗ, Dᵢ, p3, Chen2022, ρ_a, F_r, th, use_aspect_ratio) - type - a struct containing the size distribution parameters of the particle colliding with ice - Dₗ - maximum dimension of the particle colliding with ice @@ -33,11 +65,12 @@ end - ρ_a - density of air - F_r - rime mass fraction (L_rim/ L_ice) - th - P3 particle properties thresholds + - use_aspect_ratio - Bool flag set to true if we want to consider the effects of + particle aspect ratio on its terminal velocity (default: true) Returns the absolute value of the velocity difference between an ice particle and cloud or rain drop as a function of their sizes. It uses Chen 2022 velocity parameterization for ice and rain and assumes no sedimentation of cloud droplets. -Needed for P3 numerical integrals """ function velocity_difference( pdf_r::Union{ @@ -48,12 +81,22 @@ function velocity_difference( Dᵢ::FT, p3::PSP3, Chen2022::CMP.Chen2022VelType, - ρ_a::FT, + ρₐ::FT, + F_r::FT, + th, + use_aspect_ratio = true, ) where {FT} # velocity difference for rain-ice collisions return abs( - ice_particle_terminal_velocity(Dᵢ, Chen2022.snow_ice, ρ_a) - - CM2.rain_particle_terminal_velocity(Dₗ, Chen2022.rain, ρ_a), + ice_particle_terminal_velocity( + p3, + Dᵢ, + Chen2022.snow_ice, + ρₐ, + F_r, + th, + use_aspect_ratio, + ) - CM2.rain_particle_terminal_velocity(Dₗ, Chen2022.rain, ρₐ), ) end function velocity_difference( @@ -62,14 +105,27 @@ function velocity_difference( Dᵢ::FT, p3::PSP3, Chen2022::CMP.Chen2022VelType, - ρ_a::FT, + ρₐ::FT, + F_r::FT, + th, + use_aspect_ratio = true, ) where {FT} # velocity difference for cloud-ice collisions - return abs(ice_particle_terminal_velocity(Dᵢ, Chen2022.snow_ice, ρ_a)) + return abs( + ice_particle_terminal_velocity( + p3, + Dᵢ, + Chen2022.snow_ice, + ρₐ, + F_r, + th, + use_aspect_ratio, + ), + ) end """ - ice_terminal_velocity_2(p3, Chen2022, q, N, ρ_r, F_r, ρ_a) + ice_terminal_velocity(p3, Chen2022, L, N, ρ_r, F_r, ρ_a, use_aspect_ratio) - p3 - a struct with P3 scheme parameters - Chen2022 - a struch with terminal velocity parameters as in Chen(2022) @@ -78,9 +134,11 @@ end - ρ_r - rime density (L_rim/B_rim) [kg/m^3] - F_r - rime mass fraction (L_rim/q_ice) - ρ_a - density of air + - use_aspect_ratio - Bool flag set to true if we want to consider the effects + of particle aspect ratio on its terminal velocity (default: true) Returns the mass and number weighted fall speeds for ice following -eq C10 of Morrison and Milbrandt (2015). +eq C10 of Morrison and Milbrandt (2015) and using Chen 2022 terminal velocity scheme. """ function ice_terminal_velocity( p3::PSP3, @@ -90,32 +148,53 @@ function ice_terminal_velocity( ρ_r::FT, F_r::FT, ρₐ::FT, + use_aspect_ratio = true, ) where {FT} + if N < eps(FT) || L < eps(FT) + return FT(0), FT(0) + else + # get the particle properties thresholds + th = thresholds(p3, ρ_r, F_r) + # get the size distribution parameters + (λ, N₀) = distribution_parameter_solver(p3, L, N, ρ_r, F_r) + # get the integral limit + D_max = get_ice_bound(p3, λ, eps(FT)) + + # ∫N(D) m(D) v(D) dD + v_m = QGK.quadgk( + D -> + N′ice(p3, D, λ, N₀) * + p3_mass(p3, D, F_r, th) * + ice_particle_terminal_velocity( + p3, + D, + Chen2022, + ρₐ, + F_r, + th, + use_aspect_ratio, + ), + FT(0), + D_max, + rtol = FT(1e-6), + )[1] - # get the particle properties thresholds - th = thresholds(p3, ρ_r, F_r) - # get the size distribution parameters - (λ, N₀) = distribution_parameter_solver(p3, L, N, ρ_r, F_r) - # get the integral limit - D_max = get_ice_bound(p3, λ, eps(FT)) - - # ∫N(D) m(D) v(D) dD - v_m = QGK.quadgk( - D -> - N′ice(p3, D, λ, N₀) * - p3_mass(p3, D, F_r, th) * - ice_particle_terminal_velocity(D, Chen2022, ρₐ), - FT(0), - D_max, - )[1] - - # ∫N(D) v(D) dD - v_n = QGK.quadgk( - D -> - N′ice(p3, D, λ, N₀) * - ice_particle_terminal_velocity(D, Chen2022, ρₐ), - FT(0), - D_max, - )[1] - return (v_n / N, v_m / L) + # ∫N(D) v(D) dD + v_n = QGK.quadgk( + D -> + N′ice(p3, D, λ, N₀) * ice_particle_terminal_velocity( + p3, + D, + Chen2022, + ρₐ, + F_r, + th, + use_aspect_ratio, + ), + FT(0), + D_max, + rtol = FT(1e-6), + )[1] + return (v_n / N, v_m / L) + end end diff --git a/test/p3_tests.jl b/test/p3_tests.jl index c4811e863..7186ea9a0 100644 --- a/test/p3_tests.jl +++ b/test/p3_tests.jl @@ -88,7 +88,7 @@ function test_p3_thresholds(FT) end end - TT.@testset "mass and area tests" begin + TT.@testset "mass, area, density and aspect ratio tests" begin # values ρ_r = FT(500) F_r = FT(0.5) @@ -115,11 +115,22 @@ function test_p3_thresholds(FT) TT.@test P3.p3_mass(p3, D_3, F_r, th) == P3.mass_s(D_3, th.ρ_g) TT.@test P3.p3_mass(p3, D_cr, F_r, th) == P3.mass_r(p3, D_cr, F_r) + # test density + TT.@test P3.p3_density(p3, D_1, F_r, th) ≈ p3.ρ_i + TT.@test P3.p3_density(p3, D_2, F_r, th) ≈ 544.916989830 + TT.@test P3.p3_density(p3, D_3, F_r, th) ≈ th.ρ_g + TT.@test P3.p3_density(p3, D_cr, F_r, th) ≈ 383.33480937 + + # test aspect ratio + TT.@test P3.ϕᵢ(p3, D_1, F_r, th) ≈ 1 + TT.@test P3.ϕᵢ(p3, D_2, F_r, th) ≈ 0.5777887690 + TT.@test P3.ϕᵢ(p3, D_3, F_r, th) ≈ 1 + TT.@test P3.ϕᵢ(p3, D_cr, F_r, th) ≈ 0.662104776 + # test F_r = 0 and D > D_th F_r = FT(0) TT.@test P3.p3_area(p3, D_2, F_r, th) == P3.A_ns(p3, D_2) TT.@test P3.p3_mass(p3, D_2, F_r, th) == P3.mass_nl(p3, D_2) - end end @@ -193,12 +204,21 @@ function test_particle_terminal_velocities(FT) F_r = FT(0.5) ρ_r = FT(500) th = P3.thresholds(p3, ρ_r, F_r) + use_aspect_ratio = false # Allow for a D falling into every regime of the P3 Scheme Ds = range(FT(0.5e-4), stop = FT(4.5e-4), length = 5) expected = [0.08109, 0.4115, 0.7912, 1.1550, 1.4871] for i in axes(Ds, 1) D = Ds[i] - vel = P3.ice_particle_terminal_velocity(D, Chen2022.snow_ice, ρ_a) + vel = P3.ice_particle_terminal_velocity( + p3, + D, + Chen2022.snow_ice, + ρ_a, + F_r, + th, + use_aspect_ratio, + ) TT.@test vel >= 0 TT.@test vel ≈ expected[i] rtol = 1e-3 end @@ -227,6 +247,19 @@ function test_bulk_terminal_velocities(FT) [3.65, 3.37, 3.04, 2.62, 2.02], [3.64, 3.37, 3.04, 2.61, 2.01], ] + reference_vals_m_ϕ = [ + [4.23, 4.65, 4.90, 4.94, 4.89], + [4.23, 4.65, 4.88, 4.84, 4.44], + [4.23, 4.65, 4.87, 4.82, 4.33], + [4.23, 4.64, 4.87, 4.81, 4.28], + ] + reference_vals_n_ϕ = [ + [2.21, 2.34, 2.38, 2.31, 2.04], + [2.21, 2.33, 2.37, 2.27, 2.04], + [2.21, 2.33, 2.35, 2.25, 1.91], + [2.21, 2.33, 2.36, 2.24, 1.9], + ] + for i in 1:length(ρ_rs) for j in 1:length(F_rs) ρ_r = ρ_rs[i] @@ -240,15 +273,35 @@ function test_bulk_terminal_velocities(FT) ρ_r, F_r, ρ_a, + false, + ) + calculated_vel_ϕ = P3.ice_terminal_velocity( + p3, + Chen2022.snow_ice, + L, + N, + ρ_r, + F_r, + ρ_a, + true, ) # number weighted TT.@test calculated_vel[1] > 0 TT.@test reference_vals_n[i][j] ≈ calculated_vel[1] atol = 0.1 + TT.@test calculated_vel_ϕ[1] > 0 + TT.@test reference_vals_n_ϕ[i][j] ≈ calculated_vel_ϕ[1] atol = + 0.1 # mass weighted TT.@test calculated_vel[2] > 0 TT.@test reference_vals_m[i][j] ≈ calculated_vel[2] atol = 0.1 + TT.@test calculated_vel_ϕ[2] > 0 + TT.@test reference_vals_m_ϕ[i][j] ≈ calculated_vel_ϕ[2] atol = + 0.1 + + TT.@test calculated_vel_ϕ[1] <= calculated_vel[1] + TT.@test calculated_vel_ϕ[2] <= calculated_vel[2] end end end @@ -461,50 +514,59 @@ function test_integrals(FT) Chen2022 = CMP.Chen2022VelType(FT) N = FT(1e8) - qs = range(0.001, stop = 0.005, length = 5) + Ls = range(0.001, stop = 0.005, length = 5) ρ_r = FT(500) F_rs = [FT(0), FT(0.5)] ρ_a = FT(1.2) + use_aspect_ratio = false tolerance = eps(FT) TT.@testset "Gamma vs Integral Comparison" begin for F_r in F_rs - for i in axes(qs, 1) - q = qs[i] + for i in axes(Ls, 1) + L = Ls[i] # Velocity comparisons vel_N, vel_m = P3.ice_terminal_velocity( p3, Chen2022.snow_ice, - q, + L, N, ρ_r, F_r, ρ_a, + use_aspect_ratio, ) - λ, N_0 = P3.distribution_parameter_solver(p3, q, N, ρ_r, F_r) + λ, N_0 = P3.distribution_parameter_solver(p3, L, N, ρ_r, F_r) th = P3.thresholds(p3, ρ_r, F_r) ice_bound = P3.get_ice_bound(p3, λ, tolerance) - vel(d) = - P3.ice_particle_terminal_velocity(d, Chen2022.snow_ice, ρ_a) + vel(d) = P3.ice_particle_terminal_velocity( + p3, + d, + Chen2022.snow_ice, + ρ_a, + F_r, + th, + use_aspect_ratio, + ) f(d) = vel(d) * P3.N′ice(p3, d, λ, N_0) - qgk_vel_N, = QGK.quadgk(d -> f(d) / N, FT(0), 2 * ice_bound) + qgk_vel_N, = QGK.quadgk(d -> f(d) / N, FT(0), ice_bound) qgk_vel_m, = QGK.quadgk( - d -> f(d) * P3.p3_mass(p3, d, F_r, th) / q, + d -> f(d) * P3.p3_mass(p3, d, F_r, th) / L, FT(0), - 2 * ice_bound, + ice_bound, ) - TT.@test vel_N ≈ qgk_vel_N rtol = 1e-7 - TT.@test vel_m ≈ qgk_vel_m rtol = 1e-7 + TT.@test vel_N ≈ qgk_vel_N rtol = 1e-5 + TT.@test vel_m ≈ qgk_vel_m rtol = 1e-5 # Dₘ comparisons - D_m = P3.D_m(p3, q, N, ρ_r, F_r) + D_m = P3.D_m(p3, L, N, ρ_r, F_r) f_d(d) = d * P3.p3_mass(p3, d, F_r, th) * P3.N′ice(p3, d, λ, N_0) - qgk_D_m, = QGK.quadgk(d -> f_d(d) / q, FT(0), 2 * ice_bound) + qgk_D_m, = QGK.quadgk(d -> f_d(d) / L, FT(0), ice_bound) TT.@test D_m ≈ qgk_D_m rtol = 1e-8 end @@ -515,10 +577,11 @@ end println("Testing Float32") test_p3_thresholds(Float32) -test_particle_terminal_velocities(Float64) +test_particle_terminal_velocities(Float32) #TODO - only works for Float64 now. We should switch the units inside the solver # from SI base to something more managable #test_p3_shape_solver(Float32) +# test_bulk_terminal_velocities(Float32) println("Testing Float64") test_p3_thresholds(Float64) diff --git a/test/performance_tests.jl b/test/performance_tests.jl index 0b9129969..226bb13c0 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -157,10 +157,17 @@ function benchmark_test(FT) ) bench_press( P3.ice_terminal_velocity, - (p3, ch2022.snow_ice, q_ice, N, ρ_r, F_r, ρ_air), - 2.1e5, - 3e4, - 2e3, + (p3, ch2022.snow_ice, q_ice, N, ρ_r, F_r, ρ_air, false), + 2.5e5, + 800, + 4, + ) + bench_press( + P3.ice_terminal_velocity, + (p3, ch2022.snow_ice, q_ice, N, ρ_r, F_r, ρ_air, true), + 2.5e5, + 800, + 4, ) bench_press(P3.D_m, (p3, q_ice, N, ρ_r, F_r), 1e5) end