From 088aba7974c13796e01f044a4fc8bb0a5ad040e5 Mon Sep 17 00:00:00 2001 From: "inwoo.park" Date: Thu, 22 Jan 2026 14:44:50 +0900 Subject: [PATCH 01/17] CHG: basalstress.m - Fix bug in checking md.friction.coupling value. --- src/m/mech/basalstress.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/m/mech/basalstress.m b/src/m/mech/basalstress.m index 4350c7dbe..68fdc38b9 100644 --- a/src/m/mech/basalstress.m +++ b/src/m/mech/basalstress.m @@ -16,9 +16,10 @@ sealevel = 0; p_ice = g*rho_ice*md.geometry.thickness; -if isfield(md.friction, 'coupling') +if isprop(md.friction, 'coupling') coupling = md.friction.coupling; else + warning(sprintf('md.friction.coupling is not found. Default coupling is set to 0.')); coupling = 0; end From 8f3335e2b7b7b7db63d9d1dc6b7176807c2401b5 Mon Sep 17 00:00:00 2001 From: "inwoo.park" Date: Thu, 22 Jan 2026 15:18:30 +0900 Subject: [PATCH 02/17] CHG: basalstress.py - Check available md.friction classes and add Weertman type friction law (Matlab > Python). --- src/m/mech/basalstress.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/m/mech/basalstress.py b/src/m/mech/basalstress.py index a99f2eeeb..6b14d1c7a 100644 --- a/src/m/mech/basalstress.py +++ b/src/m/mech/basalstress.py @@ -1,4 +1,5 @@ from friction import friction +from frictionweertman import frictionweertman from frictionschoof import frictionschoof from averaging import averaging import numpy as np @@ -16,7 +17,7 @@ def basalstress(md): ''' #Check md.friction class - if not isinstance(md.friction,friction): + if not (isinstance(md.friction,friction) | isinstance(md.friction,frictionschoof) | isinstance(md.friction,frictionweertman)): raise Exception('Error: md.friction only supports "friction.m" class.') #Compute effective pressure @@ -74,6 +75,12 @@ def basalstress(md): Cmax=averaging(md,md.friction.Cmax,0) alpha2 = (C**2 * ub**(m-1))/(1 + (C**2/(Cmax*N))**(1/m)*ub)**m + + elif isinstance(md.friction,frictionweertman): + m = averaging(md,md.friction.m,0.0) + C = md.friction.C + alpha2 = C**2 * ub**(1/m-1) + else: raise Exception('not supported yet') From 64ffe36857ac1c0f79f36faf79f3da3bdec6e3ac Mon Sep 17 00:00:00 2001 From: NJSchlegel <95945328+NJSchlegel@users.noreply.github.com> Date: Wed, 18 Feb 2026 19:15:27 -0500 Subject: [PATCH 03/17] CHG: add pressure to GEMB precip scaling --- src/c/classes/Elements/Element.cpp | 18 +++++++++++------- test/Archives/Archive258.arch | Bin 24551 -> 24551 bytes 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/src/c/classes/Elements/Element.cpp b/src/c/classes/Elements/Element.cpp index 1f902fc7d..db7206881 100644 --- a/src/c/classes/Elements/Element.cpp +++ b/src/c/classes/Elements/Element.cpp @@ -5659,17 +5659,21 @@ void Element::SmbGemb(IssmDouble timeinputs, int count, int steps){/*{{{*/ } else pAir=pparam; - //Hold relative humidity constant and calculte new saturation vapor pressure + //Hold relative humidity constant, calculte new saturation vapor pressure, + // and new saturation specific humidity to scale precipitation //https://cran.r-project.org/web/packages/humidity/vignettes/humidity-measures.html //es over ice calculation //Ding et al., 2019 after Bolton, 1980 //ea37 = rh37*100*6.112.*exp((17.67*(t237-273.15))./(t237-29.65)); - rhparam=eaparam/6.112/exp((17.67*(taparam-273.15))/(taparam-29.65)); - eAir=fmax(rhparam*6.112*exp((17.67*(Ta-273.15))/(Ta-29.65)),0.0); - - if (isprecipmap && (eaparam>0) && (eAir>0)){ - P=prparam*eAir/eaparam; - C=C*eAir/eaparam; + IssmDouble esparam, es; + esparam=6.112*exp((17.67*(taparam-273.15))/(taparam-29.65)); + es=6.112*exp((17.67*(Ta-273.15))/(Ta-29.65)); + rhparam=eaparam/esparam; + eAir=fmax(rhparam*es,0.0); + + if ((isprecipmap) && (eaparam>0) && (pAir>0)){ + P=fmax(prparam*es/esparam*pparam/pAir,0.0); + C=fmax(C*es/esparam*pparam/pAir,0.0); } else P=prparam; diff --git a/test/Archives/Archive258.arch b/test/Archives/Archive258.arch index 48c6811b53c19d43f9276c8832b8788eeab5778d..d71423628a4ae10730136ec5b99a709a227215a1 100644 GIT binary patch literal 24551 zcmeI)c{o+u-vDq)3Q2}iq@?0rrHG4A-$l`!MvamNO=uG`>xd}2MI#~eP>Bec3(0Uz zg@kYpGMol;gnE~*rQhTJ`2FttzQ2Fo_nhbH`8?l!zI*MpziaKa+t2ZAHa51|Y;0^3 z@IQWa(<3KNS{z&co3ft4$zz5`)=e-D(8ypv}6*a=H+tQKYT#yA%Wotc?eRN>Wt9fs0RN^ zeR{`oC5YI)*k|n>5%6#GdBA)s2tlDjf|vR}LfA5yazhIq2<(!++zly=JwJE`d^UGXvPw?_f8X!t*k43Jpp#XfL1i)ol`0(D z(;EOmUtcX{&pQmEHn4qiP% zW)g+{5OCNg>ZB1H1SD47*OIva!6m^Pos`re#O(Frd@Fgl*S7e$>PLWk-}6Fb^P(X1 z)Z9r5r}QA;_f&&ohjzF-t5&kT@PBq4z-wHC(@6%5=w0MkY~K%F#2{sJW6u5H zMU2>a{$5ifcoEz0(AKr{f!DC!8VMP)H+^A5qfM=6ld|yr;7@~RamTeXieNEIC(ss%FN$SylvUOPF`J%K zD5o6?#2ofWoAzz8K+ILJRNta;28elD0rQ*v*ZBShhsV zh&Cq?h~)=UTkm|F0mL`K4@6Y?lz~|Bt-ar8E}jp>sxGu4gZmv2-+5?k z%v#L}#9A$1{nZ2cK&-o*@iw+Y8;Ep=km6{i3?MQbA8F52O#)*54gc0#!Bs$Ph~I)T zS5yMAvH6-{-g|J{zrc|#Bh<(Sz*mYj_1M#cpqm13f zl0fV)t&R0;KLNx+S)ohCuVR4sEi^;=RB;asi`eF<-{#T@#m_TiKJ>_se&|6|ScmDr zY8^zyWobyObR#PLnJUG(9Z@-`KTyjZQLorJLew{-=#SHbkO|wS0`u_u$YetCTlvTF$b`Mw z>r7WGGMN~i!Y!18OgJnKBs{r^OeXENwC0gQCY&5O?z*?ICuSellEK!N70qNIlgWo3 zaM&EemR_xX#End*?95U#b3!Iud~UBj7h^wZxosCmCfunLODyw|2~S|&s+p^i$<*77 z4S#eY6W&*%OGUep$+T{|vzsq65il6~{Qr~pZaYnd!gZNiv|J()*K5hVX`URou3u)h zxmDnC_;5?)!W=v<3k1)0y5Mo^`*0{F6H(ujgKtg5^J1{yhP^2nQT^Wv7Z&B?ac`M& z^l%!Ydfhj%)A9WI;v^$TzxzWq5w1U*>h5im0yqOA{VlLR4q3`r;i-ESzne z*NCVN{nf%In-SHnRUsGVfvC0}Y@v@t5cS#qvR=q5JkN}HRo)9A>f@dh%IemL`hXti zFQ+4lDak(LojRggrmjPSd7h~V-rn?L;n6}xJWuc4J7w?p48_Gu@cGnU?Weh3c^aHEfH1fv5-R?Uk9&u z&QX|PDwuc4ygVFH`Skvt^`6N;kW${f@8WdlNmf@pz=2o+t618&RpsPu)_a5S22y$xe3wQArhMXO(^FNkqjj4LuRxg4YGEg+cY5h>C6eV_$eEqGGb5dE#(g zABWH{86!kRUoz=4!Sf*Mj7h2CB1Apf-S%G03sI38-!}U;A?l%QxLpOV=L7MmT`^oAf4MXK*9H-FPjkw#2zf;L?a%tEH5XC7<`K15{SoEk znS47DANOv|x7lws5#?QOR+^8;-)m4JMhcI&mlXHxRr3*bXXla*iDpFIc5)AY7Kf-? zPotVcC9^68dXKPLLGY+g=>&!kpzdUSZ zl3yAi>RN7?O59AW$c|T+_95!3vFr)`Cy4s}#fr->FJtGOHcL-J)D_!|=_Po5b$_Qk zz=`|sj`Cd_@Vs^lGt1cj49maFIFdrtWxE~$tzC$^^f9BuLmzwS+N`wWh;l98z0-a+ zR$E|Ij4z^G@{-hq@bO&MuXD)9^}U!E&-=9syQM#Ebqb=K%LP^`;{BZUIs}vNBg*N+ z^?Sy+?oKwJ=Q8cEyaQ3w)DU$c@X3~6zG78)YkoCGlw;`x$(SN6<@`oE2+Mcf_33j& z(H_r#rh?C(wo9r_@;5}C@AtJnsf!J^idl9Zs~!6;7$47JaJ5n>K7WVEx$Q=H{yFHa zQdo2hD|ko3Y6_z4^A=PJvU3O4lV8aYpO@`{kkp;H9ot{I-A)3cY`!Qqhv55? z%}WW%%XhGMG6UH^V+{qe6mUOm)F)5eh|j}D_;$UD3Zkq(r}P+?V{?saw&C@`Iyi8n z&TA}f+raiSSUs~;j{q#NpYjdF&gTGjdD{JJ)hhrk| zk<;^Ja6GK#;6+)>Rz!`+oUr;Pp4XPd@`{!!Tu)14h7a4ZH8>_lCAQ~3z%kKZiBHrD z$3)K;2V?PdXi0Q&W<=sTTM`}C^VH&X(voO>{Dn>kj)|6Z%LGRp6Rl$|R;dnsKdeKY z<>8A%`w#1Cp7lHi$HW`fOO>aLLjAP=(o^^bQ=4IJh3dBO)Z^=jh@wy#UMFmc zr}nAq;e1=7N7QCxd>z{ow0A+~Tyy{i*sQKaPnuE4dFlfb*T*;}M!Yrfq~VzOC?O*g-;eBwF}JSY!uJb1Vq9Pa&#yQpKAqTVhwpcG z#KbDIPtrIhCTk|Hor7ay%07M_y#Ck`Qy2eETZ-cmbHkbVerK0DY#Mr7HVMZgdKVSr z^~o-c`0iyB9v{23VbggxGUaeQqRd2F@u8Tgp3#0~C?;o!G0{BsNzG79Jm-G={ZLG_EVcQOPqd1&oi^S> zkHTLbdd%C?ns{yAAoK3NlePOUB8!PN?)Lee$U?}JEz401S%`BPcH}XUh14SdjbcK` zV#zw44X*)Nm_KYew5A$aH0N!WIbzctbqzxj6l@0Hk&YPk0Lg>EMS}t@x7*gYnuOgRF`B{pcQ;|#H zYwfkt50LBY&FX{7n#j44ldHv6068}d@P(Y;ftr8DsJfzR;0^I(~XTCQDjv<=It0q^K54c~N7v8ELK(tu(>24D3h?Z;S zRiDO#XkR@C4|G-{N4B}Q?u1=Lj{8p61r~2ZvM393sI%A%g!U(-TSxAkA6qA zhm(p@+VK8`eGB*0Xd*|RhD6o2d&qH#W~}C`*T_+(D8%sS1>~p{C{-@xgX`07E!MUL zIWD`D`$+mAa@=}#`=|AwhRYw;=MXmq_mH18C;urPO$$H%G2DdNh$qSE=uMOo5d zkk*b|F97X*O*Fc^^4rnoYA)oyXc;{s$S^w%^>cum-u^8ZdJiFhQ5C z&)OHwc14#9HFv*xXNugFHB!42Es#5f=P?@_OI96N^MG|7u+{1^7SnB|59bl~k z|7;z=+s6E*=eN-hJrLeS*_^GM0{7ov*Ip~22an^wB}p!1gSh%{98(y)kgCvjZ|+Pz zNLniKQRJjMJk?HaJ$FM068si+weI4EgfEMi>`haG`1;Mk11UN1a^q~-RLxR&QLTMH z)I|_-3%Ad>5PKeSvxAMB3(_F>mN51EvLeX!5_3K?coMP;-tnYY9f0B+vRZP3ijYgM zuy)V6tuIRRg+D@UG$` zwLlM9eD-S9GN7MnZ!ey{7U+kcCuPl41_rm3u%gRjV1$%Ohxv*EZ3Y z!6RV2jqqLUZ3v8nT@M#E{|1b&jQpa3)xa>s(H zu(p^Aj5PsN3wnadZM`m8&eV(M|LQQ z=2ijSw}hS@eEWx2t?ohHtgouOHaG))Qph^J!g)}qs&OxRBR|w_vaM|4yAQR^jsc%V zouOuT#pzQMnxIB<|4XSjQ>Ze1uO70@1?t3=jnwyTgW4^ZPcnj4pl;gI2zTXpsM}Q( zZJ5>tb@7MZdPvNHx(8Oh8#SGvZrhhhhV;u&C-Y_LcbOwln|W@9ZciiBu1tiz3-zFO z#>OR2b}B;6g#+C&pVZ*p{3G6)>Vr^|8CUcvQViauHFNkbph1;~($V=)&bT!z*+}b>%c!* z2k@pbA9@u3^3c=vqlRqonfv`2j_(MVGdnfHw;qEg*D}jH>Q{gsyd@U zK7a;O!G=u>_d-Kuc|wotLTCu@kt|g-fX1zHExEjff40WP23MG~CvILKLSmMz3a~1` zssO73tO~Fyz^VYN0;~$KD!{71|6&EkeCSd7%R^7ypNEB-8GqDpeLD~3GVh(Y`sYI7 zf?W)KPBlnYmNAjcZ-CdI27A^Qe1LbytheQ#^MSf3*Qtspl%Q75aZ}P^XJ8~U>bE|Z zf%391>f5=cfvzhhF!u@mWoi5mv?Lc4J!|6t^6>TFbo`^>hp_))yQB9X>v*H{{?qYT z{UF=_)5rf$wNLubl2rj#1y~heRe)6iRs~oUU{!!sfiWsD=0lJ2UmkjxJefwNWm>?T zT)Wo(tu-*YZm;e?w-}gH=AN@CQ3dAoBW*8NZN^eu^e1P5$-Us{ULkj2a?YmP9O?un zU*#Uxj9nwP7Axc_%n=77bC$SHo#PT<3SQe@lI;SmEgtnDa_4?pth=*Ndk-*qB0AL5 zjzjbEd352mpGItTUA^zQZZa@=<)RHkQhz#)tKjmU6}VsfD%m+4aQhcEGY#| zUWc;qFIRn;S%ZPV6q(T5{7eIwGbEelL{PvKmP+gqQvs&X$9-(8n{hh@h0#-2;{B&+ zY8f2D-jg}*ISp@5?qClN2ByGu=663Me4g9ArDos9v0J`ju?}vpf5Pw1xGo}k>~+&| zKSjcnd!gb<4?KjKa0N3G3$+MH$_&DAA$q#Tlv)-T=YU_R~;CAM4=5}C; zr)453JWlhrEk&-lU6J+}h0+{o9e89o+4|^+O`Q?zA7q>1KOg+CV`2MsV9Njauw>O? zr~<5cz?uiFc`$Mv$PHf)V>$YI8J)Lmcz>eoa6Fcy_y4=`=>Ck3M<18;>)&;Be@3@I z`upX>=R208_y1S%Z^P{o|5g52$0ggB4!1}AyY(UM5J&fiY#+>SOyRpLXco4}r6P1R)Cz9<%#o_OXWcyf>JR)g#ER~1bC*uEa{kUHcNAE}4 zC6e>{cl}rO8|%Etc5>XY9Npj1F*#2n{tLFT|GaGl{_@BBtXX8aBgS!xY_uZ zwc)lUQ!kb%_T3&g8$Xf`w>{hM^1W%{8^_Ha8B6j8zVWh`{~XKV_oI!QJu;T$ZEE9Y zkBsHNmHl%p2NJ)CMm~=ox9yR!gkt6fhS=V+akEFp66hJhRG*gNxY;9PnSQOhs^-w1 zmT|L3#uBml%nf+6DQ(>Bk+Fn~r`u+A8f+Rjdt@x3QapXNSJ}StvgJp{5{d@T_g%Y| zJ8t&KSi(1-$%&~6FmCq9Sf-mwyC-XYA6gm4ye{Q`j^*=5?7eEGL$8h;JNxHY+SfMK hYu`LIZuZDnB92^%lO<~U<7SVHC93eiZ*(h2|1TX+Lc#z5 literal 24551 zcmeI)2{e`6+c0nyico|IrBW1SillOFO0xzH8c2#JJl#|>#z6=TXhO)GnIt4*BxxW@ z%1nvVKq+K;x92(g^k1IE^Z$R}`o6W^_dRDV>$k3R@3ZfHU)SE-x;whr*w{qb*x1J5 zAAU{CLr0HUnJB5K86G}nVrHy}-{iqRF&ffH*dp!zmNMzh&u)1q!DiXx}ui>-DX0 zQ`tC6=+|txxgtA%z1jx2WfbP4^N1g!B;UnsZ@3C!@;}CJoO&F>l%Cxb?3)j_C!Lbf znPv$wcGDwc%&j0`+j5N(V-tvw^gP{tDixxP^J7bO%OH{pTUGWk4MH!UU&`$+0il;Z zt5cF*?@0we>5M2H#ZXn=3MEM^U`l=EN zHw(gY0vgreX4S+;*_S+CHFA}nQ zUKWJE*|{!Fp$~#%Bo<7(I|G7~8z-F1-v$w9uGwX*wT9?{q1lNa!ys5w)^uUO8wg%K zUQ=|cAw=BXcq70*6#_%{mK@oC6mFt#ria28!>#bfK`(_y2=j8BYcn-s@*tw(YMsW@mk??7mA2-s350wvo}Zkv6K;)*Y4w?*1X1RL zlV2o#gi!Aw;ni0qAyAAS@Ak$WqG#HiX)wikC_Aa!&r=P-a(3J*`$Zx62vr-XYYpKd z@$~q$gAg&JF*kj+D?~LsDV*S}0ztN0?ngL&h9K{#lx1h?AapmIN#qpANi&<2PhcBJJ{<+~wuYt zqTNBcw>N;8Jj*+DUVA?`c#XQy|`ZQj}xl*a5`!J~)>iq5#B94cFHaF;zgk z?^l*6bg2f2S-%cytE=||G26n|HRMto5OaFkC66yo1!CShF+*{mCqT>>d;WH^!Zcvo z{pw1+kQ?MPET~B|xXd05#D}HRE(%F@0yZU(c3l8+ErHo zvG`Z`w)!AfAeP9QM9Y7A48$jY99|T%jRwThqRaVKyL5o~OsTqM$zcm1KF?70z2S2m zh-I?*)5mKH0kNEGUB7eOCSdx*t;tcl8_F5J>|9>>`HmwHD+c(hPEjgAd?gy#q&C|Z zh_Bb~ig>eK5QtT7@n+M*xPkbl>Q~i$;{!l^yFB39+IAHnzDxNScF!pth_%bN8mL1r z5bNITOyfS94aA1kFSG{--UG2wdO%xMP!5Pq!b=YHy!HcP^Hkpbo9`S0VvDf4MYrt* zAiiJVra|A)2*eK?%HJHC2AF>7MgH<*v=C;B-FzPL^?H-EC1DX zAbx7)@SJ$d3W#4Ueb$$aD+OW)_r!*%DZ;??0~MZHoe3QbJ728igvmZY{8prBuEn_) zh~G=Z;_`h9f!Ot~OqN%00T6!-O09ozyab5d3h&4HT=WKF&#A|;+du3EV(*i^y4=}K zK{G5>HDtl&RkVga1zC(2 zn3I(=16fRva}(h_gYC9WRZc(_?0G(~o@_%F98nuSZJ3TMIQ8TcqMl$MY!L=UWHC`y z{h(AT_Ps`$f*7)x6m%)rvkF;onOaP%Q$-fsR~ClpG$V`2?~lsNIe{#AXK(nXgOCNE z+U-LvACbirDWkhy?~#ST;o;Z+UU>^E2ObO}3!W9Rp+oy|-V94E1;cPY7p{y^ZpZl+ zpR`9d2G>I;ac_wqu9t+5*@6;84OOL^&Gy9gJRf@dqZbTGBC01; zaj0q`qPk;Nt}L%c)X(&xjXLx2yePT+@Z}yvbv+LfPI!-~@1ND~UW&x3?PQIi0aVun(6uw&nxE^F~wDg`eL#0%8438eb!&s`5Di%PpfxMt5QG| zeZi3#t$PvGKG}b|o)Mxx(mB&x3lY^?{B+O^kLyDiSN|D(M7_7E@#hIZRLd4)L$#}j zYMzt7qYwAj*fp&>1kcxohxVQ)KOw5#M_jHC_g7~yEyyMw&vR*)e69>cz3aGOYh;M1 zw;9`RU5Z6i_1TVY3GAEI-qIpUh^pdRl@ZB?s8`R_c&dvLReAZmUbHNtDz@wSr}yG@ zBUBKeDvhYJsAF!~b$ET>TiE;R1@1SE!`L7g*E_YbS4S`jQFro6=PbngK}sk+(clk6 zC7&}BR{MsiBvbuz>05|O-2N=aekY<5G@d<3>qJz%oT**^S472$-;IB*hN##n_l{4} zK~xNfgH$T6L-e3vWWWJLMfG~{`P@TPWOrkYWFMj;dSzSJxFafj;KY%GX;}8VU29e# zDr^e7hY`+GsJMyCG@O@Pa&>XMxPBqXN}+K)qJs4lQ+$;X6=Wg2Y|8~i-MlEXMg_-# ziKeAe%McY%-l7yMiYWh{D@I+d*g0pjD#H+UV@u@R;$4Wk?!+~Y7te!hsd1WxcM)~9 z-FrWG6ILpI-y9P}T{+18Ja8qV{KBFd9V`*$`@wPhNxVLM7x_B;C`FXdnL51_xDJ<# z4p<#A#R>}xhHxP2(y@L4g&;&-%;)QgLRhg4gOg_?>VkFr#4a&Jc~>t=bS}fLCcG*L~w9^8uMUEbq*x@~lpQ7J&wxtEbi6}V34^nGu>lt7ego8~5dybfHg zW4gW4vD}3j{nJYS}O*RfM*y5^FJ*iC~{JMs7(e+7uKrt(mVnK#}H-rLrYOg6I-PEm2wCs57g~wx4veslNKCev@%`xZ|_S)dwYdF6) zjv1>z&cYshQEMvuFYfp%+Hw)cL{0hJySQH);$nBP#*H{8idR*+isG2a!+YcGVH`8G z?pND%NCi<$Kd;xgF@)plJyL2v{>1{ryjYwsYhseB*UzmuCWc7QIZekg(I?YVcsh=W zZgutpxi}`CE&nCI5648C$_qYd%s?ZkvnD{~E7-u1l8QO`)YiXI|`AgiCNmIh}(2nTH z!MA$&I5ByJ@E&y>6C3zx!~4A*Lwi1rExZ@-`XQ>*n)8)#Of(KNndpUMqPO{6d3-*5 zVoJIsC5~farNJHg7aSA2)~bKA!ZAY!{zElCYbI4vvZ8N>jOp+Y@7s zS;ph_;Xq70BK{WV$ANff|B@9iaZF6JFN=1QWptEott;^V1=lVGgpCFDy!!hwpoXecym}vdW4F8^R$R*lT*BK7SL|gM9<>8p=bGr>=-SjN^+fC2u z&b)zD4vy&bkABHbi=C0xgyt~mCp(dqpskww)oaK~Qt0U6#35uQt#X)8Zzi&m+gHBI zNDx`sH)kd7sYX^=tww7NrlK<{x)VGi`;hh0oQlSs$C34s?W(hncq1FfiU_(!4zi8> z)W#8%gBImJG!j1dw-?lHk4Iq$|J_sI4v%4!_) zC_4GL&pQ#hT@b5ZS+D@P?lA~+I}n0g&fl)K%Dj$T_H&7e1zM!i9_?WRUOWvPf}F5N69d=Sx&8tU296e605 zhB!-Y{*kJi)ky_*~Wp6sY>>={FsMV_^U5Y%{-2b7p z+822*%H%1NHbq|cHx#XsCn2v!!{r9fKam%|>A_6hbmW;xaa`c=MxJ-!m6*Id@>W}U zb5pq;@}}_oWn*K>$^&a2u+{-6*Li+~@PyIc6>bG4Mzz)q=Z69vWp`*L;`02#pRo zCsjDFK;yunqo*2#p=s8q$|UjI(CBMSy;i&njTLNP_pSGW#;-ru3?)y4rg)K>_f;;? zlxntIdvhiZd`|yKdPByj{>#yS&v}aV|7v+;ev{Y--TF$hC`6rv+NH`%tesq;E~oWU{nOb{n=W4S zdWJH*yUbPiKr9{JSV>+Caj$@?zBSYPWLx2Nh>G!q$N;FAD=78`e=Yga&hlbvPar(s z(!5k?+#x8{UUh7y*fw}Pk5}X}FB?3XA=PKbEddX6+T4@-BoO=X)W(}UOI99O^MJJu zSo;8LA7JeRtbKsB53u%uG4=s`Xw;h?)xX{Jbo`l1Hu%c@C|lO%6wrlrtZVPhhgNqZ zcX^@H&=|Tabr!b(H1BLKkKDTlT23r}J3nyt@;SQ%htfRzDO23Q$jWq_3dRt8uZ_;SqwJyRWhqL5|1D&UD z>eE09+WS_lQW&(iS&6Z4vixn?a`CChe1SgM ze~*%k6?|0UJYdUh$h6&Ko~F^YPTarEtu347fzGexIRQWAjLz@Vsc%<`&#$c^x^*vh zpN!o2DxgpAEAnsh|E-YBx+`K%IF=DyAio&sf}d4W4E6E#6xH>td5Qbsmr#<+$K(4L zX|U4;_rvBRSKsw7`ZcYp&;mM$=t94mHqNdC`V7PI9C6Or=%xHR z$v~I5Y@J$y^F*AU|D#-EAN~f-JL2dW+&@no;?QmT?KD@VUEe7@F8}i>QuN<)cA9Oi zeul>vTG1}F73h3#Lf!D6$I&OdEUI|Oj@#|O$RnSQ$Gy~ZV7d$*&xN3`{MG+rR^xZ$ zsra~+|C&8GUjlzvZiMAP7e!VbK2z{`_t;V^IDkIqsPL(WnK&N;{ueyBe#>v%@?c#C z(0P}ND5y8$x|=>vZg`0EkY;rC8@^t_=d-Kt<9vuZE0@2U`&$986H^2Xfi9WTEwjlA z=rbo@=X`61=T&m@c+V?L+xu(wCi?P0>-S&I6<1}UDW`K)=;d`lSN`+Gl9h+y46x<_ zYaaX$o(D@9^LsRv7{^Cbm2rGD)fvZ$q#p&w`$YVE@y~fi#J^|%jESTjk-SeN{Sp-! zpC^u7$K?IdBP!5m<7D3dyOMs1vW)W)6=MvJewCXLw|f3wEudw!0`lA4aZFgRXLQu+C^2h^ z|1Z_-HW8koFb_l2z-vo*7pFnIHIxJt3$n5~($%(Yz9W?C9F{O>2D=3Du9 zE%TE~)2d9o$87z(mS>yCbLSKeuhG%_X4Vq17d>j_w|p{YYi2Fc8$bMb%)pNS0@xV` AW&i*H From 6bc8e4f57e1f5ee5c0648223cbf752438e534a26 Mon Sep 17 00:00:00 2001 From: Alicia Gunilla Maila Bratner Date: Thu, 19 Feb 2026 08:46:05 -0500 Subject: [PATCH 04/17] frictionjosh: cap scaled friction coefficients using user input --- src/c/classes/Loads/Friction.cpp | 7 +- src/c/shared/Enum/Enum.vim | 1 + src/c/shared/Enum/EnumDefinitions.h | 1 + src/c/shared/Enum/EnumToStringx.cpp | 1 + src/c/shared/Enum/Enumjl.vim | 1 + src/c/shared/Enum/StringToEnumx.cpp | 121 ++++++++++++++-------------- src/c/shared/Enum/issmenums.jl | 3 + src/m/classes/frictionjosh.m | 7 ++ src/m/classes/model.m | 21 ++--- 9 files changed, 88 insertions(+), 75 deletions(-) diff --git a/src/c/classes/Loads/Friction.cpp b/src/c/classes/Loads/Friction.cpp index 6a0a3cc72..76f8f8587 100644 --- a/src/c/classes/Loads/Friction.cpp +++ b/src/c/classes/Loads/Friction.cpp @@ -641,7 +641,7 @@ void Friction::GetAlpha2Josh(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ /*Intermediaries: */ IssmDouble T,Tpmp,deltaT,deltaTref,pressure,diff,drag_coefficient; - IssmDouble alpha2,time,gamma,ref,alp_new,alphascaled; + IssmDouble alpha2,time,gamma,ref,alp_new,alphascaled,max_coefficient; const IssmDouble yts = 365*24*3600.; /*Get viscous part*/ @@ -655,6 +655,8 @@ void Friction::GetAlpha2Josh(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ * element->GetInputValue(&Tpdd,gauss,TemperaturePDDEnum); * */ + element->parameters->FindParam(&max_coefficient,FrictionMaxCoefficientEnum); + /*Compute delta T*/ element->GetInputValue(&T,gauss,TemperatureEnum); element->GetInputValue(&pressure,gauss,PressureEnum); @@ -669,7 +671,7 @@ void Friction::GetAlpha2Josh(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ alp_new = ref/exp(deltaT/gamma); alphascaled = sqrt(alp_new)*drag_coefficient; - if (alphascaled > 300) alp_new = (300/drag_coefficient)*(300/drag_coefficient); + if (alphascaled > max_coefficient) alp_new = (max_coefficient/drag_coefficient)*(max_coefficient/drag_coefficient); alp_new=alp_new*alpha2; @@ -1419,6 +1421,7 @@ void FrictionUpdateParameters(Parameters* parameters,IoModel* iomodel){/*{{{*/ case 9: parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.friction.coefficient_max",FrictionMaxCoefficientEnum)); parameters->AddObject(new IntParam(FrictionCouplingEnum,2));/*comment this line to use effective pressure from Beuler and Pelt (2015)*/ break; case 10: diff --git a/src/c/shared/Enum/Enum.vim b/src/c/shared/Enum/Enum.vim index 0c6266e60..8ee80ebae 100644 --- a/src/c/shared/Enum/Enum.vim +++ b/src/c/shared/Enum/Enum.vim @@ -221,6 +221,7 @@ syn keyword cConstant FrictionFEnum syn keyword cConstant FrictionGammaEnum syn keyword cConstant FrictionLawEnum syn keyword cConstant FrictionLinearizeEnum +syn keyword cConstant FrictionMaxCoefficientEnum syn keyword cConstant FrictionPseudoplasticityExponentEnum syn keyword cConstant FrictionU0Enum syn keyword cConstant FrictionThresholdSpeedEnum diff --git a/src/c/shared/Enum/EnumDefinitions.h b/src/c/shared/Enum/EnumDefinitions.h index 70557519f..a1b15f7d0 100644 --- a/src/c/shared/Enum/EnumDefinitions.h +++ b/src/c/shared/Enum/EnumDefinitions.h @@ -215,6 +215,7 @@ enum definitions{ FrictionGammaEnum, FrictionLawEnum, FrictionLinearizeEnum, + FrictionMaxCoefficientEnum, FrictionPseudoplasticityExponentEnum, FrictionU0Enum, FrictionThresholdSpeedEnum, diff --git a/src/c/shared/Enum/EnumToStringx.cpp b/src/c/shared/Enum/EnumToStringx.cpp index 2d4d382b0..ab845e22c 100644 --- a/src/c/shared/Enum/EnumToStringx.cpp +++ b/src/c/shared/Enum/EnumToStringx.cpp @@ -223,6 +223,7 @@ const char* EnumToStringx(int en){ case FrictionGammaEnum : return "FrictionGamma"; case FrictionLawEnum : return "FrictionLaw"; case FrictionLinearizeEnum : return "FrictionLinearize"; + case FrictionMaxCoefficientEnum : return "FrictionMaxCoefficient"; case FrictionPseudoplasticityExponentEnum : return "FrictionPseudoplasticityExponent"; case FrictionU0Enum : return "FrictionU0"; case FrictionThresholdSpeedEnum : return "FrictionThresholdSpeed"; diff --git a/src/c/shared/Enum/Enumjl.vim b/src/c/shared/Enum/Enumjl.vim index 14c41d125..2245470f3 100644 --- a/src/c/shared/Enum/Enumjl.vim +++ b/src/c/shared/Enum/Enumjl.vim @@ -214,6 +214,7 @@ syn keyword juliaConstC FrictionFEnum syn keyword juliaConstC FrictionGammaEnum syn keyword juliaConstC FrictionLawEnum syn keyword juliaConstC FrictionLinearizeEnum +syn keyword juliaConstC FrictionMaxCoefficientEnum syn keyword juliaConstC FrictionPseudoplasticityExponentEnum syn keyword juliaConstC FrictionU0Enum syn keyword juliaConstC FrictionThresholdSpeedEnum diff --git a/src/c/shared/Enum/StringToEnumx.cpp b/src/c/shared/Enum/StringToEnumx.cpp index f31d143ae..ec30a99fd 100644 --- a/src/c/shared/Enum/StringToEnumx.cpp +++ b/src/c/shared/Enum/StringToEnumx.cpp @@ -226,6 +226,7 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"FrictionGamma")==0) return FrictionGammaEnum; else if (strcmp(name,"FrictionLaw")==0) return FrictionLawEnum; else if (strcmp(name,"FrictionLinearize")==0) return FrictionLinearizeEnum; + else if (strcmp(name,"FrictionMaxCoefficient")==0) return FrictionMaxCoefficientEnum; else if (strcmp(name,"FrictionPseudoplasticityExponent")==0) return FrictionPseudoplasticityExponentEnum; else if (strcmp(name,"FrictionU0")==0) return FrictionU0Enum; else if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum; @@ -258,11 +259,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"FrontalForcingsSdNumberofParams")==0) return FrontalForcingsSdNumberofParamsEnum; else if (strcmp(name,"FrontalForcingsSdpolyparams")==0) return FrontalForcingsSdpolyparamsEnum; else if (strcmp(name,"GrdModel")==0) return GrdModelEnum; - else if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"GroundinglineMeltInterpolation")==0) return GroundinglineMeltInterpolationEnum; + if (strcmp(name,"GroundinglineFrictionInterpolation")==0) return GroundinglineFrictionInterpolationEnum; + else if (strcmp(name,"GroundinglineMeltInterpolation")==0) return GroundinglineMeltInterpolationEnum; else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum; else if (strcmp(name,"GroundinglineNumRequestedOutputs")==0) return GroundinglineNumRequestedOutputsEnum; else if (strcmp(name,"GroundinglineRequestedOutputs")==0) return GroundinglineRequestedOutputsEnum; @@ -381,11 +382,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"LoveIntegrationScheme")==0) return LoveIntegrationSchemeEnum; else if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum; else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum; - else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"LoveNTemporalIterations")==0) return LoveNTemporalIterationsEnum; + if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; + else if (strcmp(name,"LoveNTemporalIterations")==0) return LoveNTemporalIterationsEnum; else if (strcmp(name,"LoveNYiEquations")==0) return LoveNYiEquationsEnum; else if (strcmp(name,"LoveR0")==0) return LoveR0Enum; else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; @@ -504,11 +505,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum; else if (strcmp(name,"TidalLoveK")==0) return TidalLoveKEnum; else if (strcmp(name,"TidalLoveL")==0) return TidalLoveLEnum; - else if (strcmp(name,"TidalLoveK2Secular")==0) return TidalLoveK2SecularEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"LoadLoveH")==0) return LoadLoveHEnum; + if (strcmp(name,"TidalLoveK2Secular")==0) return TidalLoveK2SecularEnum; + else if (strcmp(name,"LoadLoveH")==0) return LoadLoveHEnum; else if (strcmp(name,"LoadLoveK")==0) return LoadLoveKEnum; else if (strcmp(name,"LoadLoveL")==0) return LoadLoveLEnum; else if (strcmp(name,"LoveTimeFreq")==0) return LoveTimeFreqEnum; @@ -627,11 +628,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"SmbIsmungsm")==0) return SmbIsmungsmEnum; else if (strcmp(name,"SmbIsprecipforcingremapped")==0) return SmbIsprecipforcingremappedEnum; else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum; - else if (strcmp(name,"SmbIssetpddfac")==0) return SmbIssetpddfacEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum; + if (strcmp(name,"SmbIssetpddfac")==0) return SmbIssetpddfacEnum; + else if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum; else if (strcmp(name,"SmbIstemperaturescaled")==0) return SmbIstemperaturescaledEnum; else if (strcmp(name,"SmbIsthermal")==0) return SmbIsthermalEnum; else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum; @@ -750,11 +751,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; else if (strcmp(name,"TransientIsmmemasstransport")==0) return TransientIsmmemasstransportEnum; else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; - else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; + if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; + else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; else if (strcmp(name,"TransientIssampling")==0) return TransientIssamplingEnum; else if (strcmp(name,"TransientIsslc")==0) return TransientIsslcEnum; else if (strcmp(name,"TransientIssmb")==0) return TransientIssmbEnum; @@ -873,11 +874,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"StrOld")==0) return StrOldEnum; else if (strcmp(name,"Str")==0) return StrEnum; else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; - else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; + if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; + else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; @@ -996,11 +997,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; else if (strcmp(name,"Input")==0) return InputEnum; else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; - else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum; + if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; + else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum; else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; @@ -1119,11 +1120,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"SealevelchangeAzimuthIndexHydro")==0) return SealevelchangeAzimuthIndexHydroEnum; else if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum; else if (strcmp(name,"SealevelchangeViscousSG")==0) return SealevelchangeViscousSGEnum; - else if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum; + if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum; + else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum; else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum; else if (strcmp(name,"CouplingTransferCount")==0) return CouplingTransferCountEnum; else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; @@ -1242,11 +1243,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum; else if (strcmp(name,"SmbPrecipitationSubstep")==0) return SmbPrecipitationSubstepEnum; else if (strcmp(name,"SmbPrecipitationsAnomaly")==0) return SmbPrecipitationsAnomalyEnum; - else if (strcmp(name,"SmbDsradiationAnomaly")==0) return SmbDsradiationAnomalyEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"SmbDlradiationAnomaly")==0) return SmbDlradiationAnomalyEnum; + if (strcmp(name,"SmbDsradiationAnomaly")==0) return SmbDsradiationAnomalyEnum; + else if (strcmp(name,"SmbDlradiationAnomaly")==0) return SmbDlradiationAnomalyEnum; else if (strcmp(name,"SmbWindspeedAnomaly")==0) return SmbWindspeedAnomalyEnum; else if (strcmp(name,"SmbAirhumidityAnomaly")==0) return SmbAirhumidityAnomalyEnum; else if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum; @@ -1365,11 +1366,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"VxDebris")==0) return VxDebrisEnum; else if (strcmp(name,"Vx")==0) return VxEnum; else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; - else if (strcmp(name,"VxObs")==0) return VxObsEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"VxShear")==0) return VxShearEnum; + if (strcmp(name,"VxObs")==0) return VxObsEnum; + else if (strcmp(name,"VxShear")==0) return VxShearEnum; else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum; else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; else if (strcmp(name,"VyBase")==0) return VyBaseEnum; @@ -1488,11 +1489,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum; else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum; else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum; - else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; else stage=13; } if(stage==13){ - if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum; + if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; + else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum; else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum; else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum; else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum; @@ -1611,11 +1612,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition213")==0) return Outputdefinition213Enum; else if (strcmp(name,"Outputdefinition214")==0) return Outputdefinition214Enum; else if (strcmp(name,"Outputdefinition215")==0) return Outputdefinition215Enum; - else if (strcmp(name,"Outputdefinition216")==0) return Outputdefinition216Enum; else stage=14; } if(stage==14){ - if (strcmp(name,"Outputdefinition217")==0) return Outputdefinition217Enum; + if (strcmp(name,"Outputdefinition216")==0) return Outputdefinition216Enum; + else if (strcmp(name,"Outputdefinition217")==0) return Outputdefinition217Enum; else if (strcmp(name,"Outputdefinition218")==0) return Outputdefinition218Enum; else if (strcmp(name,"Outputdefinition219")==0) return Outputdefinition219Enum; else if (strcmp(name,"Outputdefinition220")==0) return Outputdefinition220Enum; @@ -1734,11 +1735,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition332")==0) return Outputdefinition332Enum; else if (strcmp(name,"Outputdefinition333")==0) return Outputdefinition333Enum; else if (strcmp(name,"Outputdefinition334")==0) return Outputdefinition334Enum; - else if (strcmp(name,"Outputdefinition335")==0) return Outputdefinition335Enum; else stage=15; } if(stage==15){ - if (strcmp(name,"Outputdefinition336")==0) return Outputdefinition336Enum; + if (strcmp(name,"Outputdefinition335")==0) return Outputdefinition335Enum; + else if (strcmp(name,"Outputdefinition336")==0) return Outputdefinition336Enum; else if (strcmp(name,"Outputdefinition337")==0) return Outputdefinition337Enum; else if (strcmp(name,"Outputdefinition338")==0) return Outputdefinition338Enum; else if (strcmp(name,"Outputdefinition339")==0) return Outputdefinition339Enum; @@ -1857,11 +1858,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition450")==0) return Outputdefinition450Enum; else if (strcmp(name,"Outputdefinition451")==0) return Outputdefinition451Enum; else if (strcmp(name,"Outputdefinition452")==0) return Outputdefinition452Enum; - else if (strcmp(name,"Outputdefinition453")==0) return Outputdefinition453Enum; else stage=16; } if(stage==16){ - if (strcmp(name,"Outputdefinition454")==0) return Outputdefinition454Enum; + if (strcmp(name,"Outputdefinition453")==0) return Outputdefinition453Enum; + else if (strcmp(name,"Outputdefinition454")==0) return Outputdefinition454Enum; else if (strcmp(name,"Outputdefinition455")==0) return Outputdefinition455Enum; else if (strcmp(name,"Outputdefinition456")==0) return Outputdefinition456Enum; else if (strcmp(name,"Outputdefinition457")==0) return Outputdefinition457Enum; @@ -1980,11 +1981,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition569")==0) return Outputdefinition569Enum; else if (strcmp(name,"Outputdefinition506")==0) return Outputdefinition506Enum; else if (strcmp(name,"Outputdefinition570")==0) return Outputdefinition570Enum; - else if (strcmp(name,"Outputdefinition571")==0) return Outputdefinition571Enum; else stage=17; } if(stage==17){ - if (strcmp(name,"Outputdefinition572")==0) return Outputdefinition572Enum; + if (strcmp(name,"Outputdefinition571")==0) return Outputdefinition571Enum; + else if (strcmp(name,"Outputdefinition572")==0) return Outputdefinition572Enum; else if (strcmp(name,"Outputdefinition573")==0) return Outputdefinition573Enum; else if (strcmp(name,"Outputdefinition574")==0) return Outputdefinition574Enum; else if (strcmp(name,"Outputdefinition575")==0) return Outputdefinition575Enum; @@ -2103,11 +2104,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition687")==0) return Outputdefinition687Enum; else if (strcmp(name,"Outputdefinition688")==0) return Outputdefinition688Enum; else if (strcmp(name,"Outputdefinition689")==0) return Outputdefinition689Enum; - else if (strcmp(name,"Outputdefinition608")==0) return Outputdefinition608Enum; else stage=18; } if(stage==18){ - if (strcmp(name,"Outputdefinition690")==0) return Outputdefinition690Enum; + if (strcmp(name,"Outputdefinition608")==0) return Outputdefinition608Enum; + else if (strcmp(name,"Outputdefinition690")==0) return Outputdefinition690Enum; else if (strcmp(name,"Outputdefinition691")==0) return Outputdefinition691Enum; else if (strcmp(name,"Outputdefinition692")==0) return Outputdefinition692Enum; else if (strcmp(name,"Outputdefinition693")==0) return Outputdefinition693Enum; @@ -2226,11 +2227,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition813")==0) return Outputdefinition813Enum; else if (strcmp(name,"Outputdefinition814")==0) return Outputdefinition814Enum; else if (strcmp(name,"Outputdefinition815")==0) return Outputdefinition815Enum; - else if (strcmp(name,"Outputdefinition816")==0) return Outputdefinition816Enum; else stage=19; } if(stage==19){ - if (strcmp(name,"Outputdefinition817")==0) return Outputdefinition817Enum; + if (strcmp(name,"Outputdefinition816")==0) return Outputdefinition816Enum; + else if (strcmp(name,"Outputdefinition817")==0) return Outputdefinition817Enum; else if (strcmp(name,"Outputdefinition818")==0) return Outputdefinition818Enum; else if (strcmp(name,"Outputdefinition819")==0) return Outputdefinition819Enum; else if (strcmp(name,"Outputdefinition820")==0) return Outputdefinition820Enum; @@ -2349,11 +2350,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition932")==0) return Outputdefinition932Enum; else if (strcmp(name,"Outputdefinition933")==0) return Outputdefinition933Enum; else if (strcmp(name,"Outputdefinition934")==0) return Outputdefinition934Enum; - else if (strcmp(name,"Outputdefinition935")==0) return Outputdefinition935Enum; else stage=20; } if(stage==20){ - if (strcmp(name,"Outputdefinition936")==0) return Outputdefinition936Enum; + if (strcmp(name,"Outputdefinition935")==0) return Outputdefinition935Enum; + else if (strcmp(name,"Outputdefinition936")==0) return Outputdefinition936Enum; else if (strcmp(name,"Outputdefinition937")==0) return Outputdefinition937Enum; else if (strcmp(name,"Outputdefinition938")==0) return Outputdefinition938Enum; else if (strcmp(name,"Outputdefinition939")==0) return Outputdefinition939Enum; @@ -2472,11 +2473,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1050")==0) return Outputdefinition1050Enum; else if (strcmp(name,"Outputdefinition1051")==0) return Outputdefinition1051Enum; else if (strcmp(name,"Outputdefinition1052")==0) return Outputdefinition1052Enum; - else if (strcmp(name,"Outputdefinition1053")==0) return Outputdefinition1053Enum; else stage=21; } if(stage==21){ - if (strcmp(name,"Outputdefinition1054")==0) return Outputdefinition1054Enum; + if (strcmp(name,"Outputdefinition1053")==0) return Outputdefinition1053Enum; + else if (strcmp(name,"Outputdefinition1054")==0) return Outputdefinition1054Enum; else if (strcmp(name,"Outputdefinition1055")==0) return Outputdefinition1055Enum; else if (strcmp(name,"Outputdefinition1056")==0) return Outputdefinition1056Enum; else if (strcmp(name,"Outputdefinition1057")==0) return Outputdefinition1057Enum; @@ -2595,11 +2596,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1169")==0) return Outputdefinition1169Enum; else if (strcmp(name,"Outputdefinition1106")==0) return Outputdefinition1106Enum; else if (strcmp(name,"Outputdefinition1170")==0) return Outputdefinition1170Enum; - else if (strcmp(name,"Outputdefinition1171")==0) return Outputdefinition1171Enum; else stage=22; } if(stage==22){ - if (strcmp(name,"Outputdefinition1172")==0) return Outputdefinition1172Enum; + if (strcmp(name,"Outputdefinition1171")==0) return Outputdefinition1171Enum; + else if (strcmp(name,"Outputdefinition1172")==0) return Outputdefinition1172Enum; else if (strcmp(name,"Outputdefinition1173")==0) return Outputdefinition1173Enum; else if (strcmp(name,"Outputdefinition1174")==0) return Outputdefinition1174Enum; else if (strcmp(name,"Outputdefinition1175")==0) return Outputdefinition1175Enum; @@ -2718,11 +2719,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1287")==0) return Outputdefinition1287Enum; else if (strcmp(name,"Outputdefinition1288")==0) return Outputdefinition1288Enum; else if (strcmp(name,"Outputdefinition1289")==0) return Outputdefinition1289Enum; - else if (strcmp(name,"Outputdefinition1208")==0) return Outputdefinition1208Enum; else stage=23; } if(stage==23){ - if (strcmp(name,"Outputdefinition1290")==0) return Outputdefinition1290Enum; + if (strcmp(name,"Outputdefinition1208")==0) return Outputdefinition1208Enum; + else if (strcmp(name,"Outputdefinition1290")==0) return Outputdefinition1290Enum; else if (strcmp(name,"Outputdefinition1291")==0) return Outputdefinition1291Enum; else if (strcmp(name,"Outputdefinition1292")==0) return Outputdefinition1292Enum; else if (strcmp(name,"Outputdefinition1293")==0) return Outputdefinition1293Enum; @@ -2841,11 +2842,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1413")==0) return Outputdefinition1413Enum; else if (strcmp(name,"Outputdefinition1414")==0) return Outputdefinition1414Enum; else if (strcmp(name,"Outputdefinition1415")==0) return Outputdefinition1415Enum; - else if (strcmp(name,"Outputdefinition1416")==0) return Outputdefinition1416Enum; else stage=24; } if(stage==24){ - if (strcmp(name,"Outputdefinition1417")==0) return Outputdefinition1417Enum; + if (strcmp(name,"Outputdefinition1416")==0) return Outputdefinition1416Enum; + else if (strcmp(name,"Outputdefinition1417")==0) return Outputdefinition1417Enum; else if (strcmp(name,"Outputdefinition1418")==0) return Outputdefinition1418Enum; else if (strcmp(name,"Outputdefinition1419")==0) return Outputdefinition1419Enum; else if (strcmp(name,"Outputdefinition1420")==0) return Outputdefinition1420Enum; @@ -2964,11 +2965,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1532")==0) return Outputdefinition1532Enum; else if (strcmp(name,"Outputdefinition1533")==0) return Outputdefinition1533Enum; else if (strcmp(name,"Outputdefinition1534")==0) return Outputdefinition1534Enum; - else if (strcmp(name,"Outputdefinition1535")==0) return Outputdefinition1535Enum; else stage=25; } if(stage==25){ - if (strcmp(name,"Outputdefinition1536")==0) return Outputdefinition1536Enum; + if (strcmp(name,"Outputdefinition1535")==0) return Outputdefinition1535Enum; + else if (strcmp(name,"Outputdefinition1536")==0) return Outputdefinition1536Enum; else if (strcmp(name,"Outputdefinition1537")==0) return Outputdefinition1537Enum; else if (strcmp(name,"Outputdefinition1538")==0) return Outputdefinition1538Enum; else if (strcmp(name,"Outputdefinition1539")==0) return Outputdefinition1539Enum; @@ -3087,11 +3088,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1650")==0) return Outputdefinition1650Enum; else if (strcmp(name,"Outputdefinition1651")==0) return Outputdefinition1651Enum; else if (strcmp(name,"Outputdefinition1652")==0) return Outputdefinition1652Enum; - else if (strcmp(name,"Outputdefinition1653")==0) return Outputdefinition1653Enum; else stage=26; } if(stage==26){ - if (strcmp(name,"Outputdefinition1654")==0) return Outputdefinition1654Enum; + if (strcmp(name,"Outputdefinition1653")==0) return Outputdefinition1653Enum; + else if (strcmp(name,"Outputdefinition1654")==0) return Outputdefinition1654Enum; else if (strcmp(name,"Outputdefinition1655")==0) return Outputdefinition1655Enum; else if (strcmp(name,"Outputdefinition1656")==0) return Outputdefinition1656Enum; else if (strcmp(name,"Outputdefinition1657")==0) return Outputdefinition1657Enum; @@ -3210,11 +3211,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1769")==0) return Outputdefinition1769Enum; else if (strcmp(name,"Outputdefinition1706")==0) return Outputdefinition1706Enum; else if (strcmp(name,"Outputdefinition1770")==0) return Outputdefinition1770Enum; - else if (strcmp(name,"Outputdefinition1771")==0) return Outputdefinition1771Enum; else stage=27; } if(stage==27){ - if (strcmp(name,"Outputdefinition1772")==0) return Outputdefinition1772Enum; + if (strcmp(name,"Outputdefinition1771")==0) return Outputdefinition1771Enum; + else if (strcmp(name,"Outputdefinition1772")==0) return Outputdefinition1772Enum; else if (strcmp(name,"Outputdefinition1773")==0) return Outputdefinition1773Enum; else if (strcmp(name,"Outputdefinition1774")==0) return Outputdefinition1774Enum; else if (strcmp(name,"Outputdefinition1775")==0) return Outputdefinition1775Enum; @@ -3333,11 +3334,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1887")==0) return Outputdefinition1887Enum; else if (strcmp(name,"Outputdefinition1888")==0) return Outputdefinition1888Enum; else if (strcmp(name,"Outputdefinition1889")==0) return Outputdefinition1889Enum; - else if (strcmp(name,"Outputdefinition1808")==0) return Outputdefinition1808Enum; else stage=28; } if(stage==28){ - if (strcmp(name,"Outputdefinition1890")==0) return Outputdefinition1890Enum; + if (strcmp(name,"Outputdefinition1808")==0) return Outputdefinition1808Enum; + else if (strcmp(name,"Outputdefinition1890")==0) return Outputdefinition1890Enum; else if (strcmp(name,"Outputdefinition1891")==0) return Outputdefinition1891Enum; else if (strcmp(name,"Outputdefinition1892")==0) return Outputdefinition1892Enum; else if (strcmp(name,"Outputdefinition1893")==0) return Outputdefinition1893Enum; @@ -3456,11 +3457,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; else if (strcmp(name,"AgeAnalysis")==0) return AgeAnalysisEnum; - else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum; else stage=29; } if(stage==29){ - if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum; + if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum; + else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum; else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum; else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; else if (strcmp(name,"AutodiffJacobian")==0) return AutodiffJacobianEnum; @@ -3579,11 +3580,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; - else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; else stage=30; } if(stage==30){ - if (strcmp(name,"GenericParam")==0) return GenericParamEnum; + if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; + else if (strcmp(name,"GenericParam")==0) return GenericParamEnum; else if (strcmp(name,"GenericExternalResult")==0) return GenericExternalResultEnum; else if (strcmp(name,"Gradient1")==0) return Gradient1Enum; else if (strcmp(name,"Gradient2")==0) return Gradient2Enum; @@ -3702,11 +3703,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"MeshY")==0) return MeshYEnum; else if (strcmp(name,"MinVel")==0) return MinVelEnum; else if (strcmp(name,"MinVx")==0) return MinVxEnum; - else if (strcmp(name,"MinVy")==0) return MinVyEnum; else stage=31; } if(stage==31){ - if (strcmp(name,"MinVz")==0) return MinVzEnum; + if (strcmp(name,"MinVy")==0) return MinVyEnum; + else if (strcmp(name,"MinVz")==0) return MinVzEnum; else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum; else if (strcmp(name,"Moulin")==0) return MoulinEnum; else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; @@ -3825,11 +3826,11 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum; else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; - else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; else stage=32; } if(stage==32){ - if (strcmp(name,"Tetra")==0) return TetraEnum; + if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; + else if (strcmp(name,"Tetra")==0) return TetraEnum; else if (strcmp(name,"TetraInput")==0) return TetraInputEnum; else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum; else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum; diff --git a/src/c/shared/Enum/issmenums.jl b/src/c/shared/Enum/issmenums.jl index 666e47e26..3088a80c9 100644 --- a/src/c/shared/Enum/issmenums.jl +++ b/src/c/shared/Enum/issmenums.jl @@ -210,6 +210,7 @@ FrictionGammaEnum FrictionLawEnum FrictionLinearizeEnum + FrictionMaxCoefficientEnum FrictionPseudoplasticityExponentEnum FrictionU0Enum FrictionThresholdSpeedEnum @@ -3986,6 +3987,7 @@ function EnumToString(enum::IssmEnum) if(enum==FrictionGammaEnum) return "FrictionGamma" end if(enum==FrictionLawEnum) return "FrictionLaw" end if(enum==FrictionLinearizeEnum) return "FrictionLinearize" end + if(enum==FrictionMaxCoefficientEnum) return "FrictionMaxCoefficient" end if(enum==FrictionPseudoplasticityExponentEnum) return "FrictionPseudoplasticityExponent" end if(enum==FrictionU0Enum) return "FrictionU0" end if(enum==FrictionThresholdSpeedEnum) return "FrictionThresholdSpeed" end @@ -7762,6 +7764,7 @@ function StringToEnum(name::String) if(name=="FrictionGamma") return FrictionGammaEnum end if(name=="FrictionLaw") return FrictionLawEnum end if(name=="FrictionLinearize") return FrictionLinearizeEnum end + if(name=="FrictionMaxCoefficient") return FrictionMaxCoefficientEnum end if(name=="FrictionPseudoplasticityExponent") return FrictionPseudoplasticityExponentEnum end if(name=="FrictionU0") return FrictionU0Enum end if(name=="FrictionThresholdSpeed") return FrictionThresholdSpeedEnum end diff --git a/src/m/classes/frictionjosh.m b/src/m/classes/frictionjosh.m index 011e468b3..2445ebc3f 100644 --- a/src/m/classes/frictionjosh.m +++ b/src/m/classes/frictionjosh.m @@ -9,6 +9,7 @@ pressure_adjusted_temperature = NaN; gamma = 0.; effective_pressure_limit = 0; + coefficient_max = 0; end methods function self = extrude(self,md) % {{{ @@ -30,6 +31,9 @@ %Default gamma: 1 self.gamma = 1.; + % Default max friction coefficient: 300 + self.coefficient_max = 300; + %Default 0 self.effective_pressure_limit = 0; @@ -43,6 +47,7 @@ md = checkfield(md,'fieldname','friction.pressure_adjusted_temperature','NaN',1,'Inf',1); md = checkfield(md,'fieldname','friction.gamma','numel',1,'NaN',1,'Inf',1,'>',0.); md = checkfield(md,'fieldname','friction.effective_pressure_limit','numel',[1],'>=',0); + md = checkfield(md,'fieldname','friction.coefficient_max','numel',1,'NaN',1,'Inf',1,'>',0.); %Check that temperature is provided md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'Inf',1,'size','universal'); @@ -53,6 +58,7 @@ function disp(self) % {{{ fielddisplay(self,'pressure_adjusted_temperature','friction pressure_adjusted_temperature (T - Tpmp) [K]'); fielddisplay(self,'gamma','(T - Tpmp)/gamma [K]'); fielddisplay(self,'effective_pressure_limit','Neff do not allow to fall below a certain limit: effective_pressure_limit*rho_ice*g*thickness (default 0)'); + fielddisplay(self,'coefficient_max', 'effective friction C = min(coefficient_max, sqrt(exp(T_b(modern) - T_b(t))/gamma) * coefficient)'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -61,6 +67,7 @@ function marshall(self,prefix,md,fid) % {{{ WriteData(fid,prefix,'class','friction','object',self,'fieldname','pressure_adjusted_temperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'class','friction','object',self,'fieldname','gamma','format','Double'); WriteData(fid,prefix,'object',self,'class','friction','fieldname','effective_pressure_limit','format','Double'); + WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient_max','format','Double'); end % }}} end end diff --git a/src/m/classes/model.m b/src/m/classes/model.m index 1586acabc..ee3953e0c 100644 --- a/src/m/classes/model.m +++ b/src/m/classes/model.m @@ -198,7 +198,9 @@ %2022 Oct 28 if ~isa(md.debris,'debris'); md.debris=debris(); end %Mmetransport: Jun 2022: - if ~isa(md.mmemasstransport,'mmemasstransport'); md.mmemasstransport=mmemasstransport(); end; + if ~isa(md.mmemasstransport,'mmemasstransport'); md.mmemasstransport=mmemasstransport(); end + % 2026 February 18 + if isa(md.friction, 'frictionjosh') && md.friction.coefficient_max==0; md.friction.coefficient_max=300; end end% }}} end methods @@ -598,7 +600,7 @@ function disp(self) % {{{ %copy model md1=md; - %recover options: + %recover optoins: options=pairoptions(varargin{:}); %some checks @@ -662,13 +664,13 @@ function disp(self) % {{{ %loop over model fields model_fields=fields(md1); - for i=1:length(model_fields) + for i=1:length(model_fields), %get field field=md1.(model_fields{i}); fieldsize=size(field); if isobject(field), %recursive call object_fields=fields(md1.(model_fields{i})); - for j=1:length(object_fields) + for j=1:length(object_fields), %get field field=md1.(model_fields{i}).(object_fields{j}); fieldsize=size(field); @@ -677,7 +679,7 @@ function disp(self) % {{{ md2.(model_fields{i}).(object_fields{j})=field(pos_node,:); elseif (fieldsize(1)==numberofvertices1+1) md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n + %size = number of elements * n elseif fieldsize(1)==numberofelements1 md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:); elseif (fieldsize(1)==numberofelements1+1) @@ -690,7 +692,7 @@ function disp(self) % {{{ md2.(model_fields{i})=field(pos_node,:); elseif (fieldsize(1)==numberofvertices1+1) md2.(model_fields{i})=[field(pos_node,:); field(end,:)]; - %size = number of elements * n + %size = number of elements * n elseif fieldsize(1)==numberofelements1 md2.(model_fields{i})=field(pos_elem,:); elseif (fieldsize(1)==numberofelements1+1) @@ -706,13 +708,6 @@ function disp(self) % {{{ md2.mesh.numberofvertices=numberofvertices2; md2.mesh.elements=elements_2; - % Extract ISMIP6 basal tf field - if isa(md1.basalforcings, 'basalforcingsismip6') - for i=1:numel(md.basalforcings.tf) - md2.basalforcings.tf{i} = [md1.basalforcings.tf{i}(pos_node); md1.basalforcings.tf{i}(end)]; - end - end - %mesh.uppervertex mesh.lowervertex if isa(md1.mesh,'mesh3dprisms'), md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node); From 38009b926a9b8131ee92f608cd7f964c387cd4a7 Mon Sep 17 00:00:00 2001 From: mmorligh Date: Thu, 19 Feb 2026 09:24:50 -0500 Subject: [PATCH 05/17] CHG: translated to python --- src/m/classes/frictionjosh.m | 24 ++++++++++---------- src/m/classes/frictionjosh.py | 41 ++++++++++++++++++++--------------- 2 files changed, 36 insertions(+), 29 deletions(-) diff --git a/src/m/classes/frictionjosh.m b/src/m/classes/frictionjosh.m index 2445ebc3f..57686ff68 100644 --- a/src/m/classes/frictionjosh.m +++ b/src/m/classes/frictionjosh.m @@ -5,11 +5,11 @@ classdef frictionjosh properties (SetAccess=public) - coefficient = NaN; + coefficient = NaN; pressure_adjusted_temperature = NaN; - gamma = 0.; - effective_pressure_limit = 0; - coefficient_max = 0; + gamma = 0.0; + effective_pressure_limit = 0.0; + coefficient_max = 0.0; end methods function self = extrude(self,md) % {{{ @@ -32,10 +32,10 @@ self.gamma = 1.; % Default max friction coefficient: 300 - self.coefficient_max = 300; + self.coefficient_max = 300.0; %Default 0 - self.effective_pressure_limit = 0; + self.effective_pressure_limit = 0.0; end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -62,12 +62,12 @@ function disp(self) % {{{ end % }}} function marshall(self,prefix,md,fid) % {{{ - WriteData(fid,prefix,'name','md.friction.law','data',9,'format','Integer'); - WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'class','friction','object',self,'fieldname','pressure_adjusted_temperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'class','friction','object',self,'fieldname','gamma','format','Double'); - WriteData(fid,prefix,'object',self,'class','friction','fieldname','effective_pressure_limit','format','Double'); - WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient_max','format','Double'); + WriteData(fid,prefix, 'name', 'md.friction.law', 'data',9, 'format', 'Integer'); + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype',1, 'timeserieslength',md.mesh.numberofvertices+1, 'yts',md.constants.yts); + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'pressure_adjusted_temperature', 'format', 'DoubleMat', 'mattype',1, 'timeserieslength',md.mesh.numberofvertices+1, 'yts',md.constants.yts); + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'gamma', 'format', 'Double'); + WriteData(fid,prefix, 'object',self, 'class', 'friction', 'fieldname', 'effective_pressure_limit', 'format', 'Double'); + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'coefficient_max', 'format', 'Double'); end % }}} end end diff --git a/src/m/classes/frictionjosh.py b/src/m/classes/frictionjosh.py index a333bc361..45a56f9cc 100644 --- a/src/m/classes/frictionjosh.py +++ b/src/m/classes/frictionjosh.py @@ -14,13 +14,13 @@ class frictionjosh(object): """ def __init__(self): # {{{ - self.coefficient = np.nan + self.coefficient = np.nan self.pressure_adjusted_temperature = np.nan - self.gamma = 0 - self.effective_pressure_limit = 0 + self.gamma = 0. + self.effective_pressure_limit = 0. + self.coefficient_max = 0. self.setdefaultparameters() - #self.requested_outputs = [] # }}} def __repr__(self): # {{{ @@ -30,7 +30,7 @@ def __repr__(self): # {{{ s += '{}\n'.format(fielddisplay(self, 'pressure_adjusted_temperature', 'friction pressure_adjusted_temperature (T - Tpmp) [K]')) s += '{}\n'.format(fielddisplay(self, 'gamma', '(T - Tpmp)/gamma [K]')) s += '{}\n'.format(fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)')) - #s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested')) + s += '{}\n'.format(fielddisplay(self, 'coefficient_max', 'effective friction C = min(coefficient_max, sqrt(exp(T_b(modern) - T_b(t))/gamma) * coefficient)')) return s # }}} @@ -41,34 +41,41 @@ def extrude(self, md): # {{{ # }}} def setdefaultparameters(self): # {{{ + + #Default gamma: 1 self.gamma = 1. - #self.requested_outputs = ['default'] - self.effective_pressure_limit = 0 - return self - # }}} - def defaultoutputs(self, md): # {{{ - list = [] - return list + #Default 0. + self.effective_pressure_limit = 0.0 + + #Default max friction coefficient: 300 + self.coefficient_max = 300. + + return self # }}} def checkconsistency(self, md, solution, analyses): # {{{ + # Early return if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses: return md + md = checkfield(md, 'fieldname', 'friction.coefficient', 'timeseries', 1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'friction.pressure_adjusted_temperature','NaN',1,'Inf',1) md = checkfield(md, 'fieldname', 'friction.gamma','numel',1,'NaN',1,'Inf',1,'>',0.) md = checkfield(md, 'fieldname', 'friction.effective_pressure_limit', 'numel', [1], '>=', 0) + md = checkfield(md,'fieldname', 'friction.coefficient_max', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0.) + # Check that temperature is provided md = checkfield(md,'fieldname','initialization.temperature','NaN',1,'Inf',1,'size','universal') return md # }}} def marshall(self, prefix, md, fid): # {{{ - WriteData(fid,prefix,'name','md.friction.law','data',9,'format','Integer') - WriteData(fid,prefix,'class','friction','object',self,'fieldname','coefficient','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - WriteData(fid,prefix,'class','friction','object',self,'fieldname','pressure_adjusted_temperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) - WriteData(fid,prefix,'class','friction','object',self,'fieldname','gamma','format','Double') - WriteData(fid,prefix,'object',self,'class','friction','fieldname','effective_pressure_limit','format','Double') + WriteData(fid,prefix, 'name', 'md.friction.law', 'data',9, 'format', 'Integer') + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'coefficient', 'format', 'DoubleMat', 'mattype',1, 'timeserieslength',md.mesh.numberofvertices+1, 'yts',md.constants.yts) + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'pressure_adjusted_temperature', 'format', 'DoubleMat', 'mattype',1, 'timeserieslength',md.mesh.numberofvertices+1, 'yts',md.constants.yts) + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'gamma', 'format', 'Double') + WriteData(fid,prefix, 'object',self, 'class', 'friction', 'fieldname', 'effective_pressure_limit', 'format', 'Double') + WriteData(fid,prefix, 'class', 'friction', 'object',self, 'fieldname', 'coefficient_max', 'format', 'Double') # }}} From e664ae1a6aabec973e4fd01fcfcb5a4fb16d48fc Mon Sep 17 00:00:00 2001 From: Alicia Gunilla Maila Bratner Date: Fri, 20 Feb 2026 10:58:28 -0500 Subject: [PATCH 06/17] smbpddsicopolis: melt factors from user input + rlaps from user input (not overwritten by sinusoidal) --- src/c/analyses/SmbAnalysis.cpp | 2 + src/c/classes/Elements/Element.cpp | 7 +- .../PddSurfaceMassBalanceSicopolis.cpp | 20 +-- src/c/shared/Elements/elements.h | 2 +- src/c/shared/Enum/Enum.vim | 2 + src/c/shared/Enum/EnumDefinitions.h | 2 + src/c/shared/Enum/EnumToStringx.cpp | 2 + src/c/shared/Enum/Enumjl.vim | 2 + src/c/shared/Enum/StringToEnumx.cpp | 170 +++++++++--------- src/c/shared/Enum/issmenums.jl | 6 + src/m/classes/SMBpddSicopolis.m | 10 ++ src/m/classes/SMBpddSicopolis.py | 12 +- src/m/classes/model.m | 3 + 13 files changed, 141 insertions(+), 99 deletions(-) diff --git a/src/c/analyses/SmbAnalysis.cpp b/src/c/analyses/SmbAnalysis.cpp index 3c17be2be..02052f30d 100644 --- a/src/c/analyses/SmbAnalysis.cpp +++ b/src/c/analyses/SmbAnalysis.cpp @@ -338,6 +338,8 @@ void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s parameters->AddObject(iomodel->CopyConstantObject("md.smb.isfirnwarming",SmbIsfirnwarmingEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.pdd_fac_ice",PddfacIceEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.pdd_fac_snow",PddfacSnowEnum)); break; case SMBdebrisEvattEnum: parameters->AddObject(iomodel->CopyConstantObject("md.smb.qlaps",SmbDesfacEnum)); diff --git a/src/c/classes/Elements/Element.cpp b/src/c/classes/Elements/Element.cpp index db7206881..b0d1be170 100644 --- a/src/c/classes/Elements/Element.cpp +++ b/src/c/classes/Elements/Element.cpp @@ -3695,6 +3695,7 @@ void Element::PositiveDegreeDaySicopolis(bool isfirnwarming){/*{{{*/ IssmDouble* p_ampl=xNew(NUM_VERTICES); // precip anomaly IssmDouble* t_ampl=xNew(NUM_VERTICES); // remperature anomaly IssmDouble rho_water,rho_ice,desfac,rlaps; + IssmDouble pdd_fac_ice,pdd_fac_snow; IssmDouble inv_twelve=1./12.; //factor for monthly average IssmDouble time,yts,time_yr; @@ -3709,6 +3710,10 @@ void Element::PositiveDegreeDaySicopolis(bool isfirnwarming){/*{{{*/ desfac=this->FindParam(SmbDesfacEnum); rlaps=this->FindParam(SmbRlapsEnum); + /*Get pdd melt factors*/ + pdd_fac_ice=this->FindParam(PddfacIceEnum); + pdd_fac_snow=this->FindParam(PddfacSnowEnum); + /* Get time */ this->parameters->FindParam(&time,TimeEnum); this->parameters->FindParam(&yts,ConstantsYtsEnum); @@ -3747,7 +3752,7 @@ void Element::PositiveDegreeDaySicopolis(bool isfirnwarming){/*{{{*/ for (int iv = 0; iv (m IE)/(a*deg C) - beta2=beta2*(0.001*365)*sconv; // (mm WE)/(d*deg C) --> (m IE)/(a*deg C) + pdd_fac_snow=pdd_fac_snow*(0.001*365)*sconv; // (mm WE)/(d*deg C) --> (m IE)/(a*deg C) + pdd_fac_ice=pdd_fac_ice*(0.001*365)*sconv; // (mm WE)/(d*deg C) --> (m IE)/(a*deg C) /* initalize fields */ precip=0.0; @@ -49,8 +46,7 @@ IssmDouble PddSurfaceMassBalanceSicopolis(IssmDouble* monthlytemperatures, IssmD /********* Surface temperature correction *******/ st=(s-s0t)/1000.; - // FIXME rlaps ?? - rlaps=-6.309e-03+(-5.426e-03-(-6.309e-03))*sin((imonth+1-4)*PI/6.0)*1000.0; + /******** Monhtly temperature correction *******/ monthlytemperatures[imonth]=monthlytemperatures[imonth]-rlaps*st;//*max(st,1e-3); tstar=monthlytemperatures[imonth]+t_ampl[0]; @@ -109,20 +105,20 @@ IssmDouble PddSurfaceMassBalanceSicopolis(IssmDouble* monthlytemperatures, IssmD if(rainfall<0.0) rainfall=0.0; // negative values if(rainfall<=(Pmax*snowfall)){ - if((rainfall+beta1*pdd)<=(Pmax*snowfall)) { - smelt_star = rainfall+beta1*pdd; + if((rainfall+pdd_fac_snow*pdd)<=(Pmax*snowfall)) { + smelt_star = rainfall+pdd_fac_snow*pdd; smelt = 0.0; runoff = smelt; } else{ smelt_star = Pmax*snowfall; - smelt = beta2*(pdd-(smelt_star-rainfall)/beta1); + smelt = pdd_fac_ice*(pdd-(smelt_star-rainfall)/pdd_fac_snow); runoff = smelt; } } else{ smelt_star = Pmax*snowfall; - smelt = beta2*pdd; + smelt = pdd_fac_ice*pdd; runoff = smelt+rainfall-Pmax*snowfall; } diff --git a/src/c/shared/Elements/elements.h b/src/c/shared/Elements/elements.h index 5f66ef92a..f11a8677b 100644 --- a/src/c/shared/Elements/elements.h +++ b/src/c/shared/Elements/elements.h @@ -28,7 +28,7 @@ IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* m IssmDouble PddSurfaceMassBalanceSicopolis(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, IssmDouble* melt, IssmDouble* accu, IssmDouble* melt_star, IssmDouble* t_ampl, IssmDouble* p_ampl, IssmDouble yts, IssmDouble s, IssmDouble desfac,IssmDouble s0t, - IssmDouble s0p, IssmDouble rlaps, IssmDouble rho_water, IssmDouble rho_ice); + IssmDouble s0p, IssmDouble rlaps, IssmDouble rho_water, IssmDouble rho_ice, IssmDouble pdd_fac_ice, IssmDouble pdd_fac_snow); void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime, IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, IssmDouble* PrecipitationsPresentday, diff --git a/src/c/shared/Enum/Enum.vim b/src/c/shared/Enum/Enum.vim index 8ee80ebae..bb95a37d1 100644 --- a/src/c/shared/Enum/Enum.vim +++ b/src/c/shared/Enum/Enum.vim @@ -434,6 +434,8 @@ syn keyword cConstant OutputBufferSizePointerEnum syn keyword cConstant OutputFileNameEnum syn keyword cConstant OutputFilePointerEnum syn keyword cConstant OutputdefinitionEnum +syn keyword cConstant PddfacIceEnum +syn keyword cConstant PddfacSnowEnum syn keyword cConstant QmuErrNameEnum syn keyword cConstant QmuInNameEnum syn keyword cConstant QmuIsdakotaEnum diff --git a/src/c/shared/Enum/EnumDefinitions.h b/src/c/shared/Enum/EnumDefinitions.h index a1b15f7d0..f941e36b5 100644 --- a/src/c/shared/Enum/EnumDefinitions.h +++ b/src/c/shared/Enum/EnumDefinitions.h @@ -428,6 +428,8 @@ enum definitions{ OutputFileNameEnum, OutputFilePointerEnum, OutputdefinitionEnum, + PddfacIceEnum, + PddfacSnowEnum, QmuErrNameEnum, QmuInNameEnum, QmuIsdakotaEnum, diff --git a/src/c/shared/Enum/EnumToStringx.cpp b/src/c/shared/Enum/EnumToStringx.cpp index ab845e22c..14e78bb32 100644 --- a/src/c/shared/Enum/EnumToStringx.cpp +++ b/src/c/shared/Enum/EnumToStringx.cpp @@ -436,6 +436,8 @@ const char* EnumToStringx(int en){ case OutputFileNameEnum : return "OutputFileName"; case OutputFilePointerEnum : return "OutputFilePointer"; case OutputdefinitionEnum : return "Outputdefinition"; + case PddfacIceEnum : return "PddfacIce"; + case PddfacSnowEnum : return "PddfacSnow"; case QmuErrNameEnum : return "QmuErrName"; case QmuInNameEnum : return "QmuInName"; case QmuIsdakotaEnum : return "QmuIsdakota"; diff --git a/src/c/shared/Enum/Enumjl.vim b/src/c/shared/Enum/Enumjl.vim index 2245470f3..ce3c81591 100644 --- a/src/c/shared/Enum/Enumjl.vim +++ b/src/c/shared/Enum/Enumjl.vim @@ -427,6 +427,8 @@ syn keyword juliaConstC OutputBufferSizePointerEnum syn keyword juliaConstC OutputFileNameEnum syn keyword juliaConstC OutputFilePointerEnum syn keyword juliaConstC OutputdefinitionEnum +syn keyword juliaConstC PddfacIceEnum +syn keyword juliaConstC PddfacSnowEnum syn keyword juliaConstC QmuErrNameEnum syn keyword juliaConstC QmuInNameEnum syn keyword juliaConstC QmuIsdakotaEnum diff --git a/src/c/shared/Enum/StringToEnumx.cpp b/src/c/shared/Enum/StringToEnumx.cpp index ec30a99fd..b75f68d5c 100644 --- a/src/c/shared/Enum/StringToEnumx.cpp +++ b/src/c/shared/Enum/StringToEnumx.cpp @@ -445,6 +445,8 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum; else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; + else if (strcmp(name,"PddfacIce")==0) return PddfacIceEnum; + else if (strcmp(name,"PddfacSnow")==0) return PddfacSnowEnum; else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum; else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum; else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum; @@ -503,12 +505,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"LovePolarMotionTransferFunctionColinear")==0) return LovePolarMotionTransferFunctionColinearEnum; else if (strcmp(name,"LovePolarMotionTransferFunctionOrthogonal")==0) return LovePolarMotionTransferFunctionOrthogonalEnum; else if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum; - else if (strcmp(name,"TidalLoveK")==0) return TidalLoveKEnum; - else if (strcmp(name,"TidalLoveL")==0) return TidalLoveLEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"TidalLoveK2Secular")==0) return TidalLoveK2SecularEnum; + if (strcmp(name,"TidalLoveK")==0) return TidalLoveKEnum; + else if (strcmp(name,"TidalLoveL")==0) return TidalLoveLEnum; + else if (strcmp(name,"TidalLoveK2Secular")==0) return TidalLoveK2SecularEnum; else if (strcmp(name,"LoadLoveH")==0) return LoadLoveHEnum; else if (strcmp(name,"LoadLoveK")==0) return LoadLoveKEnum; else if (strcmp(name,"LoadLoveL")==0) return LoadLoveLEnum; @@ -626,12 +628,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"SmbIsmappedforcing")==0) return SmbIsmappedforcingEnum; else if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum; else if (strcmp(name,"SmbIsmungsm")==0) return SmbIsmungsmEnum; - else if (strcmp(name,"SmbIsprecipforcingremapped")==0) return SmbIsprecipforcingremappedEnum; - else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"SmbIssetpddfac")==0) return SmbIssetpddfacEnum; + if (strcmp(name,"SmbIsprecipforcingremapped")==0) return SmbIsprecipforcingremappedEnum; + else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum; + else if (strcmp(name,"SmbIssetpddfac")==0) return SmbIssetpddfacEnum; else if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum; else if (strcmp(name,"SmbIstemperaturescaled")==0) return SmbIstemperaturescaledEnum; else if (strcmp(name,"SmbIsthermal")==0) return SmbIsthermalEnum; @@ -749,12 +751,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; - else if (strcmp(name,"TransientIsmmemasstransport")==0) return TransientIsmmemasstransportEnum; - else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; + if (strcmp(name,"TransientIsmmemasstransport")==0) return TransientIsmmemasstransportEnum; + else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; + else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; else if (strcmp(name,"TransientIssampling")==0) return TransientIssamplingEnum; else if (strcmp(name,"TransientIsslc")==0) return TransientIsslcEnum; @@ -872,12 +874,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Dsl")==0) return DslEnum; else if (strcmp(name,"DeltaStr")==0) return DeltaStrEnum; else if (strcmp(name,"StrOld")==0) return StrOldEnum; - else if (strcmp(name,"Str")==0) return StrEnum; - else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; + if (strcmp(name,"Str")==0) return StrEnum; + else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; + else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; @@ -995,12 +997,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"DebrisMaskNodeActivation")==0) return DebrisMaskNodeActivationEnum; else if (strcmp(name,"Ice")==0) return IceEnum; else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; - else if (strcmp(name,"Input")==0) return InputEnum; - else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; + if (strcmp(name,"Input")==0) return InputEnum; + else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; + else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum; else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; @@ -1118,12 +1120,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"SealevelchangeAzimuthIndexOcean")==0) return SealevelchangeAzimuthIndexOceanEnum; else if (strcmp(name,"SealevelchangeAzimuthIndexIce")==0) return SealevelchangeAzimuthIndexIceEnum; else if (strcmp(name,"SealevelchangeAzimuthIndexHydro")==0) return SealevelchangeAzimuthIndexHydroEnum; - else if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum; - else if (strcmp(name,"SealevelchangeViscousSG")==0) return SealevelchangeViscousSGEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum; + if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum; + else if (strcmp(name,"SealevelchangeViscousSG")==0) return SealevelchangeViscousSGEnum; + else if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum; else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum; else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum; else if (strcmp(name,"CouplingTransferCount")==0) return CouplingTransferCountEnum; @@ -1241,12 +1243,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"SmbPddfacIce")==0) return SmbPddfacIceEnum; else if (strcmp(name,"SmbPddfacSnow")==0) return SmbPddfacSnowEnum; else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum; - else if (strcmp(name,"SmbPrecipitationSubstep")==0) return SmbPrecipitationSubstepEnum; - else if (strcmp(name,"SmbPrecipitationsAnomaly")==0) return SmbPrecipitationsAnomalyEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"SmbDsradiationAnomaly")==0) return SmbDsradiationAnomalyEnum; + if (strcmp(name,"SmbPrecipitationSubstep")==0) return SmbPrecipitationSubstepEnum; + else if (strcmp(name,"SmbPrecipitationsAnomaly")==0) return SmbPrecipitationsAnomalyEnum; + else if (strcmp(name,"SmbDsradiationAnomaly")==0) return SmbDsradiationAnomalyEnum; else if (strcmp(name,"SmbDlradiationAnomaly")==0) return SmbDlradiationAnomalyEnum; else if (strcmp(name,"SmbWindspeedAnomaly")==0) return SmbWindspeedAnomalyEnum; else if (strcmp(name,"SmbAirhumidityAnomaly")==0) return SmbAirhumidityAnomalyEnum; @@ -1364,12 +1366,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; else if (strcmp(name,"VxBase")==0) return VxBaseEnum; else if (strcmp(name,"VxDebris")==0) return VxDebrisEnum; - else if (strcmp(name,"Vx")==0) return VxEnum; - else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"VxObs")==0) return VxObsEnum; + if (strcmp(name,"Vx")==0) return VxEnum; + else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; + else if (strcmp(name,"VxObs")==0) return VxObsEnum; else if (strcmp(name,"VxShear")==0) return VxShearEnum; else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum; else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; @@ -1487,12 +1489,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum; else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum; else if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum; - else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum; - else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum; else stage=13; } if(stage==13){ - if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; + if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum; + else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum; + else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum; else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum; else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum; @@ -1610,12 +1612,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition211")==0) return Outputdefinition211Enum; else if (strcmp(name,"Outputdefinition212")==0) return Outputdefinition212Enum; else if (strcmp(name,"Outputdefinition213")==0) return Outputdefinition213Enum; - else if (strcmp(name,"Outputdefinition214")==0) return Outputdefinition214Enum; - else if (strcmp(name,"Outputdefinition215")==0) return Outputdefinition215Enum; else stage=14; } if(stage==14){ - if (strcmp(name,"Outputdefinition216")==0) return Outputdefinition216Enum; + if (strcmp(name,"Outputdefinition214")==0) return Outputdefinition214Enum; + else if (strcmp(name,"Outputdefinition215")==0) return Outputdefinition215Enum; + else if (strcmp(name,"Outputdefinition216")==0) return Outputdefinition216Enum; else if (strcmp(name,"Outputdefinition217")==0) return Outputdefinition217Enum; else if (strcmp(name,"Outputdefinition218")==0) return Outputdefinition218Enum; else if (strcmp(name,"Outputdefinition219")==0) return Outputdefinition219Enum; @@ -1733,12 +1735,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition330")==0) return Outputdefinition330Enum; else if (strcmp(name,"Outputdefinition331")==0) return Outputdefinition331Enum; else if (strcmp(name,"Outputdefinition332")==0) return Outputdefinition332Enum; - else if (strcmp(name,"Outputdefinition333")==0) return Outputdefinition333Enum; - else if (strcmp(name,"Outputdefinition334")==0) return Outputdefinition334Enum; else stage=15; } if(stage==15){ - if (strcmp(name,"Outputdefinition335")==0) return Outputdefinition335Enum; + if (strcmp(name,"Outputdefinition333")==0) return Outputdefinition333Enum; + else if (strcmp(name,"Outputdefinition334")==0) return Outputdefinition334Enum; + else if (strcmp(name,"Outputdefinition335")==0) return Outputdefinition335Enum; else if (strcmp(name,"Outputdefinition336")==0) return Outputdefinition336Enum; else if (strcmp(name,"Outputdefinition337")==0) return Outputdefinition337Enum; else if (strcmp(name,"Outputdefinition338")==0) return Outputdefinition338Enum; @@ -1856,12 +1858,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition449")==0) return Outputdefinition449Enum; else if (strcmp(name,"Outputdefinition404")==0) return Outputdefinition404Enum; else if (strcmp(name,"Outputdefinition450")==0) return Outputdefinition450Enum; - else if (strcmp(name,"Outputdefinition451")==0) return Outputdefinition451Enum; - else if (strcmp(name,"Outputdefinition452")==0) return Outputdefinition452Enum; else stage=16; } if(stage==16){ - if (strcmp(name,"Outputdefinition453")==0) return Outputdefinition453Enum; + if (strcmp(name,"Outputdefinition451")==0) return Outputdefinition451Enum; + else if (strcmp(name,"Outputdefinition452")==0) return Outputdefinition452Enum; + else if (strcmp(name,"Outputdefinition453")==0) return Outputdefinition453Enum; else if (strcmp(name,"Outputdefinition454")==0) return Outputdefinition454Enum; else if (strcmp(name,"Outputdefinition455")==0) return Outputdefinition455Enum; else if (strcmp(name,"Outputdefinition456")==0) return Outputdefinition456Enum; @@ -1979,12 +1981,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition567")==0) return Outputdefinition567Enum; else if (strcmp(name,"Outputdefinition568")==0) return Outputdefinition568Enum; else if (strcmp(name,"Outputdefinition569")==0) return Outputdefinition569Enum; - else if (strcmp(name,"Outputdefinition506")==0) return Outputdefinition506Enum; - else if (strcmp(name,"Outputdefinition570")==0) return Outputdefinition570Enum; else stage=17; } if(stage==17){ - if (strcmp(name,"Outputdefinition571")==0) return Outputdefinition571Enum; + if (strcmp(name,"Outputdefinition506")==0) return Outputdefinition506Enum; + else if (strcmp(name,"Outputdefinition570")==0) return Outputdefinition570Enum; + else if (strcmp(name,"Outputdefinition571")==0) return Outputdefinition571Enum; else if (strcmp(name,"Outputdefinition572")==0) return Outputdefinition572Enum; else if (strcmp(name,"Outputdefinition573")==0) return Outputdefinition573Enum; else if (strcmp(name,"Outputdefinition574")==0) return Outputdefinition574Enum; @@ -2102,12 +2104,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition685")==0) return Outputdefinition685Enum; else if (strcmp(name,"Outputdefinition686")==0) return Outputdefinition686Enum; else if (strcmp(name,"Outputdefinition687")==0) return Outputdefinition687Enum; - else if (strcmp(name,"Outputdefinition688")==0) return Outputdefinition688Enum; - else if (strcmp(name,"Outputdefinition689")==0) return Outputdefinition689Enum; else stage=18; } if(stage==18){ - if (strcmp(name,"Outputdefinition608")==0) return Outputdefinition608Enum; + if (strcmp(name,"Outputdefinition688")==0) return Outputdefinition688Enum; + else if (strcmp(name,"Outputdefinition689")==0) return Outputdefinition689Enum; + else if (strcmp(name,"Outputdefinition608")==0) return Outputdefinition608Enum; else if (strcmp(name,"Outputdefinition690")==0) return Outputdefinition690Enum; else if (strcmp(name,"Outputdefinition691")==0) return Outputdefinition691Enum; else if (strcmp(name,"Outputdefinition692")==0) return Outputdefinition692Enum; @@ -2225,12 +2227,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition811")==0) return Outputdefinition811Enum; else if (strcmp(name,"Outputdefinition812")==0) return Outputdefinition812Enum; else if (strcmp(name,"Outputdefinition813")==0) return Outputdefinition813Enum; - else if (strcmp(name,"Outputdefinition814")==0) return Outputdefinition814Enum; - else if (strcmp(name,"Outputdefinition815")==0) return Outputdefinition815Enum; else stage=19; } if(stage==19){ - if (strcmp(name,"Outputdefinition816")==0) return Outputdefinition816Enum; + if (strcmp(name,"Outputdefinition814")==0) return Outputdefinition814Enum; + else if (strcmp(name,"Outputdefinition815")==0) return Outputdefinition815Enum; + else if (strcmp(name,"Outputdefinition816")==0) return Outputdefinition816Enum; else if (strcmp(name,"Outputdefinition817")==0) return Outputdefinition817Enum; else if (strcmp(name,"Outputdefinition818")==0) return Outputdefinition818Enum; else if (strcmp(name,"Outputdefinition819")==0) return Outputdefinition819Enum; @@ -2348,12 +2350,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition930")==0) return Outputdefinition930Enum; else if (strcmp(name,"Outputdefinition931")==0) return Outputdefinition931Enum; else if (strcmp(name,"Outputdefinition932")==0) return Outputdefinition932Enum; - else if (strcmp(name,"Outputdefinition933")==0) return Outputdefinition933Enum; - else if (strcmp(name,"Outputdefinition934")==0) return Outputdefinition934Enum; else stage=20; } if(stage==20){ - if (strcmp(name,"Outputdefinition935")==0) return Outputdefinition935Enum; + if (strcmp(name,"Outputdefinition933")==0) return Outputdefinition933Enum; + else if (strcmp(name,"Outputdefinition934")==0) return Outputdefinition934Enum; + else if (strcmp(name,"Outputdefinition935")==0) return Outputdefinition935Enum; else if (strcmp(name,"Outputdefinition936")==0) return Outputdefinition936Enum; else if (strcmp(name,"Outputdefinition937")==0) return Outputdefinition937Enum; else if (strcmp(name,"Outputdefinition938")==0) return Outputdefinition938Enum; @@ -2471,12 +2473,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1049")==0) return Outputdefinition1049Enum; else if (strcmp(name,"Outputdefinition1004")==0) return Outputdefinition1004Enum; else if (strcmp(name,"Outputdefinition1050")==0) return Outputdefinition1050Enum; - else if (strcmp(name,"Outputdefinition1051")==0) return Outputdefinition1051Enum; - else if (strcmp(name,"Outputdefinition1052")==0) return Outputdefinition1052Enum; else stage=21; } if(stage==21){ - if (strcmp(name,"Outputdefinition1053")==0) return Outputdefinition1053Enum; + if (strcmp(name,"Outputdefinition1051")==0) return Outputdefinition1051Enum; + else if (strcmp(name,"Outputdefinition1052")==0) return Outputdefinition1052Enum; + else if (strcmp(name,"Outputdefinition1053")==0) return Outputdefinition1053Enum; else if (strcmp(name,"Outputdefinition1054")==0) return Outputdefinition1054Enum; else if (strcmp(name,"Outputdefinition1055")==0) return Outputdefinition1055Enum; else if (strcmp(name,"Outputdefinition1056")==0) return Outputdefinition1056Enum; @@ -2594,12 +2596,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1167")==0) return Outputdefinition1167Enum; else if (strcmp(name,"Outputdefinition1168")==0) return Outputdefinition1168Enum; else if (strcmp(name,"Outputdefinition1169")==0) return Outputdefinition1169Enum; - else if (strcmp(name,"Outputdefinition1106")==0) return Outputdefinition1106Enum; - else if (strcmp(name,"Outputdefinition1170")==0) return Outputdefinition1170Enum; else stage=22; } if(stage==22){ - if (strcmp(name,"Outputdefinition1171")==0) return Outputdefinition1171Enum; + if (strcmp(name,"Outputdefinition1106")==0) return Outputdefinition1106Enum; + else if (strcmp(name,"Outputdefinition1170")==0) return Outputdefinition1170Enum; + else if (strcmp(name,"Outputdefinition1171")==0) return Outputdefinition1171Enum; else if (strcmp(name,"Outputdefinition1172")==0) return Outputdefinition1172Enum; else if (strcmp(name,"Outputdefinition1173")==0) return Outputdefinition1173Enum; else if (strcmp(name,"Outputdefinition1174")==0) return Outputdefinition1174Enum; @@ -2717,12 +2719,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1285")==0) return Outputdefinition1285Enum; else if (strcmp(name,"Outputdefinition1286")==0) return Outputdefinition1286Enum; else if (strcmp(name,"Outputdefinition1287")==0) return Outputdefinition1287Enum; - else if (strcmp(name,"Outputdefinition1288")==0) return Outputdefinition1288Enum; - else if (strcmp(name,"Outputdefinition1289")==0) return Outputdefinition1289Enum; else stage=23; } if(stage==23){ - if (strcmp(name,"Outputdefinition1208")==0) return Outputdefinition1208Enum; + if (strcmp(name,"Outputdefinition1288")==0) return Outputdefinition1288Enum; + else if (strcmp(name,"Outputdefinition1289")==0) return Outputdefinition1289Enum; + else if (strcmp(name,"Outputdefinition1208")==0) return Outputdefinition1208Enum; else if (strcmp(name,"Outputdefinition1290")==0) return Outputdefinition1290Enum; else if (strcmp(name,"Outputdefinition1291")==0) return Outputdefinition1291Enum; else if (strcmp(name,"Outputdefinition1292")==0) return Outputdefinition1292Enum; @@ -2840,12 +2842,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1411")==0) return Outputdefinition1411Enum; else if (strcmp(name,"Outputdefinition1412")==0) return Outputdefinition1412Enum; else if (strcmp(name,"Outputdefinition1413")==0) return Outputdefinition1413Enum; - else if (strcmp(name,"Outputdefinition1414")==0) return Outputdefinition1414Enum; - else if (strcmp(name,"Outputdefinition1415")==0) return Outputdefinition1415Enum; else stage=24; } if(stage==24){ - if (strcmp(name,"Outputdefinition1416")==0) return Outputdefinition1416Enum; + if (strcmp(name,"Outputdefinition1414")==0) return Outputdefinition1414Enum; + else if (strcmp(name,"Outputdefinition1415")==0) return Outputdefinition1415Enum; + else if (strcmp(name,"Outputdefinition1416")==0) return Outputdefinition1416Enum; else if (strcmp(name,"Outputdefinition1417")==0) return Outputdefinition1417Enum; else if (strcmp(name,"Outputdefinition1418")==0) return Outputdefinition1418Enum; else if (strcmp(name,"Outputdefinition1419")==0) return Outputdefinition1419Enum; @@ -2963,12 +2965,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1530")==0) return Outputdefinition1530Enum; else if (strcmp(name,"Outputdefinition1531")==0) return Outputdefinition1531Enum; else if (strcmp(name,"Outputdefinition1532")==0) return Outputdefinition1532Enum; - else if (strcmp(name,"Outputdefinition1533")==0) return Outputdefinition1533Enum; - else if (strcmp(name,"Outputdefinition1534")==0) return Outputdefinition1534Enum; else stage=25; } if(stage==25){ - if (strcmp(name,"Outputdefinition1535")==0) return Outputdefinition1535Enum; + if (strcmp(name,"Outputdefinition1533")==0) return Outputdefinition1533Enum; + else if (strcmp(name,"Outputdefinition1534")==0) return Outputdefinition1534Enum; + else if (strcmp(name,"Outputdefinition1535")==0) return Outputdefinition1535Enum; else if (strcmp(name,"Outputdefinition1536")==0) return Outputdefinition1536Enum; else if (strcmp(name,"Outputdefinition1537")==0) return Outputdefinition1537Enum; else if (strcmp(name,"Outputdefinition1538")==0) return Outputdefinition1538Enum; @@ -3086,12 +3088,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1649")==0) return Outputdefinition1649Enum; else if (strcmp(name,"Outputdefinition1604")==0) return Outputdefinition1604Enum; else if (strcmp(name,"Outputdefinition1650")==0) return Outputdefinition1650Enum; - else if (strcmp(name,"Outputdefinition1651")==0) return Outputdefinition1651Enum; - else if (strcmp(name,"Outputdefinition1652")==0) return Outputdefinition1652Enum; else stage=26; } if(stage==26){ - if (strcmp(name,"Outputdefinition1653")==0) return Outputdefinition1653Enum; + if (strcmp(name,"Outputdefinition1651")==0) return Outputdefinition1651Enum; + else if (strcmp(name,"Outputdefinition1652")==0) return Outputdefinition1652Enum; + else if (strcmp(name,"Outputdefinition1653")==0) return Outputdefinition1653Enum; else if (strcmp(name,"Outputdefinition1654")==0) return Outputdefinition1654Enum; else if (strcmp(name,"Outputdefinition1655")==0) return Outputdefinition1655Enum; else if (strcmp(name,"Outputdefinition1656")==0) return Outputdefinition1656Enum; @@ -3209,12 +3211,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1767")==0) return Outputdefinition1767Enum; else if (strcmp(name,"Outputdefinition1768")==0) return Outputdefinition1768Enum; else if (strcmp(name,"Outputdefinition1769")==0) return Outputdefinition1769Enum; - else if (strcmp(name,"Outputdefinition1706")==0) return Outputdefinition1706Enum; - else if (strcmp(name,"Outputdefinition1770")==0) return Outputdefinition1770Enum; else stage=27; } if(stage==27){ - if (strcmp(name,"Outputdefinition1771")==0) return Outputdefinition1771Enum; + if (strcmp(name,"Outputdefinition1706")==0) return Outputdefinition1706Enum; + else if (strcmp(name,"Outputdefinition1770")==0) return Outputdefinition1770Enum; + else if (strcmp(name,"Outputdefinition1771")==0) return Outputdefinition1771Enum; else if (strcmp(name,"Outputdefinition1772")==0) return Outputdefinition1772Enum; else if (strcmp(name,"Outputdefinition1773")==0) return Outputdefinition1773Enum; else if (strcmp(name,"Outputdefinition1774")==0) return Outputdefinition1774Enum; @@ -3332,12 +3334,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"Outputdefinition1885")==0) return Outputdefinition1885Enum; else if (strcmp(name,"Outputdefinition1886")==0) return Outputdefinition1886Enum; else if (strcmp(name,"Outputdefinition1887")==0) return Outputdefinition1887Enum; - else if (strcmp(name,"Outputdefinition1888")==0) return Outputdefinition1888Enum; - else if (strcmp(name,"Outputdefinition1889")==0) return Outputdefinition1889Enum; else stage=28; } if(stage==28){ - if (strcmp(name,"Outputdefinition1808")==0) return Outputdefinition1808Enum; + if (strcmp(name,"Outputdefinition1888")==0) return Outputdefinition1888Enum; + else if (strcmp(name,"Outputdefinition1889")==0) return Outputdefinition1889Enum; + else if (strcmp(name,"Outputdefinition1808")==0) return Outputdefinition1808Enum; else if (strcmp(name,"Outputdefinition1890")==0) return Outputdefinition1890Enum; else if (strcmp(name,"Outputdefinition1891")==0) return Outputdefinition1891Enum; else if (strcmp(name,"Outputdefinition1892")==0) return Outputdefinition1892Enum; @@ -3455,12 +3457,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"AdaptiveTimestepping")==0) return AdaptiveTimesteppingEnum; else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum; else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; - else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; - else if (strcmp(name,"AgeAnalysis")==0) return AgeAnalysisEnum; else stage=29; } if(stage==29){ - if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum; + if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; + else if (strcmp(name,"AgeAnalysis")==0) return AgeAnalysisEnum; + else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum; else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum; else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum; else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; @@ -3578,12 +3580,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"GLheightadvectionAnalysis")==0) return GLheightadvectionAnalysisEnum; else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum; else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; - else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; - else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; else stage=30; } if(stage==30){ - if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; + if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; + else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; + else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; else if (strcmp(name,"GenericParam")==0) return GenericParamEnum; else if (strcmp(name,"GenericExternalResult")==0) return GenericExternalResultEnum; else if (strcmp(name,"Gradient1")==0) return Gradient1Enum; @@ -3701,12 +3703,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum; else if (strcmp(name,"MeshX")==0) return MeshXEnum; else if (strcmp(name,"MeshY")==0) return MeshYEnum; - else if (strcmp(name,"MinVel")==0) return MinVelEnum; - else if (strcmp(name,"MinVx")==0) return MinVxEnum; else stage=31; } if(stage==31){ - if (strcmp(name,"MinVy")==0) return MinVyEnum; + if (strcmp(name,"MinVel")==0) return MinVelEnum; + else if (strcmp(name,"MinVx")==0) return MinVxEnum; + else if (strcmp(name,"MinVy")==0) return MinVyEnum; else if (strcmp(name,"MinVz")==0) return MinVzEnum; else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum; else if (strcmp(name,"Moulin")==0) return MoulinEnum; @@ -3824,12 +3826,12 @@ int StringToEnumx(const char* name,bool notfounderror){ else if (strcmp(name,"IntrusionMelt")==0) return IntrusionMeltEnum; else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum; else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum; - else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; - else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; else stage=32; } if(stage==32){ - if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; + if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; + else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; + else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; else if (strcmp(name,"Tetra")==0) return TetraEnum; else if (strcmp(name,"TetraInput")==0) return TetraInputEnum; else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum; diff --git a/src/c/shared/Enum/issmenums.jl b/src/c/shared/Enum/issmenums.jl index 3088a80c9..98ba8ce2d 100644 --- a/src/c/shared/Enum/issmenums.jl +++ b/src/c/shared/Enum/issmenums.jl @@ -423,6 +423,8 @@ OutputFileNameEnum OutputFilePointerEnum OutputdefinitionEnum + PddfacIceEnum + PddfacSnowEnum QmuErrNameEnum QmuInNameEnum QmuIsdakotaEnum @@ -4200,6 +4202,8 @@ function EnumToString(enum::IssmEnum) if(enum==OutputFileNameEnum) return "OutputFileName" end if(enum==OutputFilePointerEnum) return "OutputFilePointer" end if(enum==OutputdefinitionEnum) return "Outputdefinition" end + if(enum==PddfacIceEnum) return "PddfacIce" end + if(enum==PddfacSnowEnum) return "PddfacSnow" end if(enum==QmuErrNameEnum) return "QmuErrName" end if(enum==QmuInNameEnum) return "QmuInName" end if(enum==QmuIsdakotaEnum) return "QmuIsdakota" end @@ -7977,6 +7981,8 @@ function StringToEnum(name::String) if(name=="OutputFileName") return OutputFileNameEnum end if(name=="OutputFilePointer") return OutputFilePointerEnum end if(name=="Outputdefinition") return OutputdefinitionEnum end + if(name=="PddfacIce") return PddfacIceEnum end + if(name=="PddfacSnow") return PddfacSnowEnum end if(name=="QmuErrName") return QmuErrNameEnum end if(name=="QmuInName") return QmuInNameEnum end if(name=="QmuIsdakota") return QmuIsdakotaEnum end diff --git a/src/m/classes/SMBpddSicopolis.m b/src/m/classes/SMBpddSicopolis.m index 153264197..a94cee8a4 100644 --- a/src/m/classes/SMBpddSicopolis.m +++ b/src/m/classes/SMBpddSicopolis.m @@ -16,6 +16,8 @@ s0t = NaN; rlaps = 0; isfirnwarming = 0; + pdd_fac_ice = 0; + pdd_fac_snow = 0; steps_per_step = 1 averaging = 0 requested_outputs = {}; @@ -71,6 +73,8 @@ self.isfirnwarming = 1; self.desfac = -log(2.0)/1000; self.rlaps = 7.4; + self.pdd_fac_ice = 7.28; + self.pdd_fac_snow = 2.73; self.requested_outputs={'default'}; end % }}} @@ -85,6 +89,8 @@ md = checkfield(md,'fieldname','smb.rlaps','>=',0,'numel',1); md = checkfield(md,'fieldname','smb.monthlytemperatures','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 12]); md = checkfield(md,'fieldname','smb.precipitation','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 12]); + md = checkfield(md,'fieldname','smb.pdd_fac_ice','>',0,'numel',1); + md = checkfield(md,'fieldname','smb.pdd_fac_snow','>',0,'numel',1); end md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]); @@ -108,6 +114,8 @@ function disp(self) % {{{ fielddisplay(self,'isfirnwarming','is firnwarming (Reeh 1991) activated (0 or 1, default is 1)'); fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'); fielddisplay(self,'averaging','averaging methods from short to long steps'); + fielddisplay(self,'pdd_fac_ice','Pdd factor for ice for all the domain [mm ice equiv/day/degree C]'); + fielddisplay(self,'pdd_fac_snow','Pdd factor for snow for all the domain [mm ice equiv/day/degree C]'); disp(sprintf('%51s 0: Arithmetic (default)',' ')); disp(sprintf('%51s 1: Geometric',' ')); disp(sprintf('%51s 2: Harmonic',' ')); @@ -124,6 +132,8 @@ function marshall(self,prefix,md,fid) % {{{ WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0p','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0t','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'class','smb','fieldname','rlaps','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','pdd_fac_ice','format','Double'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','pdd_fac_snow','format','Double'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); diff --git a/src/m/classes/SMBpddSicopolis.py b/src/m/classes/SMBpddSicopolis.py index 9d627a5e6..502f1a785 100644 --- a/src/m/classes/SMBpddSicopolis.py +++ b/src/m/classes/SMBpddSicopolis.py @@ -26,6 +26,8 @@ def __init__(self, *args): # {{{ self.s0t = np.nan self.rlaps = 0 self.isfirnwarming = 0 + self.pdd_fac_ice = 0 + self.pdd_fac_snow = 0 self.steps_per_step = 1 self.averaging = 0 self.requested_outputs = [] @@ -50,7 +52,9 @@ def __repr__(self): # {{{ s += '{}\n'.format(fielddisplay(self, 'desfac', 'desertification elevation factor (default is -log(2.0)/1000)')) s += '{}\n'.format(fielddisplay(self, 'isfirnwarming', 'is firnwarming (Reeh 1991) activated (0 or 1, default is 1)')) s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step')) - s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) + s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) + s += '{}\n'.format(fielddisplay(self, 'pdd_fac_ice', 'Pdd factor for ice for all the domain [mm ice equiv/day/degree C]')) + s += '{}\n'.format(fielddisplay(self, 'pdd_fac_snow', 'Pdd factor for snow for all the domain [mm ice equiv/day/degree C]')) s += '\t\t{}\n'.format('0: Arithmetic (default)') s += '\t\t{}\n'.format('1: Geometric') s += '\t\t{}\n'.format('2: Harmonic') @@ -100,6 +104,8 @@ def setdefaultparameters(self): # {{{ self.isfirnwarming = 1 self.desfac = -np.log(2.0) / 1000 self.rlaps = 7.4 + self.pdd_fac_ice = 7.28 + self.pdd_fac_snow = 2.73 self.requested_outputs = ['default'] return self # }}} @@ -114,6 +120,8 @@ def checkconsistency(self, md, solution, analyses): # {{{ md = checkfield(md, 'fieldname', 'smb.rlaps', '>=', 0, 'numel', 1) md = checkfield(md, 'fieldname', 'smb.monthlytemperatures', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 12]) md = checkfield(md, 'fieldname', 'smb.precipitation', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 12]) + md = checkfield(md, 'fieldname', 'smb.pdd_fac_ice', '>', 0, 'numel', 1) + md = checkfield(md, 'fieldname', 'smb.pdd_fac_snow', '>', 0, 'numel', 1) md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1]) md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2]) md = checkfield(md, 'fieldname', 'smb.requested_outputs', 'stringrow', 1) @@ -128,6 +136,8 @@ def marshall(self, prefix, md, fid): # {{{ WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 's0p', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 's0t', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'rlaps', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'pdd_fac_ice', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'pdd_fac_snow', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'monthlytemperatures', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'precipitation', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'temperature_anomaly', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) diff --git a/src/m/classes/model.m b/src/m/classes/model.m index ee3953e0c..8623f0eaf 100644 --- a/src/m/classes/model.m +++ b/src/m/classes/model.m @@ -201,6 +201,9 @@ if ~isa(md.mmemasstransport,'mmemasstransport'); md.mmemasstransport=mmemasstransport(); end % 2026 February 18 if isa(md.friction, 'frictionjosh') && md.friction.coefficient_max==0; md.friction.coefficient_max=300; end + % 2026 February 20 + if isa(md.smb, 'SMBpddSicopolis') && md.smb.pdd_fac_ice==0; md.smb.pdd_fac_ice=7.28; end + if isa(md.smb, 'SMBpddSicopolis') && md.smb.pdd_fac_snow==0; md.smb.pdd_fac_snow=2.73; end end% }}} end methods From 372c47c52437db48ef4fcdde0a7ac1066fb094c0 Mon Sep 17 00:00:00 2001 From: NJSchlegel <95945328+NJSchlegel@users.noreply.github.com> Date: Fri, 20 Feb 2026 11:28:51 -0500 Subject: [PATCH 07/17] CHG: GEMB precipitation mapping unit valideation --- src/c/classes/Elements/Element.cpp | 11 +++++++---- test/Archives/Archive258.arch | Bin 24551 -> 24551 bytes 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/c/classes/Elements/Element.cpp b/src/c/classes/Elements/Element.cpp index db7206881..b4fab4806 100644 --- a/src/c/classes/Elements/Element.cpp +++ b/src/c/classes/Elements/Element.cpp @@ -5665,15 +5665,18 @@ void Element::SmbGemb(IssmDouble timeinputs, int count, int steps){/*{{{*/ //es over ice calculation //Ding et al., 2019 after Bolton, 1980 //ea37 = rh37*100*6.112.*exp((17.67*(t237-273.15))./(t237-29.65)); - IssmDouble esparam, es; + //ea37s (hPa) = 6.1121*exp(22.587*(t-273.15)/(t-273.15+273.86)); (with respect to ice) + IssmDouble esparam, es, qsparam, qs; esparam=6.112*exp((17.67*(taparam-273.15))/(taparam-29.65)); es=6.112*exp((17.67*(Ta-273.15))/(Ta-29.65)); rhparam=eaparam/esparam; eAir=fmax(rhparam*es,0.0); + qs=fmax(0.622*es/(pAir/100 - 0.378*es),0); + qsparam=fmax(0.622*esparam/(pparam/100 - 0.378*esparam),0); - if ((isprecipmap) && (eaparam>0) && (pAir>0)){ - P=fmax(prparam*es/esparam*pparam/pAir,0.0); - C=fmax(C*es/esparam*pparam/pAir,0.0); + if ((isprecipmap) && (qsparam>0)){ + P=fmax(prparam*qs/qsparam,0.0); + C=fmax(C*qs/qsparam,0.0); } else P=prparam; diff --git a/test/Archives/Archive258.arch b/test/Archives/Archive258.arch index d71423628a4ae10730136ec5b99a709a227215a1..92461413d602a8c684828f66fee6a2c792e6c9b1 100644 GIT binary patch delta 6350 zcmZvgc{r3^8^94+vS%;ZLLyQ~_39jjii#{Hdz7~lovCCu%~%UX_W3+;gLqtjm6`EDUJ#T({^qG$d1H z$pMMqg|xLS7r;$d-@KzrL@sc%yi;T9P9w-{O*4RuCm`#5`|1QKNsu%xxk@VsXiXH@T(?(e^C3y|FX(eRS zFxv?71of1*mOzliqEEby2t*~%tW_l+KvY$I zq@VN{qUVOH*&MtOwQ=k`MxKMHOMQnby%(ZSB87G5PeJtM^+fddK8U8aY6y1-Lo~B| z{&mhih~@;+Rrzm0^vkR1@6m6Nz;-}jy-f@xu(yG>$}mXa7*-nSFM)(*MaNH#6+i+f z(N~Rkq7@Rjjv()CN|3O8d7sjOJD8_TJl|#ssRamAg&kgyAUOSF zMUg8ch`KX>Ac1m3a1j3r69uT@{YCcp0f7KspD{d80nyI|PJa$Lydy(vAM8_xXgb=qF_{md zA5`(8>SBn#=hspLe&U_0@)6MphiJ0N`$MKRB%+CF%}+0QAesQKJCyJZ<3~3Yt+IgV zn@Qu>s?!jSS^i+RZ-wYMQ?>2HnGy;o_XQ>b|RF;RC20%2- zAXhcWFQx;sSg%T&HxPk`uM z@Kn>MLWnwzs`NZ@&AjE!Xso6{^ycN$x!YY3y^e#zwLuWITkA{6;T>v|Xri(EK-BV1 zW+9!3BYKr=^th!5qGrbvVpK(lUdmGT%dIF-x%#Q~giXdt@AyzlN6&I=c zN3f6sQQesN;dp6?YD+F}@c=@P9-xPo;VvqAR=2x}3!=iTAMd9mAS$>N>2|;vqI+i| zZRh`n=&mhy#4a3o{w{I;$72xP*;r-&DGj3A3Sz|V&+H+(B_{v8wHl&4-ftJgu;+4I z6LiB9X8#f>dA9|ktT`{nIs8VM^HrDIu3~cVG^;Tn%HXxK*y93G`pO)^=SLvAsj%X- z-6%viv?Dfat3i}zET0hU08y%a*EUf$h*G?3)Ok1|qF!qsiFSf0xtJ*7U7!I`Qb&xg z`5uT8=eLV9WfAHFWVWm7gzs3U)LNIe1@r}fA+vn zh~nL2L#}EKK6N+VJ;hog6QWC_?@qVl1uwZ|dkh6Y6xrcPys{S4d>t&nIge;;5tG0>5wX`vWBe>c z;f+mK&l4abhMSA5+&u?TScj@R556(Xgf+DLfUxO$l z{rP1_+$|y76ox2th{&(S+1}$6kx?QsH3?HFdqKniB2u9)rvy$B$!0sd7T)TMKMQOR z2w-MLJ7gPSn%7yMYR3IfoVW2ODdR};wb@iZU^-g}#@b-6%zN&s08wzGuGu-<1;ON@ z!WMiK2OCtbypB^8G@BPGfe)ad@(fkIZcP7lz6yOzgB0^%+)Y6o(Gt11{{tJ2hc`2M zP5z}=2Bt&Fn!j-M2dYVz=0*}B3ix46JJg9;zj5d#J{SW2ZQD5Og?U;*jyn+3P<;Jb zrb}jZ>^TUL|JTNnK3t^!4FRW);y3iqyqCE>4wK}pxH}!wnN;v*91{qspH0A&TZp`Z z7xw2coq1IUk>9|usyQxNzs3%qFybzJQOKuLf^*^b??rQ~HJGH$=~r>YejY{&`eOgX zLk>}w@ImCq;?9JGdR%mVENU?}?B~GF;+ky??r`iZ3V1#Q`Pi4}J0E+HX#Nf&=EWiz zE7n)xF7RcsgC3nZjGe`rtV8+Rv9oxu=DQ#Bc@}Sg+CW@R1X z__9c}AU-)$f(M8e|E2#KMZyer7JUbUhH&5evglu~6Xe1Ce2K(Q%@1ps&o41D;F;)N z>@3EI`nUeV&SECNZ)zWQ7E6h%ZeKZq2TGUtvba3}}gdTt|l5ac+N8 zW$RS82zC}tjgJu8u(L=odDVw&$dARO`NB_l$NX3%*2`I6@xcQuPT9UZfsatXCHgB{ zjGE#U{N|xsyPx8A^JnoAmogD|jX#Tz*i9v^u(S9T1TWy?^6*l-aAuvq$}C$j-Ni=S2I-lkw*VxSo9iQZA%PAtNEX{IZ&vv|tx3I{Ik zKo;ZfT^`0=68K*dQQ(76-A6dW-|T!4TIYO&`yQNPLSYT(d=sIOi2|JTS=0!<*r!CI$ZXFoMSJmEpP0^~ zTR|lLTM$&`whUiT?n!Uqr7_jJ_&$yDRa$=Wzx2&A;>4kpsUGM(HRE|iR1scS)=e?H z9|$jqon0~bq7z<_=MCc4^@A5wl}id9{)87a_MkzBIC#N3U*vXZ0=zKT7#(Tw3tkvH zW|^z!0)3ApghmJnLthtI7)rN?e&?Q;1XmECKV{&`#=8@cxM%u+GZKYFgMq=x4SevT z7|~GIkM}XWm%uVeI@Mj2a2S&1l2zN?2FZgb7o7X`A(@^hnF{P6`Ht1!Neb1F?31o8BP#{T+6j~6 z_^&Z}z0R#TyB0`Z&OhQdGyzG((yGEP_0N!W)hy-yKno=Kn8_rB>qFA9WQpOZFi1MO z=VT1K4J0|r8*)KCNOFQc#{#M$$;w%jOI-|-PA8dV1-*x)6FaD{^wRKi+wIdft3i^- zm$SnsEv$IsA*b1Cq)Kj+Q2>kW{P{_L*ZGk{Vnu8ChO|q*sHD3Ek#RlbrtUby+kG()NulKFWL+j&!o zkgOr*uc7-Ek~gfN?zcS$$@)Q}NMjU|)kTB6)hZxa%WK8<+)a>dk$!q%557MhSdn`& zAlWl~YPnb*3|V`xW#i8x7@BwTq3iWQ7_sH+XLGhu7)dVg)pGm-FYj3GF*PLyqm@;% zz=hK=%2WD#zTI9J^>A0xtc=1B7?rT>9;dqwyzGCHTj%X3c)8i3OgA_WM(giAaZY55Q(yoHSQa1C8!G@rupFMn}>mqop^PPInzz00kTsII;H3DVkHM_6N zQb1{{daI;gHYoj}*PinH5O~mcxKDRoET}h=){m@D0kt13u1!jsfrjRTElr+9prJ0) ztH-ebG+dWF8RuaP8WLswd{3_^H0h6a5$u?n5=c?I@f!-3Up{DMSpqHBYQTV7P=*?0euOog3y=Nw# zSW_lJuWR*J)!cGG=hH}D7&Hg;?Bz52J)Q%4A@!4}b1tBdn1`Iz&jEDeNA;E{_YHvl zj=yAO;x<4p%j!8$AO+|X^gCXATL691N6y(o0MLC6UT6XmpjWo@X(Vh1^v!XOQkJ2B zZgjLWB{2igEq8DlS>exx6H6D$v_NmtU&beWt}T)}W(|7Fw@q3UxdOV|LA`ZS`16z= zg*~m;0g-tjQfUon(rXJ3+th(xE|2Ov@k)TU?XzWuMhBqn%X+D&k`8(XyXVgDQUu+H zn*}GN(m}UMyJFEX4e;UY1)Y7p4uB@_wBGnr6zJI>o=K-|1~kDvdnP>31DZ{hx2xk< zKzndfeXo!-pxurnp4nFbXa|QK4_k2&0d2#Z>#?DKgPx}}wuyiQ&|}aP;(I(5^oUzT zt{rg)-67ZPbh(W{mwfwH|0W60{q%Wtreiqhs^ax4NSFa1;=KLs2y&qFR7gpcYBhK} z!X8mP7XUh_+V009xB#ed zp^8uO7Ar+4w+;WLR2a$^+wY88LHV*@lX!AqE*>zrdyE|4n*auH0=V-c2rziooKO3v z2pGa{WA!`sVpJxqFl5;Oia1*W z7~(35$mk7VNUBz%M9DL&+yo z4`qyTVA0R7tj6DfGWc&*-aG0Fesg8$>Pc!WV2CB20_k|GSzP|jOg`5Z2XfPM>ZUwk z2zEvrY&1y;(%hvHx^$14<9Bho6pCeUqJKXLO5`* zW8Iv`iGU#ui7grcM|{LSA+rfE6i-NKKfya7ENbFr#`T9H$HK_I&u|L#8cGn}5k3#< z)Ibd0iTizC6BBr$siUQyB*2h*3#-=_|DkBq8P`}9z>u$e>hkkAV8|WRbk+C6>wPQS zQ?0NI4A0$-(!AmFo1Z80$W#GSe9Vm;dSFviTwsBFkon!?oE<+`Ut#X4aEgM1EM>ht zKUZ~KN~-&!+1j7rB5s!EI>9s%9o%L7_?H6;nZzo~DkUx6BkaJ>HFz^f)EnD+PQ7wd z)y3h{c=*y*@3hiIOL7q4EOB`%yC~1GFu(1j#5u}}Z2Udx(idxsquT#>bkqMErCfC6 z=b9_X1i)b)fxLPumnZ^;O(ti zCV}o!Dz{hN{DdMCjcuGXY*3;>or&f@DkxF3$}@C)97Xx1`Sd9Epp^0A^z_A3=+3u` zBLAd5KzF{j^5?&EM3J?3+mu%Lpvc$T(^5iR(VYQzxrL4gQ8Y1iVdHv36s0edbm-qo z6tzoU$R_F*ilZL?`(c(>1WHkUy0}D*f)duTLL9i)qlCW$quzNMpwtOo%gTsx6uHJr z#rx28Fhs)g6%|D?Wi3|h+KCcYhh@cY>Ozs@BbV%c^rE{wKWwj^m_aesD-?p`+)-3} zq2fx-Je0~UIjkLW8Aa6#2k(Bl0>#<98W#vNM$|aDBP*(Tub`AW!|4VbZ7A-;5ysFu zZIpgcEo$3sX%sgQb)WT23?-(Bi(Q=fic;68H#xZSp@bi*7h@DEP(tW2ZgV+Hl&Bi_ zQ>jG&CALPG$fWK_4(lOlpPliDhmm$$Q8^OlUFOkclswcMgi#Fmn7Y#}z7y z*;wCPXoOHwoTRjuh&8&KT$|-LeGWw#4{>@F7ooW588`Nq87LvpW8%bze3aOrdE8<= z9wq*KyPCb~C`!4$YHw1i7)r<&-MaI65KKXQ(zr(wigdKu7_qkxrTHIpV()o@?u?4N zDArD*_@iDKryUod_`KG8M(Wf6l=LQPSD>~&N_KgtR_(DK#SW;Q)ccB1>`Ybi+Nunc za%L%~<{4WQ|4-qGdg>s$yJTE)-J~!|l`ihLQRP61!6|FC#&)CVb6>Ovx9g$kbZtYO z`nUvmG`Whvt{j%54ItXhcrCu)Y42m~XqH0Q8^6^;h1>K&Y_{Fjfsz{pBs9@F!5OgXKa6giEj-@*D!K1F{0z*)})Dv2`-g{ zhEhyS@!#iebHv0i$IP$7iI}orcg=;gdzivj|9<`B986*Fy>sqIKc+0qEZ`Ndz!VOu zt99VWtY0!dtlTje9kw@Fi4jme_(R{I{4p|NVg}0&itj%Zy-(;0Plq zEE-cpPt5*c3f)q0ka~BJrc1zQ`Rmpx=0P8;^=FInEX|%R96kf=;O9OT7-%5@ZIbTDC$_0bbe_+ zCPp7N+lVzn(TXn#1o30yNA*RvM;9?c9Uj+L+t0$pFgpK9{u3sKj%<;n7-M44s6`_+ z0uuxK7o98uOSZ{xq56&jx}QZc@zvs#j=m=*zTn5z>lm0|DYGyBsE>($Zs(3( zTTJx!$~#L#d(TTH`&uoSz1Y|@Cw^h#lUs(2DHkTXjh$clp2t+8YZ>ph)J#k;hHi8f z!-YFbwbWL@9d%sgW9%-+M7yhlysMkejJDFwt0(`bO0q6Ackc915^J-W`zDTmJx-l|-h|DOF6=x6EX< zRZuZecXwr_w;LvEZH#?W6);gFfAK<78YZe4lcPK2Fj19M`ssN)CMr)(mg+&tD%8rN zRFyIDYCv1qW(X56vz-k>DVQkty4u_ni-|IWhl33+m?-7#UKz>6L`j2j*tQv1-?0lL zUcrtxafp2k$I6aL4mc9agNg9 zumTg04i0>lyMu{`2EUA>dob}}ZJJLDRP?^W^kyoCQr#2hxu3NV6X`r3HblT3r?Gox zaXMfkWqNd@YbPd>rw*E}fH_Z^PHp)n2(k^b>%agLdGyZCz#J!tT@mMkIgeKuwvmGh z#c7-qxjK!BSVQjP>FY5ObEy2M(NauAUr6t`5{HSXYx%eG;JSC8{Gv*|H^f9_lS@N2 ztbseziceHv&hM!3N^M?-iQ5NO4Jr0w;#OdIT1hr0Za!08GYE5i<2x5eTrId#VQs%P zCa#}wH!-||33{r)<@+}=ajm~SvxEcG=$Du{kBJCx_59Z-FmbgqRX1A#l-~dLq7^2t zI8oQ0I`R|~|GeCA>GdU0{;W&!b4*R4;SgK3up%PGJ zjnhK{6PJ8OMUBibaq(;EP{a|??&^}Flb8r?I(WcW3N#Vj{3IF^Ayv;eio*p%c5J6s zLxqE@as+<1g48{eMOz9m5!58QSqn}KvKnEry(7w_M6aEcs%1`P*BB)t45fq|_K6`-)+T&#p$qWPMD6>rtuGTeMT40@Q~~2OVJ%Qf)VDt~d^g-B$GA#)NOx@(+^G?(5%j zK@}?DyXhg{^DIpG{JfpL7`CBL-nEOeuqb^zBjPVF2RB`o*#wpGVfU8Y0|WK0G%Gl; zAG{`_?Ep9CZJk_r0Q&J>$s2YW!T$Fe)9Ouz$F|pNMdeGk!P{jC?B78L(Q-`~h?oAN zg}dNRyd-aR>FQ#_^LxRlQxjO}*uD>*6P`&4yUgE#w0%>i=RjMR!iabfv7d=f0Of}! zx+6jM&Cxut);vCaoLaUDtg7;f8HN4tk@@Rv@m=t?OU_s5*n{?$uBkkDR;M@B^FNyJ z^WPc;nbaV>@FcvKJV@7c(uiY_Ntdp5^cR9mDsb$qygKB0-KT@sx>NfxF?TR4i*3sd zSlz#EYUyi*lDd3Q;{c`{^DcZPeeKV(wRC&8_HvwxqL+Ub5Yc=r5xLobxScmkR9y6;=< z50FWd7Afyo44E`ng8HZxR-=2a#Q&_xSgUauGHFlckK8wq=k?$;>v?zsN=T}#wUwI> zGU-nJkwRGg9;B2lxvysDnl$n!(`$B$NncZCtKwmQdc2-@#5~B?6!N^DEB16{X~5zk zHBWPosenv+abb9^9c0q{nXYPh4tSD&-@I*6D&%>+6pz$(!n@|TR5Rhh(@D@k8a%pP z2ukWjn)_k{DQ~=C!(fQ-w^L zZ|$v^{;m2#Z-v=Rs$V*Ic{Y>o{p*(GY$iRlIQ`*lCOwuJ$~&8nWwM9z z(sQM>dD>^=g{-IT^f4&~cg$I3hTidW}Z7*xPuP>*UDq|%imBVjtPjwpF?Cb-=N>;7>?cme z_l&mjV80|(nB#N>hor1uqBYEqLlWMZY*W3DL*E(ePwN=spif*pecqxtsC!B%*>67% zGATnS54t(BUtk)A?UVWu2~4wE{=IYQDoi`h zneOM*hG~cI4wMKbW7_^x{g<?~x2X)gN<3No0OW-X%qng0N$d07S#`;KBN97IXgzy)5HMBflSj%g8Z zdqQ=;U|ODGzDS!$!2V$2>*rMmF6~9WK4oA;S^84Mci=*5a zr62ufyeM}}ZIxw_Hp=NTPMRvHK(BX6tt~WcKrh=&?xlo?p~_m*#Q|AaSDR8fi>xp-M)JXL#9nu=!nBN*H?C)ans; z&>u0>Qa?U%TRv<1#K?JYT=wjSV#LU{r99Hxh!{~xX7$H95F=Ub{FT-AGXbMSABSz$34$6xT&3UZ+FT}{ZptyjgjF`vvYs*x&B1ZHZ#;c@Tb6)X? zMbyqEKlRLZ1|bG#@^;%=dDN+A5SzJ61aJ(g^9? z^S;OUF2g|7xo?cqfpH0Ss*kOnQ9njS9c6ACY(_t!j!k*UQppx|EZ(*1=>aX&9$-E4 zadl?{wh@WF3@~vJOfo#h&VPOBdCbJJ8Pa%MUfFe-7FkSBYwun=y}YZAAGx>J;VbZuIWk^yrS7FX-cO z&wbTyQK&N`lwa$VHtNvu-~Ie(5Mt&ryY{|NM@@}m`lh_9h+(55y7cMa^axD_Zq%1w zQ&%#;LBFqmo4(Ls37w_CB+JX-JRdzylbim{cs13rq0tDj7Ikd%eea1_Jh!$?x~U-+ z_fj|4H+qOAd~D$LW@C`xVLUyLSiH;aEycqTi%W{(WjBmiLLV$bOU-`UU$0rEDXV}; zSxXemJN;K7me^I(H?Klae_upbvWD9q>!~&eYE3K(5o)h{b=eSRuI$QON_ASfK1Ldl|$MzjJ1%k@p|D26s-(+(s;sct_VC zL5L;T5-y}V02K?6ne2gqh#&Tf*<1?)*0#{uDh>m>>>qZG>pwn^df|Ko{_aZj8XbUv zJJ~m1-G*3FnA&#G8cHW^c_>0>6tNUfNoWMX0E9)2628U#QN%vziP?TAy;jwpw@@kG zvoTA;Pe3J}z9~8V3NAEqB>%oOVl8RLFE?)gqiEQKRVoasQPIJnFveR~n>emF$(|q|cZwp?8=GXN{TlanJ`W z7qMy0{ev@Q#n0K+IvjtLzRX&Vu6l@{v+#)Yh{{~DICP5BYsGiX9kWS96E+?{{IX*q zqv(tce{Y}rjIcF7XGh;2Cksfb9LwDfjLbDVKy&aLt%+N6DbiYiBUUhN*IcN8F#F{v zJ;Fv#WbN-x=f7;5Yt{IFTi5-6t#m&de$FwOhc7bW%$1xwnR1qmvwhIppgoD0Fmjkn zyodQYTifj{`sR`jq^0rAyNhO%hz{8=Dt)$Z$?%EYbFCYGw|-D4-g2kWiq2;HpYfu) zX}^i9S1V^*=>yQ;qP)Uf1ZgE_c)sBb+g$r@TYgUV@r@1etuK7a6H-r?$&E8}DGie1 W={FmX%q9^n@@lvx0=`13DgHl^U6K0$ From feeb041b0ec5950d1954b9f32b1faa14a268f756 Mon Sep 17 00:00:00 2001 From: "inwoo.park" Date: Tue, 24 Feb 2026 01:18:09 +0900 Subject: [PATCH 08/17] CHG: processdata.py - Fix bug in processing projection2d. --- src/m/plot/processdata.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/m/plot/processdata.py b/src/m/plot/processdata.py index e770adc23..dfb455890 100644 --- a/src/m/plot/processdata.py +++ b/src/m/plot/processdata.py @@ -158,9 +158,9 @@ def processdata(md, data, options): # }}} - # layer procesiong? {{{ + # layer processing? {{{ if options.getfieldvalue('layer',0)>=1: - procdata=project2d(md,data,options.getfieldvalue('layer')) + procdata=project2d(md,procdata,options.getfieldvalue('layer')) if options.getfieldvalue('depthaverage',0): raise Exception('Error: "depthaverage" option is not supported in Python.') From a90c4bd54f6318af4e3020f0ebe0c850a78fd149 Mon Sep 17 00:00:00 2001 From: "inwoo.park" Date: Tue, 24 Feb 2026 01:20:27 +0900 Subject: [PATCH 09/17] CHG: plot_unit.py - Fix bugs processing np.ma.core.MaskedArray dataset. --- src/m/plot/plot_unit.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/m/plot/plot_unit.py b/src/m/plot/plot_unit.py index 85c11b73b..0706d18a9 100644 --- a/src/m/plot/plot_unit.py +++ b/src/m/plot/plot_unit.py @@ -174,10 +174,10 @@ def plot_unit(x, y, z, elements, data, is2d, isplanet, datatype, options, fig, a if is2d: if np.ma.is_masked(data): if hasattr(np, 'isin'): #Numpy 2017+ - tmp = np.isin(index, np.where(data.mask)) + tmp = np.isin(range(len(data)), np.where(data.mask)) else: #For backward compatibility - tmp = np.in1d(index, np.where(data.mask)) - EltMask = np.asarray([np.any(tmp) for index in elements]) + tmp = np.in1d(range(len(data)), np.where(data.mask)) + EltMask = np.asarray([np.any(tmp[index]) for index in elements]) triangles = mpl.tri.Triangulation(x, y, elements, EltMask) else: triangles = mpl.tri.Triangulation(x, y, elements) From 819aec143e6479182ce2779d47ec3bcdd56f7c90 Mon Sep 17 00:00:00 2001 From: "inwoo.park" Date: Tue, 24 Feb 2026 13:20:32 +0900 Subject: [PATCH 10/17] CHG: plot_landsat - Check increasing order for md.radarpower.x and md.radarpower.y values. Synchronize Python > Matlab. --- src/m/plot/plot_landsat.m | 296 ++++++++++++++++++++----------------- src/m/plot/plot_landsat.py | 15 +- 2 files changed, 169 insertions(+), 142 deletions(-) diff --git a/src/m/plot/plot_landsat.m b/src/m/plot/plot_landsat.m index 56dba6ce7..147f6fabd 100644 --- a/src/m/plot/plot_landsat.m +++ b/src/m/plot/plot_landsat.m @@ -1,152 +1,170 @@ -% Explain -% This funtion loads Landsat Image Mosaic Antarctica (LIMA) for background image. -% -% Usage -% plot_landsat(md,data,options,plotlines,plotcols,i), -% function plot_landsat(md,data,options,plotlines,plotcols,i), + % Explain + % This funtion loads Landsat Image Mosaic Antarctica (LIMA) for background image. + % + % Usage + % plot_landsat(md,data,options,plotlines,plotcols,i), + % + + %process mesh and data + [x2d y2d z2d elements2d is2d isplanet]=processmesh(md,[],options); + [data datatype]=processdata(md,data,options); + + %check is2d + if ~is2d, + error('buildgridded error message: gridded not supported for 3d meshes, project on a layer'); + end -%process mesh and data -[x2d y2d z2d elements2d is2d isplanet]=processmesh(md,[],options); -[data datatype]=processdata(md,data,options); - -%check is2d -if ~is2d, - error('buildgridded error message: gridded not supported for 3d meshes, project on a layer'); -end - -%Get some options -transparency = getfieldvalue(options,'transparency',.3); -highres = getfieldvalue(options,'highres',0); - -% get xlim, and ylim -xlim=getfieldvalue(options,'xlim',[min(x2d) max(x2d)])/getfieldvalue(options,'unit',1); -ylim=getfieldvalue(options,'ylim',[min(y2d) max(y2d)])/getfieldvalue(options,'unit',1); - -if md.mesh.epsg == 3031 & (isempty(md.radaroverlay.pwr) | isempty(md.radaroverlay.x) | isempty(md.radaroverlay.y) | length(size(md.radaroverlay.pwr)) < 3), % Antarctica region {{{ - if highres, - disp(' LIMA with geotiff'), % {{{ - disp('WARNING : this image shoud be collected with geocoded tif file'); - % find merged mosaic landsat image {{{ - limapath = {'/drive/project_inwoo/issm/Data/LIMA/AntarcticaLandsat.tif'}; - pos = zeros(length(limapath),1); - for ii = 1:length(limapath) - if exist(limapath{ii}), pos(ii) = 1; end + %Get some options + transparency = getfieldvalue(options,'transparency',.3); + highres = getfieldvalue(options,'highres',0); + + %Get xlim, and ylim + xlim=getfieldvalue(options,'xlim',[min(x2d) max(x2d)])/getfieldvalue(options,'unit',1); + ylim=getfieldvalue(options,'ylim',[min(y2d) max(y2d)])/getfieldvalue(options,'unit',1); + + pwr = md.radaroverlay.pwr; + xm = md.radaroverlay.x; + ym = md.radaroverlay.y; + nx = numel(xm); + ny = numel(ym); + + if md.mesh.epsg == 3031 & (isempty(md.radaroverlay.pwr) | isempty(md.radaroverlay.x) | isempty(md.radaroverlay.y) | length(size(md.radaroverlay.pwr)) < 3), % Antarctica region {{{ + if highres, + disp(' LIMA with geotiff'), % {{{ + disp('WARNING : this image shoud be collected with geocoded tif file'); + % find merged mosaic landsat image {{{ + limapath = {'/drive/project_inwoo/issm/Data/LIMA/AntarcticaLandsat.tif'}; + pos = zeros(length(limapath),1); + for ii = 1:length(limapath) + if exist(limapath{ii}), pos(ii) = 1; end + end + limapath = limapath{find(pos)}; + fprintf(' LIMA path is %s\n', limapath); + % }}} + + % read image + im = imread(limapath); + + % Region of LIMA data set + info = gdalinfo(limapath); % get geotiff info + xm = info.xmin + info.dx*[0:info.nx-1]; + ym = info.ymax - info.dy*[0:info.ny-1]; + + % find region of model at LIMA + offset = 1e+4; + posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset)); + posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset)); + % }}} + else + disp(' LIMA with reduced tiff'), + % find merged mosaic landsat image {{{ + limapath = {'/drive/project_inwoo/issm/Data/LIMA/tiff_90pct/00000-20080319-092059124.tif'}; + pos = zeros(length(limapath),1); + for ii = 1:length(limapath) + if exist(limapath{ii}), pos(ii) = 1; end + end + + if sum(pos) == 0, + fprintf('download website : https://lima.usgs.gov/fullcontinent.php\n'); + error('Landsat image at Antarctic region should be downloaded at above website'); + end + limapath = limapath{find(pos)}; + fprintf(' LIMA path is %s\n', limapath); + % }}} + + % read image + im = imread(limapath); + + % Region of LIMA data set + info = gdalinfo(limapath); % get geotiff info + xm = info.xmin + info.dx*[0:info.nx-1]; + ym = info.ymax - info.dy*[0:info.ny-1]; + + % find region of model at LIMA + offset = 1e+4; + posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset)); + posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset)); end - limapath = limapath{find(pos)}; - fprintf(' LIMA path is %s\n', limapath); + + % update region of radaroverlay + md.radaroverlay.x = xm(posx); + md.radaroverlay.y = ym(posy); + md.radaroverlay.pwr = im(posy, posx,:); % }}} + elseif length(size(md.radaroverlay.pwr)) == 3 + % it already contains LIMA image. + elseif md.mesh.epsg == 3431 & (isempty(md.radaroverlay.pwr) | isempty(md.radaroverlay.x) | isempty(md.radaroverlay.y) | length(size(md.radaroverlay.pwr)) < 3), % Greenladn region + error('Greenland region is not yet available.'); + else + error('Check md.mesh.epsg, available Landsat regeion is at Antarctica (EPSG:3031)'); + end + + %Check dataset + assert(ndims(pwr) ~= 3, 'Error: Check ndims(md.radaroverlay.pwr) should be equal to 3.'); + assert(size(pwr) == [nx, ny, 3], 'Error: Given md.radaroverlay.pwr shoule be equal to (nx, ny, 3)'); - % read image - im = imread(limapath); + if any(diff(xm)) < 0 + disp('WARNING: md.radaroverlay.x should be increasing order.') + xm = flip(xm); + pwr= flip(pwr, 1); + end + if any(diff(ym)) < 0 + disp('WARNING: md.radaroverlay.y should be increasing order.') + ym = flip(ym); + pwr= flip(pwr, 2); + end - % Region of LIMA data set - info = gdalinfo(limapath); % get geotiff info - xm = info.xmin + info.dx*[0:info.nx-1]; - ym = info.ymax - info.dy*[0:info.ny-1]; + %Process image from model + final = double(pwr)/double(max(md.radaroverlay.pwr(:))); %rescale between 0 and 1 - % find region of model at LIMA - offset = 1e+4; - posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset)); - posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset)); - % }}} + %Prepare grid + if size(md.radaroverlay.x,1)==1 | size(md.radaroverlay.x,2)==1, + data_grid=InterpFromMeshToGrid(elements2d,x2d/getfieldvalue(options,'unit',1),y2d/getfieldvalue(options,'unit',1),data,xm,ym,NaN); + %data_grid=InterpFromMeshToGrid(md.mesh.elements,md.mesh.x/getfieldvalue(options,'unit',1),md.mesh.y/getfieldvalue(options,'unit',1),data,x_m,y_m,NaN); else - disp(' LIMA with reduced tiff'), - % find merged mosaic landsat image {{{ - limapath = {'/drive/project_inwoo/issm/Data/LIMA/tiff_90pct/00000-20080319-092059124.tif'}; - pos = zeros(length(limapath),1); - for ii = 1:length(limapath) - if exist(limapath{ii}), pos(ii) = 1; end - end - - if sum(pos) == 0, - fprintf('download website : https://lima.usgs.gov/fullcontinent.php\n'); - error('Landsat image at Antarctic region should be downloaded at above website'); - end - limapath = limapath{find(pos)}; - fprintf(' LIMA path is %s\n', limapath); - % }}} + X = md.radaroverlay.x; + Y = md.radaroverlay.y; + data_grid=InterpFromMeshToMesh2d(elements2d,x2d,y2d,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X)); + %data_grid=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X)); + xm=X(1,:); ym=Y(:,1); + end - % read image - im = imread(limapath); + data_nan=isnan(data_grid); + if exist(options,'caxis'), + caxis_opt=getfieldvalue(options,'caxis'); + data_grid(find(data_gridcaxis_opt(2)))=caxis_opt(2); + data_min=caxis_opt(1); + data_max=caxis_opt(2); + else + data_min=min(data_grid(:)); + data_max=max(data_grid(:)); + end + colorm = getcolormap(options); + image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(colorm)/(data_max-data_min))),colorm); - % Region of LIMA data set - info = gdalinfo(limapath); % get geotiff info - xm = info.xmin + info.dx*[0:info.nx-1]; - ym = info.ymax - info.dy*[0:info.ny-1]; + alpha=ones(size(data_grid)); + alpha(find(~data_nan))=transparency; + alpha=repmat(alpha,[1 1 3]); - % find region of model at LIMA - offset = 1e+4; - posx = find((xm > xlim(1)-offset).* (xm < xlim(2)+offset)); - posy = find((ym > ylim(1)-offset).* (ym < ylim(2)+offset)); - % }}} + final=alpha.*final+(1-alpha).*image_rgb; + + %Select plot area + subplotmodel(plotlines,plotcols,i,options); + + h=imagesc(xm*getfieldvalue(options,'unit',1),ym*getfieldvalue(options,'unit',1),final); + + %last step: mesh gridded? + if exist(options,'edgecolor'), + A=elements(:,1); B=elements(:,2); C=elements(:,3); + patch('Faces',[A B C],'Vertices', [x y z],'FaceVertexCData',data_grid(1)*ones(size(x)),'FaceColor','none','EdgeColor',getfieldvalue(options,'edgecolor')); end - % update region of radaroverlay - md.radaroverlay.x = xm(posx); - md.radaroverlay.y = ym(posy); - md.radaroverlay.pwr = im(posy, posx,:); - % }}} -elseif length(size(md.radaroverlay.pwr)) == 3 - % it already contains LIMA image. -elseif md.mesh.epsg == 3431 & (isempty(md.radaroverlay.pwr) | isempty(md.radaroverlay.x) | isempty(md.radaroverlay.y) | length(size(md.radaroverlay.pwr)) < 3), % Greenladn region - error('Greenland region is not yet available.'); -else - error('Check md.mesh.epsg, available Landsat regeion is at Antarctica (EPSG:3031)'); -end - -%Process image from model -final = double(md.radaroverlay.pwr)/double(max(md.radaroverlay.pwr(:))); %rescale between 0 and 1 - -%Prepare grid -if size(md.radaroverlay.x,1)==1 | size(md.radaroverlay.x,2)==1, - x_m = md.radaroverlay.x; - y_m = md.radaroverlay.y; - data_grid=InterpFromMeshToGrid(elements2d,x2d/getfieldvalue(options,'unit',1),y2d/getfieldvalue(options,'unit',1),data,x_m,y_m,NaN); - %data_grid=InterpFromMeshToGrid(md.mesh.elements,md.mesh.x/getfieldvalue(options,'unit',1),md.mesh.y/getfieldvalue(options,'unit',1),data,x_m,y_m,NaN); -else - X = md.radaroverlay.x; - Y = md.radaroverlay.y; - data_grid=InterpFromMeshToMesh2d(elements2d,x2d,y2d,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X)); - %data_grid=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X)); - x_m=X(1,:); y_m=Y(:,1); -end - -data_nan=isnan(data_grid); -if exist(options,'caxis'), - caxis_opt=getfieldvalue(options,'caxis'); - data_grid(find(data_gridcaxis_opt(2)))=caxis_opt(2); - data_min=caxis_opt(1); - data_max=caxis_opt(2); -else - data_min=min(data_grid(:)); - data_max=max(data_grid(:)); -end -colorm = getcolormap(options); -image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(colorm)/(data_max-data_min))),colorm); - -alpha=ones(size(data_grid)); -alpha(find(~data_nan))=transparency; -alpha=repmat(alpha,[1 1 3]); - -final=alpha.*final+(1-alpha).*image_rgb; - -%Select plot area -subplotmodel(plotlines,plotcols,i,options); - -h=imagesc(x_m*getfieldvalue(options,'unit',1),y_m*getfieldvalue(options,'unit',1),final); - -%last step: mesh gridded? -if exist(options,'edgecolor'), - A=elements(:,1); B=elements(:,2); C=elements(:,3); - patch('Faces',[A B C],'Vertices', [x y z],'FaceVertexCData',data_grid(1)*ones(size(x)),'FaceColor','none','EdgeColor',getfieldvalue(options,'edgecolor')); -end - -%Apply options -if ~isnan(data_min), - options=changefieldvalue(options,'caxis',[data_min data_max]); % force caxis so that the colorbar is ready -end -options=addfielddefault(options,'axis','xy equal off'); % default axis -applyoptions(md,data,options); -end + %Apply options + if ~isnan(data_min), + options=changefieldvalue(options,'caxis',[data_min data_max]); % force caxis so that the colorbar is ready + end + options=addfielddefault(options,'axis','xy equal off'); % default axis + applyoptions(md,data,options); + end diff --git a/src/m/plot/plot_landsat.py b/src/m/plot/plot_landsat.py index c54fe4959..e6b2af3a4 100644 --- a/src/m/plot/plot_landsat.py +++ b/src/m/plot/plot_landsat.py @@ -45,8 +45,8 @@ def plot_landsat(md,data,options,fig,axgrid,gridindex): isunit = options.getfieldvalue('unit',1) #Get xlim, and ylim - xlim=options.getfieldvalue('xlim',np.array([min(x2d),max(x2d)]))/isunit - ylim=options.getfieldvalue('ylim',np.array([min(y2d),max(y2d)]))/isunit + xlim=np.array(options.getfieldvalue('xlim',[min(x2d),max(x2d)]))/isunit + ylim=np.array(options.getfieldvalue('ylim',[min(y2d),max(y2d)]))/isunit pwr = md.radaroverlay.pwr xm = md.radaroverlay.x @@ -113,8 +113,17 @@ def plot_landsat(md,data,options,fig,axgrid,gridindex): if np.ndim(pwr) != 3: raise Exception('Error: Check np.ndim(md.radaroverlay.pwr) should be equal to 3.') + if np.any(np.diff(xm) < 0): + print('WARNING: md.radaroverlay.x should be increasing order.') + xm = np.flip(xm) + pwr= np.flip(pwr,axis=0) + if np.any(np.diff(md.radaroverlay.y) < 0): + print('WARNING: md.radaroverlay.y should be increasing order.') + ym = np.flip(ym) + pwr= np.flip(pwr,axis=1) + #Check image size - #shape of image should be (ny, nx, band) + #shape of image should be (nx, ny, band) nx = len(xm) ny = len(ym) if (np.shape(pwr)[0]==nx) & (np.shape(pwr)[1]==ny): From 2f8b97b89a546a756df0fd35b8caf254a6b6eb56 Mon Sep 17 00:00:00 2001 From: "inwoo.park" Date: Tue, 24 Feb 2026 01:18:09 +0900 Subject: [PATCH 11/17] CHG: processdata.py - Fix bug in processing projection2d. --- src/m/plot/processdata.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/m/plot/processdata.py b/src/m/plot/processdata.py index e770adc23..dfb455890 100644 --- a/src/m/plot/processdata.py +++ b/src/m/plot/processdata.py @@ -158,9 +158,9 @@ def processdata(md, data, options): # }}} - # layer procesiong? {{{ + # layer processing? {{{ if options.getfieldvalue('layer',0)>=1: - procdata=project2d(md,data,options.getfieldvalue('layer')) + procdata=project2d(md,procdata,options.getfieldvalue('layer')) if options.getfieldvalue('depthaverage',0): raise Exception('Error: "depthaverage" option is not supported in Python.') From c1ce7e4a757b2602c79fcf837fbe93bc3c1bebe7 Mon Sep 17 00:00:00 2001 From: "inwoo.park" Date: Tue, 24 Feb 2026 01:20:27 +0900 Subject: [PATCH 12/17] CHG: plot_unit.py - Fix bugs processing np.ma.core.MaskedArray dataset. --- src/m/plot/plot_unit.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/m/plot/plot_unit.py b/src/m/plot/plot_unit.py index 85c11b73b..0706d18a9 100644 --- a/src/m/plot/plot_unit.py +++ b/src/m/plot/plot_unit.py @@ -174,10 +174,10 @@ def plot_unit(x, y, z, elements, data, is2d, isplanet, datatype, options, fig, a if is2d: if np.ma.is_masked(data): if hasattr(np, 'isin'): #Numpy 2017+ - tmp = np.isin(index, np.where(data.mask)) + tmp = np.isin(range(len(data)), np.where(data.mask)) else: #For backward compatibility - tmp = np.in1d(index, np.where(data.mask)) - EltMask = np.asarray([np.any(tmp) for index in elements]) + tmp = np.in1d(range(len(data)), np.where(data.mask)) + EltMask = np.asarray([np.any(tmp[index]) for index in elements]) triangles = mpl.tri.Triangulation(x, y, elements, EltMask) else: triangles = mpl.tri.Triangulation(x, y, elements) From 4f55e8eab39c3a51e3e1fca3337f806609912341 Mon Sep 17 00:00:00 2001 From: "inwoo.park" Date: Tue, 24 Feb 2026 13:56:20 +0900 Subject: [PATCH 13/17] CHG: plot_landsat.py - Enable "log" scale plot. --- src/m/plot/plot_landsat.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/m/plot/plot_landsat.py b/src/m/plot/plot_landsat.py index e6b2af3a4..5a0ed494b 100644 --- a/src/m/plot/plot_landsat.py +++ b/src/m/plot/plot_landsat.py @@ -155,11 +155,19 @@ def plot_landsat(md,data,options,fig,axgrid,gridindex): else: data_min=np.nanmin(data_grid) data_max=np.nanmax(data_grid) + caxis_opt=[data_min, data_max] # Back-up caxis for 'applyoptions' options.addfielddefault('colormap',plt.cm.viridis) cmap = getcolormap(copy.deepcopy(options)) #TODO: Matlab version #image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(colorm)/(data_max-data_min))),colorm); + #Set log scale + if options.exist('log'): + #NOTE: Tricy part for log scale dataset. "log" scale option does not rely on "processdata.py" function. + data_grid=np.log(data_grid)/np.log(options.getfieldvalue('log')) + data_min =np.log(data_min)/np.log(options.getfieldvalue('log')) + data_max =np.log(data_max)/np.log(options.getfieldvalue('log')) + #NOTE: Python version for ind2rgb image_rgb = cmap((data_grid-data_min)/(data_max-data_min)) @@ -176,6 +184,7 @@ def plot_landsat(md,data,options,fig,axgrid,gridindex): xmax = max(xm)/isunit ymin = min(ym)/isunit ymax = max(ym)/isunit + #Draw RGB image h=ax.imshow(final, extent=[xmin, xmax, ymin, ymax], origin='lower') @@ -188,6 +197,6 @@ def plot_landsat(md,data,options,fig,axgrid,gridindex): #Apply options if ~np.isnan(data_min): - options.changefieldvalue('caxis',[data_min, data_max]) # force caxis so that the colorbar is ready + options.changefieldvalue('caxis',caxis_opt) # force caxis so that the colorbar is ready options.addfielddefault('axis','xy equal off') # default axis applyoptions(md,data,options,fig,axgrid,gridindex) From 6de6bed437771a9dd4b2bb2a8edfa379002b526f Mon Sep 17 00:00:00 2001 From: "inwoo.park" Date: Tue, 24 Feb 2026 15:35:35 +0900 Subject: [PATCH 14/17] CHG: Fix minor: Treaing colorbar axis. --- src/m/plot/applyoptions.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/m/plot/applyoptions.py b/src/m/plot/applyoptions.py index 80e6b0862..e4074a73d 100644 --- a/src/m/plot/applyoptions.py +++ b/src/m/plot/applyoptions.py @@ -231,12 +231,12 @@ def applyoptions(md, data, options, fig, axgrid, gridindex): if options.exist('log'): formatter = mpl.ticker.LogFormatterSciNotation(base=options.getfieldvalue('log')) - try: - #NOTE: axis made by ImageGrid contains "cax". If axis generated by "plt.subplots()", we have to manually generate cax... + #NOTE: axis made by ImageGrid contains "cax". If axis generated by "plt.subplots()", we have to manually generate cax... + if hasattr(ax, 'cax'): # type(ax) == mpl_toolkits.axes_grid1.mpl_axes.Axes cb = mpl.colorbar.ColorbarBase(ax.cax, cmap=cmap, norm=norm, extend=extend, format=formatter) - except: + else: from mpl_toolkits.axes_grid1 import make_axes_locatable driver = make_axes_locatable(ax) cax = driver.append_axes(options.getfieldvalue('colorbarpos','right'), From d7786f5abc67cd8d9ae7db2d10821f2a3b1f0643 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alicia=20Br=C3=A5tner?= Date: Tue, 24 Feb 2026 14:50:42 -0500 Subject: [PATCH 15/17] updated archive test 245 --- test/Archives/Archive245.arch | Bin 1588 -> 1588 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/test/Archives/Archive245.arch b/test/Archives/Archive245.arch index afe5124371f5b6051ec74e43bbdea82b3022ef7a..4dafb9726d7726756a7ce1ec2bde1699e9ae2d93 100644 GIT binary patch delta 1167 zcmWkseN>Zm6n^Id!GVRO8zUKp9`Qg53prBOeGw!g`4maZkCA*Y=9{0RfGs_mVByAK z!a-swaKyn81@VyMyf2LJfcTUT1x8p1OpHX1v?g8udd~gb=RWs$?|rIms%;cISt-ho zxYiFyD(>=)6*&&dLtNHr<*8#1Jd#QZSemo&CaFXYvq|WJ=px)5_qoz6U`?pjvjh$S zNV)#>3a6B$#3suBTsMj||19kORJD(_+3qDt#fMpFjl%dA68(aQw=8)y&vHi>B6c=o=#&k34ovYK&=CTOk5 zE|qIR?`$`5HMQJ_a3IJVprk(xwBy{Ww&;n;1`k$~3 z4Lu6zkLyT|*`k5I(4btJ9ZKJ`j?$;wu^?|e_Nd@4L07NW=c^+~7ku7fw) zL){H9F%w%UJh&e^{_8$$S$PCncW;w^R=cQyD(4XC`W6xlwfY*~+}}`F^~%-DuO4n~ z=e^0CVL`9QY0Ijs26*)NWX9Sy0!FUHmGCRZVeH%5mX^s%7+Lx2laR%F+^9GEZP|Oz zyFJ(8wyhN1ybh8&mYGq0!HT3LIy?JQYw$)RGllBoc$=YMqZoPOki@XTV9YP3wlw5m zKFQg@ohGT|Z5+vJFT7T6R!z8uq`uar9SF`u{@UB=OF@`V`E#^48TH*i{PLRttoL|e zDmg+@GRbk1hyN@Rs4=sC9t(@kNF!{JaQ|Yi_7w`u`@Fw=7LlD*CmBS%WYJpZhK&Q$ zWBw^eU-A3=)|n7DTMM$Y5M$LlPMT5gS=XCMBWT^p0PkXaAfaH`h4h@q>tD9vIl~CPyam=uWrB zF%z1tQ6g~jyZqYBIfEe$3>Ypp%Z%sW0%mb~%&=2VH%|?FI@*h%onf!JJlW5pH!cYhJ(*rqdy3I}| zO2!5Vl*x9M)BuT+dqD&t zl@kA+ax8hQc)inu}RA!I(!{$k{0M-P^83>TM6pma3tIHTH*s_2HY8;3|3TAz1$mP3PX z?wyVUL3$Xoq*buF3D7lfJnB}VgtqAI!noRLJ7}5lLV9Zhm}?C+zLT?1ueGs4qV({6 zsPA&UkbuX*AGm36PQmb#?7S82%fVE_zrE_a1QmmDa&Ov-5(|B`xuj9a_0cUmiRDQ zl5js`*7aB%No9uTR($*+zShP4%Aa2&sf_mPhkjbXdaog-=?v~?3Gz6ecrJyR)|ZR@ zw{?;=eIzB$N<*|bcPtQX^SE#m58hQ}PUT_2^AJV)9c(b`#67$cKV+*v+!BfN{Dnp5 zQ*j=}^;}4zK%P~wk z1K;RZGF zs0KfDm8kH=_4KlgS8l7Zu?unjvC6?M(fBal!7Ytn$--MUxOj`jQz4t&gTcIZE+w6P z1^sw+!@EZ|VVL)Ro!k~C!7vZWWWfp;Bj+*)XPMBjo%z+Javo}HNY+mFD4;`?%^L`i zqH=m{olbyZp4d0reoXVHy~&r8MqaYewofU9Qn1X{gh?9=-JEA>DZNle!u` z66=k&HF2Z@X~L{_QK$q8=GW<_f-%e|hkNpO(m;FekFC$GG-%>3tH`qupu4eACweJ? zE^<^}a`zbYj}g_LgJ$SYbJqOg!Gyk8*ODg%3FyJ=FVt`6NTD<hN>7HKIgVz znCah0pQ$P6PU3Q{gp>zkzn+T7*xm+RFOHn^3{pT_xLS0{i3IIdj$1yD2Ig<)MJIz; zP_KMzAaEImIZSt07Emo!wd|s F{{gfsE$#pS From 01161cc2f95107082f2aa8aafcf3e4d8ed1971e3 Mon Sep 17 00:00:00 2001 From: "inwoo.park" Date: Wed, 25 Feb 2026 15:41:30 +0900 Subject: [PATCH 16/17] CHG: Synchronize Matlab > Python for "ProfileValues.py". --- src/m/interp/ProfileValues.py | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 src/m/interp/ProfileValues.py diff --git a/src/m/interp/ProfileValues.py b/src/m/interp/ProfileValues.py new file mode 100644 index 000000000..5eee099fa --- /dev/null +++ b/src/m/interp/ProfileValues.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python3 +from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d +from InterpFromMeshToMesh3d import InterpFromMeshToMesh3d +from project2d import project2d +import numpy as np +from numpy import ceil, mean, linspace, arange, transpose, reshape, ones, shape + +def ProfileValues(md,data,xprof,yprof,resolution): + """ + PROFILEVALUES - compute the value of a field on a vertical profile + + This routine gets the value of a given field of the model on + a point given by its coordinates + + Usage: + z,data=ProfileValues(md,data,xcoord,ycoord,resolution) + """ + + #Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system + offset=10e-3 + bed=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.base,1),xprof,yprof)+offset + surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),xprof,yprof)-offset + + #Some useful parameters + layers=ceil(mean(md.geometry.thickness)/resolution) + Z=arange(bed,surface+resolution,resolution) + Z=transpose(Z) + X=xprof*ones(shape(Z)) + Y=yprof*ones(shape(Z)) + data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X,Y,Z,np.nan) + + #Get np.ndarray in list + data_interp=data_interp[0] + + return Z, data_interp From c8d298f496661e1eb3bb3621c942b41a43f63f0b Mon Sep 17 00:00:00 2001 From: "inwoo.park" Date: Wed, 25 Feb 2026 15:42:01 +0900 Subject: [PATCH 17/17] CHG: InterpFromMeshToMesh3d.py - Fix bugs in checking number of arguments as input. This method is only for matlab, not python. This section is changed to annotation. --- src/m/modules/InterpFromMeshToMesh3d.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/m/modules/InterpFromMeshToMesh3d.py b/src/m/modules/InterpFromMeshToMesh3d.py index 743ded6cf..b5805b06c 100644 --- a/src/m/modules/InterpFromMeshToMesh3d.py +++ b/src/m/modules/InterpFromMeshToMesh3d.py @@ -18,10 +18,10 @@ def InterpFromMeshToMesh3d(index, x, y, z, data, x_prime, y_prime, z_prime, defa """ # Check usage - nargs = len(args) - if nargs != 9: - print(InterpFromMeshToMesh3d.__doc__) - raise Exception('Wrong usage (see above)') + #nargs = len(args) + #if nargs != 9: + # print(InterpFromMeshToMesh3d.__doc__) + # raise Exception('Wrong usage (see above)') # Call Python module data_prime = InterpFromMeshToMesh3d_python(index, x, y, z, data, x_prime, y_prime, z_prime, default_value)