From 10d974283ec7cc9d4c0bda9178e2b52e647dd27a Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Wed, 21 Aug 2024 09:10:21 -0400 Subject: [PATCH] fix time frequency impl (add/use arange instead of linspace) --- autotest/TestMathUtil.f90 | 41 +++++++++++++++++- .../test_mf6model[prtreltfreq].npy | Bin 108384 -> 108384 bytes src/Model/ParticleTracking/prt-prp.f90 | 6 +-- src/Utilities/MathUtil.f90 | 35 ++++++++++++++- 4 files changed, 76 insertions(+), 6 deletions(-) diff --git a/autotest/TestMathUtil.f90 b/autotest/TestMathUtil.f90 index 199b6d9654c..331aa0d5c61 100644 --- a/autotest/TestMathUtil.f90 +++ b/autotest/TestMathUtil.f90 @@ -5,7 +5,8 @@ module TestMathUtil to_string, unittest_type use MathUtilModule, only: f1d, is_close, mod_offset, & zero_ch, zero_br, & - get_perturbation, linspace + get_perturbation, & + arange, linspace implicit none private public :: collect_mathutil @@ -26,6 +27,7 @@ subroutine collect_mathutil(testsuite) test_zero_br), & new_unittest("get_perturbation", & test_get_perturbation), & + new_unittest("arange", test_arange), & new_unittest("linspace", test_linspace) & ] end subroutine collect_mathutil @@ -190,14 +192,17 @@ subroutine test_zero_ch(error) z = zero_ch(-1.0_DP, 1.0_DP, f, 0.001_DP) call check(error, is_close(z, 0.0_DP, atol=1d-6), & 'expected 0, got: '//to_string(z)) + if (allocated(error)) return z = zero_ch(-4.0_DP, -1.0_DP, f, 0.001_DP) call check(error, is_close(z, -pi, atol=1d-6), & 'expected -pi, got: '//to_string(z)) + if (allocated(error)) return z = zero_ch(1.0_DP, 4.0_DP, f, 0.001_DP) call check(error, is_close(z, pi, atol=1d-6), & 'expected pi, got: '//to_string(z)) + if (allocated(error)) return end subroutine test_zero_ch subroutine test_zero_br(error) @@ -211,14 +216,17 @@ subroutine test_zero_br(error) z = zero_br(-1.0_DP, 1.0_DP, f, 0.001_DP) call check(error, is_close(z, 0.0_DP, atol=1d-6), & 'expected 0, got: '//to_string(z)) + if (allocated(error)) return z = zero_br(-4.0_DP, -1.0_DP, f, 0.001_DP) call check(error, is_close(z, -pi, atol=1d-6), & 'expected -pi, got: '//to_string(z)) + if (allocated(error)) return z = zero_br(1.0_DP, 4.0_DP, f, 0.001_DP) call check(error, is_close(z, pi, atol=1d-6), & 'expected pi, got: '//to_string(z)) + if (allocated(error)) return end subroutine test_zero_br subroutine test_get_perturbation(error) @@ -232,6 +240,7 @@ subroutine test_get_perturbation(error) call check(error, & is_close(v1, v2, atol=1d-12), & 'expected '//to_string(v1)//' got: '//to_string(v2)) + if (allocated(error)) return ! test derivative calculation for sin(x) where x=1 x = 1.d0 @@ -241,6 +250,7 @@ subroutine test_get_perturbation(error) call check(error, & is_close(v1, v2, atol=1d-5), & 'expected '//to_string(v1)//' got: '//to_string(v2)) + if (allocated(error)) return ! test derivative calculation for sin(x) where x=0 x = 0.d0 @@ -250,6 +260,7 @@ subroutine test_get_perturbation(error) call check(error, & is_close(v1, v2, atol=1d-5), & 'expected '//to_string(v1)//' got: '//to_string(v2)) + if (allocated(error)) return ! test derivative calculation for x ** 2 x = 1.d6 @@ -259,9 +270,35 @@ subroutine test_get_perturbation(error) call check(error, & is_close(v1, v2, atol=1d-1), & 'expected '//to_string(v1)//' got: '//to_string(v2)) + if (allocated(error)) return end subroutine test_get_perturbation + subroutine test_arange(error) + type(error_type), allocatable, intent(out) :: error + real(DP), allocatable :: a(:) + integer(I4B) :: i + + a = arange(DZERO, 10.0_DP, DONE) + call check(error, size(a) == 10, "wrong size: "//to_string(size(a))) + if (allocated(error)) return + do i = 1, 10 + call check(error, is_close(a(i), real(i - 1, DP))) + if (allocated(error)) return + end do + + deallocate (a) + + a = arange(DZERO, DONE, 0.1_DP) + call check(error, size(a) == 10, "wrong size: "//to_string(size(a))) + if (allocated(error)) return + do i = 1, 10 + call check(error, is_close(a(i), real(i - 1, DP) / 10.0_DP)) + if (allocated(error)) return + end do + + end subroutine test_arange + subroutine test_linspace(error) type(error_type), allocatable, intent(out) :: error real(DP), allocatable :: a(:) @@ -269,8 +306,10 @@ subroutine test_linspace(error) a = linspace(DONE, 10.0_DP, 10) call check(error, size(a) == 10) + if (allocated(error)) return do i = 1, 10 call check(error, is_close(a(i), real(i, DP))) + if (allocated(error)) return end do end subroutine test_linspace diff --git a/autotest/__snapshots__/test_prt_release_timing/test_mf6model[prtreltfreq].npy b/autotest/__snapshots__/test_prt_release_timing/test_mf6model[prtreltfreq].npy index 1bb7d8859c5a45082d5e199b83d56214022863bb..ffa70cc253c8ab407987febecb7fa3062f181440 100644 GIT binary patch literal 108384 zcmc(odz2kTnZPeId5k=V5W*x(lHpE567ql$ATOYcHwXwpWg#LG2p+(YKoTPYf)8Nz z00%W9yNgjld~p;#=mC6zhysc$QQ=_2vVu9l4uK_vhhPNv-tVj5z1}*9)75uxeas)J z^z=;sruueQ{krO_uh+ioy~iGZ@>Z2KTU|6|>6xp}SUF|E!2707`OtZ1uADMsV9Ney z&zmxB#z5^Gt5%=a@{P04UAnCG8|SSIeq!Z$t6P6!*^-a8e&dWYmjws8aK-X7TYqNN z>Lsf`ysGtkE6-fAYQ^%F@2y^W=Cay%^_M^B|7zv4b%Fc>Tl%N4bzx)Q`)k$y;EI*2 zS1wupffXy4*3QF%fy0+9TXkmb`>W1ca^9J>?Wyzk+G|GbynJln|M|7`Ru%m!Tl!zC z@Ox~B<$wHp{p&~B9Qmp`-Ty?!@0Be%e)xNQ|2-SNe*G^8zMuWR`q=B2ob%I6zh8A) zmgjDxzu5Tqt+?0x;+3hVT=Q(Uc)D_ZSKTK?ut>L`uT5e*v|d%{xx@gW!xj4+t)UCS|s<&FMRXfqi^5EjlSiG zCodfR-<~_VZSI^%?)EcpUH8H_cXN9_K54`!uHE3dy=`+3h~y6KKD*~Tr%ZG6K% z#$xr_d7e`V_~qM;cSmxinl~-hwyxv)D$hMNlKb4E zJ0@Jd#d&UTa6gsjejt)tzpU0QKiG3eJ2~&eXg|txPmAPEuC*42TzQD+_PMPe`Pz90 zY)yMvp1UNHTfcp`J^fJ6?RW3owrl#oj-h`j&pkbod+=#bt*>r@JK$RT!}8oSBDunO z|6#tqRoB}8m**}GxWoN+ckumAbLr>7gxl&5YfhoB0*oi=@Nf zWad|J>vv85yRFRp3U1GL`fHcg4rYD@SK7tP`poU**}{&m45xQ^$`r_NzRTc4M`*U*%c-Y2kutdzg)1asSHW?$|wjf5dnS z<16K#2aCq9;C`guA+59at30bec0cyIwO{%BXzTs`dFK3bzRREIqVX%vbNx12|4D1V z%Cm`e!J1|5S9w;~_HW-{?N@mw<7RhS`&FLFyvDWGewAnP96AfnvGS`tt3QT*>a%dI z{HpWZ(5J#}od?|D=L0u*&r0uCaOHV$?o(!d1-CX?Ra2~2&HM_kj01NKH}fmFdY*SQ z^DDTuL8O_tjvr^{S8!$A`!PGef;+xG0DX;}U%_2fZ=f%=^DDSFzgQb+Sa0W7aAkOM z=6vrt>E{RG%JqDJonOJN&%O@3&Cajj)@BDcE;`cAui(l&-B0ZN3a*Ts{mIU+;L5zl zEp~ndSK=s+z;n|372GMcFpB=0;imViRNT;~!fky%8uDI6nclD9N_(|#qnTg9)p(%4 znfVo583(?8gqdH#)%b-`W_|@%=GC9I^DDUDy7%bX(Al8z{NPYlk8(TgRXe|;zWMKJ z6AVw;`4wCVZ|j~D#6P8rzk+*T{o_5&&adFgcHP}s-L+=&wtL{F_p45EyDIvHM^XB3 zN5;4RI{Tn}Rh|A`@wj&^H~f1v&iig37rWk(}+coa%xlE=RE_u{*2imT0Y<#=M zjZ~cUONS75-QPBMqJBV)e>uOD>`;H--pO_S(vYujU)$W+_;zV0Cynu(O__|d+!Y_+ z9$r^_P$6GZPfB)+Z%poUI*-KsP5az^IB)@1{A$0dtIclIMs{BH=Q%#^E3EHl31POm}|U>2j43e!G0TaTEvkC7&xLyJcME zj=i1E^TVsO&5ez3mvsx5&heZc_s3e8&`I+--{tH39xCKZ>Pg9NdEUQpuIKi;kI8}t zuCMak_j2F@?nmnjHm+LWxue~ywN*jCc!m3^Jom&%u3XRCEY$1yXWX~#b5D!p{<^;W z=Ry3rFK92z*LO)ISJp3Fj6e4W{X==~>5<&}>aOc<$(;W9fNSj!%X80&5bLhjqKj8dwzRREIBEN#G<7O|K`4wEPZ&z^t zr1LAd|M}bc*A78@khqwTU%{1l@6C381y|M&e#y?S;Qp!}*7rF(zk(~vCMM4C z?PDdmU%{1iZ5!?U3hudealX~gui(nMr|;VN6~ z+1Rg$%k4Jyt6byje$Z>|SHz3%;@hkA<5yhII!-p$Z2StYv>V-nX5&|IWm)Q9t^F#` zB<;uL)_&#Vt|!bf8^7W_f3>!}Y|_)#e&yqGpRo3;JgX0HEZ=PHS9vDS3+G$=Rh~&c zmwm1MD$iv7!e6ZYiuD?JPWtgH?jv~)T?Uu&6vkJ|kE;}oU!~%PJ{4{g_eH6HM|!`4 zSAXh${^DoM{0gqdBfVzkS8z3N&u}xpf?KJ7eSAkVzk;jzWEdw)?^keT`sAu{W_|@% z#-Go&^DDUWynpI+v-v@|GJo`(onOI~xS&OLeg#+N1IO;;QltEQ>OPTxVuQ8%=2b`1y|S8Z87sJxVr9Xq?upA zmHEMAMw$5)Tv?BE(pZ20r1vYhlCR+8UCsOouB;z?(9W;mN`LfAJHLV};~OW;^q;TQ z{NPX~d6_P^^DF8r=_ik|^DDR#r+T=ZU%{35x|8hu3a*UD*7A&5`4!wR9#V^FxZcjM z;L7@iqwM?&uC9m0bJF`2Txlo20yn*1b&A_v(JwrT(to=!zFj@$!rUY6&cREB%^ zJa?pgz_iP>)8)CRM{>nGa@OGvxYqu#Jok)9?gO>ugD^EgAB&HM_kj+=Ek{k+rr6&w*P$7`xRUn&$-Udui)zXeLKH`E9(>=w(~2v@2WrWK4j-tE!>yw z{0gqjYix_>r1vYhZ_4Z8;HLMhR9yNyI>K;Uxx}^mIOsX)$FJb(`sxkFewAkuU)TSt zv0o8a-|gFHE5FJ!%||iJ*st9#Y+5}HL=W-&>Bq0Q zuIg_eTsD&Nob=;YE!-*d&Bm|biWfLM$I7pKzc7q>XDh$zJU8^IaGSWV)-t_c!P7jS z&zkuaT$%qp?sYT2f~)Vb9B$@Ua3#;l@)2f!1y`O^X76O?S8z2?{&+LLf-CD6c4D46 zy9_|=eqPn?}!QD0d<_#->Nf-Ct7X76Jbe+5_GGquIeui#4D$O=2Z zf-8A}&$aWb7Vgz{eg#+3h=0S*ui#2t?m2dT)hX^U8StyTWjt!=SN_Kvh5j9RruitV zPInZ?=frYDzO2){?T?v0k3QwYYTM_=#GpBb$AF+vkqqz`le>$!?8v8{s5x!M1<*-1zvb@bewuzy(~*_q?6+Kjx&~+l=$u z<ieGd46Y|W=knZRBDt~-cmCeKzJ0;{RGxcm zB=^wz1n_6(1not|{~hJIvGMJak8#aDp40E-{hGbBm*u%9Mb=lIL&pfte1S>eop$sM*kNdS0CQLvGMKV8Lv~{s%!24%hz{lP~YKxn{_Pg?ld=SAHr?* zsb+k;m0!WtcH;>%zk(~*)s`2%XQcNlxb?S}Z2ib9W_|@%;u)s)nE4f4^*Z|9O3+@G zzmIm~M)4)rGbLR0xSfOkAe~>qmAv8~BK|79U%^$MC?QEwWzk;Lf#?!`rm1lBY?fHtaU*(y+ z9&YQ`jr}UmWWDn*-eK%lc_#67Utm1L%CB<0p8jO*S9vCJBmIm=TKQF;Nxp)cto_Qz z^FM0sSALzy@l&0)YrVnvS$S_ihz8P^w3C#ySDLTh+OP6VUN1RvFSGG0>MQR-_=>e( z<(Z5dZL;>OJd=Fa|7Go0KJNN5Yro1f$sh9zYro31-#yLRuR71AeW3XWw~70rE)$zLBg88a+<5zGc-0=T~qg-uqcQzk)0It*6d4i@$;^arO7w`BkU5#C73Sl=xJ*&G%c&C%ktI`gWxF zIW(3V@-;5^(T3;f-%)Lw+ocC+{0ong9i8UJ#d<7M}>SzJt^7U zY3?o@xPYtgsbTzzxS;sB`tZ8yi{xq^s?Kq%kz9G+-#p&wk28)FANLj3cQBGGaZoRH zjyp4wE8{sccK7vVyo|E^^*rls;?9fY>U$n}9}u4_pL6;8&VQS@heUGaIrL%HOLE_G zAD6H1p>GrSU6EXQ@54go2WgjSr_0y(xVMQL8{aPZ>hH!Ut4{yi`1tm)y;>Yu-%fL5 z#=8yw({38_bvBv`TMuXuTpV~ z{0go-|4wH7iguZHx_o_${0gp&tL)x6ZndD@fGhDTcd>quewTi@e0_`j3a-AFr*qsQ zzk)08J^H+zU!~#}`4wD!&!e4RrQ#O(6R-0=t5n<~ze>doeJb2mZZp2!)~|eA?xbgp{i^faIKRp>xt@R6W$ag- z=f?R}p6Pf_=R7CQukuXCb2`V3^Q%0Q^;TE0e$dLV@=WtGbdDS6S9vD+3SP7JtIl)d z{L1H1pFYdqzj!mBb9vt=8o#1_mi+A<aXi`<5yhIo#w{*Rp;v)`c$~x68ANSaisL)S8!!L&S89>lzh(R+pDfh z>lgK zo#Phy6fxvGjfgSDt?_c*%JDYDk~EJIF2aE9xuW{tWNOsrY`0ezAOiRODBw zxJ7;iSI2YcZ`1n~TpiEp9Jk616x>JLPjKb^bf4$*Oz&54<-I&>JI5{ZE4cD}^ea2R zO2sYmE4cDLpugJrRVr?gU%{1i);HSuRVr?gU%{34!uK#=#{E?B?R5F}xhVcB6}QN* z;L7v<)pmZBid*DYo#Kvoi}s)Hf&PH)u*CDU9CW(B__~VahJ5n3D`!(Cc>`~*w$1I* z12z7IN68NIw?Fp!CAC!0Z+O3}I^AEszD0fom;CL1+sS(YZ`rMFZjoQXC4aka_c)1r zyM*?me0{g)z`o@CQnH);?YF)3M(vI1+c^Dx``n`OE4Z3(6Tj+pk{|hF#3z-n?^q7( zOX^9-AK1`gfH7e2e@FF8SM)Q%UYuaLM1U-1L40m;CLm?PEH>f=m8(t#5k2f=mAP z*8VS@U%^$M}_;y>r$}{q} z>*s6bSIzwGwtm&j->&tw@~dY4cI8_6RWpCPtzR|sw`+Z^{HmG1-PW&~`P;R=R({pY z->zIMziQ@hxAm($BY(TCUp4c$>v^{Ft2`rryRBa}^S5h#t^BH)zunfan)%yp{i>P2 zUC*&C5yaQOp`^@Fo}Y7q{T+UIsP zYNG#vN6BtkN5AX^&#gFp4{_D${$igD>)XwNeF=|}-O?|t^W#>_M>_p})tBYDvGMII zWMX40&*^azC$up>u0E`9Y<#=48^0kQw%2LD*gkg$JzYTO|@~+y!d$j@bld@ zk}LT(FPr8hU;d+<@A!DLkQ*D{F5{DvX88K{2mJD%@9vTHm3HH~y@K-_A6FmNcdtmU zoad#yr)I#(I*JKgU*+pNJCZB&3$yld^8SZ|xsK!G`NR4yjO5BV@QwR9$(KKY_JVe# ze0^i%+aFuUp{sKBnD@k8cmpb8LKjr@3MKA8xDHHRIc@{0gqdKR4fxliIJ~>Uhd#Grxi> z^BQxzbll3oui)0Tb?L%xf8V6{E4VW5eI4(qN$*#1Wu5g|te;EoS8z2i!~`?Hf~(in zB>(xQ_ba%P_j041U%{36&@1iy3U2MKcx&#QINQvx;7VSI>+Jjrt~{44vGXgqGT#2M zonOI~{Lc^B`Be+|WjnuuD{>exDwC*dpo~M#SMKb+*WQgzTMWZ@=WHBZXvJf zrp(5#@@#1LI_Xtozw&Xp>${EpD$i~_tZtOUjQz^TZ7v<=ZS9zx6qSk(uXOib+(H_D05$C=9^Ka4k744PGhrVd-S9vD+%}%iPE1x(10c*d? zGl`%0g0)}enas2I%rzUo;<}Rjgd^vhjbF8Jr_A^5Wcu+dxRv^sv-T_BFAQ5~Hhx8Y z#VFb+*eh}_g$JHV?udwqgxH3=nrk!8GmG#~C+4&V*S&y^7onOI~eC7-6 z{Hle!*v_xuO1a9;ui(nO##}qUO2rL*D%@7DYsRuc)uawcGg>Tp1^vy^r_F^zm1453XO(TkQM_uEdS3u=6Xp zlF$5HJHKk-UTx=BaHYThhMix*mAKq#%v36GN9>hIe--R^bosV_u6t?IYSbBF7J8vnwhWQWEtj&Zt< zAJ3{f-Cv#?8{aPLo+cA_JX-7f1aY$Exv}x>@;-w7s9&Gc?e@7l>WOOn%lV~bhra)J z52xGxK0b(e+1I7dgA&nU*Fl0Txl1tTHqwl%9C*8X?YdqubFKX!{Zn~vY<#=KKhJmi`8N7{##zd9W8>R9%?;bba9cg98Q*T@S8(q; zy1w%1VKcvitMMu?n)wx6t#6lJSC#aB1y|x2f4hU{ruQqj8Yg4tS8ye7&u03c^nL|b z;`+W+_2)UgU%}P)J=yscT+L@b!_2SX>iSnZzk)0ASGU^v6r^DDTruKY%4RKps`=&nZ9?}+ONoGZtYijruhlx`1728{EF*J;;){z z_NzRT@!uz`{VLD?d|ZvZ+1juCx{33x{VLBSFXFz|ewAmEUga;=ewFL{ez$gJ<5%2A z@*KL%+OPb&whnUfnTr1kw~70r)cU9QD|p{uSzqb(jG14-)i|iv%=`+j#6Pbc?)!!G zeg#)P{*O2FB&41ngsbc5`^@|buDmzxD)KR<_ba$E{(QEbU%{2<{Zps=>ni>HAY9F7 zZs%8UPml>hJHLV}@mFK_F`FNREAM&SX`b2qAY4hCH+iAi{2*M(_jQwSNp^k(SK=nFu=A@_+|Z}O?N*<{t0?IwQjK2?Ws)yDCBNt~c76p{;#3c}^DDR# zXLgdEU%{2}*hlRA3hoyVsYNtgZ|7HV<-Gt$+4&V*iIX|N&adD~JNXqmzv>kC9WqeS zd=F0OQ{lG%@kXI>I?rWV(oYooRV+8;%Q(uw_Zxoo}>d6_Csw|n$?bWQu**!Xs>-*%qc>*PEi(LOgezFprFw~ObF zc5>c#XrDVmPh8_)&Mzf9G~Swc+&-u43aU=`m%pAza$sM=qhz;?d!I7R>2|-9=h73a ze!D!Ecv(H(Ti@3*UU$qa&mnJn``p<0c8#}g#zE^1S~~<0{JIe3$1=phCW+o|NpC_~-ruoy1c#uBZ68`mlYT%7F{G66dzX z1@{kejpggRcO+NVhb}+ZNk80ZFXH3s!}`X?x9fQOA?`>C=pS_Y=Mq;@zP_>X?UGOC z^h2Grmmh%}AKxC|iYp+*c76p{`oDh~SJD1A zonOI~_~+y7{0gp&6YhMFnP0(`aqmyt`4wD!FT9;!!Ik%UU1aB1a3z22=)=tX3a*UD zuCep0RNT;~!foX?HZ zf8WHl)_#>|`hH1kzw+xfK5Ok)o#%!=6>byvMQP<%@H7tUDKo!++>CEeJ$?mO=J&@pPFV+&Nueb9nxEjYe-)w#muEwj_`4wF8^xN$G3a%_8IMU9q;L1GRPwe~(u8f=g z$doeJb46=Og0ThkV>s`uHoj(q66GXy#XN z^}RfQGakPh%48h)`Vs#5*J^%nDAVt z^QxU+!IgZfPucku+{Y)?BEGwuaptY@Q02#SisG-}Y92W|zk(~{vFF+O6EG{nPu;Gfru%Z{ z>;2>Q^gVaZ_x$>t)BU->{`K!{`0PdJp0}_^w5G@U(VfdyUAc1f%$7xCM*q{b%T|tV zYZ*Q3OEX4~ZEGoiuxjgQ5pReMhCVi=0Xa7NJbSP3=SJ`%1+3c&X zS-E=Uk`-ULW@Tr&9cH$CX36qZ%gUdx`r?vnmzC>d+9yqFE4RxHEr(31Zx zlK+#sTmI*r+3Nu{#lA|Kw#Q2TFSL~XyZG!o?|L(41^v@q#^UEI${uApq)vP<#Q}_FC-2KSB`v>a&ThDxP zUH`wZZvUEfr+Di29sjdUZ`?gxPuxAC*G=DfnROd$);-lzch9KF4fie{tJ`fuzs9-^ zHS3<~sjE)8{K2#5Phg#tb?a-^JjyRKHahoHxb7D`b+gNA&5Aj!+h5yuub2KKT=xo3-I3+q;`CcjXWb@!^q+oo?Wspg ze;Ka3#8WrBefumugLRwrIeqHJePe)(58=8?J$2_?@i$XCD7r1WdOQr*z0yRdLL z?Jeo*@jqO*)6qSG>ph*%i`JEKUS!-leQvm4Gz!aCMm<}mjL8kpfceU(v&y6W)-&G( z^OaG=S9M^%GK%=>STJ81MSRr+<}0HtUyUv~<3V@xl~IJ z`N}BES5Kq)$|yU|zk=o~qby%7M)Q?Xmam4))-+C__{u2DSKpe$x+uOf%JS7cXudM` z_U-dk;T#lS8D;ruh;R;yuMCa5H;68ZuZq`odCINReG@ZZ#bCaQslTd!3Cve9^VPdx zzKWTz>cM;!tLeAQF8@+#1oKtQd^HHnS26L`-DtjwiLVx)1omIW#8>kt@^~cvlzOg* z{~q-7U&X{%+t7R!6JH(Q&i0}FSKMEXn8mu%zS7>|?d<2jiixkjhvuu8_^JcVS26L` z);ZkHD87n`uXdvODki?#jOMGD_^LxV2gO%0jk}vf7sXe_>$*JU)+75$n6DBrUnR^} zuYvh0VZPcA=BtGHsu!5A66UMpz4r*h(C`%TsQhzg{UQ{_1gcm3@^oZRomx_SJRuDb9OG zCH{-_YS;CRZ{Iw3+uI*Kvsaa)`^u4{e`fubn(g(CZ?|FBk3`yrU7s#-*UdHSo=6Q= zd6)K;QZ~?fg}!>n@9rvFzK8Agk8gL+uWx*NIhejY@DE`7Qne}zK>L+{e7mcAk{HCL z*eRuKe*Ve#FL-%mlg`37ckP$o54Gz~l^u7Xo2??4JFTSue(=BF{NuUraXb6ReYy6Y z=Bb;7$JUKINta`B%K^C`v#x(!y{qfTS5fqt>^9myPV-@ETCUQI>*4*z*`D^6qpn{+ zfBVVWE}PgYdQr`~zVYpL-NKDiSf@cheeTk72%5BaxV;y8+B>?uY@&JHbk=RuH`rgM za(#vCUgW8}HVY`cZ6@pX*YA~A1wH((+)v@U7kld3>$!IaUC+Oe`?hx7D?D`{&zAo@ zCH&mv^q1lGF7ed0@padWf3n#bAHsE)dg^AYyYAR3+}@(A$HQ>lD?N3yZsWRkt?BC} zE)Dlv;kv%@?M3UlM)4J|Uw8=3 zS4LUB+K1*Vqb%nvMDtbgx-L(-^~k>5_;wUu#fpB1y#nT|nEC2GFki*YSM^}N;&lsy zjsx>m%zPyrhvKW4`D!qjuVUh>X{})YRZM&}ZXDQu6%$`=NAp!od^KYV*nbrhUwwk+ ztC;xe0Gh92;;TNgUJ~WM;_LZVG+)KUR}0X5#p@S7MDrD|KN={UgYsX+#8+!Yx5PXv z{*`dQil6_gcwLvL+&b}<8{dxNtAzP#FPN_q=Boo>zDk&{dV%?h*DVa{2j;7U`Dzf~ zKPbLRn6HL_`6^-m6`HRS_Ftj-Dq;WCcs{=<|5ZYKbsd_o65=a0mF-2%4<^J{GbJAd zYJM;wzM7BbE54q0qWLNzzPcUFSG;~ zE>F4j$i7nM8;Ng6@m0#>?^|HLN|~=d2J=x%Vt^ z{8dVP^%R<~c>TilXujh0N1qqYLGe{eeD#3nmUw*VZoWEv-JX`0qMGB|DGvH+UPmZ) zckQq9)phwUi&7o-^_|LQZ(qHOIjeSEzxb;tvT1skzrp?VtE}6?=Q-To6S&=_5Id!m z&+9MFGG9illd3E__UC8$$6vYaJV18LrRYg1n<>9&54|PJM*Jh;>DqN$WXE0Tes@Kd zIKEM5VU#QSz0d6|^M&E<>>Jl@#mS^2eJ57pc(`ok>D z{iS|f57#|~cO-?JS1ILl%5T=H>A3_wS`y+(n)&9jtk1_GNrJXYVhwJ*y4@ObZy6%1I)~Pv!8x21f2gO%L+570> z-C(|wIOMm%d}S2zRh_1GNAZ>9Z|~{!7jhqmU(YRMK;>PoD=BG%QN&lrf%(cP;;SKG zzA}pVYMADR0r8cg{2S+uaK`)Y<}0J@xLvjWXmJTAP1<}0HtU%iOtE7tAN4(2Oa_avNy;wz&p=j;|;6kipuTRq<8el70* zxya+&vHq)A(eJTW!FNaagG4xt@9vVs^P~J% zG4a)K@#moYSFAf5%~!ncX`yfqimzhgtD!RQjN+@}b^Z7%k8j8FRl&D6;}Gm5V`KmU0n!Tzg+{a2Fz z3gy4zbrW0Ae3cMi-HYZcZs!qG!SjO&@zp*wUnRs>OVNCl5MM1p^A+pfgyt)GUJIJ9 z65^}N(0oi17%uac(iUB~D6>bl26`52#8Z~ZY#g#0%StF2wvH@^Mm3w|$WXE0T zQoiS7G=ELo=hOEd%=Peg9?mK!v$(>j%7*-rns zFSng}2NiW`Unyn7)eo)S&~)?Yie9yTcstLQ9djvqQp#q^S3iMu`#bF&t~*C|+=cE> z9?2Y|6P@cx&U3i#d{5o-Ljk@qp#5aFw@KS}ukXqAaNTn~b#~kTHW#hy_785Ibo2OjEMFOV{?V=%!F*+C z-QHuuHz>X`iume1Fkcx(eANKvE2D_7WZaFM{xbYN8ZI{qmvTKzDc3Xc)d|jc(A|7x z6!DeBU!nNQ@ckpt50i1XyZOrS^E=UeW%#-9V>F%bMu+dhudCU7JxL+gla%ziL`N}BESC6Cl$|%cM+t7T)ao3j#Z=m?fD9cwH(R{^m_3h%H zMDbPey4C$d*YWLGzKWTz_JH{+X1S!=u#cF;3@)P0?h(80x zS1~;w_OtV7HjuPPmw{UeWW$MRL8hTCbn6wFr%^VI?H{9wX-^)Z;Q66UKuV7^M2 zua47vy*d45`1m`#`ocZL_?s|a4FvNQxARkIzDoH18Sczib+`X2VgD7HuM+lON!~4# z|BCZ$H&5d08Rfr9h_9|d^HoB8wI9t_3Gvn0Xue8_uck@78ESrz z;Kds$^HpCkU!}}f{dhb=`L9yuD>q)GyXObFofiyedr|%?{(bTsny*s!UyTL(uTu73 zq4|pQY%fOhRZ4vILo{Ed#8>U|yhZnZ3h%G{;;&NTt6R`~l@edALGu;IT@MxBK+O-P z#8->Ze8rrxU+hKA4<5d*#B~W*Nh!myTbIwPOUZ}T%Hu@W_mQuzt50#cf3EnE=r~F> z>(FY92=FfGmy{(?QJARiXL@rKu zyr=HkJ)Q_`K`sr&7JTwh)yBF_Vo^D5^#+}`$2N%wS5-7KB? zP0z}DNx5(3J`T6{j894Tv!1#(FLB3IO<$LOS^DX4dl!64y1w!4HlOMa;mH!)E&i5p zdoT60w`g78`1Yc8-Er2fQ*-3;?NmqcmBi)BxqzGK4DTE z>bfGnlJ$c!?#ehEZm%C-8AW_moUR{V8G1h5y=cBFUe}MWjIw<77c^fLuj|KGMp?c( zLFNbLzLonp{Cf7|tKxP2_==zVupZ4<#q0X<6^{!&WL~nwetcEDZZ*$z9p4_=mvB}q_aou`#-VhpKTr7gZ=NTC^h-7C`teoq zx_*2md3ujv`%rvUysjT#aeVAT8E;X1l`vlwr|ZX8iF!IqYi&##QTG7g8|M}GdR zL|ySWQ`%Lat{-2?eT3$#;&uJ_Dq;T>ny-r2_2a9A`06`ozA9eVkFWT-4+qhFRlKer zU&*+D=Bwg${rIYQ-D;jW6kmxyP)ZbErOa1(o`lH0N<2P=k4JudRlKerU!{~E_@*P- zJ``USuj|KG93Q*z1hyB&S1I#Vak_qdl@ed=NAp!md{v;XA7Al()Ha6Cukfst;rFkf z|EhRhKfX%Ye>IWq6>b)e4!76Oe^tD$A77=!SI?sPs(4*LzT)RTbezKHSMFo!r^C;$ zU;I_^x_*4c zw(UE%4zF4F7~YW-VyBeyx#VxBdV{w4H@_nNN4ULxWyf5Ko|Lj#^0#k$`=e*B+tXXq z=e6tl%@0PA!m*Y)G8D3bi`)#HIb*SYQN$5&A#`P->^O1gSH47c}E&vuqP z>eb_Y?%$Dq=1`?;uODAUQPH|`d}6j+r_T-Zg<-Bj@l~F`osI*=S9$*S>g%bb={Ul# zD?h%<^S6^u1mdeae>>@-_$treUfn--H(wdf&mjIQ6kp}}+pEXF?&d35uOaP<;;TG= zJ6%^OzRL5rWBDr2-;U+0Jbyd2Gm5YB{Owr2%Ja9Ay(qrQ^S5L9D$n1J<*Ph@JGC>4 zukt+VSiUM=*X1d|94NlZ^S5L9D$n0e_M-SI&)-hED89<`w`2J# z&)-h=qWCJ$-;U+0Jbydci{h(1e>>@-_$trej^(R7e>;}1^8D@8&M3ah^S5L9D$n0e z_M-SI&)<&at2}=@map>s?bOaFzRL5cWBIChU6-faI<;$_za2Y2nCEY&<3RCMp1&Q- zS9$(+vKPfydH#0NMe$Xhza7h0dH!~?7sXe3{&p;1<@wvmUKC&D`P)es#aDU$b}V1z z`P;F4mFI7#c1H14p1&Q-S9$(+vKPfydH!}RU*-ARv3!;1Z>M%f@l~Ej9m`k6>$*JU z)~Q`(-ZKUJuk!ruavas|9X?N2M?0##i@j2cPI;AA`_Iey+p&C==Wpi&1kDfT`P)es z#aDU$b}V1z`P<1}6kp}}+p&C==Wi!_QGAu>Zzo+8U*-ARv3!;1Z^!agp1+;i8O2w5 z{&p;1<@wvmUKC&D`P;F4mFI8A@>QO{o!S}2S9u=w>hWGUBz)ZWi@!R2-6N@?Dsk-9 zTyd!NgZxkB;wt-kRZ-chX-n7k`s%vJ^Yz!yFF!G6C&z!K%BI~OQM;}mUqw-`o4#}U zt;^qF-AHGTAzt_IrJC*SNexwb7dxet4fPlI9sjdU?6+DW>j+DlwuiT~Z+!bT{j!w= zJy@qf@6DS1b^o|}_x$?Cw|{pXcZ!fPND5t(`V|$w&{qXZU%F|xkZ|t4qw6A|$ zy=(6zPhH#2o${QT7H!v243q0C+}_Eax?kSVRDbb=soFmOVUAqK{_*^-y&ayqzj=4e z;@dY*(>9;$FzGL(9|^bDH@@BG_r7a}(?9yhx4ZWG#m@fn($r=7{c zJlomEdq3aqjDHayN4TzUe7oiKO*1v!$I@^4^Of7qetgC22ny77`+v7iry-AT$MTg? z_C9)eH<+)CqIHXJgZauR;;TCHTfz9sDB`Q0eBYq>$|&M1c}@+AuQ(s_Mp-|H;wz(w zuZDs7$|%cM=Z)a=i{dMz?6_ZllEXLM%~wX*et8R;uZ*&M)ixQ-SKQ89(0pZ-<*OBF zzT$D=B{W|d%AfNhny*;5M?2rYD8Ax&??E%ce8utnyU~1AyspbrZk_llk8j8FRjgiI zmwAp}1@l$ReDywc_e$zT$ky8)ckF@m0)xbt0IrV&W@< z=Bt?a>Si=w#l%-5#`ATB@?XWoSIg0S#qInGny+HwtF36his^TEX*$?{#reI5w}bsx zl2;tfR~+xX5Y1N{&p)&S?7u2r*X1dL8e}66UL3V7^M2 zuf*Ss@?Ry)SI2|-it{0F9Ll;VzDm@i^D<9T3z)AG_Fs);dr|(Yg#B0J!2YX*{a0wd zN{FxSMe`N6^N6W@eo^y-3GvlFG+!meS4+`+#reILp!rJjilh08;;&NXt53jul`>y7fcYwAzB-EgXB1zh%vTbx zf%0E*KIDxOe}v+zl=vz}^Hs|JD>PrF?7u?uRm%P=c}@+=f0Ytn?M3qyw{s_&uTtWx z38#SLuTtWx&1k;j{N8t=`AYJNqxp*Cy`M+(702^$Li5$(>oyQD)g0e$)3uNL{CcjeDmCFMcdvAU#;CUYH~wc?Ycc> z$6Sh@l(M-T-2TShk4XGrgSOS3B~9DI+qqtL%%$i_DVr(pMqf?!MtxzH@cW7q*TZ#> zpdD4-MNdlEK=F$Mw2e#B;>Rj!+8(a!8{hu?lkZ>f@<@q0?oal{5+@t3>l@#074OMxU_LWjLT>jwM-qrPuZ>RN&N>`76GCqa3vu}L6 zT~FSw>HJp4d-3yx>-xsG7p?2|hi;ute;(hCTAvb~fDrUZtd;uuF;&l_#TETo36JL!R$GRy0mER_;xH`CCpbcPsjejioY=2k5oqnRo=y3DaBr! z*X-&8;Q7IX`Kp(u<3RBh$1x7-$GRy06~`fu7zF04g!wAZ&tQ47+vf*m-9#(5Gs=IJ zu>T6pSBau|vhiU5RYH7q9h$Ea;wv>3JU^HaU(K8Wo*ztzujZrqim&ILXue8_uWm>4 z6+ahXmIBWYCd60M(R`HY_;^**O{V|xY zQsyg}S4Z&`pVvVWzkuQ^IlqI!e3dd^Nn8xdf5q!2rsa7OF!5I@`>*nRP?TS^g*h+$ zcZy&9RZ4ud9nDuM@zr)TU!}xXH=_9}CB9mR=Bt$W>NYfA@%7w$7GJ-p_^Xuo>M1l| z@pA#LNAp!meD!%WU!}xX51{$#@O6)h=p|ewC4XFrdzZRf{{D&VwO3Iozr(xkpM7=R ztt0?-SO-LMRYWD{EJ-ViL zUElb2vhNtyZPd1%&#Yb7H@@BapYI>Yy8X3n_hW059} zl{9S+zn+hj9djvqQp)CRnZ&(|$7-rK>vC{mIRSQw>*2Z*w@UlVe;>KJcA9;@tf!N8 zleO#m@l_O2ymcN2b>GVD3YT^dxA%B4luNNwO4&?tn=^E^oq&k;gzJv+)V0rfIH;Zb$3L##ZRe9cb?x(M zR?N{uFYWC3Yo)*NkE?feedF6l+61zv>$7Y?f2*d$mAKXLcJ_^LxB2mxo}q1j`LCiY zzr(|IedF8bl%M)gI(UYz9{*%~3fJ|GZ_iHag3IR8I8+($#m^J2>l@!*w65DPxplgJ z40DEIu0rvZQO}lF4^D1)2FzDR(f#u#n6Hc?zN+K<7R6VRzx`O&Me&uP-$(aLe~jWQ zqby&IF0s8RzB0=4)eC68GRpGR6KK9N%JS73G+*($n5WTvWt8QsSI~UL@z0CVd}Wm7 zt0A+&d}Wm7t8bzC$|%cM_n`U8*xR?ySJ8ZBl;x`-bHRLNl;xZo(0o9_mH>#Onh~BiQ5^) zS26L`;(Q$+#(xzPU(KJ$_KH76><|Av=;yzRiLbVy`6{M$$H%vW{a4&yjmYP%BKwl| z4sT~a{}snSe-F)9G4WLgny+HwtF36hiixjwqWLN&zS@lDtC;wz1I<@4@zo|YUlp(G z@|0Vr{vnTV$MRLed?oXAslv<;Cd^m+!F-i4U*+-b^!vY6(|HK*pZ(?sIls5b*AHUm z2RVN7{(K!DhOZL#UrD>7_$p!l6`HRS_Ftj-Dj~jl56xEz@zn)rzT*C>6U|o~|NH`) zuM*;`MQFZCh_9YR^HoB8^*Wlb65^|e(R`Hed<~>vJ z{9sDgRh`V!q4+9gzWNBvS1I$A%&VjPSCZemA9Ebaf5q{W_vh>QFnpD={|e1lDf_R+ zaXX{puTu5H!4r3n7&VD?#h+5*dN{ZH`L9yqt7& @brief Generate evenly spaced reals over the specified interval. + !> @brief Return reals separated by the given step over the given interval. + !! + !! This function is designed to behave like NumPy's arange. + !! Adapted from https://stackoverflow.com/a/65839162/6514033. + !< + pure function arange(start, stop, step) result(array) + ! dummy + real(DP), intent(in) :: start, stop + real(DP), intent(in), optional :: step + real(DP), allocatable :: array(:) + ! local + real(DP) :: lstep + integer(I4B) :: i, n + + if (present(step)) then + lstep = step + else + lstep = DONE + end if + + if (step <= 0) then + allocate (array(0)) + else + n = ceiling((stop - start) / step) + allocate (array(n)) + do i = 0, n - 1 + array(i + 1) = start + step * i + end do + end if + end function arange + + !> @brief Return evenly spaced reals over the given interval. !! !! This function is designed to behave like NumPy's linspace. !! Adapted from https://stackoverflow.com/a/57211848/6514033.